// ESAF : Euso Simulation and Analysis Framework
// $Id: FastRadiativeProcessesCalculator.cc,v 1.18 2005/04/14 15:57:46 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"
ClassImp(FastRadiativeProcessesCalculator)
using EConst::Avogadro;
//_____________________________________________________________________________
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) const {
//
// total Transmission between SinglePhoton position and finalpos
//
// only rayleigh
const Atmosphere* atmo = Atmosphere::Get();
Double_t X = atmo->Grammage(p.Pos(),finalpos);
//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());
return To * exp(-OD);
//return exp(-X/(2970*gram/cm2)*pow(400*nm/p.Wl(),4)); //DEL these kevin Reil, HiRes
}
//_____________________________________________________________________________
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;
}
}