// $Id: SinglePhotonPropagator.cc,v 1.15 2005/11/03 12:01:33 moreggia Exp $
// Author: S. Moreggia 2005/08/16
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: SinglePhotonPropagator *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// SinglePhotonPropagator
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "SinglePhotonPropagator.hh"
#include "SinglePhoton.hh"
#include "FastRadiativeProcessesCalculator.hh"
#include "Clouds.hh"
#include "Atmosphere.hh"
#include "EsafRandom.hh"
#include "Ground.hh"
#include "EConst.hh"
using namespace sou;
using namespace TMath;
using namespace EConst;
ClassImp(SinglePhotonPropagator)
//_____________________________________________________________________________
SinglePhotonPropagator::SinglePhotonPropagator() : EsafMsgSource() {
//
// Constructor
//
Msg(EsafMsg::Panic) <<"This constructor SHOULD not be used : ground must be given in argument" << MsgDispatch;
}
//_____________________________________________________________________________
SinglePhotonPropagator::SinglePhotonPropagator(const Ground* g, Double_t maxphfunc, Double_t maxtof) {
//
// Constructor which copy RT ground description
//
CopyGround(g);
fCalcul = new FastRadiativeProcessesCalculator();
fDetGeom = 0;
fMaxPhaseFunction = maxphfunc;
fTofCut = maxtof;
}
//_____________________________________________________________________________
SinglePhotonPropagator::~SinglePhotonPropagator() {
//
// Destructor
//
SafeDelete(fCalcul);
}
//_____________________________________________________________________________
void SinglePhotonPropagator::Go(const vector<SinglePhoton*>& list,Int_t scatorder) const {
//
// MC code for photon scattering simulation
//
size_t nbPh = list.size();
Msg(EsafMsg::Info) <<"Nb of SinglePhoton to propagate = " << nbPh << MsgDispatch;
// for future info dumping
Bool_t next(0),subnext(0);
Double_t nbTracked(0), fraction(0), left(0),right(10),subleft(0),subright(2.5);
Msg(EsafMsg::Info) << "Scattering Simulation: 0%" << MsgFlush;
// loop over the list
SinglePhoton* s = 0;
for(size_t i=0; i<nbPh; i++) {
s = list[i];
DoOnePhoton(*s,scatorder);
nbTracked++;
// dump info
fraction = 100*nbTracked/Double_t(nbPh);
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;
}
//_____________________________________________________________________________
void SinglePhotonPropagator::DoOnePhoton(SinglePhoton& s,Int_t scatorder) const {
//
// MC propagation : handle rayleigh, clouds and ground for fluo and cerenkov
// Each scattering order simulation is independent from another
//
const Clouds* clouds = Atmosphere::Get()->GetClouds();
const Atmosphere* atmo = Atmosphere::Get();
Int_t nbInter = 0;
Double_t TOF = 0.;
Int_t invert_status(0); // for Atmosphere::InvertGrammage() method
EarthVector dir(1.);
Int_t count_debug(0); // to stop simu if infinite loop
Bool_t IsCloudInter = false;
// init
TRandom* rndm = EsafRandom::Get();
EarthVector nextimpact(1), scatpos_rayl(0.,0,HUGE), scatpos_clouds(0,0,HUGE);
Medium nextmedium = CLEARSKY;
// propagation stops when :
// - photon is absorbed
// - photon is out of atmosphere
// - relevant scattering order is reached
while(true) {
// if weird situations appear
if(nbInter == scatorder) Msg(EsafMsg::Warning) <<"<DoOnePhoton> Loop should have stopped before"<<MsgDispatch;
if(count_debug++ > 100) {
Msg(EsafMsg::Warning) <<"<DoOnePhoton> infinite loop broken by hand "<<MsgDispatch;
break;
}
// cut in time of flight defined here
if(TOF > fTofCut) {
s.SetStatus(LostByCut);
#ifdef DEBUG
//Msg(EsafMsg::Debug) <<"<DoOnePhoton> TOF cut reached : photon is lost"<<MsgDispatch;
#endif
break;
}
// if within clouds
// - single scattering albedo ~ 1 -> no absorption by water droplets or ice crystals (cf. Liou)
// - scattering cross-section is ~ lambda independent in UV band (cf. Liou)
// so far works if clouds is a SINGLE LAYER medium
if(clouds->IsInClouds(s.Pos()) || (nextmedium == CLOUDY) ) {
// 1. get ground impact
nextmedium = GetNextImpact(s.Pos(),s.Dir(),nextimpact,true);
// 2. get random clouds scattering position,
// check it is within clouds. If not -> scatpos_clouds = (0,0,HUGE)
// check it is before ground impact. If not -> scatpos_clouds = (0,0,HUGE)
scatpos_clouds = clouds->RandomScatPos(s.Pos(),s.Dir());
if(!clouds->IsInClouds(scatpos_clouds)) scatpos_clouds.SetXYZ(0,0,HUGE);
if(fGround->IsUnderGround(scatpos_clouds)) scatpos_clouds.SetXYZ(0,0,HUGE);
// 3. get random rayleigh scattering position
// check it is before ground impact. If not -> grammage_rayl = -1
Double_t grammage_rayl = fCalcul->RandomGrammage_RaylScat(s.Wl());
if(atmo->Grammage(s.Pos(),nextimpact) < grammage_rayl) grammage_rayl = -1;
// 4. if ground reached before any scattering,
// position sets to impact and goes to next step
if((scatpos_clouds.Z() == HUGE) && (grammage_rayl == -1) ) {
TOF += (nextimpact - s.Pos()).Mag() / Clight();
s.AddToPosTof(nextimpact - s.Pos());
continue;
}
// else photon interacts with clouds or air
s.AddInteraction();
nbInter++;
// 5. compare the two scattering positions and choose the nearest
IsCloudInter = false;
if(scatpos_clouds.Z() == HUGE) IsCloudInter = false;
else if((atmo->Grammage(s.Pos(),scatpos_clouds) < grammage_rayl) || (grammage_rayl == -1)) IsCloudInter = true;
// then the corresponding type of scattering is applied
// 5.1 clouds scattering
if(IsCloudInter) {
TOF += (scatpos_clouds - s.Pos()).Mag() / Clight();
s.AddToPosTof(scatpos_clouds - s.Pos());
s.AddHistory(CloudScat);
// 5.10 if photon has come to its order of scattering
// look if it points toward detector and stop its simu
if(scatorder == nbInter) {
if(rndm->Rndm() < (clouds->PhaseFunction(s.Dir(),(EUSO() - s.Pos()).Unit()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
s.SetStatus(CloudScat);
s.SetDir(EUSO() - s.Pos());
}
else s.SetStatus(LostByCut);
break;
}
// 5.11 else chose a new direction, and goes to next step
else {
dir = s.Dir();
clouds->RandomDir(dir);
s.SetDir(dir);
s.SetStatus(LostByCut);
nextmedium = CLOUDY;
}
}
// 5.2 or rayleigh scattering
else {
if(grammage_rayl > -1) invert_status = atmo->InvertGrammage(s.Pos(),s.Dir(),grammage_rayl,scatpos_rayl,fTofCut - TOF);
else Msg(EsafMsg::Warning) <<"<DoOnePhoton> This situation SHOULD NOT occur"<<MsgDispatch;
if(invert_status > 2) Msg(EsafMsg::Warning) <<"<DoOnePhoton> Atmo::InvertGrammage() returns status = " <<invert_status<<MsgDispatch;
if(invert_status == 1) {
s.SetStatus(LostByCut);
break;
}
TOF += (scatpos_rayl - s.Pos()).Mag() / Clight();
s.AddToPosTof(scatpos_rayl - s.Pos());
s.AddHistory(RaylScat);
// 5.20 if photon has come to its order of scattering
// look if it points toward detector and stop its simu
if(scatorder == nbInter) {
if(rndm->Rndm() < (fCalcul->RayleighPhaseFunction(s.Dir(),(EUSO() - s.Pos()).Unit()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
s.SetStatus(RaylScat);
s.SetDir(EUSO() - s.Pos());
}
else s.SetStatus(LostByCut);
break;
}
// 5.21 else chose a new direction, and goes to next step
else {
dir = s.Dir();
fCalcul->RandomDir("rayleigh",dir);
s.SetDir(dir);
s.SetStatus(LostByCut);
nextmedium = CLEARSKY;
}
}
}
// if clear sky
else if(nextmedium == CLEARSKY) {
// 1. get next medium type, and impact position upon it
nextmedium = GetNextImpact(s.Pos(),s.Dir(),nextimpact);
// 2. get random rayleigh scattering position, expressed in grammage
Double_t grammage_rayl = fCalcul->RandomGrammage_RaylScat(s.Wl());
// 3. if photon does not rayleigh-interact before next impact
// position sets to impact and goes to next step
if(atmo->Grammage(s.Pos(),nextimpact) < grammage_rayl) {
TOF += (nextimpact - s.Pos()).Mag() / Clight();
s.AddToPosTof(nextimpact - s.Pos());
continue;
}
// 4. else photon is rayleigh scattered at corresponding position
s.AddInteraction();
nbInter++;
invert_status = atmo->InvertGrammage(s.Pos(),s.Dir(),grammage_rayl,scatpos_rayl,fTofCut - TOF);
if(invert_status > 1) Msg(EsafMsg::Warning) <<"<DoOnePhoton> Atmo::InvertGrammage() returns status = " <<invert_status<<MsgDispatch;
if(invert_status == 1) {
s.SetStatus(LostByCut);
break;
}
TOF += (scatpos_rayl - s.Pos()).Mag() / Clight();
s.AddToPosTof(scatpos_rayl - s.Pos());
s.AddHistory(RaylScat);
// 5. if photon has come to its order of scattering
// look if it points toward detector and stop its simu
if(scatorder == nbInter) {
if(rndm->Rndm() < (fCalcul->RayleighPhaseFunction(s.Dir(),(EUSO() - s.Pos()).Unit()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
s.SetStatus(RaylScat);
s.SetDir(EUSO() - s.Pos());
}
else s.SetStatus(LostByCut);
break;
}
// 6. else chose a new direction, and goes to next step
else {
dir = s.Dir();
fCalcul->RandomDir("rayleigh",dir);
s.SetDir(dir);
s.SetStatus(LostByCut);
nextmedium = CLEARSKY;
}
}
// if on ground
else if(nextmedium == GROUND) {
s.AddInteraction();
s.AddHistory(Reflected);
nbInter++;
// 1. if not absorbed at ground
if(rndm->Rndm() < fGround->Albedo(s.Pos())) {
// 1.1. if photon has come to its order of scattering
// look if it points toward detector and stop its simu
if(scatorder == nbInter) {
if(rndm->Rndm() < (fGround->Outgoing_phase_function(s.Pos(),EUSO()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
s.SetStatus(Reflected);
s.SetDir(EUSO() - s.Pos());
}
else s.SetStatus(LostByCut);
break;
}
// 1.2. else get a new direction and go to next step
else {
fGround->RandomDir(s.Pos(),dir);
s.SetDir(dir);
s.SetStatus(LostByCut);
nextmedium = CLEARSKY;
}
}
// 2. else photon is absorbed at ground
else {
s.SetAbsorbed();
s.SetStatus(ScatAbsorb);
break;
}
}
else if(nextmedium == TOA) {
s.AddHistory(OutAtmo);
s.SetStatus(OutAtmo);
break;
}
else if(nextmedium == NONE) Msg(EsafMsg::Panic) <<"<DoOnePhoton> next medium is NONE -> SHOULD NOT" <<MsgDispatch;
}
}
//_____________________________________________________________________________
Medium SinglePhotonPropagator::GetNextImpact(const EarthVector& pos,const EarthVector& dir,EarthVector& res,Bool_t without_clouds) const {
//
// Calculate the impact upon the next medium according to the given track
// Return type of the next encountered medium
//
// - If pos is underground, NONE medium is returned and impact set to (0,0,-HUGE)
// OUT OF TOA is handled
//
Medium rtn = GROUND;
EarthVector gd_imp, cloud_imp;
// get impact on ground
gd_imp = fGround->GetImpact(pos,dir);
res = gd_imp;
// if pos is underground
if(gd_imp.Z() == -HUGE) return NONE;
// get impact on clouds, then make appropriate changes
cloud_imp = Atmosphere::Get()->GetClouds()->GetCloudImpact(pos,dir);
if( ((cloud_imp - pos).Mag() < (gd_imp - pos).Mag()) && !without_clouds ) {
res = cloud_imp;
rtn = CLOUDY;
}
// if no impact, means top of atmosphere reached
if(res.Z() == HUGE) {
rtn = TOA;
res = Atmosphere::Get()->ImpactAtTOA(pos,dir,50*km); // (should be same as in atmosphere::InvertGrammage() method) TOA at 50km to avoid pbic cases at high altitudes
}
return rtn;
}
//_____________________________________________________________________________
void SinglePhotonPropagator::Reset() {
//
// ready for next event
//
fCalcul->Reset();
}