Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TestLightSource - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TestLightSource.cc,v 1.31 2005/04/14 12:44:17 moreggia Exp $
// Anne Stutz created Dec,  2 2003
//
// This class brings test photons in atmosphere to Euso
// The type of events supported are the following:
//     SPOT  light from a single point at the same time
//     RANDSPOT light from a random point at the same time
//     TRACK light making a track in atmosphere
//     RANDTRACK 
//
// The parameters needed are the following :
//     TestLightSource.Type         see above can be SPOT or TRACK
//     TestLightSource.Photons      number of photons generated per event (if<1000 SinglePhoton else BunchOfPhotons)
//     TestLightSource.Theta1       lower zenith angle                                 
//     TestLightSource.Theta2       higher zenith angle                                 
//     TestLightSource.Phi1         lower azimuth angle
//     TestLightSource.Phi2         higher azimuth angle
//     TestLightSource.FirstPointX  X first point source in km 
//     TestLightSource.FirstPointY  Y first point source in km 
//     TestLightSource.FirstPointZ  Z first point source in km 
//     TestLightSource.LongExtension = 0          longitudinal extension of the bunch in g/cm2 can be 0
//     TestLightSource.LatDistribution = NKG2      lateral distribution of the bunch can be NKG NULL
//     TestLightSource.AngDistribution = baltru   angular distribution of the bunch can be baltru NULL
//     TestLightSource.SpectrumType   can be FLUO,MONO,CERENKOV
//     TestLightSource.Lambda         Wavelenght if MONO
//     TestLightSource.DirectionType  can be EUSO,ISOTROPIC,MONO

#include "TestLightSource.hh"
#include "ListPhotonsInAtmosphere.hh"
#include "EsafRandom.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "FluoCalculator.hh"
#include "LightSourceFactory.hh"
#include "TFormula.h"
#include "TMath.h"
#include "EConst.hh"
#include "EarthVector.hh"
#include "RadiativeFactory.hh"
#include "Ground.hh"
#include "Atmosphere.hh"

ClassImp(TestLightSource)

using EConst::C;

 /* dNe/dr derived from 'standard' NGK formula for Lateral distribution. r is the distance to shower axis 
 Rm=moliere radius is a parameter . The integration of NGK2 over r from 0. to infinity is normalized to 1 . The integral
 from radius r1 to radius r2 give the fraction of total number of electrons which lie inside such interval.*/
Double_t TestLightSourceNKG2(Double_t *x, Double_t *par) {
  Double_t R = x[0], age = x[1], Rm=par[0];
  Double_t e1=2.0;
  Double_t e2=4.5;
  Double_t D=R/Rm;
  return TMath::Gamma(e2-age)/TMath::Gamma(age)/TMath::Gamma(e2-2.*age)*TMath::Power(D,age-(e1-1.0))
        *TMath::Power((1.0+D),age-e2)/Rm;
}

/* dNe/dtheta(theta,Et) from Baltrusaitis et al. J.Phys.G:Nucl. Phys. 13 (1987) where theta is the
angle between the electrons and the shower axis and Et (MeV) is the energy thr for electrons considered.
The integration of dNe/dtheta(theta,Et) over dtheta from 0 to pi is normalized to 1. 
The integral from theta1 to theta2 give the fraction of total number of electrons which lie inside 
such angular interval (distribution in phi is supposed uniform)
.*/
Double_t TestLightSourceBaltru(Double_t *x, Double_t *par) {
  Double_t theta = x[0], Et = x[1];
  Double_t a=0.85;
  Double_t b=-0.66;
  Double_t theta0=a*TMath::Power(Et,-b);
  return TMath::Exp(-theta/theta0)/theta0/(1.-TMath::Exp(-180/theta0) );
}

//_________________________________________________________________________________________________
TestLightSource::TestLightSource() : LightSource("TEST"), EsafMsgSource(), fFluocalcul(0) {
    // ctor

    Msg(EsafMsg::Info) << "Start Building TestLightSource" << MsgDispatch;

    fPh_in_atmo = new ListPhotonsInAtmosphere;
    if(!fPh_in_atmo) Msg(EsafMsg::Panic) << "Pb for memory allocation of fPh_in_atmo" << MsgDispatch;
  
    Configure();
    Msg( EsafMsg::Info) << "TestLightSource built"<<MsgDispatch;
}

//_________________________________________________________________________________________________
TestLightSource::~TestLightSource() {
    //
    // dtor 
    //
    SafeDelete(fPh_in_atmo);
    SafeDelete(fFluocalcul);
    SafeDelete(fLateralDistribution);
    SafeDelete(fAngularDistribution);
}

