// ESAF : Euso Simulation and Analysis Framework
// $Id: O2_ClearSkyPropagator.cc,v 1.40 2005/05/10 15:34:25 moreggia Exp $
// Sylvain Moreggia created Jun, 2 2004
#include "O2_ClearSkyPropagator.hh"
#include "SinglePhoton.hh"
#include "BunchOfPhotons.hh"
#include <math.h>
#include "EsafRandom.hh"
#include "ListPhotonsInAtmosphere.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "RadiativeFactory.hh"
#include "EConst.hh"
#include "Atmosphere.hh"
#include "Clouds.hh"
#include "LowtranRadiativeProcessesCalculator.hh"
#include "TestGround.hh"
#include "TBenchmark.h"
ClassImp(O2_ClearSkyPropagator)
using EConst::C;
using EConst::EarthRadius;
//_____________________________________________________________________________
O2_ClearSkyPropagator::O2_ClearSkyPropagator() : ClearSkyPropagator() {
//
// ctor (should not be used)
//
fOrder = 2;
fFinal_medium = NONE;
}
//_____________________________________________________________________________
O2_ClearSkyPropagator::O2_ClearSkyPropagator(const Ground* g) : ClearSkyPropagator(g) {
//
// ctor, copy RadiatvieTransfer ground description
//
fOrder = 2;
fFinal_medium = NONE;
}
//_____________________________________________________________________________
O2_ClearSkyPropagator::~O2_ClearSkyPropagator() {
//
// dtor
//
}
//_____________________________________________________________________________
Medium O2_ClearSkyPropagator::Go(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
//
// call methods for propagation
// NB: bunch can be modified, its copy b cannot
//
// to assess backscattering algorithm time consumption
#ifdef DEBUG
char prt[100];
sprintf(prt,"Bunch order2 propagation");
TBenchmark bench;
bench.Start(prt);
#endif
// initializations
static UInt_t bunch_id = 0; // to know if bunch has already come here (for clouds etc..)
Medium rtn = NONE;
const BunchOfPhotons& b = bunch;
EarthVector impact;
// if fluo, not propagated //TOFIX: fluo propagation should be handled in the future
if(b.GetType() == Fluo) {
bunch.SetFate(1);
return NONE;
}
// if first treatment for this bunch,
// get next impact and next medium encountered
//TOFIX : valid algorithm if a SINGLE intermediate medium exists
if(bunch_id != b.GetId()) {
// bunch preprocessing -> generate a lowtran table along the bunch track
BunchPreProcess(bunch);
if(fFinal_medium == NONE) return NONE;
bunch_id = b.GetId();
rtn = GetNextImpact(b,impact);
bunch.SetNextImpact(impact);
if(impact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an "impact"" << MsgDispatch;
if(b.GetFinalImpact().Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an "impact"" << MsgDispatch;
// if no intermediate medium, go to final impact
if((b.GetNextImpact() - b.GetFinalImpact()).Mag() == 0) rtn = fFinal_medium;
}
// otherwise go to final impact
else {
bunch.SetNextImpact(b.GetFinalImpact());
rtn = fFinal_medium;
}
// if bunch showerpos is within clouds, switch off the present clear sky propagation
if(rtn == CLOUDY && (b.GetNextImpact() - b.GetPos()).Mag() < TOLERANCE) {
return rtn;
}
// run propagation algorithms
Propagate(bunch,list);
#ifdef DEBUG
bench.Stop(prt);
MsgForm(EsafMsg::Debug,"REAL=%6.2f s CPU=%6.2f s",bench.GetRealTime(prt),bench.GetCpuTime(prt));
#endif
if(rtn == GROUND) bunch.SetFate(4);
// in the case impact is out of FoV, rtn set to NONE to stop the propagation (no ground reflexion)
if(rtn == CLEARSKY) {
rtn = NONE;
bunch.SetFate(5);
}
return rtn;
}
//_____________________________________________________________________________
void O2_ClearSkyPropagator::Propagate(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
//
// Propagate a BunchOfPhotons in clear sky conditions,
// generating SinglePhotons all along the travel -> ListPhotonsInAtmosphere filled
// NB: bunch can be modified, its copy b cannot
//
const BunchOfPhotons& b = bunch;
const EarthVector& impact = b.GetNextImpact();
const Atmosphere* atmo = Atmosphere::Get();
// initializations
Double_t conv = 1.;
Int_t nb = b.GetWlSpectrum().GetAxis().GetNbins();
vector<Double_t*> coeff;
Double_t* TotalTrans = new Double_t[nb];
Double_t* RayleighTrans = new Double_t[nb];
Double_t* OzoneTrans = new Double_t[nb];
Double_t* AerosolTrans = new Double_t[nb];
if(!TotalTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for TotalTrans" << MsgDispatch;
if(!RayleighTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for RayleighTrans" << MsgDispatch;
if(!OzoneTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for OzoneTrans" << MsgDispatch;
if(!AerosolTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for AerosolTrans" << MsgDispatch;
coeff.push_back(TotalTrans);
coeff.push_back(RayleighTrans);
coeff.push_back(OzoneTrans);
coeff.push_back(AerosolTrans);
Double_t raylScat,phfunc;
Int_t NbRaylBackscat = 0;
EarthVector pos, step, entry, exit, dir;
Double_t depthstep = Conf()->GetNum("BunchRadiativeTransfer.DepthStep")*g/cm2;
entry = b.GetPos() - 0.5*(b.GetParent()->GetShowerPosf() - b.GetParent()->GetShowerPosi());
dir = (impact - entry).Unit();
exit = atmo->InvertGrammage(entry,dir,depthstep);
#ifdef DEBUG
// to get nb of iteration for propagation of this bunch
Float_t count = 0;
#endif
while((impact - b.GetPos()).Mag() > TOLERANCE) {
#ifdef DEBUG
count++;
#endif
step = exit - entry;
if((impact - b.GetPos()).Mag() < step.Mag()) step = impact - b.GetPos();
fCalcul->Trans(b,b.GetPos()+step,coeff);
phfunc = fCalcul->PhaseFunction("rayleigh",b,EUSO() - b.GetPos() - 0.5*step,EusoOmega(b.GetPos() + 0.5*step)); //TOFIX
// position and time of flight incrementation (bunch position is at the step end)
bunch.AddToPosTof(step);
// cerenkov backscattering (only rayleigh so far) //TOFIX
for(Int_t i=0; i<nb; i++) {
//raylScat = (1 - coeff[1][i])*coeff[2][i]*coeff[3][i];
if(coeff[1][i]*coeff[2][i] + coeff[1][i]*coeff[3][i] + coeff[2][i]*coeff[3][i] - 3*coeff[0][i]) {
raylScat = (1 - coeff[1][i])*coeff[2][i]*coeff[3][i]*(1-coeff[0][i]) /
(coeff[1][i]*coeff[2][i] + coeff[1][i]*coeff[3][i] + coeff[2][i]*coeff[3][i] - 3*coeff[0][i]);
}
NbRaylBackscat = EsafRandom::Get()->Poisson(b.GetWeight(i) * raylScat*phfunc);
// Singles generation (bunch position is at the step end)
if(NbRaylBackscat) GenerateSingles(i,NbRaylBackscat,raylScat,step,b,list);
}
// convolutions with transmission coefficients
conv = bunch.ConvoluteWl(coeff[0]);
bunch.Convolute(conv);
entry = exit;
exit = atmo->InvertGrammage(entry,dir,depthstep);
}
// to have position precisely, for later processes
bunch.SetPos(impact);
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"Nb iteration = " <<count << MsgDispatch;
Msg(EsafMsg::Debug) <<"after backscat, SinglePhoton list size = "<<list.GetListOfSingle().size() << MsgDispatch;
#endif
// cleaning of coeff
for(UInt_t m=0; m<coeff.size(); m++) {
if(coeff[m]) delete[] coeff[m];
coeff[m] = 0;
}
coeff.clear();
}
//_____________________________________________________________________________
void O2_ClearSkyPropagator::GenerateSingles(Int_t specbin, Int_t nb, Double_t scattrate, const EarthVector& step,
const BunchOfPhotons& b, ListPhotonsInAtmosphere& list) const {
//
// Generate a SinglePhoton list, coming from backscattering process applied to the given BunchOfPhotons
// Tuning of scattering position done thanks to additional data : scattrate (rate of scattering),
// step (step used for backscattering calculations)
// NB : bunch position is at the end of the considered propagation step
//
Double_t wl, date, tof;
EarthVector showerpos, pos, dir, diff;
SinglePhoton* s = 0;
TRandom* rndm = EsafRandom::Get();
UInt_t bid = b.GetId();
PhotonType type = b.GetType();
for(Int_t i=0; i<nb; i++) {
// corrections for date and showerpos, from BunchOfPhotons mean values, longitudinal and lateral distributions at creation
showerpos = b.RandomPosInShower();
if(fGround->IsUnderGround(showerpos)) continue;
diff = showerpos - b.GetShowerPos();
date = b.GetDate() + diff.Dot(b.GetDir().Unit())/C();
// calculation of the scattering position within the step
pos = b.GetPos() - step; // as if bunch was at the beginning of the step
pos += (log(1 - rndm->Rndm()*scattrate)/log(1 - scattrate))*step; //TOFIX: only rayleigh
// position correction, due to angular distributions at creation
pos = b.PosRandomAngCorrec(showerpos,pos);
// if underground, photon is kept for another try
if(fGround->IsUnderGround(pos)) {
EarthVector axis, impact;
impact = b.GetNextImpact();
if(impact.Z() == HUGE) {
Msg(EsafMsg::Warning) <<"nno bunch impact, but singles underground" <<MsgDispatch;
continue;
}
else {
EarthVector localvertical(impact);
localvertical(2) += EarthRadius();
axis = localvertical.Cross(b.GetDir());
if(axis.Mag()) { // if localvertical and bunch dir not colinear
EarthVector v1 = pos - impact;
v1.Rotate(TMath::Pi(),axis);
pos = v1 + impact;
}
if(fGround->IsUnderGround(pos)) {
Msg(EsafMsg::Warning) <<"GenerateSingles() -> pos is still underground" <<MsgDispatch;
continue;
}
}
}
tof = b.GetTof() * (showerpos - pos).Mag()/(b.GetPos() - b.GetShowerPos()).Mag();
dir = EUSO() - pos;
wl = b.GetWlSpectrum().GetLambda(specbin);
s = new SinglePhoton(type,date,tof,wl,showerpos,pos,dir,RaylScat,bid);
list.Add(s);
}
}
//_____________________________________________________________________________
Medium O2_ClearSkyPropagator::BunchPreProcess(BunchOfPhotons& b) const {
//
// Preprocess run once for each bunch before propagation
// Call RadiativeProcessesCalculator::BunchPreProcess
// Get the Medium of final impact
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"IN O2_ClearSkyPropagator::bunchPreProcess" << MsgDispatch;
#endif
const BunchOfPhotons& bunch = b;
EarthVector finalimpact;
// assess the next impact
fFinal_medium = GetFinalImpact(bunch,finalimpact);
b.SetFinalImpact(finalimpact);
if(finalimpact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an "impact"" << MsgDispatch;
// if bunch is underground, its simulation is stopped now
if(fFinal_medium == NONE) {
b.SetFate(3);
return fFinal_medium;
}
// run Calculator BunchPreProcessing
if(fCalcul->GetName() == "lowtran") {
LowtranRadiativeProcessesCalculator* calcul = (LowtranRadiativeProcessesCalculator*) fCalcul;
calcul->PreProcess(bunch);
}
return fFinal_medium;
}
//_____________________________________________________________________________
void O2_ClearSkyPropagator::Reset() {
//
// reset imbricated objects
//
fCalcul->Reset();
fFinal_medium = NONE;
}