Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

BunchRadiativeTransfer - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: BunchRadiativeTransfer.cc,v 1.93 2005/11/14 11:24:52 moreggia Exp $
// Sylvain Moreggia created Jan, 15 2004

#include "BunchRadiativeTransfer.hh"
#include "Atmosphere.hh"  
#include "RadiativeFactory.hh"
#include "ListPhotonsInAtmosphere.hh"
#include "ListPhotonsOnPupil.hh"
#include "EsafRandom.hh"
#include <TROOT.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TProfile.h>
#include <math.h>  
#include "ParentPhoton.hh" 
#include "SinglePhoton.hh" 
#include "Config.hh"
#include "LowtranRadiativeProcessesCalculator.hh"
#include "EConst.hh"
#include "BunchOfPhotons.hh"
#include "Ground.hh"

#include "EEvent.hh"
#include "EAtmosphere.hh"
#include "EAtmosphereBunchAdder.hh"
#include "EAtmosphereSingleAdder.hh"

#include "TBenchmark.h"

ClassImp(BunchRadiativeTransfer)

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

//_______________________________________________________________________________________________________________________
 BunchRadiativeTransfer::BunchRadiativeTransfer() : RadiativeTransfer(), fPhotons(0), fTotalList(0),
                                                   fCSpropag(0), fICpropag(0), fNbBunch(0), fNbTot(0) {
    //
    // ctor
    //
    Msg(EsafMsg::Info) << "Enabled" << MsgDispatch;
    fGround = RadiativeFactory::Get()->GetGround();
    if(!fGround) Msg(EsafMsg::Panic) << "Pb of memory allocation for fGround" << MsgDispatch;
    
    fCSpropag = RadiativeFactory::Get()->GetClearSkyPropagator(fGround);
    fICpropag = RadiativeFactory::Get()->GetInCloudsPropagator(fGround);
    if(!fCSpropag) Msg(EsafMsg::Panic) << "Problem of memory allocation for CSpropag" << MsgDispatch;
    if(!fICpropag) Msg(EsafMsg::Panic) << "Problem of memory allocation for ICpropag" << MsgDispatch;
    fLowToDetec = new LowtranRadiativeProcessesCalculator();
    
    // detector geometry : decoupled or optimized
    // decoupled : detector seems to be a sphere. Allows to study several detector geometries with the same RadiativeTransfer simulation
    // optimized : true detector geometry is used in RadiativeTransfer simulation
    string name = Conf()->GetStr("BunchRadiativeTransfer.fDecoupled");
    if(name == "decoupled") fDecoupled = true;
    else if(name == "optimized") fDecoupled = false;
    else Msg(EsafMsg::Panic) << "Wrong option for BunchRadiativeTransfer.fDecoupled" << MsgDispatch;
    
}

//_______________________________________________________________________________________________________________________
 BunchRadiativeTransfer::~BunchRadiativeTransfer() {
    //
    // dtor
    //    
    SafeDelete(fPhotons);
    SafeDelete(fCSpropag);
    SafeDelete(fICpropag);
    SafeDelete(fLowToDetec);
}

