Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

DetectorPhotonTransporter - source file

// $Id: DetectorPhotonTransporter.cc,v 1.15 2005/04/17 15:50:08 thea Exp $
// Author: D.Demarco, M.Pallavicini

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: DetectorPhotontransporter                                            *
 *  Package: Optics                                                          *
 *  Coordinator: Alessandro.Thea                                             *
 *                                                                           *
 *****************************************************************************/

//______________________________________________________________________________
//  
//  DetectorPhotonTransporter
//  =========================
//
//  Abstarct base class providing the interface of a generic photon
//  transporter. Every element in the detector that is able to interact
//  with photons is a DetectorPhotonTransporter. 
//  The DetectorPhotonTransporters are represented as cylinders with position
//  and height (fDZdown+fDZup).
//

#include "DetectorPhotonTransporter.hh"
#include "Etypes.hh"

ClassImp(DetectorPhotonTransporter)

using namespace TMath;

//______________________________________________________________________________
DetectorPhotonTransporter::DetectorPhotonTransporter(): EsafConfigurable(),
    fPos(0,0,0), fDZdown(0), fDZup(0), fR(0) {
    //
    // Constructor
    //

//    fLastCylInt.first=0;
//    fLastCylInt.second=new Double_t[2];
}

//______________________________________________________________________________
 DetectorPhotonTransporter::~DetectorPhotonTransporter(){
    //
    // Destructor
    // 
    
//    delete [] fLastCylInt.second;
}


//______________________________________________________________________________
 Bool_t DetectorPhotonTransporter::IsInside( Photon *p ) const {
    //
    // True if the photon is inside the boundaries of the transporter
    //
    
    return ((p->pos[Z]-Bottom()) > -kTolerance && (p->pos[Z]-Top()) < kTolerance &&
            (p->pos.Perp()-Radius()) < kTolerance);
}

//______________________________________________________________________________
 Double_t DetectorPhotonTransporter::CylinderIntersection( const TVector3 &pos,
       const TVector3 &dir, Bool_t zero) const {
    //
    // Calculate the next intersection of p with this transporter boundaries
    //
    
    return CylinderIntersection( pos, dir, fR, fPos[Z]+fDZup, fPos[Z]-fDZdown, zero);
}

//______________________________________________________________________________
 Double_t DetectorPhotonTransporter::CylinderIntersection( const Photon *p,
       Bool_t zero) const {
    //
    // Calculate the next intersection of p with this transporter boundaries
    //
    
    return CylinderIntersection( p, fR, fPos[Z]+fDZup, fPos[Z]-fDZdown, zero );
}

//______________________________________________________________________________
 Double_t DetectorPhotonTransporter::CylinderIntersection( const Photon *p,
    Double_t radius, Double_t  zup, Double_t zdown, Bool_t zero ) const {
    // 
    // Finds, if exists, the next interaction of p over cyl.  In any case saves
    // the distances between p->pos and the two int points in fLastCylInt
    // 

    return CylinderIntersection(p->pos, p->dir, radius, zup, zdown, zero );
}

