Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

MCRadiativeTransfer - source file

// $Id: MCRadiativeTransfer.cc,v 1.25 2005/11/14 11:24:52 moreggia Exp $
// Author: Sylvain Moreggia   2005/08/16

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: MCRadiativeTransfer                                                  *
 *  Package: radiativetransfer                                               *
 *  Coordinator: Sylvain Moreggia                                            *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// MCRadiativeTransfer
//
// Detection from a ground detector no possible in this mode
//
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "MCRadiativeTransfer.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(MCRadiativeTransfer)

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

//_____________________________________________________________________________
 MCRadiativeTransfer::MCRadiativeTransfer() : RadiativeTransfer() , fPhotons(0), fTotalList(0), fPropagator(0), fLowToDetec(0),
                                             fNbBunch(0), fNbTot(0), fNbSingles(0), fOmegaMax(0) {
    //
    // ctor
    //
    Msg(EsafMsg::Info) << "Enabled" << MsgDispatch;
    fGround = RadiativeFactory::Get()->GetGround();
    if(!fGround) Msg(EsafMsg::Panic) << "Pb of memory allocation for fGround" << MsgDispatch;
    
    fMaxPhaseFunction = Conf()->GetNum("MCRadiativeTransfer.fMaxPhaseFunc");
    Double_t maxtof = Conf()->GetNum("MCRadiativeTransfer.fTofCut")*microsecond;
    fPropagator = new SinglePhotonPropagator(fGround,fMaxPhaseFunction,maxtof);
    fLowToDetec = new LowtranRadiativeProcessesCalculator();
    
    fScatOrder = Int_t(Conf()->GetNum("MCRadiativeTransfer.fScatOrder"));
    Msg(EsafMsg::Info) << "Max Scattering order to be simulated = " << fScatOrder << MsgDispatch;
    
    // 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("MCRadiativeTransfer.fDecoupled");
    if(name == "decoupled") fDecoupled = true;
    else if(name == "optimized") fDecoupled = false;
    else Msg(EsafMsg::Panic) << "Wrong option for MCRadiativeTransfer.fDecoupled" << MsgDispatch;
}

//_____________________________________________________________________________
 MCRadiativeTransfer::~MCRadiativeTransfer() {
    //
    // Destructor
    //
    SafeDelete(fPhotons);
    SafeDelete(fPropagator);
    SafeDelete(fLowToDetec);
}

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



    // copy DetectorGeometry
    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(fDetAtGrnd) Msg(EsafMsg::Panic) << "Ground detector mode not possible with MCRadiativeTransfer" << MsgDispatch; //TOFIX
    
    // if given list of photons is empty  
    if(!photons) return NULL;

    if(photons->GetType() == "empty") Msg(EsafMsg::Panic) <<"When NoLightSource used, NoRadiativeTransfer should be used"<< MsgDispatch;
    if(photons->GetType() != "list") Msg(EsafMsg::Panic) <<"Wrong PhotonsInAtmosphere format. ListPhotonsInAtmosphere expected."<< MsgDispatch;
    fTotalList = (ListPhotonsInAtmosphere*) photons;




    // 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;
    }
    else Msg(EsafMsg::Info) <<"NO LOWTRAN datacard in decoupled mode"<< MsgDispatch;





    // find position of lighttrack max
    // number of bunches
    fNbBunch = fTotalList->GetListOfBunch().size();
    fNbTot = 0;
    fNbSingles = 0;
    Double_t nf(0), nc(0), maxweight(0);
    EarthVector showermaxpos(0,0,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();
	    if((fTotalList->GetListOfBunch()[i])->GetWeight() > maxweight) {
	        maxweight = (fTotalList->GetListOfBunch()[i])->GetWeight();
		showermaxpos = (fTotalList->GetListOfBunch()[i])->GetShowerPos();
	    }
	}
        else nc += fTotalList->GetListOfBunch()[i]->GetWeight();
    }
    
    // set omega max
    SetOmegaMax(showermaxpos);
    fPropagator->SetOmegaMax(fOmegaMax);