//________________________________________________________________________
 void TestLightSource::Configure() {
    //
    // Configure TestLightSource
    //

    // Get fluorescence calculator
    string fluoname = Conf()->GetStr("TestLightSource.FluoCalculator");
    fFluocalcul = LightSourceFactory::Get()->GetFluoCalculator( fluoname );
    Msg(EsafMsg::Info) << "Fluo calculator name " <<fFluocalcul->GetName()<< MsgDispatch;

    // Lateral and Angular distributions
    fLateralDistribution = NULL;
    fAngularDistribution = NULL;
    string LateralDistributionName = Conf()->GetStr("TestLightSource.LatDistribution");
    string AngularDistributionName = Conf()->GetStr("TestLightSource.AngDistribution");
 
    if (LateralDistributionName == "NKG2" ) 
       fLateralDistribution = new TF2("LatDist",TestLightSourceNKG2,0.001,5000,0,2,1);
    if (AngularDistributionName == "baltru" )
       fAngularDistribution = new TF2("AngDist",TestLightSourceBaltru,0.,TMath::Pi(),.5,1000.);
}

//_________________________________________________________________________________________________
 void TestLightSource::Reset() {
    //
    // reset internal list of photons
    //

    if(fPh_in_atmo) fPh_in_atmo->Reset();
    if(fFluocalcul) fFluocalcul->Reset();
}

//_________________________________________________________________________________________________
 PhotonsInAtmosphere *TestLightSource::Get( const PhysicsData *dummy) {
    // 
    // return the list of photons in atmosphere
    //
    Reset();
    string type = Conf()->GetStr("TestLightSource.Type");
    Double_t nbPhotons = Conf()->GetNum("TestLightSource.Photons");
    Double_t fFirstPointX = (Conf()->GetNum("TestLightSource.FirstPointX"))*km;
    Double_t fFirstPointY = (Conf()->GetNum("TestLightSource.FirstPointY"))*km;
    Double_t fFirstPointZ = (Conf()->GetNum("TestLightSource.FirstPointZ"))*km;

    TRandom* rndm = EsafRandom::Get();

    EarthVector Posi;
    if (fFirstPointZ > 100*km ) {
        Msg(EsafMsg::Warning) << "In TestLightSource FirstPointZ set at 100 km" << MsgDispatch; 
        fFirstPointZ = 100*km;
    }
    Posi.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);

    // spot in fixed position
    if ( type == "SPOT" ) {
       MakeSpot(Posi,nbPhotons);
    }

    else if ( type == "RANDSPOT" ) {    // spot in random position // FIXME not uniform
       return NULL;
    }

   // track in fixed position
    else if ( type == "TRACK" ) {    
       // random direction between theta1 and theta2, and between phi1 and phi2
       Double_t theta1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
       Double_t theta2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
       Double_t phi1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
       Double_t phi2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi2");

       Double_t thmax = theta2;
       Double_t thmin = theta1;
       if (theta1>theta2) {
            thmax = thmin;
            thmin = theta2;
       }
       Double_t theta = thmin + rndm->Rndm()*(thmax-thmin);

       Double_t phmax = phi2;
       Double_t phmin = phi1;
       if (phi1>phi2) {
            phmax = phmin;
            phmin = phi2;
       }
       Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);

       MakeTrack(theta,phi,Posi,nbPhotons); 
    }

    else if ( type == "RANDTRACK" ) {   // track in random position
        return NULL;
    }
    else {
        Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.Type = " << type <<MsgDispatch;
        return (PhotonsInAtmosphere*)0;
    }
    return fPh_in_atmo;
}

