// ESAF : Euso Simulation and Analysis Framework
// $Id: LowtranRadiativeProcessesCalculator.cc,v 1.47 2005/05/11 13:42:33 moreggia Exp $
// Sylvain Moreggia created Jan, 29 2004
#include "LowtranRadiativeProcessesCalculator.hh"
#include <vector>
#include "Config.hh"
#include "SinglePhoton.hh"
#include "BunchOfPhotons.hh"
#include "EarthVector.hh"
#include "EConst.hh"
#include "Atmosphere.hh"
#include "Ground.hh"
#include "DetectorGeometry.hh"
#include "ListPhotonsInAtmosphere.hh"
ClassImp(LowtranRadiativeProcessesCalculator)
using EConst::EarthRadius;
//____________________________________________________________________________________
//
// for lowtran interface
//
// lowtran output container
// places of declarations are of importance, don't change them
struct LowtranTransferOutput {
Float_t R,Altitude,Angle;
Float_t WaveLength[100];
Float_t Transmission[100];
Float_t Ozone[100];
Float_t Rayl[100];
Float_t Aerosol[100];
};
// Fortran interfaces
extern "C"{
#define lwtrans lwtrans_
#define lwtrans_v2 lwtrans_v2_
// calculate transmission between a point and EUSO
void lwtrans(Int_t* init, Float_t* altitude, Float_t* angle, Float_t* trans);
// calculate transmission between two points, expressed in lowtran frame
void lwtrans_v2(Int_t* init, Float_t* h1, Float_t* h2, Float_t* theta);
}
extern "C"{
#define lowtrantransfer tuple_
struct LowtranTransferOutput tuple_;
}
void InitLOWTRAN() {
// lowtran initialisation
Float_t trans[100];
Int_t flag(0);
Float_t h(0),theta(0);
lwtrans(&flag,&h,&theta,trans);
}
//________________________________________________________________________________________
LowtranRadiativeProcessesCalculator::LowtranRadiativeProcessesCalculator() : RadiativeProcessesCalculator() {
//
// ctor
//
// Configure
fStep = 0;
fFlag = false;
fInitPos.SetXYZ(0,0,HUGE);
fNbwl = 0;
fName = "lowtran";
}
//________________________________________________________________________________________
LowtranRadiativeProcessesCalculator::~LowtranRadiativeProcessesCalculator() {
//
// dtor
//
Reset();
}
//________________________________________________________________________________________
Double_t LowtranRadiativeProcessesCalculator::Trans(const SinglePhoton& s, const EarthVector& finalpos) const {
//
// Total transmission between photon position and finalpos
//
// Principle of interpolation from datacard :
// - find a reference transmission (calcuated from ground for reflected -- interpolated within a step for scattered)
// - then get the interaction length and apply it to the real travel
//
// init
Double_t rtn, trans1, trans2, coeff1, coeff2;
UInt_t Wli, Stepi;
Float_t wl = s.Wl();
Double_t OD = Atmosphere::Get()->GetClouds()->TotalOD(s.Pos(),s.Dir());
Double_t travel, entrystep, endstep, ontrack;
travel = (finalpos - s.Pos()).Mag();
const ListPhotonsInAtmosphere& track = *fTrack;
UInt_t tr_size = track.GetNbTrackSteps();
// if fToDetectorCoeff card doesn't exist
if(fToDetectorCoeff.size() == 0) {
LowtranTransferOutput output;
CallLowtran100(s.Pos(),finalpos,output);
Double_t rtn(0);
// lambda research (array is in decreasing order)
UInt_t nabove, nbelow, middle;
nabove = 101;
nbelow = 0;
while(nabove-nbelow > 1) {
middle = (nabove+nbelow)/2;
if (wl == output.WaveLength[middle-1]) Wli = middle-1;
if (wl > output.WaveLength[middle-1]) nabove = middle;
else nbelow = middle;
}
Wli = nbelow-1;
rtn = InterpolWl("linear",wl,output.WaveLength[Wli],output.WaveLength[Wli+1],output.Transmission[Wli],output.Transmission[Wli+1]);
return rtn * exp(-OD);
}
// if transmission card to detector exists :
// lambda research (array is in decreasing order)
UInt_t nabove, nbelow, middle;
nabove = 101;
nbelow = 0;
while(nabove-nbelow > 1) {
middle = (nabove+nbelow)/2;
if (wl == fWlArray[middle-1]) Wli = middle-1;
if (wl > fWlArray[middle-1]) nabove = middle;
else nbelow = middle;
}
Wli = nbelow-1;
// step research
EarthVector projec;
Double_t mag_projec = 0;
Stepi = track.TrackStepIndex(s.Pos(),projec);
mag_projec = (projec - track.GetTrackStep(0)).Mag();
// reflected photons are treated differently
// scattered photons with position further than the end of last fTrack step are treated here as reflected ones
if(s.Status() == Reflected || Stepi == (tr_size-1)) {
if(s.Status() != Reflected) Msg(EsafMsg::Debug) <<"scattered photons position further than track end" <<MsgDispatch;
if(fImpact.Z() == HUGE) Msg(EsafMsg::Warning) <<"SinglePhoton reflected but No track impact" <<MsgDispatch;
coeff1 = fToDetectorCoeff[fToDetectorCoeff.size()-1][Wli];
if(Wli == 99) coeff2 = coeff1;
else coeff2 = fToDetectorCoeff[fToDetectorCoeff.size()-1][Wli+1];
rtn = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1,coeff2); // reference transmission
endstep = (finalpos - fImpact).Mag();
return pow(rtn,travel/endstep) * exp(-OD);
}
// for backscattered photons : get reference transmission
endstep = (track.GetTrackStep(0) - track.GetTrackStep(Stepi+1)).Mag();
entrystep = (track.GetTrackStep(0) - track.GetTrackStep(Stepi)).Mag();
// interpolation in lambda for the two step bounds
coeff1 = fToDetectorCoeff[Stepi][Wli];
if(Wli == 99) coeff2 = coeff1;
else coeff2 = fToDetectorCoeff[Stepi][Wli+1];
trans1 = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1,coeff2);
//
coeff1 = fToDetectorCoeff[Stepi+1][Wli];
if(Wli == 99) coeff2 = coeff1;
else coeff2 = fToDetectorCoeff[Stepi+1][Wli+1];
trans2 = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1,coeff2);
rtn = (trans2 - trans1)/(endstep - entrystep) * (mag_projec - entrystep) + trans1; // reference transmission
ontrack = (finalpos - projec).Mag();
return pow(rtn,travel/ontrack) * exp(-OD);
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::PreProcess(const BunchOfPhotons& b) {
//
// Preprocessing for bunch transmission along the track (for ckov backscattering simulation)
// Initialization of tables (fAlongTrackCoeff) for a step by step transfer process of a given bunchofphotons
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"Bunch Preprocessing.." << MsgDispatch;
#endif
// initializations
ConfigFileParser* pConfig = Config::Get()->GetCF("RadiativeTransfer","BunchRadiativeTransfer");
fStep = pConfig->GetNum("BunchRadiativeTransfer.fStep")*km;
fNbwl = b.GetWlSpectrum().GetAxis().GetNbins();
LowtranTransferOutput output;
Double_t wl;
Int_t j;
fInitPos = b.GetPos();
CleanBuffer("alongtrack");
vector<vector<Double_t> >* Transmis;
// calculation of lowtran number of steps, and step length tunning
EarthVector pos, impact, increm;
impact = b.GetFinalImpact();
increm = impact - fInitPos;
Double_t Nb = ceil(increm.Mag()/fStep);
if(Nb == 0) Msg(EsafMsg::Warning) << "BunchPreProcess() : Nb should not be equal to zero" << MsgDispatch;
if(Nb > 1000) Nb = 1000;
fStep = increm.Mag() / Nb;
// loop makes calculations then fills each vector (one vector for one step)
for(Int_t k=0; k<Nb; k++) {
increm.SetMag((k+1)*fStep);
pos = fInitPos + increm;
Transmis = new vector<vector<Double_t> >(4);
fAlongTrackCoeff.push_back(Transmis);
CallLowtran100(fInitPos,pos,output);
Int_t lim = 99;
// fill the fAlongTrackCoeff[k]
for(Int_t i=0; i<fNbwl; i++) {
wl = b.GetWlSpectrum().GetLambdaEntry(i);
for(j=0; j<lim; j++){
if(output.WaveLength[j+1] <= wl && wl <= output.WaveLength[j]){
(*fAlongTrackCoeff[k])[0].push_back(InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Transmission[j],output.Transmission[j+1]));
(*fAlongTrackCoeff[k])[1].push_back(InterpolWl("rayleigh",wl,output.WaveLength[j],output.WaveLength[j+1],output.Rayl[j],output.Rayl[j+1]));
(*fAlongTrackCoeff[k])[2].push_back(InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Ozone[j],output.Ozone[j+1]));
(*fAlongTrackCoeff[k])[3].push_back(InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Aerosol[j],output.Aerosol[j+1]));
break;
}
}
if(j == lim) Msg(EsafMsg::Panic) << "PB in searching lambda in LowtranRadiativeProcessCalculator" << MsgDispatch;
lim = j+1; // to optimize following lambda researches
}
}
fFlag = true;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::PreProcessSingles(const Ground* ground, const EarthVector& detector, const ListPhotonsInAtmosphere& track) {
//
// Copy the light track and
// build a transmission map (fToDetectorCoeff), from positions along the track until detector
// used in BunchRadiativeTransfer::PropagationOfSingles()
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"Preprocessing Singles propagation" << MsgDispatch;
#endif
// initializations
LowtranTransferOutput output;
CleanBuffer("todetector");
fTrack = &track;
UInt_t tr_size = track.GetNbTrackSteps();
fToDetectorCoeff.reserve(tr_size + 1); // +1 for photons on ground
EarthVector start(1);
Double_t* Transmis = 0;
// loop makes calculations then fills the map
// an array for each step
for(UInt_t k=0; k < tr_size; k++) {
Transmis = new Double_t[100];
fToDetectorCoeff.push_back(Transmis);
start = track.GetTrackStep(k);
CallLowtran100(start,detector,output);
for(UInt_t i=0; i<100; i++) {
if(!k) fWlArray[i] = output.WaveLength[i];
fToDetectorCoeff[k][i] = output.Transmission[i];
}
}
// last step with init point at the track impact
Transmis = new Double_t[100];
fToDetectorCoeff.push_back(Transmis);
fImpact = ground->GetImpact(track.GetTrackStep(0),track.GetTrackStep(tr_size-1) - track.GetTrackStep(0));
start = fImpact;
if(start.Z() != HUGE) {
CallLowtran100(start,detector,output);
for(UInt_t i=0; i<100; i++) fToDetectorCoeff[tr_size][i] = output.Transmission[i];
}
else {
for(UInt_t i=0; i<100; i++) fToDetectorCoeff[tr_size][i] = 0;
Msg(EsafMsg::Info) << "PreProcessSingles: NO IMPACT on ground -> Trans=0" << MsgDispatch;
}
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::Trans(const BunchOfPhotons& b,const EarthVector& nextpos,vector<Double_t*>& coeff) const {
//
// Total transmission between bunch position and a final position
// NB: for a considered bunch, step for running lowtran is independant of step defined for bunch propagation
// (lowtran step is larger.. cf computing time consumption)
//
Int_t j,k,m;
LowtranTransferOutput output;
// if single step process
if(!fFlag) {
// lowtran called for transfer from bunch position to nextpos (one step process)
if((nextpos - b.GetNextImpact()).Mag() < TOLERANCE) {
CallLowtran100(b.GetPos(),nextpos,output);
// then fill coeff
Int_t lim = 99;
Double_t wl;
for(Int_t i=0; i<b.GetWlSpectrum().GetAxis().GetNbins(); i++) {
wl = b.GetWlSpectrum().GetLambdaEntry(i);
for(j=0; j<lim; j++){
if(output.WaveLength[j+1] <= wl && wl <= output.WaveLength[j]){
coeff[0][i] = InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Transmission[j],output.Transmission[j+1]);
coeff[1][i] = InterpolWl("rayleigh",wl,output.WaveLength[j],output.WaveLength[j+1],output.Rayl[j],output.Rayl[j+1]);
coeff[2][i] = InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Ozone[j],output.Ozone[j+1]);
coeff[3][i] = InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Aerosol[j],output.Aerosol[j+1]);
break;
}
}
if(j == lim) Msg(EsafMsg::Panic) << "PB in searching lambda in LowtranRadiativeProcessCalculator" << MsgDispatch;
lim = j+1; // to optimize following lambda researches
}
}
else {
Msg(EsafMsg::Warning) << "<PreProcess()> method has to be called before step by step transmission" << MsgDispatch;
Msg(EsafMsg::Warning) << "May rarely occur if clouds exist" <<MsgDispatch;
}
}
// return transmission between bunch position and nextpos, using interpolation
for(k=0; k<4; k++) {
for(m=0; m<fNbwl; m++)
coeff[k][m] = Interpolate(k,m,nextpos)/Interpolate(k,m,b.GetPos());
}
}
//________________________________________________________________________________________
Double_t LowtranRadiativeProcessesCalculator::Interpolate(Int_t i,Int_t j,const EarthVector& pos) const {
//
// Interpolate the specified transmission coefficient fAlongTrackCoeff[][i][j]
// to get the value at pos (on the transfer track)
//
Double_t rtn, Ti, Tf;
Double_t mag = (pos - fInitPos).Mag();
if(mag/fStep == Double_t(fAlongTrackCoeff.size())) return (*fAlongTrackCoeff[fAlongTrackCoeff.size() - 1])[i][j];
if(mag > fAlongTrackCoeff.size()*fStep) {
if(fabs(mag - fAlongTrackCoeff.size()*fStep) > 1*mm) {
#ifdef DEBUG
//Msg(EsafMsg::Debug) << "<Interpolate> mag is too large -> " <<fabs(mag - fAlongTrackCoeff.size()*fStep) << MsgDispatch;
if(fabs(mag - fAlongTrackCoeff.size()*fStep) > 5e4*mm)
Msg(EsafMsg::Warning) << "<Interpolate> mag is TOO large (>50m) -> " <<fabs(mag - fAlongTrackCoeff.size()*fStep) << MsgDispatch;
#endif
mag = fAlongTrackCoeff.size()*fStep;
}
else return (*fAlongTrackCoeff[fAlongTrackCoeff.size() - 1])[i][j];
}
Double_t iter = floor(mag/fStep);
UInt_t n = UInt_t(iter);
if(n > (fAlongTrackCoeff.size() - 1)) n = fAlongTrackCoeff.size() - 1;
if(n == 0) {
Ti = 1;
Tf = (*fAlongTrackCoeff[n])[i][j];
}
else {
Ti = (*fAlongTrackCoeff[n-1])[i][j];
Tf = (*fAlongTrackCoeff[n])[i][j];
}
// exponential interpolation
rtn = Ti*exp((mag - n*fStep)/fStep * log(Tf/Ti));
// linear interpolation
//rtn = Ti + (Tf - Ti) * (mag - n*fStep)/fStep;
return rtn;
}
//________________________________________________________________________________________
EarthVector LowtranRadiativeProcessesCalculator::InLowtranFrame(const EarthVector& v1,const EarthVector& v2) const {
//
// Change of frame to match lowtran (v1 is the reference vector)
// the frame moves, then v2 expression in the new frame is returned
//
EarthVector rtn(v2);
EarthVector new1(v1);
EarthVector vrot;
if(v1.Theta() == 0) return rtn;
// translation -> frame is earth-centered
rtn(2) += EarthRadius();
new1(2) += EarthRadius();
// rotation defined to get earth-centered frame z-axis aligned with position new1
EarthVector Uz(0,0,1);
vrot = new1.Cross(Uz);
rtn.Rotate(new1.Theta(),vrot);
rtn(2) -= EarthRadius();
return rtn;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::CallLowtran100(const EarthVector& pos, const EarthVector& nextpos, LowtranTransferOutput& output) const {
//
// Call lowtran (100 wavelengths processed) to calculate transmission between pos and nextpos, then fill output
//
EarthVector pos1(1);
EarthVector pos2(1);
pos1 = InLowtranFrame(pos,pos);
pos2 = InLowtranFrame(pos,nextpos);
if(pos1.Zv() > pos2.Zv()) { // to get h1 < h2 (then good theta)
pos1 = InLowtranFrame(nextpos,nextpos);
pos2 = InLowtranFrame(nextpos,pos);
}
Float_t h1 = (Float_t)pos1.Zv()/km;
Float_t h2 = (Float_t)pos2.Zv()/km;
h1 = floor(h1*1000 + 0.5)/1000; // to avoid lowtran round-off bugs
h2 = floor(h2*1000 + 0.5)/1000;
Float_t theta = (Float_t)(pos2 - pos1).Theta()/deg;
Int_t flag(1);
lwtrans_v2(&flag,&h1,&h2,&theta);
// copy lowtran tuple_ in output
output.R = lowtrantransfer.R * km;
output.Altitude = lowtrantransfer.Altitude * km;
output.Angle = lowtrantransfer.Angle;
for(Int_t i=0; i<100; i++) {
output.WaveLength[i] = lowtrantransfer.WaveLength[i] * micrometer;
output.Transmission[i] = lowtrantransfer.Transmission[i];
output.Rayl[i] = lowtrantransfer.Rayl[i];
output.Ozone[i] = lowtrantransfer.Ozone[i];
output.Aerosol[i] = lowtrantransfer.Aerosol[i];
}
}
//________________________________________________________________________________________
Double_t LowtranRadiativeProcessesCalculator::InterpolWl(string type,Double_t w,Double_t w1,Double_t w2,Double_t T1,Double_t T2) const {
//
// Interpolation of transmission coefficients as function of wavelength
//
Double_t rtn(0);
if(type == "rayleigh") {
Double_t expon = (w*w*w*w - w1*w1*w1*w1) /
(w2*w2*w2*w2 - w1*w1*w1*w1)* w2*w2*w2*w2;
rtn = T1*exp(expon * log(T2/T1));
}
else if(type == "linear") {
rtn = T1 + (w - w1)/(w2 - w1) * (T2 - T1);
}
else Msg(EsafMsg::Panic) << "Unvalid type of interpolation (LowtranRadiativeProcessesCalculator)" << MsgDispatch;
return rtn;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::CleanBuffer(string DM) {
//
// Clean the indicated internal trans map (fAlong.. or fToDetect..)
//
if(DM == "alongtrack") {
if(fAlongTrackCoeff.size()) {
for(size_t i=0; i<fAlongTrackCoeff.size(); i++) {
SafeDelete(fAlongTrackCoeff[i]);
}
fAlongTrackCoeff.clear();
}
}
else if(DM == "todetector") {
if(fToDetectorCoeff.size()) {
for(size_t i=0; i<fToDetectorCoeff.size(); i++) {
if(fToDetectorCoeff[i]) delete[] fToDetectorCoeff[i];
fToDetectorCoeff[i] = 0;
}
fToDetectorCoeff.clear();
}
}
else Msg(EsafMsg::Warning) << "<CleanBuffer> Wrong argument -> POTENTIAL MEMORY LEAK" << MsgDispatch;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::Reset() {
//
// get ready for next event
//
CleanBuffer("alongtrack");
CleanBuffer("todetector");
fFlag = false;
fInitPos.SetXYZ(0,0,HUGE);
fStep = 0;
fNbwl = 0;
}