#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 namespace sou;
using namespace EConst;
O2_ClearSkyPropagator::O2_ClearSkyPropagator() : ClearSkyPropagator() {
fOrder = 2;
fFinal_medium = NONE;
}
O2_ClearSkyPropagator::O2_ClearSkyPropagator(const Ground* g) : ClearSkyPropagator(g) {
fOrder = 2;
fFinal_medium = NONE;
fNbSinglesCorrected = 0;
ConfigFileParser *pConfig = Config::Get()->GetCF("Atmosphere", "LowtranAtmosphere");
fAerosolModel = (Int_t)pConfig->GetNum("LowtranAtmosphere.Ihaze");
if(fAerosolModel != 0 && fAerosolModel != 1 && fAerosolModel != 2 && fAerosolModel != 4 && fAerosolModel != 5) {
fAerosolModel = 0;
Msg(EsafMsg::Warning) <<"Aerosol model is such that scattering simu cannot be done using BUNCH algorithm" << MsgDispatch;
Msg(EsafMsg::Warning) <<"Aerosol model will be used only for last transmission to detector" << MsgDispatch;
}
if(fAerosolModel == 0) fAerosol_correction_factor = 0.;
else if(fAerosolModel == 1) fAerosol_correction_factor = 0.95;
else if(fAerosolModel == 2) fAerosol_correction_factor = 0.95;
else if(fAerosolModel == 4) fAerosol_correction_factor = 0.98;
else if(fAerosolModel == 5) fAerosol_correction_factor = 0.70;
}
O2_ClearSkyPropagator::~O2_ClearSkyPropagator() {
}
Medium O2_ClearSkyPropagator::Go(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
#ifdef DEBUG
char prt[100];
sprintf(prt,"Bunch order2 propagation");
//root.cern.ch/root/html/TBenchmark.html">TBenchmark bench;
bench.Start(prt);
#endif
static UInt_t bunch_id = 0;
Medium rtn = NONE;
const BunchOfPhotons& b = bunch;
EarthVector impact;
if(b.GetType() == Fluo) {
bunch.SetFate(1);
return NONE;
}
if(bunch_id != b.GetId()) {
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((b.GetNextImpact() - b.GetFinalImpact()).Mag() == 0) rtn = fFinal_medium;
}
else {
bunch.SetNextImpact(b.GetFinalImpact());
rtn = fFinal_medium;
}
if(rtn == CLOUDY && (b.GetNextImpact() - b.GetPos()).Mag() < TOLERANCE) {
return rtn;
}
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);
if(rtn == CLEARSKY) {
rtn = NONE;
bunch.SetFate(5);
}
return rtn;
}
void O2_ClearSkyPropagator::Propagate(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
Int_t status = -100;
const BunchOfPhotons& b = bunch;
const EarthVector& impact = b.GetNextImpact();
const Atmosphere* atmo = Atmosphere::Get();
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 = 0;
Double_t phfunc = 0;
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();
status = atmo->InvertGrammage(entry,dir,depthstep,exit);
#ifdef DEBUG
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->RayleighPhaseFunction(b.GetDir(),EUSO() - b.GetPos() - 0.5*step);
bunch.AddToPosTof(step);
for(Int_t i=0; i<nb; i++) {
raylScat = (1 - coeff[1][i]);
NbRaylBackscat = EsafRandom::Get()->Poisson(b.GetWeight(i) * EusoOmega(b.GetPos() - 0.5*step) * raylScat*phfunc);
if(NbRaylBackscat) GenerateRayleighSingles(i,NbRaylBackscat,raylScat,step,b,list,coeff[3][i],coeff[2][i]);
}
conv = bunch.ConvoluteWl(coeff[0]);
bunch.Convolute(conv);
entry = exit;
status = atmo->InvertGrammage(entry,dir,depthstep,exit);
}
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
for(UInt_t m=0; m<coeff.size(); m++) {
if(coeff[m]) delete[] coeff[m];
coeff[m] = 0;
}
coeff.clear();
}
void O2_ClearSkyPropagator::GenerateRayleighSingles(Int_t specbin, Int_t nb, Double_t scattrate, const EarthVector& step,
const BunchOfPhotons& b, ListPhotonsInAtmosphere& list, Double_t Ta,Double_t To) const {
Double_t wl, date, tof, scatlength;
EarthVector showerpos, pos, dir, diff;
SinglePhoton* s = 0;
//root.cern.ch/root/html/TRandom.html">TRandom* rndm = EsafRandom::Get();
UInt_t bid = b.GetId();
PhotonType type = b.GetType();
for(Int_t i=0; i<nb; i++) {
showerpos = b.RandomPosInShower();
if(fGround->IsUnderGround(showerpos)) continue;
diff = showerpos - b.GetShowerPos();
date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
pos = b.GetPos() - step;
scatlength = (log(1 - rndm->Rndm()*scattrate)/log(1 - scattrate));
pos += scatlength *step;
if(rndm->Rndm() > exp(scatlength*log(Ta))) continue;
if(rndm->Rndm() > exp(scatlength*log(To))) continue;
pos = b.PosRandomAngCorrec(showerpos,pos);
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()) {
EarthVector v1 = pos - impact;
v1.Rotate(//root.cern.ch/root/html/TMath.html">TMath::Pi(),axis);
pos = v1 + impact;
}
else Msg(EsafMsg::Warning) <<"<GenerateSingles> OK if vertical track" <<MsgDispatch;
if(fGround->IsUnderGround(pos)) {
Msg(EsafMsg::Warning) <<"GenerateSingles() -> pos is still underground" <<MsgDispatch;
continue;
}
fNbSinglesCorrected++;
}
}
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,b.GetShowerAge());
s->AddInteraction();
s->AddHistory(RaylScat);
list.Add(s);
}
}
void O2_ClearSkyPropagator::GenerateAerosolSingles(Int_t specbin, Int_t nb, Double_t scattrate, const EarthVector& step,
const BunchOfPhotons& b, ListPhotonsInAtmosphere& list, Double_t Tr,Double_t To) const {
Double_t wl, date, tof, scatlength;
EarthVector showerpos, pos, dir, diff;
SinglePhoton* s = 0;
//root.cern.ch/root/html/TRandom.html">TRandom* rndm = EsafRandom::Get();
UInt_t bid = b.GetId();
PhotonType type = b.GetType();
for(Int_t i=0; i<nb; i++) {
showerpos = b.RandomPosInShower();
if(fGround->IsUnderGround(showerpos)) continue;
diff = showerpos - b.GetShowerPos();
date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
pos = b.GetPos() - step;
scatlength = (log(1 - rndm->Rndm()*scattrate)/log(1 - scattrate));
pos += scatlength *step;
if(rndm->Rndm() > exp(scatlength*log(Tr))) continue;
if(rndm->Rndm() > exp(scatlength*log(To))) continue;
pos = b.PosRandomAngCorrec(showerpos,pos);
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()) {
EarthVector v1 = pos - impact;
v1.Rotate(//root.cern.ch/root/html/TMath.html">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,AeroScat,bid,b.GetShowerAge());
s->AddInteraction();
s->AddHistory(AeroScat);
list.Add(s);
}
}
Medium O2_ClearSkyPropagator::BunchPreProcess(BunchOfPhotons& b) const {
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"IN O2_ClearSkyPropagator::bunchPreProcess" << MsgDispatch;
#endif
const BunchOfPhotons& bunch = b;
EarthVector finalimpact;
fFinal_medium = GetFinalImpact(bunch,finalimpact);
b.SetFinalImpact(finalimpact);
if(finalimpact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an \"impact\"" << MsgDispatch;
if(fFinal_medium == NONE) {
b.SetFate(3);
return fFinal_medium;
}
if(fCalcul->GetName() == "lowtran") {
LowtranRadiativeProcessesCalculator* calcul = (LowtranRadiativeProcessesCalculator*) fCalcul;
calcul->PreProcess(bunch);
}
return fFinal_medium;
}
void O2_ClearSkyPropagator::Reset() {
fCalcul->Reset();
fFinal_medium = NONE;
fNbSinglesCorrected = 0;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.