//______________________________________________________________________________
 Double_t DetectorPhotonTransporter::CylinderIntersection( const TVector3 &pos,
    const TVector3 &dir, Double_t radius, Double_t  zup, Double_t zdown,
    Bool_t zero ) const {
    // 
    // Finds, if exists, the next interaction of p over cyl.  In any case saves
    // the distances between p->pos and the two int points in fLastCylInt
    //

    Double_t dist[4] = {-kHuge,-kHuge,-kHuge,-kHuge};
    // array of intersections, 0, 1 basis, 3, 4 lateral surface

    // int with the basis
    if ( Abs(dir[Z]) > kTolerance ) {
        dist[0] = (zup-pos[Z])/dir[Z];
        dist[1] = (-zdown-pos[Z])/dir[Z];
    }

    Double_t a, b, c;
    a = dir[X]*dir[X]+dir[Y]*dir[Y];
    b = 2*(dir[X]*pos[X]+dir[Y]*pos[Y]);
    c = pos[X]*pos[X]+pos[Y]*pos[Y]-radius*radius;


    Double_t delta = (b*b)-4*a*c;

    if ( a > 0 ) 
        if ( Abs(delta) < kTolerance ) {
            dist[2] = -b/(2*a);
        } else if ( delta >= kTolerance ) {
            dist[2] = (-b-Sqrt(delta))/(2*a);
            dist[3] = (-b+Sqrt(delta))/(2*a);
        }


    Double_t threshold = zero ? -kTolerance : kTolerance;
    Double_t dt = MaxElement(4,dist); 
    for(Int_t i(0); i<4; i++) 
        if ( dist[i] > threshold && dist[i] < dt ) dt = dist[i];

    // inbetween +/- ktolerance dt is basically 0
    if ( zero && dt > -kTolerance && dt < kTolerance) dt = 0;
    //MsgForm(EsafMsg::Info,"dt %3ft,d(%8f,%8f,%8f,%8f)",dt,dist[0],dist[1],dist[2],dist[3]); 
    return dt;

/*************************************************************************** 
 *
 * Old code, kept for the time being as reference
 *
 * *************************************************************************
 
    // intesection point
    TVector3 intPoint(0,0,0), uDir;
    uDir = dir.Unit();

    Double_t a, b, c, dup, ddown, dummy[2];

    fLastCylInt.first=0;
    fLastCylInt.second[0]=fLastCylInt.second[1]=0;

    // distances to the upper and lower planes
    if ( TMath::Abs(uDir[Z]) >= kTolerance ){
        dup = (zUp-pos[Z])/uDir[Z];
        ddown = (zDown-pos[Z])/uDir[Z];

    }

    if((TMath::Abs(uDir.Theta()) < kTolerance) || 
            (TMath::Abs(uDir.Theta() - TMath::Pi()) < kTolerance)) {
        // photon directed along Z: special case
        // doesn't hit the cyl.
        if(pos.Perp()-radius > kTolerance) {
            return -1*kHuge;
        }
        // intersection with the bases
        if(uDir[Z] > 0 && pos[Z] < zDown || 
                uDir[Z] < 0 && pos[Z] > zDown && pos[Z] < zUp ){
            // lower base
            fLastCylInt.second[0]= ddown;
            fLastCylInt.second[1]= dup;
        }
        else if(uDir[Z] < 0 && pos[Z] > zUp || 
                uDir[Z] > 0 && pos[Z] > zDown && pos[Z] < zUp ){
            // upper base
            fLastCylInt.second[0]= dup;
            fLastCylInt.second[1]= ddown;
        } else {
            MsgForm(EsafMsg::Warning,"CylinderIntersection(): photon lost");
            return -2*kHuge;
        }
        return fLastCylInt.second[0];
    }

    // fill the coefficients of the equation
    a=uDir[X]*uDir[X]+uDir[Y]*uDir[Y];
    b=2*(uDir[X]*pos[X]+uDir[Y]*pos[Y]);
    c=pos[X]*pos[X]+pos[Y]*pos[Y]-radius*radius;
    // find the first intersection with the cylinder if it exists
    pair<int, double* > res;
    res=findRoots(a,b,c);

    if (res.first == 0) {
        // photon doesn't cross the cylinder
        // return something;
        return -3*kHuge;

    } else if (res.first == 1) {
        // one solution, tangent photon
        intPoint=pos+uDir*res.second[0];
        // check the sol to be inside cylinder
        if ( intPoint[Z] > zUp || intPoint[Z] < zDown )
            return -4*kHuge;
        fLastCylInt.first=1;
        fLastCylInt.second[0]=fLastCylInt.second[1]=res.second[0];
        return fLastCylInt.second[0];

    } else if ( TMath::Abs(uDir[Z]) < kTolerance ){ 
        // two crossings but no intersection with the bases
        // horizontal photon
        intPoint=pos+uDir*res.second[0];
        // check the sol to be inside cylinder
        if ( intPoint[Z] > zUp || intPoint[Z] < zDown )
            return -5*kHuge;

    }

    // 4 candidates left, just 2 of them on the borders of the cyl; in fact 3 if
    // the photons crosses the intersection between the side surface and one of
    // the bases
    fLastCylInt.first=2;

    // distances of the 4 intpoints from pos
    // 0: upper base plane
    // 1: lower base plane
    // 2: 1st lateral surface
    // 3: 2nd lateral surface
    Double_t intDist[4];
    intDist[0]=dup;
    intDist[1]=ddown;
    intDist[2]=res.second[0];
    intDist[3]=res.second[1];

    Int_t nTrueInts(0);
    Bool_t onSide, onBase;
    for(Int_t i(0); i<4; i++ ){
        intPoint = pos+uDir*intDist[i];
        onSide = TMath::Abs(intPoint.Perp()-radius) < kTolerance &&
            (intPoint[Z] < zUp && intPoint[Z] > zDown);
        onBase = (TMath::Abs(intPoint[Z]-zUp) < kTolerance ||
                TMath::Abs(intPoint[Z]-zDown) < kTolerance) &&
            intPoint.Perp() <= radius;
        if ( onSide || onBase )
            // if intDist[i] < kTolerance, save just 0
            dummy[nTrueInts++]=( TMath::Abs(intDist[i]) > kTolerance ? intDist[i] : 0 );
    }

    if ( nTrueInts == 3 ) {
        // check if two of the ints are the same point (one with the bases and
        // one with the side)
        if ( TMath::Abs(intDist[0]-intDist[2]) < kTolerance ||
                TMath::Abs(intDist[0]-intDist[3]) < kTolerance || 
                TMath::Abs(intDist[1]-intDist[2]) < kTolerance ||
                TMath::Abs(intDist[1]-intDist[3]) < kTolerance) {
            nTrueInts--;
        }

    }

    if ( nTrueInts == 2) {
        // last thing to do, sort the sols
        if ( dummy[0]*dummy[1] < 0 || dummy[0]+dummy[1] < 0 ){
            // opposite sign sols or sols < 0
            fLastCylInt.second[0]=TMath::Max(dummy[0],dummy[1]);
            fLastCylInt.second[1]=TMath::Min(dummy[0],dummy[1]);
        } else {
            // sols >= 0
            fLastCylInt.second[0]=TMath::Min(dummy[0],dummy[1]);
            fLastCylInt.second[1]=TMath::Max(dummy[0],dummy[1]);
        }
        if ( fLastCylInt.second[0] < 0 )
            return -6*kHuge;
        else 
            return fLastCylInt.second[0];

    } else if ( nTrueInts==0 ) {
        // non of the sols was on the cyl
        return -7*kHuge;
    } else {
        // dump to screen all the intersections
        Msg(EsafMsg::Warning)  << "pos" << pos << "   dir" << uDir << endl;
        for(Int_t i(0); i<4; i++ ){
            MsgForm(EsafMsg::Warning, "Intersection %d:",i);
            intPoint = pos+uDir*intDist[i];
            MsgForm(EsafMsg::Warning,"intPoint.PerP() = %f intPoint[Z] = %f",
                    intPoint.Perp(),intPoint[Z]);
            MsgForm(EsafMsg::Warning,"Radius = %f Zup = %f Zdown = %f",
                    radius, zUp, zDown);

            onSide = TMath::Abs(intPoint.Perp()-radius) < kTolerance &&
                (intPoint[Z] < zUp && intPoint[Z] > zDown);

            onBase = (TMath::Abs(intPoint[Z]-zUp) < kTolerance ||
                    TMath::Abs(intPoint[Z]-zDown) < kTolerance) &&
                intPoint.Perp()-radius <= kTolerance;
            MsgForm(EsafMsg::Warning," side = %f base = %f", onSide, onBase);
        }
        Msg(EsafMsg::Warning) << " nTrueInts = " << nTrueInts << MsgDispatch;
        MsgForm(EsafMsg::Panic,"CylinderIntersection(): something gone wrong!");
        return -8*kHuge;
    }
*******************************************************************************/
} 
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005