Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

RadiativeProcessesCalculator - source file

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