// some printing
#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
    }
    
    



    // root output init
    EEvent* ev = EEvent::GetCurrent();
    ev->GetAtmosphere()->SetMaxScatOrder(fScatOrder);

    // 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();
        }
    }





    // ALGO STARTS HERE
    BunchOfPhotons* bunch = 0;
    Bool_t next(0),subnext(0);
    Double_t nbTracked(0), fraction(0), left(0),right(10),subleft(0),subright(2.5);

    // 1. Direct photons simulation
    while(true) {
        bunch = fTotalList->GetBunch();
        if(!bunch) break;

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

        // photons going directly toward detector
        DirectToEuso(*bunch);
    }

    // 2. SCATTERING SIMULATION loop
    for(Int_t ord=1; ord <= fScatOrder; ord++) {
        Msg(EsafMsg::Info) << "nSCATTERING ORDER = " << ord<< MsgDispatch;
#ifdef DEBUG
        TString jobMCRT = "Time spent for Scattering simu:";
        TBenchmark gB;
        gB.Start(jobMCRT);
#endif
        nbTracked = 0;
        fraction = 0;
        left = 0;
        right = 10;
        subleft = 0;
        subright = 2.5;
        Msg(EsafMsg::Info) << "Bunches Transformation: 0%" << MsgFlush;
        fTotalList->ResetBunchCounter();  // for below extraction from list
        while(true) {
            bunch = fTotalList->GetBunch();
            if(!bunch) break;
            nbTracked++;

            // 2.1 generate reduce nb of photons
            PhotonsInOmega(*bunch);

            // dump info
            fraction = 100*nbTracked/fNbBunch;
            if (left<fraction && fraction <=right) 
                next = kFALSE;
            else if (fraction>right) {
                next = kTRUE;
                left = right;
                right += 10;
            }
            if (subleft<fraction  && fraction <=subright) {
                subnext = kFALSE;
            }
            if (fraction > subright) {
                subnext = kTRUE;
                subleft = subright;
                subright +=2;
            }
            if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
            if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;	 
        }
        Msg(EsafMsg::Info) << " -done" << MsgDispatch;
        Msg(EsafMsg::Info) << "NB of SinglePhoton to PROPAGATE = "<<fNbSingles<< MsgDispatch;
        fNbSingles = 0;  // reset


        // 2.2 scattering simu for this bunch at the relevant scattering order
        fPropagator->Go(fTempList,ord);

        // 2.3 Move photons from fTempList to fTotalList : fTempList does not contain photons anymore
        //     Only photons which have scattered are written in fTotalList, others are lost
        fTotalList->Add(fTempList);

#ifdef DEBUG
        gB.Stop(jobMCRT);
        MsgForm(EsafMsg::Info,"Time spent for this scattering order :  REAL=%6.2f s   CPU=%6.2f s",gB.GetRealTime(jobMCRT),gB.GetCpuTime(jobMCRT));
#endif
    }





    // 3. saves single photon parameters after scattering simulation
    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();
        }
    }

#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;




    // 4. propagates direct and scattered photons until detector (Lowtran is used, c.f. Ozone)
    PropagateToDetector();
    
    
    

    // saves single photons after propagation until detector
    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);
        }
    }

    return fPhotons;
}

//_______________________________________________________________________________________________________________________
 void MCRadiativeTransfer::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(fPropagator) fPropagator->Reset();
    for(size_t i=0; i<fTempList.size(); i++) {
        SafeDelete(fTempList[i]);
        Msg(EsafMsg::Panic) <<"fTempList SHOULD BE EMPTY, remains : "<<fTempList.size()<<" photons inside"<<MsgDispatch;
	
    }
    if(fLowToDetec) fLowToDetec->Reset();
    fTempList.clear();
    fNbBunch = 0;
    fNbTot = 0;
    fNbSingles = 0;
    fOmegaMax = 0;
}

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

    // SinglePhotons created and added to fTotalList (THEY WILL NOT UNDERGO SCATTERING SIMU)
    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);
    }

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

//_______________________________________________________________________________________________________________________
 void MCRadiativeTransfer::PhotonsInOmega(const BunchOfPhotons& b) {
    //
    // Detector solid angle is applied here -> the resulting photons will undergo scattering simu
    // Photon features are set here (lambda, position, direction)
    // SinglePhoton are stored into fTempList for scattering simu
    //

    Int_t nb = EsafRandom::Get()->Poisson(b.GetWeight() * fOmegaMax * fMaxPhaseFunction);
    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++) {
        // get features at creation
        showerpos = b.RandomPosInShower();
        if(fGround->IsUnderGround(showerpos)) continue;
        diff = showerpos - b.GetShowerPos();
        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
        dir  = b.RandomDirection();
        wl   = b.GetWlSpectrum().GetLambda();
        s    = new SinglePhoton(type,date,tof,wl,showerpos,showerpos,dir,None,bid);
        fNbSingles++;
        fTempList.push_back(s);
    }
