Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

FastRadiativeProcessesCalculator - source file

// 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;
    }
}





About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005