//_______________________________________________________________________________________________________________________
 PhotonsOnPupil* BunchRadiativeTransfer::Get(PhotonsInAtmosphere* photons, const DetectorGeometry* dg) {
    //
    // Transport photons from source to Euso pupil
    //


    // set DetectorGeometry for all propagators
    CopyDetectorGeometry(dg);
    
    
    // in case ground detector simulation    //GRNDetec
    if(EUSO().Zv() == 0) fDetAtGrnd = true;
    else fDetAtGrnd = false;
    if(fDetAtGrnd && !fDecoupled) Msg(EsafMsg::Panic) << "In Ground detector mode, decoupled mode must be switched on" << MsgDispatch;


    // if given list of photons is empty  
    if(!photons) return NULL;
    

#ifdef DEBUG
    TString jobRT = "Time spent for transfer process:";
    TBenchmark gB;
    gB.Start(jobRT);
#endif

    if(photons->GetType() == "empty") Msg(EsafMsg::Panic) <<"When NoLightSource used, NoRadiativeTransfer must be used"<< MsgDispatch;
    if(photons->GetType() != "list") Msg(EsafMsg::Panic) <<"Wrong PhotonsInAtmosphere format. ListPhotonsInAtmosphere expected."<< MsgDispatch;
    fTotalList = (ListPhotonsInAtmosphere*) photons;
    fCSpropag->SetEndFoV(fTotalList->GetTrackStep(fTotalList->GetNbTrackSteps()-1)); //TOFIX when general FoV handling
    fICpropag->SetEndFoV(fTotalList->GetTrackStep(fTotalList->GetNbTrackSteps()-1)); //TOFIX when general FoV handling



    // Build a datacard of VERTICAL Lowtran transmission for PropagateToDetector() method
    if(!fDecoupled || fDetAtGrnd) {
        Msg(EsafMsg::Info) <<"Building LOWTRAN vertical transmission datacard"<< MsgDispatch;
        fLowToDetec->MakeVerticalDatacard(fGround,EUSO().Zv());
        Msg(EsafMsg::Info) <<"LOWTRAN datacard built"<< MsgDispatch;
    }
    // Build a datacard of ALONG TRACK Lowtran transmission for PropagateToDetector() method
    else if(fDecoupled) {
        Msg(EsafMsg::Info) <<"Building LOWTRAN datacard ALONG TRACK (not vertical)"<< MsgDispatch;
        fLowToDetec->MakeTrackToDetectorDatacard(fGround,EUSO(),*fTotalList);
	Msg(EsafMsg::Info) <<"LOWTRAN datacard built"<< MsgDispatch;
    }


    // number of bunches
    fNbBunch = fTotalList->GetListOfBunch().size();
    fNbTot = 0;
    Double_t nf = 0;
    Double_t nc = 0;
    for(size_t i=0; i<fNbBunch; i++) {
        fNbTot += (fTotalList->GetListOfBunch()[i])->GetWeight();
        if(fTotalList->GetListOfBunch()[i]->GetType() == Fluo) nf += fTotalList->GetListOfBunch()[i]->GetWeight();
        else nc += fTotalList->GetListOfBunch()[i]->GetWeight();
    }
#ifdef DEBUG
    Msg(EsafMsg::Debug) << "SIZE OF LIST OF BUNCHES = "<<fNbBunch << MsgDispatch;
    Msg(EsafMsg::Debug) << "Nb fluo = " <<nf << MsgDispatch;
    Msg(EsafMsg::Debug) << "Nb ckov = " <<nc << MsgDispatch;
#endif
    Msg(EsafMsg::Info) << "TOTAL number of photons ="<<fNbTot << MsgDispatch;

    if ( fNbBunch >=1 ) {
        EarthVector track = fTotalList->GetListOfBunch()[fNbBunch - 1]->GetPos() 
            - fTotalList->GetListOfBunch()[0]->GetPos();
#ifdef DEBUG
        Msg(EsafMsg::Debug) << "size of the photon track = "<<track.Mag()/km <<" km"<< MsgDispatch;
#endif
    }


    EEvent* ev = EEvent::GetCurrent();
    ev->GetAtmosphere()->SetMaxScatOrder(1);


    // if some SinglePhotons created directly by LightSource module, they are saved into root
    if ( ev ) {
        const vector<SinglePhoton*>& list_single_2 = fTotalList->GetListOfSingle();
        size_t nbnotsaved = fTotalList->GetNbNotSaved();
        for (size_t i = list_single_2.size() - nbnotsaved; i<list_single_2.size(); i++ ) {
            EAtmosphereSingleAdder sa( list_single_2[i],true,false,0 );
            ev->Fill(sa);
            fTotalList->OneSingleSaved();
        }
    }


    // loop to transport all the bunches
    BunchOfPhotons* bunch = 0;
    Double_t nbTracked(0);
    Int_t progress = 0;

    while(true) {
        bunch = fTotalList->GetBunch();
        if(!bunch) break;
        nbTracked++;

        // saves bunch parameters at creation
        if ( ev ) {
            EAtmosphereBunchAdder ba( bunch, true );
            ev->Fill(ba);
        }

        // propagation of bunches and creation of single photons
        BunchPropagation(*bunch);
#ifdef DEBUG
        Msg(EsafMsg::Debug) <<"after the process of this bunch, size of list of single = " <<fTotalList->GetListOfSingle().size() <<MsgDispatch;
#endif

        // saves single photon parameters at creation
        if ( ev ) {
            const vector<SinglePhoton*>& list_single = fTotalList->GetListOfSingle();
            size_t nbnotsaved = fTotalList->GetNbNotSaved();
            for (size_t i = list_single.size() - nbnotsaved; i<list_single.size(); i++ ) {
                EAtmosphereSingleAdder sa( list_single[i],true,true );
                ev->Fill(sa);
                fTotalList->OneSingleSaved();
            }
        }

        if (100.*nbTracked/fNbBunch >= progress) {
            Msg(EsafMsg::Info).SetProgress(progress);
            Msg(EsafMsg::Info) << "Bunches Processing:" << MsgCount;
            progress+=10;
        }
    }

#ifdef DEBUG
    Msg(EsafMsg::Debug) <<"final list of single, size = " <<fTotalList->GetListOfSingle().size() <<MsgDispatch;
#endif
    if(fTotalList->GetNbNotSaved()) Msg(EsafMsg::Warning) <<fTotalList->GetNbNotSaved()<<" SinglePhoton not saved into root" <<MsgDispatch;

    // propagates single photons
    PropagationOfSingles();

    // saves single photons after propagation in atmosphere
    if ( ev ) {
        const vector<SinglePhoton*>& list_single_3 = fTotalList->GetListOfSingle();
        for ( size_t i=0; i<list_single_3.size(); i++ ) {
            EAtmosphereSingleAdder sa( list_single_3[i],false,false,i );
            ev->Fill(sa);
        }
    }


#ifdef DEBUG
    gB.Stop(jobRT);
    MsgForm(EsafMsg::Debug,"Time spent for transfer process:  REAL=%6.2f s   CPU=%6.2f s",gB.GetRealTime(jobRT),gB.GetCpuTime(jobRT));
#endif

    return fPhotons;
}