//___________________________________________________________________________________________
 PhotonsInAtmosphere *TestLightSource::MakeSpot(const EarthVector& Posi, Double_t nbPhotons) {
    // Produce light in a single spot 

    Msg(EsafMsg::Info)<<"TestLightSource SPOT CALLED at (km) "<<Posi.X()/km<<" "
		      <<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;

    string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
    Double_t theta1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
    Double_t theta2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
    Double_t phi1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
    Double_t phi2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi2");

    Double_t thmax = theta2;
    Double_t thmin = theta1;
    Double_t phmax = phi2;
    Double_t phmin = phi1;
    if (theta1>theta2) {
       thmax = thmin;
       thmin = theta2;
    }
    if (phi1>phi2) {
        phmax = phmin;
        phmin = phi2;
    }

    if (directiontype == "MONO") {
       Msg(EsafMsg::Info)<<"                with theta between "<<thmin<<" and "<<thmax<<MsgDispatch;
       Msg(EsafMsg::Info)<<"                with phi between   "<<phmin<<" and "<<phmax<<MsgDispatch;
    }
    else if (directiontype == "ISOTROPIC") {
       Msg(EsafMsg::Info)<<"                with isotropic direction" << MsgDispatch;
    }
    else if (directiontype == "EUSO") {
       Msg(EsafMsg::Info)<<"                direct to EUSO" << MsgDispatch;
    }
    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.DirectionType = " << directiontype <<MsgDispatch;

    // wavelenght spectrum
    EsafSpectrum spectrum(357*nm);
    Double_t TotalYield = MakeSpectrum(Posi,&spectrum);

    //generation of SinglePhoton or BunchOfPhotons
    TRandom* rndm = EsafRandom::Get();

    if (nbPhotons<1000) {
       // only single photons direct to Euso
       if (directiontype == "ISOTROPIC" || directiontype == "MONO") 
           Msg(EsafMsg::Info)<<"In TestLightSource SinglePhoton are only emitted direct to Euso"<< MsgDispatch;
        
       for (int i=0; i<(int)nbPhotons; i++) {          
           Double_t wl = spectrum.GetLambda();
           EarthVector dir(0,0,0);
           dir = EUSO()-Posi;
           SinglePhoton *p = new SinglePhoton(0,wl,Posi,dir,Fluo,Direct);
           fPh_in_atmo->Add(p);
       }
    } 

    else {
    // one bunch of nbPhotons
        EarthVector dir(1);
        dir.SetTheta(thmin + rndm->Rndm()*(thmax-thmin));
        dir.SetPhi(phmin + rndm->Rndm()*(phmax-phmin));

        EarthVector Posf = GetLongitudinalExtension(Posi,dir);
        Double_t tf = (Posf-Posi).Mag()/C();
        PhotonType type=Fluo;
        if (directiontype == "MONO") type = Cerenkov;

        BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,type);

        if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
        if( fAngularDistribution ) b->SetParentAngular(GetAngularDistribution());

        fPh_in_atmo->Add(b);
    }
    return fPh_in_atmo; 
}

//_____________________________________________________________________________________________________
 PhotonsInAtmosphere *TestLightSource::MakeTrack(const Double_t& theta, const Double_t& phi, const EarthVector& Posi, Double_t nbPhotons) {
     // 
    // Produce light along a track 
    //

    Msg(EsafMsg::Info)<<"TestLightSource TRACK CALLED theta="<<theta<<" phi ="<<phi<<MsgDispatch;
    Msg(EsafMsg::Info)<<"               Starting point at (km) "<<Posi.X()/km<<" "<<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
    string directiontype = Conf()->GetStr("TestLightSource.DirectionType");

    EarthVector DirTrack(1);
    DirTrack.SetTheta(theta);
    DirTrack.SetPhi(phi);

    EsafSpectrum spectrum(357*nm);
    Double_t TotalYield;

    //generation of SinglePhoton or BunchOfPhotons
    // impact on ground
    Ground* ground = RadiativeFactory::Get()->GetGround();
    EarthVector impact = ground->GetImpact(Posi,DirTrack);
    Msg(EsafMsg::Info)<<"               Impact on ground at "<<impact.X()/km<<" "<<impact.Y()/km<<" "<<impact.Z()/km<< MsgDispatch;
    Double_t magmax = (Posi - (ground->GetImpact(Posi,DirTrack))).Mag();

    // only single photons
    if (nbPhotons<1000) {
        Msg(EsafMsg::Info)<<"                direct to EUSO" << MsgDispatch;
        if (directiontype == "ISOTROPIC" || directiontype == "MONO") 
           Msg(EsafMsg::Info)<<"In TestLightSource SinglePhoton are only emitted direct to Euso"<< MsgDispatch;
        
        for (int i=0; i<(int)nbPhotons; i++) {
           Double_t mag = magmax*(i+1)/nbPhotons;
           EarthVector pos = Posi + mag*DirTrack;
           Double_t t = mag/C();

           TotalYield = MakeSpectrum(pos,&spectrum);

           Double_t wl = spectrum.GetLambda();
           EarthVector dir = EUSO()-pos;
           SinglePhoton *p = new SinglePhoton(t,wl,pos,dir,Fluo,Direct); 
  
           fPh_in_atmo->Add(p);
        }
    }
    // 100 bunchs of photons
    else {  
        if (directiontype == "EUSO")
           Msg(EsafMsg::Info)<<"In TestLightSource BunchPhoton are not emitted directly to Euso but isotropically"<< MsgDispatch;

        Double_t nbPhotonsInBunch = nbPhotons/100.;
        PhotonType type=Fluo;
        if (directiontype == "MONO") type = Cerenkov;

        for (Float_t i=0; i<100; i++) {
	    Double_t mag = magmax*(i+1)/100;
            EarthVector pos = Posi + mag*DirTrack;
            Double_t t = mag/C();

            TotalYield = MakeSpectrum(pos,&spectrum);

            EarthVector posf = GetLongitudinalExtension(pos,DirTrack);
            if ( posf == pos ) break;
            Double_t tf = t + (posf-pos).Mag()/C();

	    BunchOfPhotons*  b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,pos,posf,t,tf,spectrum,DirTrack,type);

            if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
	    if( fAngularDistribution ) b->SetParentAngular(GetAngularDistribution());

           fPh_in_atmo->Add(b);
       }
    }
    return fPh_in_atmo; 
}

