// ESAF : Euso Simulation and Analysis Framework
// $Id: RadiativeProcessesCalculator.cc,v 1.12 2005/10/06 08:09:55 moreggia Exp $
// Sylvain Moreggia created Jan, 9 2004
#include "RadiativeProcessesCalculator.hh"
#include "EarthVector.hh"
#include "Atmosphere.hh"
#include <math.h>
#include "BunchOfPhotons.hh"
#include "SinglePhoton.hh"
#include "EsafMsgSource.hh"
#include "EsafRandom.hh"
#include <TF1.h>
ClassImp(RadiativeProcessesCalculator)
using namespace TMath;
//
// C-function used for TF1 root objects definition
//
////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////
//______________________________________________________________________________
Double_t RadiativeCalculator_HG_Mie_phase_function(Double_t *x, Double_t *par) {
//
// cos(theta) distribution for rayleigh phase function
// it is dP / d(theta) = 2pi * Phase_function(theta) * sin(theta) (integrated over phi)
// in which the variable change x = cos(theta) is done
// dx = sin(theta)d(theta) is done
// its integration over theta within [-1,1] is normalized to 1
//
Double_t g = par[0];
return 0.5 * (1 - g*g) / pow(1 - 2*g*x[0] + g*g,1.5);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////
//_____________________________________________________________________________
RadiativeProcessesCalculator::RadiativeProcessesCalculator() : EsafMsgSource() {
//
// ctor
//
}
//_____________________________________________________________________________
RadiativeProcessesCalculator::~RadiativeProcessesCalculator() {
//
// dtor
//
}
//_____________________________________________________________________________
void RadiativeProcessesCalculator::RandomDir(string type,EarthVector& dir, Double_t g) const {
//
// 'dir' is the incoming direction -> this method set dir to the scattering outgoing direction
// get a random direction, sampling relevant scattering phase functions
//
Double_t theta(0);
Double_t u1, u2;
TRandom* rndm = EsafRandom::Get();
// get random theta, assuming isotropic distribution in phi
// theta is the angle between incoming and outgoing directions
if(type == "rayleigh") {
// von neumann sampling -> much faster than using TF1
do {
u1 = -1. + rndm->Rndm()*2.;
u2 = rndm->Rndm() * 3./4.;
} while(u2 >= (3./8. * (1. + u1*u1)));
theta = acos(u1);
}
else if(type == "HG_mie") {
Msg(EsafMsg::Panic) << "<RandomTheta> Mie scattering not implemented yet"<<MsgDispatch;
if(g == -10) Msg(EsafMsg::Panic) << "<RandomTheta> Asymmetry parameter should be defined" <<MsgDispatch;
}
else Msg(EsafMsg::Panic) << "<RandomTheta> Wrong type of distribution : "<<type<<MsgDispatch;
// then get a random direction
EarthVector v = dir.Orthogonal();
v.Rotate(rndm->Rndm()*2.*Pi(),dir);
EarthVector axis = v.Cross(dir);
dir.Rotate(theta,axis);
}
//_____________________________________________________________________________
Double_t RadiativeProcessesCalculator::RayleighPhaseFunction(const EarthVector& bunchdir,const EarthVector& dir) const {
//
// rayleigh scattering phase function - Normalized when integrated over full solid angle
//
Double_t rtn;
Double_t theta = dir.Angle(bunchdir);
rtn = 3 * (1. + cos(theta)*cos(theta)) / (16.*Pi());
return rtn;
}
//_____________________________________________________________________________
Double_t RadiativeProcessesCalculator::Mie_HG_PhaseFunction(const EarthVector& bunchdir,const EarthVector& dir,Double_t g) const {
//
// Mie scattering phase function - Henyey-Greenstein parametrization
// Normalized when integrated over full solid angle
//
Double_t rtn;
Double_t theta = dir.Angle(bunchdir);
rtn = 1./(4*Pi()) * (1 - g*g) / pow(1 - 2*g*cos(theta) + g*g,1.5);
return rtn;
}