// ESAF : Euso Simulation and Analysis Framework
// $Id: FastRadiativeProcessesCalculator.cc,v 1.23 2005/10/31 13:01:09 moreggia Exp $
// Sylvain Moreggia created Jan, 29 2004
#include "FastRadiativeProcessesCalculator.hh"
#include "EarthVector.hh"
#include "Atmosphere.hh"
#include <math.h>
#include "BunchOfPhotons.hh"
#include "SinglePhoton.hh"
#include "Config.hh"
#include "EConst.hh"
#include "EsafRandom.hh"
ClassImp(FastRadiativeProcessesCalculator)
using EConst::Avogadro;
using namespace sou;
//_____________________________________________________________________________
FastRadiativeProcessesCalculator::FastRadiativeProcessesCalculator() : RadiativeProcessesCalculator() {
//
// ctor
//
Msg(EsafMsg::Info) <<"*** FastRadiativeProcessesCalculator built ***" << MsgDispatch;
fName = "fast";
}
//_____________________________________________________________________________
FastRadiativeProcessesCalculator::~FastRadiativeProcessesCalculator() {
//
// dtor
//
}
//_____________________________________________________________________________
Double_t FastRadiativeProcessesCalculator::Trans(const SinglePhoton& p, const EarthVector& finalpos, Double_t* dummy) const {
//
// total Transmission between SinglePhoton position and finalpos
//
// only rayleigh
const Atmosphere* atmo = Atmosphere::Get();
Double_t X = atmo->Grammage(p.Pos(),finalpos,"pos",50*km);
//Double_t ro_0 = atmo->Air_Density(0.);
//Double_t norm = 8*pow(pi,3)*4/3;
//norm *= (6 + 3*0.0279) / (6 - 7*0.0279);
//norm *= pow(2.83e-4,2)*28.966/Avogadro();
//Double_t To = exp( -(norm*X/(gram/cm2)) / (pow(p.Wl()/cm,4)*pow(ro_0/(gram/cm3),2)) );
Double_t wlterm = pow(p.Wl()/cm,-4) / (1 - 1.08426/9.38076 * 1e-9 * pow(p.Wl()/cm,-2));
Double_t To = exp( -8.24689e-22*X*cm2/gram * wlterm);
Double_t OD = atmo->GetClouds()->TotalOD(p.Pos(),p.Dir());
if(dummy) {
dummy[0] = To * exp(-OD);
dummy[1] = To;
dummy[2] = 1.;
dummy[3] = 1.;
}
return To * exp(-OD);
//return exp(-X/(2970*gram/cm2)*pow(400*nm/p.Wl(),4)); //DELETE these kevin Reil, HiRes
}
//_____________________________________________________________________________
Double_t FastRadiativeProcessesCalculator::RandomGrammage_RaylScat(Double_t Wl) const {
//
// sample the rayleigh exponential distribution in grammage -> get grammage at which scattering occurs
//
TRandom* rndm = EsafRandom::Get();
Double_t wlterm = pow(Wl/cm,-4) / (1 - 1.08426/9.38076 * 1e-9 * pow(Wl/cm,-2));
Double_t grammage = log(1 - rndm->Rndm()) / (-8.24689e-22*cm2/gram * wlterm);
return grammage; // in esaf units
}
//_____________________________________________________________________________
void FastRadiativeProcessesCalculator::Trans(const BunchOfPhotons& b,const EarthVector& nextpos,vector<Double_t*>& coeff) const {
//
// total transmission between bunch and a final position.
//
// Only analytical rayleigh
Double_t wl, wlterm;
const Atmosphere* atmo = Atmosphere::Get();
//Double_t ro_0 = atmo->Air_Density(0.);
Double_t X = atmo->Grammage(b.GetPos(),nextpos);
Int_t nb = b.GetWlSpectrum().GetAxis().GetNbins();
for(Int_t i=0; i<nb; i++) {
wl = b.GetWlSpectrum().GetLambdaEntry(i);
wlterm = pow(wl/cm,-4) / (1 - 1.08426/9.38076 * 1e-9 * pow(wl/cm,-2));
//norm = 8*pow(pi,3)*4/3;
//norm *= (6 + 3*0.0279) / (6 - 7*0.0279);
//norm *= pow(2.83e-4,2)*28.966/Avogadro();
//coeff[0][i] = exp( -(norm*X/(gram/cm2)) / (pow(wl/cm,4)*pow(ro_0/(gram/cm3),2)) );
coeff[0][i] = exp( -8.24689e-22*X*cm2/gram * wlterm);
coeff[1][i] = coeff[0][i];
coeff[2][i] = 1;
coeff[3][i] = 1;
}
}