// ESAF : Euso Simulation and Analysis Framework
// $Id: BunchRadiativeTransfer.cc,v 1.65 2005/05/10 15:34:25 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 "ClearSkyPropagator.hh"
#include "InCloudsPropagator.hh"
#include "BunchOfPhotons.hh"
#include "Ground.hh"
#include "EEvent.hh"
#include "EAtmosphereBunchAdder.hh"
#include "EAtmosphereSingleAdder.hh"
#include "TBenchmark.h"
ClassImp(BunchRadiativeTransfer)
using EConst::EarthRadius;
using EConst::C;
//_______________________________________________________________________________________________________________________
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;
}
//_______________________________________________________________________________________________________________________
BunchRadiativeTransfer::~BunchRadiativeTransfer() {
//
// dtor
//
SafeDelete(fPhotons);
SafeDelete(fCSpropag);
SafeDelete(fICpropag);
}
//_______________________________________________________________________________________________________________________
PhotonsOnPupil* BunchRadiativeTransfer::Get(PhotonsInAtmosphere* photons, const DetectorGeometry* dg) {
//
// Transport photons from source to Euso pupil
//
// set DetectorGeometry for all propagators
fDetGeom = dg;
fCSpropag->CopyDetectorGeometry(fDetGeom);
fICpropag->CopyDetectorGeometry(fDetGeom);
// 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() != "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
// 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();
// 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;
Bool_t next(0),subnext(0);
Double_t nbTracked(0), fraction(0), left(0),right(10),subleft(0),subright(2.5);
Msg(EsafMsg::Info) << "Bunches Processing: 0%" << MsgFlush;
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();
}
}
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;
#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
//
if(fGround) fGround->Reset();
if(fPhotons) fPhotons->Clear();
if(fCSpropag) fCSpropag->Reset();
if(fICpropag) fICpropag->Reset();
fNbBunch = 0;
fNbTot = 0;
}
//_______________________________________________________________________________________________________________________
void BunchRadiativeTransfer::DirectToEuso(const BunchOfPhotons& bigbunch) const {
//
// Creates a list of SinglePhoton considering the EUSO solid angle
//
if(bigbunch.GetType() != Fluo) return;
Int_t nb = 0;
nb = EsafRandom::Get()->Poisson(bigbunch.GetWeight() * EusoOmega(bigbunch.GetPos())/(4*pi)); //TOFIX : here only isotropic source considered so far
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.01*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->Brdf(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;
EarthVector dirtest(1);
RadiativeProcessesCalculator* calcul = RadiativeFactory::Get()->GetRadiativeProcessesCalculator();
Double_t TotTrans;
TRandom* rndm = EsafRandom::Get();
SinglePhoton* p = 0;
Double_t r, phi; //, radius;
TVector3 pos;
// ConfigFileParser* pConfig = Config::Get()->GetCF("General","Euso");
// radius = pConfig->GetNum("Euso.fRadius")*mm;
size_t nTotal(0),nTracked(0);
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) << "Singles Processing: 0%" << MsgFlush;
// if lowtran used and if a light track has been defined,
// build a card for the following transmission process
if(calcul->GetName() == "lowtran" && fTotalList->GetNbTrackSteps()) {
LowtranRadiativeProcessesCalculator* cal = (LowtranRadiativeProcessesCalculator*) calcul;
cal->PreProcessSingles(fGround,EUSO(),*fTotalList);
}
// loop over the list of SinglePhoton
while(true) {
p = fTotalList->GetSingle();
if(!p) break;
nTracked++;
dirtest = (EUSO() - p->Pos()).Unit();
dirtest -= p->Dir();
if(dirtest.Mag() > TOLERANCE) Msg(EsafMsg::Warning) << "PropagationOfSingles : SinglePhoton must be directed toward EUSO" << MsgDispatch;
// calculate transmission from photon position to EUSO
TotTrans = calcul->Trans(*p,EUSO());
if(TotTrans > rndm->Rndm()) {
p->AddToPosTof(EUSO() - p->Pos());
r = EusoRadius() * sqrt(rndm->Rndm());
phi = TMath::TwoPi() * rndm->Rndm();
pos.SetXYZ(r*cos(phi),r*sin(phi),0);
fPhotons->Add(*p,pos);
}
else p->Absorbed();
// 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
SafeDelete(calcul);
}
//_______________________________________________________________________________________________________________________
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())/C();
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())/C();
// 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);
fTotalList->Add(s);
}
}