Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

DetectorTransportManager - source file

// $Id: DetectorTransportManager.cc,v 1.83 2005/10/31 14:43:56 thea Exp $
// Author: D.DeMarco

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

//_____________________________________________________________________________
//  
//   DetectorTransportManager
//   ========================
//
//   The raytracer of the inner detector.
//

#include <math.h>
#include <algorithm>
#include <iomanip>

#include "DetectorTransportManager.hh"
#include "utils.hh"
#include "EEvent.hh"
#include "Etypes.hh"
#include "PhotonsOnPupil.hh"
#include "ParentPhoton.hh"
#include "EDetectorPhotonAdder.hh"
#include "EConst.hh"

ClassImp(DetectorTransportManager)

using namespace TMath;
using namespace sou;
using namespace EConst;

//______________________________________________________________________________
 DetectorTransportManager::DetectorTransportManager() : DetectorPhotonTransporter(),
    fDebug(kFALSE), fScaleFactor(1.), fEnableFocalPlane(kTRUE) {
    //
    //  Constructor
    //

    // size of the cyl that includes detector and the baffle
    fDZup = 10000*mm;
    fDZdown = 7000*mm;
    fR = 2000*mm;
    fPos = EVector(0,0,0);
    fInnerRadius=Conf()->GetNum("DetectorTransportManager.fInnerRadius")*mm;
    
    of = OpticsFactory::Get();
    of->Build();
    
    fOptics     = of->GetOpticalSystem();
    fFocalPlane = of->GetFocalPlane();
    fIdealFocalSurface=of->GetIdealFocalSurface();
    fWalls      = of->GetWalls();
    fBaffle     = of->GetBaffle();

    if ( fFocalPlane->Radius() > fInnerRadius ||
         fOptics->Radius() > fInnerRadius ) 
        FatalError("FocalPlane radius or Optics radius is bigger than fInnerRadius. Aborting");
    else if ( fBaffle->Radius() < fInnerRadius ) 
        FatalError("Baffle cannot be smaller than fInnerRadius. Aborting");
    else
        fOuterRadius = fBaffle->Radius();

    fMaxIterations=(int)Conf()->GetNum("DetectorTransportManager.fMaxIterations");

    Reset();

}
 
//______________________________________________________________________________
 void DetectorTransportManager::SetScaleFactor( Double_t scale ) {
    //
    // Set the detector size scaling factor
    //
    
    if ( scale > 0 )
        fScaleFactor = scale;
}

//______________________________________________________________________________
 void DetectorTransportManager::Go( PhotonsOnPupil* pupil ) {
    //
    // Get and trace all the photons in PhotonsOnPupil
    //

    if ( !pupil ) { 
        MsgForm(EsafMsg::Warning,"PhotonsOnPupil is NULL. Photon tracing disabled.");
        return;
    }
    
    EEvent *ev = EEvent::GetCurrent();

    Photon *p=0;
    Photon *p_old=0;
    int nTracked=0, nTotal(pupil->GetNphotons());
    Int_t progress = 0;
    
    fTimeFirstPhoton = kHuge;
    fTimeLastPhoton = -kHuge;

    if (!fEnableFocalPlane) Msg(EsafMsg::Info) << "FocalPlane disabled" << MsgDispatch;

    while(1) {
        if ( p_old ) {

            //
            // Update the statistics
            //
            switch ( p_old->history ) {
                case kNotInside:
                    fNout++;
                    break;
                case kBaffle:
                    fNbaffle++;
                    break;
                case kOptics:
                    fNoptics++;
                    break;
                case kWalls:
                    fNwalls++;
                    break;
                case kFocalPlane:
                    fNfocalplane++;
                    break;
                case kEarth:
                    fNreflected++;
                    break;
                default:
                    fNlost++;
            }
            
            // 
            if( ev ) {
                EDetectorPhotonAdder a(p_old,0,kFALSE);  
                ev->Fill(a);
            }
        }

        // get a new photon
        p = pupil->GetPhoton();

        // build an empty EDetPhoton in the event
        if ( !p )  break;

        // scale ph position in respect of the pupil center
        // simulate greater or smaller experiments
        if ( Abs(fScaleFactor-1) > kTolerance )
            p->pos *= 1./fScaleFactor;

        if ( ev ) {
            EDetectorPhotonAdder a(p,0,kTRUE);  
            ev->Fill(a);
        }


        p_old=p;  // keep a copy for later use

        // keep track of first and last photon times
        if ( p->time > fTimeLastPhoton )
            fTimeLastPhoton = p->time;
        if ( p->time < fTimeFirstPhoton )
            fTimeFirstPhoton = p->time;

        // update display

        nTracked++;

        if (100.*nTracked/nTotal >= progress) {
            Msg(EsafMsg::Info).SetProgress(progress);
            Msg(EsafMsg::Info) << "Tracking Photons:  " << MsgCount;
            progress+=10;
        }

        Transport(p);

    }
    
    // fill a summary of photons history
    if ( ev ) {
        TArrayI photonHistory(7);
        photonHistory.AddAt(fNlost, kUndefined);
        photonHistory.AddAt(fNout, kNotInside);
        photonHistory.AddAt(fNbaffle, kBaffle);
        photonHistory.AddAt(fNoptics, kOptics);
        photonHistory.AddAt(fNfocalplane, kFocalPlane);
        photonHistory.AddAt(fNwalls, kWalls);
        photonHistory.AddAt(fNreflected, kEarth);
    
        EDetectorPhotonAdder a(photonHistory);
        ev->Fill(a);
    }
    
    // debug staistics
    Msg(EsafMsg::Debug) << "TransportManager Photon's transport statistics" << MsgDispatch;
    Msg(EsafMsg::Debug) << "NotInside    : " << fNout << MsgDispatch;
    Msg(EsafMsg::Debug) << "Baffle       : " << fNbaffle << MsgDispatch;
    Msg(EsafMsg::Debug) << "Optics       : " << fNoptics << MsgDispatch;
    Msg(EsafMsg::Debug) << "Walls        : " << fNwalls << MsgDispatch;
    Msg(EsafMsg::Debug) << "FocalPlane   : " << fNfocalplane << MsgDispatch;
    Msg(EsafMsg::Debug) << "Reflected    : " << fNreflected << MsgDispatch;
    Msg(EsafMsg::Debug) << "Lost         : " << fNlost << MsgDispatch;

    Msg(EsafMsg::Info) << "Tracking complete. " << nTracked <<" traced." << MsgDispatch;
}