//_________________________________________________________________________________________________
 Double_t TestLightSource::MakeSpectrum(const EarthVector& pos,EsafSpectrum* spectrum) {
// Calculation of the wavelenght spectrum

    Double_t TotalYield=0;
    string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");

    if (spectrumtype == "MONO") {
       double lambda = (Conf()->GetNum("TestLightSource.Lambda"))*nm;
       if ( spectrum!=0) spectrum -> Reset(lambda);
       TotalYield = 1.;
    }
    else if (spectrumtype == "FLUO") {
       Double_t energy=80.*MeV;
       TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),energy,spectrum);
    }
    else if (spectrumtype == "CERENKOV") {
         TFormula cerenkov("cerenkov","1 /(x*x)");
         if ( spectrum !=0 ) spectrum -> Reset(&cerenkov,150,300*nm,450*nm);
	 TotalYield = 1.;
    }
    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.SpectrumType = " << spectrumtype <<MsgDispatch;

    if(!spectrum) Msg(EsafMsg::Panic)<<"Pb of memory allocation in TestLightSource"<<MsgDispatch;
    return TotalYield;
}


//_________________________________________________________________________________________________
 const TF12 TestLightSource::GetLateralDistribution(const EarthVector& pos) {
    //
    // return TF12 the lateral distribution at pos for age = 1
    //   

    TF12 LatDist;
    Double_t age = 1.;

    // if NKG in meter gives moliere radius as parameter in meter
    if ( fLateralDistribution) {
        if ( fLateralDistribution->GetNpar() == 1 ) { 
           // Get Atmosphere for Density
           const Atmosphere* atmo = Atmosphere::Get();
           Double_t Rm = 9.6 * gram/cm2 / atmo->Air_Density( pos.Zv() ) / m; // in meters
	   Msg(EsafMsg::Debug)<<"rayon moliere =  " << Rm <<MsgDispatch;
           fLateralDistribution->SetParameters(Rm,0);
        }    
        else Msg(EsafMsg::Panic)<<"Lateral Distribution should have 1 parameter" <<MsgDispatch;
        LatDist = TF12("LDS_name",fLateralDistribution,age,"X"); 
    }
    else Msg(EsafMsg::Panic)<<"No Lateral Distribution" <<MsgDispatch;

    return LatDist;
}

//_________________________________________________________________________________________________
 const TF12 TestLightSource::GetAngularDistribution() {
//
// return Pointer to the angular distribution 
//

   TF12 AngDist;

   Double_t EthCer = 20.;
   if ( fAngularDistribution ) AngDist = TF12("ADS_name",fAngularDistribution,EthCer,"X");
   else Msg(EsafMsg::Panic)<<"No Angular Distribution "<<MsgDispatch;

   return AngDist;
}

//__________________________________________________________________________________________________
 EarthVector TestLightSource::GetLongitudinalExtension(const EarthVector& Pos,const EarthVector& Dir) {
//
// return last point of a bunch with first point Pos and mean direction Dir
//
    Double_t LgExt = (Conf()->GetNum("TestLightSource.LongExtension"))*g/cm2;
    Double_t L,h(0),depth(0);
    EarthVector Posf = Pos;

    L = LgExt / Atmosphere::Get()->Air_Density(h)/100.;
    Int_t cycle=0;
    while(1) {
        Posf += Dir*L;
        if ( Posf.IsUnderSeaLevel() ) return Pos;
        depth += Atmosphere::Get()->Air_Density(Posf.Zv())*L;
        if ( depth >= LgExt ) break;
        if ( cycle>100000 ) {
	  MsgForm( EsafMsg::Debug,"Next point not found : cycle > 100000 with Lgext %f",LgExt/g*cm2 );
          return Pos;
	}
    }
    return Posf;
}

//__________________________________________________________________________________________________
 EarthVector TestLightSource::EUSO() const {
    //
    // TOOL FUNCTION : Returns the EUSO position
    //
    static Int_t counter = 0;
    static EarthVector pos(0,0,0);
    if(!counter) {
        ConfigFileParser* pConf = Config::Get()->GetCF("General","Euso");
        pos.SetZ(pConf->GetNum("Euso.fAltitude")*km);
    }
    if(!counter) counter++;
    return pos;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31