#ifdef DEBUG
    //Msg(EsafMsg::Debug)  <<"  InOmega gives -> " <<nb <<" photons" << MsgDispatch;
#endif
}

//_______________________________________________________________________________________________________________________
 void MCRadiativeTransfer::PropagateToDetector() {
    //
    // Final phase of transfer. Propagate all the SinglePhoton til EUSO pupil
    // USE LOWTRAN (because of Ozone transmission)
    //

    // initializations
    if(!fPhotons) fPhotons = new ListPhotonsOnPupil((vector<ParentPhoton*>*)NULL);
    if(!fPhotons) Msg(EsafMsg::Panic) << "MCRadiativeTransfer::PropagateToDetector, 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);
    TVector3 local_dir, pos;
    Float_t fraction(0); 
    Bool_t next(0),subnext(0);
    Float_t left(0),right(10),subleft(0),subright(2.5);
    nTotal = fTotalList->GetSingleEntries();
    Msg(EsafMsg::Info) << "Propagation to Detector: 0%" << MsgFlush;

    // loop over the list of SinglePhoton
    while(true) {
        p = fTotalList->GetSingle();
        if(!p) break;
        nTracked++;
        dirtest = (EUSO() - p->Pos()).Unit();

        // test if photons really point toward detector
        if((dirtest - p->Dir()).Mag() > TOLERANCE)
            Msg(EsafMsg::Warning) << "<PropagateToDetector()> 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);

        // Dump CPU commentaries
        fraction = 100*Float_t(nTracked)/Float_t(nTotal);
        if (left<fraction && fraction <=right) 
            next = kFALSE;
        else if (fraction>right) {
            next = kTRUE;
            left = right;
            right += 10;
        }
        if (subleft<fraction  && fraction <=subright) {
            subnext = kFALSE;
        }
        if (fraction > subright) {
            subnext = kTRUE;
            subleft = subright;
            subright +=2;
        }
        if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
        if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;	 
    }
    Msg(EsafMsg::Info) << " -done" << MsgDispatch;
#ifdef DEBUG
    Msg(EsafMsg::Debug) << "Number of ParentPhoton = " << fPhotons->GetNphotons() << MsgDispatch;
#endif
}

//_______________________________________________________________________________________________________________________
 void MCRadiativeTransfer::SetOmegaMax(const EarthVector& showermaxpos) {
    //
    // calculate Omega max according to the position of light track maximum
    //
    
    EarthVector pos_omegamax(1);
    Double_t limit_alt = 30*km;
    if(showermaxpos.Zv() > limit_alt) pos_omegamax.SetXYZ(0,0,limit_alt);
    else {    
	// find point of intersection between direction (showermaxpos - euso) and a sphere of radius (EarthRadius + 30km)
        Double_t mag(0);
	EarthVector direc = (showermaxpos - EUSO()).Unit();

	// spherical earth
	Double_t b = EUSO()*direc + direc.Z()*EarthRadius();
	Double_t c = EUSO().Mag2() + 2*EarthRadius()*(EUSO().Z() - limit_alt) - limit_alt*limit_alt;
	pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
	if(p.first == 0) {
	    Msg(EsafMsg::Warning) << "<SetOmegaMax> This case SHOULD not occur with space telescope" << MsgDispatch;
            pos_omegamax.SetXYZ(0,0,limit_alt);
	}
	else if(p.first == 1) mag = p.second[0];
	else if(p.first == 2) mag = min(p.second[0],p.second[1]);
        pos_omegamax = EUSO() + mag*direc;
    }
    fOmegaMax = EusoOmega(pos_omegamax);
#ifdef DEBUG
    Msg(EsafMsg::Debug) << "PosOmegaMax = " <<pos_omegamax<< MsgDispatch;
    Msg(EsafMsg::Debug) << "OmegaMax = " <<fOmegaMax<< MsgDispatch;
#endif
}


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