//______________________________________________________________________________
 void DetectorTransportManager::Reset() {
    //
    // Reset the manager.
    //

    fNout        = 0;
    fNbaffle     = 0;
    fNoptics     = 0;
    fNwalls      = 0;
    fNfocalplane = 0;
    fNreflected  = 0;
    fNlost       = 0;
}

//______________________________________________________________________________
 Photon* DetectorTransportManager::Transport( Photon* p) const {
    //
    // Trace one photon in the detector.
    //

    // check if the photon has been already rejected by the pupil 
    if ( p->history == kNotInside ) return 0;
    
    
    EVector intPoint;
    
    if ( fDebug ) Msg(EsafMsg::Debug) << "@ new photon" << MsgDispatch;

    // check if the photon is in front of the detector
    // should be improved
    if ( p->pos.Perp()-fBaffle->Radius() > -kTolerance ) {
        p->history=kNotInside;
        return 0;
    }

    while( p ){
        if((p->iterations)++ > fMaxIterations) {
            MsgForm(EsafMsg::Warning,"Photon reached maximum number of iterations. killed.");
            Msg(EsafMsg::Warning) << " p->id         = " << p->parent->Id()  << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->pos        = " << p->pos << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->pos.Perp() = " << p->pos.Perp() << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->dir        = " << p->dir << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->history    = " << p->history << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->iterations = " << p->iterations << MsgDispatch;
            p = 0;
            continue;
        }

        if(p->fate != 0) {
            FatalError("Shouldn't happen p->fate != 0 while photon is still alive.");
        }
        // reality check on the photon position
        // FIXME these are debug conditions
        Bool_t isInEuso;
        isInEuso = p->pos[Z]-fBaffle->Bottom() > -kTolerance &&
            p->pos[Z]-fFocalPlane->Top() < kTolerance &&
            p->pos.Perp()-fInnerRadius < kTolerance;


        if ( !(isInEuso || fBaffle->IsInside( p )) ) {
            Msg(EsafMsg::Warning) << "This photon seems to be where it shouldn't:" << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->pos        = " << p->pos << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->pos.Perp() = " << p->pos.Perp() << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->dir        = " << p->dir << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->history    = " << p->history << MsgDispatch;
            Msg(EsafMsg::Warning) << " p->iterations = " << p->iterations << MsgDispatch;
            Msg(EsafMsg::Warning) << " isInEuso      = " << isInEuso << " isInsideBaffle = "
                                  << fBaffle->IsInside( p ) << MsgDispatch;
            p->history=kNotInside;
            p = 0;
            continue;
        }

        // photon before the optics
        if ( p->pos[Z] <= fOptics->FirstLensTop() ){

            if ( fDebug ) Msg(EsafMsg::Debug) << "--first zone" << MsgDispatch;

            if ( Abs(p->pos[Z]-fBaffle->Bottom()) < kTolerance && p->dir[Z] < 0) {

                if ( fDebug ) Msg(EsafMsg::Debug) << "  ~ back to the earth" << MsgDispatch;

                p->history=kEarth;
                p = 0;
                continue;

            } else 

                if ( p->pos.Perp() > fInnerRadius ) {
                    // photon outside the inner cylinder and inside the baffle

                    if ( fDebug ) Msg(EsafMsg::Debug) << " - OUTSIDE" << MsgDispatch;

                    // check if gets throught inner cylinder
                    Double_t DL = CylinderIntersection( p , fInnerRadius, fBaffle->Top(), fBaffle->Bottom(), kTRUE );

                    if ( DL > 0 ) {
                        // bring the photon to the borders of the inner cyl

                        if ( fDebug ) Msg(EsafMsg::Debug) << "  ~ going in" << MsgDispatch;

                        intPoint = p->pos+p->dir.Unit()*DL;
                        p->time += DL/Clight();
                        p->pos   = intPoint;

                    } else {

                        if ( fDebug ) Msg(EsafMsg::Debug) << "  ~ baffle" << MsgDispatch;

                        p->history = kBaffle;
                        p = fBaffle->Transport( p );         
                    }
                } else {//if ( (p->pos.Perp()-fInnerRadius) < kTolerance )

                    if ( Abs((p->pos-fPos).Perp()-fInnerRadius) < kTolerance){

                        if ( fDebug ) Msg(EsafMsg::Debug) << " - BORDERS ";

                        if ((p->pos-fPos).XYvector()*p->dir.XYvector() <0) {
                            // inward, nothing to do

                            if ( fDebug ) Msg(EsafMsg::Debug) << "inward" << MsgDispatch;

                        } else {
                            // going outward

                            if ( fDebug ) Msg(EsafMsg::Debug) << "outward" << MsgDispatch;

                            if (p->pos[Z] < fBaffle->Top() ){
                                p->history = kBaffle;

                                if ( fDebug ) Msg(EsafMsg::Debug) << "  ~ baffle" << MsgDispatch;

                                p = fBaffle->Transport( p );
                            } else {
                                p->history = kWalls;

                                if ( fDebug ) Msg(EsafMsg::Debug) << "  ~ walls" << MsgDispatch;

                                p = fWalls->Transport( p );
                            }
                            continue;
                        }
                    }
                    // photon inside the inner cylinder

                    if ( fDebug ) Msg(EsafMsg::Debug) << " - INSIDE" << MsgDispatch;

                    intPoint=InteractionPoint(p);

                    if ( Abs((p->pos-intPoint).Mag()) < kTolerance ){
                        //                    if (fLastCylInt.second[1] < 0)
                        Msg(EsafMsg::Warning) << "TraceOnePhoton: got confused in bringing ph in the inner cyl" << MsgDispatch;
                        // try with the second option
                        //                    intPoint = p->pos+p->dir.Unit()*fLastCylInt.second[1];
                    }

                    if ( intPoint[Z] >= fOptics->Bottom() ){

                        if ( fDebug ) Msg(EsafMsg::Debug) << "  ~ optics" << MsgDispatch;

                        Double_t DZ = fOptics->Bottom()-p->pos[Z];
                        intPoint=NewInteractionPoint(p,DZ);
                        p->time+=(intPoint-p->pos).Mag()/Clight();
                        p->pos=intPoint;
                        p->history=kOptics;
                        p=fOptics->Transport( p );
                    } else if ( intPoint[Z] < fBaffle->Bottom()) {  

                        if ( fDebug ) Msg(EsafMsg::Debug) << "  ~ back to the earth" << MsgDispatch;

                        Double_t DL = fBaffle->CylinderIntersection( p );
                        p->time+=DL/Clight();
                        p->pos+=p->dir.Unit()*DL;
                        p->history=kEarth;
                        p=0;

                    } else {

                        if ( fDebug ) Msg(EsafMsg::Debug) << "  ~out of the inner cyl" << MsgDispatch;

                        p->time+=(intPoint-p->pos).Mag()/Clight();
                        p->pos=intPoint;
                    }
                }

        } else if (p->pos[Z] > Max(fOptics->FirstLensTop(), fBaffle->Top())
                && p->pos[Z] < fOptics->SecondLensBottom()) {
            Msg(EsafMsg::Debug) << "--optics forbidden zone" << MsgDispatch;
            Msg(EsafMsg::Debug) << "- p->pos="<< p->pos << "   - p->dir=" << p->dir << MsgDispatch;
            Msg(EsafMsg::Debug) << "- p->pos.Perp()=" << p->pos.Perp() << MsgDispatch;
            p->history=kNotInside;
            p=0;
        } else {

            // calculate possible impact point of the photon on the cylinder

            if ( fDebug ) Msg(EsafMsg::Debug) <<  "--second zone" << MsgDispatch;

            intPoint=InteractionPoint(p);
            // photon is in the second zone:
            // between the optics and the focal plane

            if(intPoint[Z] > fFocalPlane->Bottom()) {
                // photon on focal plane

                if ( fDebug ) Msg(EsafMsg::Debug) << "  ~focalplane" << MsgDispatch;

                fIdealFocalSurface->HitPosition(p);
                p->history=kFocalPlane;

                if (!fEnableFocalPlane) {
                    p=0;
                    continue;
                }

                // DZ distance along Z from photon to focal plane
                double DZ = fFocalPlane->Bottom() - p->pos[Z];

                // calculate impact point on the focal plane
                intPoint=NewInteractionPoint(p, DZ);

                // calculate interaction of photon with focal plane 
                p->time+=(p->pos - intPoint).Mag()/Clight();
                p->pos=intPoint;
                p=fFocalPlane->Transport(p);
            } else if(intPoint[Z] < fOptics->Top() && p->dir[Z] < 0) {
                // photon on optics

                if ( fDebug ) Msg(EsafMsg::Debug) << "  ~optics from above" << MsgDispatch;

                if ( p->pos[Z] > fOptics->Top() ) { 
                    // DZ distance along Z from photon to optics
                    double DZ = fOptics->Top() - p->pos[Z];

                    // calculate impact point on the optics surface
                    intPoint=NewInteractionPoint(p, DZ);
                    // calculate interaction of photon with optics
                    p->time+=(p->pos - intPoint).Mag()/Clight();
                    p->pos=intPoint;
                    p->history=kOptics;
                    p=fOptics->Transport(p);
                } else if (p->history == kOptics ){
                    p->history=kWalls;
                    p=fWalls->Transport(p);
                } else {
                    p->history=kOptics;
                    p=fOptics->Transport(p);
                }

            } else {       
                // photon on wall: reflection or absorption

                if ( fDebug ) Msg(EsafMsg::Debug) << "  ~walls" << MsgDispatch;

                p->time+=(p->pos - intPoint).Mag()/Clight();
                p->pos=intPoint;
                p->history=kWalls;
                p=fWalls->Transport(p);
            }
        }
    }

    return 0;

}