//_______________________________________________________________________________________________________________________
 void BunchRadiativeTransfer::Reset() {
    //
    // get ready for next event
    // NB : fLowToDetec not reset, its datacard is used for all the events of a run
    //
    if(fGround)     fGround->Reset();
    if(fPhotons)    fPhotons->Clear();
    if(fCSpropag)   fCSpropag->Reset();
    if(fICpropag)   fICpropag->Reset();
    if(fLowToDetec) fLowToDetec->Reset();
    fNbBunch = 0;
    fNbTot = 0;
}

//_______________________________________________________________________________________________________________________
 void BunchRadiativeTransfer::DirectToEuso(const BunchOfPhotons& bigbunch) const {
    //
    // Creates a list of SinglePhoton considering the EUSO solid angle
    // handles fluo and cerenkov (angular distrib used for the later)
    //
    
    Int_t nb = 0;
    EarthVector towardEUSO = (EUSO() - bigbunch.GetPos()).Unit();
    Double_t theta = fabs(towardEUSO.Angle(bigbunch.GetDir())); // to keep theta within 0-Pi()
    if(theta > Pi()) Msg(EsafMsg::Warning) << "<DirectToEuso> Pb with angle definition" << MsgDispatch;
    
    Double_t angular_distrib_value = bigbunch.AngularDist_OverTwoPi(theta);
    nb = EsafRandom::Get()->Poisson(bigbunch.GetWeight() * EusoOmega(bigbunch.GetPos()) * angular_distrib_value);
    GenerateDirectSingles(bigbunch,nb);

#ifdef DEBUG
    Msg(EsafMsg::Debug)  <<"  Direct gives -> " <<nb <<" photons" << MsgDispatch;
#endif
}

