// $Id: AlongTrack_CSPropagator.cc,v 1.18 2005/05/10 15:34:25 moreggia Exp $
// Author: Sylvain Moreggia 2004/09/29
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: AlongTrack_CSPropagator *
* Package: RadiativeTransfer *
* Coordinator: Sylvain Moreggia *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// AlongTrack_CSPropagator
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "AlongTrack_CSPropagator.hh"
#include "TBenchmark.h"
#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"
ClassImp(AlongTrack_CSPropagator)
using EConst::C;
//_____________________________________________________________________________
AlongTrack_CSPropagator::AlongTrack_CSPropagator() : O2_ClearSkyPropagator() {
//
// Constructor (should not be used)
//
fOrder = 2;
fFlag = false;
fFinal_medium = CLEARSKY;
fFirstBunch = 0;
}
//_____________________________________________________________________________
AlongTrack_CSPropagator::AlongTrack_CSPropagator(const Ground* g) : O2_ClearSkyPropagator(g) {
//
// Constructor, copy RadiatvieTransfer ground description
//
fOrder = 2;
fFlag = false;
fFinal_medium = CLEARSKY;
fFirstBunch = 0;
}
//_____________________________________________________________________________
AlongTrack_CSPropagator::~AlongTrack_CSPropagator() {
//
// Destructor
//
SafeDelete(fFirstBunch);
}
//_____________________________________________________________________________
Medium AlongTrack_CSPropagator::Go(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
//
// call preprocessing, then propagation method
//
// to assess backscattering algorithm time consumption
#ifdef DEBUG
char prt[100];
sprintf(prt,"Bunch ALONG TRACK 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..)
const BunchOfPhotons& b = bunch;
EarthVector impact;
Medium rtn = NONE;
// if fluo, not propagated //TOFIX: fluo propagation should be handled in the future
if(b.GetType() == Fluo) {
bunch.SetFate(1);
return NONE;
}
// run preprocess -> generation of lowtran table along the track
if(!fFlag) PreProcess(b);
// if PreProcess failed because track beginning is under ground
if(!fFlag) {
bunch.SetFate(3);
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_id = b.GetId();
bunch.SetFinalImpact(fFirstBunch->GetFinalImpact());
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;
}
// test if propagator is really dealing with a track
if(!IsAlongTrack(b)) Msg(EsafMsg::Warning) << "Go() : AlongTrack used, but it does not seem to be a track" << MsgDispatch;
// if bunch is under ground, propagation is stopped
if((b.GetFinalImpact() - fFirstBunch->GetPos()).Mag() < (b.GetPos() - fFirstBunch->GetPos()).Mag()) {
bunch.SetFate(3);
return NONE;
}
// 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
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 the impact is out of FoV, NONE is returned to stop the propagation (no ground reflexion)
if(rtn == CLEARSKY) {
rtn = NONE;
bunch.SetFate(5);
}
return rtn;
}
//_____________________________________________________________________________
void AlongTrack_CSPropagator::Propagate(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
//
// Propagate along a track a BunchOfPhotons in clear sky conditions,
// generating SinglePhotons all along the travel -> ListPhotonsInAtmosphere filled
// NB: bunch can be modified, its copy b cannot
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"Propagate()" << MsgDispatch;
#endif
// initializations
const BunchOfPhotons& b = bunch;
const EarthVector& impact = b.GetNextImpact();
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),phfunc(0);
Int_t NbRaylBackscat = 0;
EarthVector pos, step, entry, exit;
if(list.GetNbTrackSteps() == 0) Msg(EsafMsg::Panic) << "this propagator needs ListPhotonsInAtmosphere::fTrack" << MsgDispatch;
// find the matching track step
UInt_t index = list.TrackStepIndex(b.GetPos());
entry = b.GetPos();
if(index < (list.GetNbTrackSteps()-1)) exit = list.GetTrackStep(++index);
else {
Msg(EsafMsg::Warning) << "pb with array index : a bunch position is OUTSIDE light track"<<MsgDispatch;
exit = impact;
}
#ifdef DEBUG
// to get nb of iteration for propagation of this bunch
Float_t count = 0;
#endif
// loop of step by step propagation for backscattering simulation
while((impact - b.GetPos()).Dot(impact - b.GetShowerPos()) > TOLERANCE) {
#ifdef DEBUG
count++;
#endif
step = exit - entry;
if((impact - (b.GetPos()+step)).Dot(impact - b.GetPos()) < TOLERANCE) 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 only rayleigh here
// position and time of flight incrementation (bunch position is at the step end)
bunch.AddToPosTof(step);
bunch.SetPos(exit); // often useless, but SAFER : handle misaligned bunches
// 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;
if(index != list.GetNbTrackSteps()-1) {
index++;
exit = list.GetTrackStep(index);
}
else exit = impact;
}
// 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 AlongTrack_CSPropagator::PreProcess(const BunchOfPhotons& b) const {
//
// Set fFlag and fImpact, Call RadiativeProcessesCalculator::PreProcess()
// Called once (for the first bunch) and used for propagation of all the bunches
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"PreProcess()" << MsgDispatch;
#endif
// get final impact
EarthVector finalImpact;
fFinal_medium = GetFinalImpact(b,finalImpact);
fFirstBunch = new BunchOfPhotons(b);
fFirstBunch->SetFinalImpact(finalImpact);
if(finalImpact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an "impact"" << MsgDispatch;
if(fFinal_medium != NONE) {
// create lowtran trans table along the whole track (until final track impact)
if(fCalcul->GetName() == "lowtran") {
LowtranRadiativeProcessesCalculator* calcul = (LowtranRadiativeProcessesCalculator*) fCalcul;
calcul->PreProcess(*fFirstBunch);
}
fFlag = true;
}
// if first bunch is under ground, it is not used to initialize the propagator. The next bunch will be
else {
Msg(EsafMsg::Info) <<"PreProcess() : track starts underground" << MsgDispatch;
SafeDelete(fFirstBunch);
}
}
//_____________________________________________________________________________
Bool_t AlongTrack_CSPropagator::IsAlongTrack(const BunchOfPhotons& b) const {
//
// Check if bunches are along a track, with first at the track beginning
//
if(b.GetId() == fFirstBunch->GetId()) return true;
Bool_t rtn = false;
EarthVector diff = (b.GetPos() - fFirstBunch->GetPos()).Unit();
// check if bunch is on the track (defined by firstbunch pos and dir)
// and check if firstbunch is really the first one
if((diff - fFirstBunch->GetDir()).Mag() < 1.e5*TOLERANCE) rtn = true;
if((diff - fFirstBunch->GetDir()).Mag() > TOLERANCE)
Msg(EsafMsg::Debug) << "Bunch nb "<<b.GetId()<<" not perfectly along track -> "<<(diff - fFirstBunch->GetDir()).Mag() <<MsgDispatch;
return rtn;
}
//_____________________________________________________________________________
void AlongTrack_CSPropagator::Reset() {
//
// Reset method
//
fCalcul->Reset();
fFlag = false;
fFinal_medium = CLEARSKY;
SafeDelete(fFirstBunch);
}