//______________________________________________________________________________
 Double_t DetectorTransportManager::EusoCylIntersection( Photon *p ) const {
    //
    // Returns the distance (if any) of p from the next intersection with the
    // EUSO shell.
    //
    
    return CylinderIntersection( p , fInnerRadius, fPos[Z]+fDZup, fPos[Z]-fDZdown);
}

//______________________________________________________________________________
 EVector DetectorTransportManager::InteractionPoint(Photon *p, Double_t radius) const {
    //
    // Returns the next interaction point of p with a cylinder with a given
    // radius. The cylinder axis is the optical axis.
    //
    
    // intesection point
    EVector intPoint(p->pos);
    if ( radius < 0 ) radius = fInnerRadius;

    Double_t DL;
    DL = CylinderIntersection( p , radius, fPos[Z]+fDZup, fPos[Z]-fDZdown, kTRUE);
    if ( DL < 0 ){
        Msg(EsafMsg::Warning) << "InteractionPoint" << MsgDispatch;
//        MsgForm(EsafMsg::Warning,"DL = %f  - LastCyInt[0] = %f  - LastCylInt[1] = %f",
//                DL, fLastCylInt.second[0], fLastCylInt.second[1]); 
        Msg(EsafMsg::Warning) << "p->pos: " << p->pos << " - p->pos.Perp(): " << p->pos.Perp() << MsgDispatch;
        MsgForm(EsafMsg::Warning,"R: - %f  pos*dir: %f  - pos*dir:",
                radius, p->pos.XYvector()*p->dir.XYvector());
        p->Dump();
        FatalError("should't happen DL<0. Photon killed");
    }
    intPoint=p->pos+p->dir.Unit()*DL;
    return intPoint;
}

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