//_______________________________________________________________________________________________________________________
 void BunchRadiativeTransfer::BunchPropagation( BunchOfPhotons& bigbunch ) {
    //
    // Propagation of a bunch through the atmosphere
    //
    // temp remark : bunch splitting not used so far, thus don't handle for root filling //TOFIX
    //

#ifdef DEBUG
    Msg(EsafMsg::Debug)  << "\nPropagation of the BUNCH nb " << bigbunch.GetId() <<MsgDispatch;
    Msg(EsafMsg::Debug)  << "Weight = " << bigbunch.GetWeight() <<MsgDispatch;
    const ParentBunch* parentb = bigbunch.GetParent();
    Double_t length = (parentb->GetShowerPosi() - parentb->GetShowerPosf()).Mag()/km;
    Msg(EsafMsg::Debug)  << "Bunch longitudinal extension = "<<length <<" km" << MsgDispatch;
#endif    

    // if bunch is underground, its simulation stops here
    if(fGround->IsUnderGround(bigbunch.GetParent()->GetShowerPosf()))   //TOFIX : a part can be above ground (also see related point in ListPinAtmo::fTrack, O2_CSPropag)
        bigbunch.SetFate(3);
    else {     
	// photons directly emitted within the Euso solid angle
        DirectToEuso(bigbunch);
        Medium state = CLEARSKY;

        // if bunch weight < 1% "mean bunch weight" -> not propagated, because it doesn't contribute
        if( (bigbunch.GetWeight() < 0.001*fNbTot/fNbBunch)) {
            bigbunch.SetFate(2);
            state = NONE;
        }

        // relevant propagator called
        while(state > NONE) {
            switch(state) {

                case CLEARSKY :       state = fCSpropag->Go(bigbunch,*fTotalList);
                                      break;

                case CLOUDY   :       state = fICpropag->Go(bigbunch,*fTotalList);
                                      break;

                case GROUND   :       GroundReflection(bigbunch);
                                      state = NONE;
                                      bigbunch.SetFate(4);
                                      break;

                case AEROSOLS :       Msg(EsafMsg::Panic) << "AerosolsPropagator (for MS) not implemented" << MsgDispatch;
                                      break;

                default       :       Msg(EsafMsg::Panic) << "Invalid type of medium in bunch propagation" << MsgDispatch;
            }
        }
    }
#ifdef DEBUG
    switch(bigbunch.GetFate()) {
        case 1    : Msg(EsafMsg::Debug)  << "NOT PROPAGATED -> FLUO" << MsgDispatch; break;
        case 2    : Msg(EsafMsg::Debug)  << "NOT PROPAGATED -> TOO SMALL WEIGHT" << MsgDispatch; break;
        case 3    : Msg(EsafMsg::Debug)  << "NOT PROPAGATED -> CREATED UNDERGROUND" << MsgDispatch; break;
        case 4    : Msg(EsafMsg::Debug)  << "HAS REACHED GROUND" << MsgDispatch; break;
        case 5    : Msg(EsafMsg::Debug)  << "GONE OUT Euso FoV" << MsgDispatch; break;
        default   : Msg(EsafMsg::Debug)  << "PB with STATE" << MsgDispatch;
    }
#endif


    // saves propagated bunch parameters  //TOFIX : if one day splitting is relevant
    EEvent* ev = EEvent::GetCurrent();
    if ( ev ) {
        EAtmosphereBunchAdder ba(&bigbunch,false);
        ev->Fill(ba);
    }
}

//_____________________________________________________________________________________________________
 void BunchRadiativeTransfer::GroundReflection(const BunchOfPhotons& b) const {
    //
    // Process of a bunch reaching ground
    //
    
    // initializations
#ifdef DEBUG
    Msg(EsafMsg::Debug) << "GroundReflection" << MsgDispatch; 
#endif
    if(fabs(b.GetPos().Zv() - fGround->Altitude(b.GetPos())) > 1)
        Msg(EsafMsg::Warning) << "GroundReflection() must be called if bunch has reached the ground" << MsgDispatch;
    
    // determination of number of singlephotons produced
    Int_t nb = EsafRandom::Get()->Poisson(b.GetWeight() * fGround->Albedo(b.GetPos()) * EusoOmega(b.GetPos()) * fGround->Outgoing_phase_function(b.GetPos(),EUSO()));
    
    // creation of singlephotons
    GenerateReflectedSingles(nb,b);
}
    
//_______________________________________________________________________________________________________________________
 void BunchRadiativeTransfer::PropagationOfSingles() {
    //
    // Final phase of transfer. Propagate all the SinglePhoton til EUSO pupil
    //

    // initializations
    if(!fPhotons) fPhotons = new ListPhotonsOnPupil((vector<ParentPhoton*>*)NULL);
    if(!fPhotons) Msg(EsafMsg::Panic) << "BunchRadiativeTransfer::PropagationOfSingles, NULL fPhotons : Memory pb" << MsgDispatch;

    // build PhotonsOnPupil's frame
    BuildPupilFrame(fPhotons);
    
    EarthVector dirtest(1);
    Double_t Trans[4];
    Double_t TotTrans(0.);
    TRandom* rndm = EsafRandom::Get();
    SinglePhoton* p = 0;
    size_t nTotal(0),nTracked(0);
    Int_t progress = 0;
    TVector3 local_dir, pos;

    nTotal = fTotalList->GetSingleEntries();

    // loop over the list of SinglePhoton
    while(true) {
        p = fTotalList->GetSingle();
        if(!p) break;
        nTracked++;
        dirtest = (EUSO() - p->Pos()).Unit();
        if((dirtest - p->Dir()).Mag() > TOLERANCE) Msg(EsafMsg::Warning) << "PropagationOfSingles : SinglePhoton must be directed toward EUSO" << MsgDispatch;

        // calculate transmission from photon position to EUSO
        TotTrans = fLowToDetec->Trans(*p,EUSO(),Trans);
	p->SetLastTrans(TotTrans,"tot");
	p->SetLastTrans(Trans[1],"rayl");
	p->SetLastTrans(Trans[2],"ozone");
	p->SetLastTrans(Trans[3],"aero");
	p->SetLastTrans(TotTrans/Trans[0],"cloud");

        // photon is transmitted or absorbed
        if(TotTrans < rndm->Rndm())  p->SetAbsorbed();
	
        // Sample a random photon position on pupil
	pos =  p->Pos();
	local_dir = p->Dir();
        p->AddToPosTof(EUSO() - p->Pos());
        RamdomPosOnPupil(fPhotons,pos,local_dir);
	
	// check if photon is within the FoV
	if(!GetDetGeometry()->IsInFoV(local_dir)) p->SetOutFoV();
	else p->SetOutFoV(false);
	
	// if photon transmitted, becomes a photon on pupil
	// FoV 'status' not relevant here (too simply treated, ONLY a flag in RT part)
	// -->  will be considered in detector part
        if(!p->IsAbsorbed()) fPhotons->Add(*p,pos,local_dir);

        //GRNDetec
        // for ground detector simulation, cos(theta) effect is taken into account here
	if(fDetAtGrnd) {
	    EarthVector enterdir = -dirtest;
            if((rndm->Rndm() > cos(enterdir.Theta())) || (enterdir.Theta() > PiOver2())) p->SetAbsorbed();
	}

        // Dump CPU commentaries
        if (100.*nTracked/nTotal >= progress) {
            Msg(EsafMsg::Info).SetProgress(progress);
            Msg(EsafMsg::Info) << "Singles Processing:" << MsgCount;
            progress+=10;
        }
    }

#ifdef DEBUG
    Msg(EsafMsg::Debug) << "Number of ParentPhoton = " << fPhotons->GetNphotons() << MsgDispatch;
#endif
}

//_______________________________________________________________________________________________________________________
 void BunchRadiativeTransfer::GenerateDirectSingles(const BunchOfPhotons& b,Int_t nb) const {
    //
    // SinglePhoton generation (BunchOfPhotons portion which is created within EUSO solid angle)
    // SinglePhotons created added to fTotalList
    //
    Double_t wl, date, tof;
    EarthVector showerpos, dir, diff;
    SinglePhoton* s = 0;
    
    UInt_t bid = b.GetId();
    PhotonType type = b.GetType();
    tof = 0;
    
    for(Int_t i=0; i<nb; i++) {
        // corrections for date and showerpos from BunchOfPhotons mean values and longitudinal dispersion
        showerpos = b.RandomPosInShower();
	if(fGround->IsUnderGround(showerpos)) continue;
	diff = showerpos - b.GetShowerPos();
        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
        dir = EUSO() - showerpos;
	wl = b.GetWlSpectrum().GetLambda();
	s = new SinglePhoton(type,date,tof,wl,showerpos,showerpos,dir,Direct,bid);
	fTotalList->Add(s);
    }
}

//_____________________________________________________________________________________________________
 void BunchRadiativeTransfer::GenerateReflectedSingles(Int_t nb, const BunchOfPhotons& b) const {
    //
    // Creates SinglePhoton objects coming from a BunchOfPhotons reflection on ground
    //    
#ifdef DEBUG
    Msg(EsafMsg::Debug)<<"nb of single produced in reflexion process = "<<nb <<MsgDispatch;
#endif
    Double_t wl, date, tof;
    EarthVector showerpos, pos, dir, diff;
    SinglePhoton* s = 0;
    
    UInt_t bid = b.GetId();
    PhotonType type = b.GetType();
    date = b.GetDate();
    
    for(Int_t i=0; i<nb; i++) {
        // corrections for date and showerpos, from BunchOfPhotons mean values and longitudinal+lateral distributions at creation
        showerpos = b.RandomPosInShower();
	
	// if showerpos is underground
	if(fGround->IsUnderGround(showerpos)) continue;
	
	diff = showerpos - b.GetShowerPos();
	date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
	// tof and position corrections, due to angular distributions at creation
	pos = b.PosRandomAngCorrec(showerpos,b.GetPos());   // here is a fake correction in position, used to get direction for new impact calculation
	pos = fGround->GetImpact(showerpos,(pos - showerpos).Unit());
	// if no impact
	if(pos.Z() == HUGE) continue;
	tof = b.GetTof() * (showerpos - pos).Mag()/(b.GetPos() - b.GetShowerPos()).Mag();
        dir = EUSO() - pos;
	wl = b.GetWlSpectrum().GetLambda();
	s = new SinglePhoton(type,date,tof,wl,showerpos,pos,dir,Reflected,bid);
	s->AddInteraction();
	s->AddHistory(Reflected);
	fTotalList->Add(s);
    } 
}

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