Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TestLightSource - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TestLightSource.cc,v 1.40 2005/10/27 13:54:35 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.Option         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 namespace TMath;
using namespace sou;
using namespace EConst;

 /* 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 Gamma(e2-age)/Gamma(age)/Gamma(e2-2.*age)*Power(D,age-(e1-1.0))
        *Power((1.0+D),age-e2)/Rm;
}

Double_t TestLightSourceBaltru(Double_t *x, Double_t *par) {
    //
    // 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 theta = x[0], Et = x[1];
    Double_t a = 0.85;
    Double_t b = 0.66;
    Double_t theta0 = a*Power(Et,-b);
    Double_t pigreco = Pi();
    
    return Exp(-theta/theta0) / theta0 / (1.- Exp(-pigreco/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 detector position
    ConfigFileParser* pConf = Config::Get()->GetCF("General","Euso");
    fEUSO.SetZ(pConf->GetNum("Euso.fAltitude")*km);
    
    // 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.,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 option = Conf()->GetStr("TestLightSource.Option");
    Double_t nbPhotons = Conf()->GetNum("TestLightSource.Nbph");
    Double_t fFirstPointX = (Conf()->GetNum("TestLightSource.FirstPointX"))*km;
    Double_t fFirstPointY = (Conf()->GetNum("TestLightSource.FirstPointY"))*km;
    Double_t fFirstPointZ = (Conf()->GetNum("TestLightSource.FirstPointZ"))*km;
    
    string impactMode = Conf()->GetStr("TestLightSource.ImpactMode");
    Double_t ImpactXmin = Conf()->GetNum("TestLightSource.ImpactXmin")*km;
    Double_t ImpactXmax = Conf()->GetNum("TestLightSource.ImpactXmax")*km;
    Double_t ImpactYmin = Conf()->GetNum("TestLightSource.ImpactYmin")*km;
    Double_t ImpactYmax = Conf()->GetNum("TestLightSource.ImpactYmax")*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 ( option == "SPOT" ) MakeSpot(Posi,nbPhotons);

    // spot mimicking shower maximum, depends on theta
    else if ( option == "HmaxSPOT" ) {
	Double_t x = rndm->Rndm()*(ImpactXmax - ImpactXmin) + ImpactXmin;
	Double_t y = rndm->Rndm()*(ImpactYmax - ImpactYmin) + ImpactYmin;
	Double_t z = - (x*x + y*y)/(2*EarthRadius());   // approched value
	MakeHmaxSpot(EarthVector(x,y,z),nbPhotons);
    }

    // track in fixed position
    else if ( option == "TRACK" ) {    
       // random direction between theta1 and theta2, and between phi1 and phi2
       Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
       Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
       Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
       Double_t phi2 = 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 {
        Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.Option = " << option <<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");
    string photontype = Conf()->GetStr("TestLightSource.PhotonType");
    string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
    Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
    Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
    Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
    Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
    PhotonType phtype = Fluo;
    if(photontype == "cerenkov") phtype = Cerenkov;
    
    TRandom* rndm = EsafRandom::Get();

    string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
    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;
    }
    Double_t theta = thmin + rndm->Rndm()*(thmax-thmin);
    Double_t phi   = phmin + rndm->Rndm()*(phmax-phmin);
    EarthVector Omega(1);
    if(thetaMode == "local") {
        EarthVector v1 = Posi.Unit();
        EarthVector vrot;
        Omega.SetXYZ(1,0,0);
        EarthVector Uz(0,0,1);
        vrot = Uz.Cross(v1);
        if(vrot.Mag() > 0) {
            Omega.Rotate(v1.Theta(),vrot);
            Omega.Rotate(phi,v1);
            vrot = v1.Cross(Omega);
            Omega.Rotate(Pi()/2+theta,vrot);
        }
	else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
	theta = Omega.Theta();
	phi   = Omega.Phi();
    }
    else if(thetaMode == "mes") {
        theta = Pi() - theta;
	phi += Pi();
    }
    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;

    
    if (directiontype == "UNIQUE") {
       Msg(EsafMsg::Info)<<"                with theta between "<<thmin/deg<<" and "<<thmax/deg<<MsgDispatch;
       Msg(EsafMsg::Info)<<"                with phi between   "<<phmin/deg<<" and "<<phmax/deg<<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,photontype);

    //generation of SinglePhoton or BunchOfPhotons
    if(photonDescription == "single") {
       // only single photons direct to Euso
       if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE") 
           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,phtype,Direct);
           fPh_in_atmo->Add(p);
       }
    } 

    // one bunch of nbPhotons
    else if(photonDescription == "bunch"){  
        EarthVector dir(1);
        dir.SetTheta(theta);
        dir.SetPhi(phi);
        EarthVector Posf = GetLongitudinalExtension(Posi,dir);
        Double_t tf = (Posf-Posi).Mag()/Clight();
        BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,phtype);
        if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
        if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));

        fPh_in_atmo->Add(b);
	
	// track needed for ckov AND fluo (for last transfer until detector)
        BuildLightTrack(Posi,dir);
	//DELETE
	// for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
	//if(phtype == Cerenkov) BuildLightTrack(Posi,dir);  //DELETE
    }
    else Msg(EsafMsg::Panic)<<"<MakeSpot> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
   
    return fPh_in_atmo; 
}

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

    Msg(EsafMsg::Info)<<"TestLightSource HmaxSPOT with impact at (km) "<<impact.X()/km<<" "
		      <<impact.Y()/km<<" "<<impact.Z()/km<<MsgDispatch;

    string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
    string photontype = Conf()->GetStr("TestLightSource.PhotonType");
    string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
    Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
    Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
    Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
    Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
    PhotonType phtype = Fluo;
    if(photontype == "cerenkov") phtype = Cerenkov;

    TRandom* rndm = EsafRandom::Get();

    // bunch direction (for cerenkov bunches only)
    string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
    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(thmax > 80*DegToRad()) {
        thmax = 80*DegToRad();
        Msg(EsafMsg::Warning)<<"With HmaxSPOT, theta limited to 80 deg -> theta=80 applied"<<MsgDispatch;
    }
    Double_t theta_true = thmin + rndm->Rndm()*(thmax-thmin);
    Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);
    Double_t theta = theta_true;
    
    EarthVector Omega(1);
    // if theta given is local -> theta in MES is calculated
    if(thetaMode == "local") {
        EarthVector v1 = impact.Unit();
        EarthVector vrot;
        Omega.SetXYZ(1,0,0);
        EarthVector Uz(0,0,1);
        if(vrot.Mag() > 0) {
            Omega.Rotate(v1.Theta(),vrot);
            Omega.Rotate(phi,v1);
            vrot = v1.Cross(Omega);
            Omega.Rotate(Pi()/2+theta,vrot);
        }
	else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
	theta = Omega.Theta();
	phi   = Omega.Phi();
    }
    else if(thetaMode == "mes") Msg(EsafMsg::Panic)<<"mes ThetaMode NOT possible with HmaxSPOT option"<<MsgDispatch;
    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;
    
    // get Hmax from theta_true
    EarthVector Posi = GetHmaxPos(impact,Omega,theta_true);
    Msg(EsafMsg::Info)<<"TestLightSource HmaxSPOT with INITIAL pos at (km) "<<Posi.X()/km<<" " <<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
    
    if (directiontype == "UNIQUE") {
       Msg(EsafMsg::Info)<<"                with theta between "<<thmin/deg<<" and "<<thmax/deg<<MsgDispatch;
       Msg(EsafMsg::Info)<<"                with phi between   "<<phmin/deg<<" and "<<phmax/deg<<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,photontype);

    // only single photons direct to Euso
    if(photonDescription == "single") {
       if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE") 
           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,phtype,Direct);
           fPh_in_atmo->Add(p);
       }
    } 

    // one bunch of nbPhotons
    else if(photonDescription == "bunch"){  
        EarthVector dir(1);
        dir.SetTheta(Omega.Theta());
        dir.SetPhi(Omega.Phi());
        EarthVector Posf = GetLongitudinalExtension(Posi,dir);
        Double_t tf = (Posf - Posi).Mag()/Clight();
        BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,phtype);
        if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
        if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));

        fPh_in_atmo->Add(b);

	// track needed for ckov AND fluo (for last transfer until detector)
        BuildLightTrack(Posi,dir);
	//DELETE
	// for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
	//if(phtype == Cerenkov) BuildLightTrack(Posi,dir);  //DELETE
    }
    else Msg(EsafMsg::Panic)<<"<MakeHmaxSpot> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;

    return fPh_in_atmo; 
}

//_____________________________________________________________________________________________________
 PhotonsInAtmosphere *TestLightSource::MakeTrack(Double_t theta, 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");
    string photontype = Conf()->GetStr("TestLightSource.PhotonType");
    string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
    PhotonType phtype = Fluo;
    if(photontype == "cerenkov") phtype = Cerenkov;

    string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
    EarthVector Omega(1);
    if(thetaMode == "local") {
        EarthVector v1 = Posi.Unit();
        EarthVector vrot;
        Omega.SetXYZ(1,0,0);
        EarthVector Uz(0,0,1);
        if(vrot.Mag() > 0) {
            Omega.Rotate(v1.Theta(),vrot);
            Omega.Rotate(phi,v1);
            vrot = v1.Cross(Omega);
            Omega.Rotate(Pi()/2+theta,vrot);
        }
	else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
	theta = Omega.Theta();
	phi   = Omega.Phi();
    }
    else if(thetaMode == "mes") {
        theta = Pi() - theta;
	phi += Pi();
    }
    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;    
    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(photonDescription == "single") {
        Msg(EsafMsg::Info)<<"                direct to EUSO" << MsgDispatch;
        if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE") 
           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/Clight();

           TotalYield = MakeSpectrum(pos,&spectrum,photontype);

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

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

            TotalYield = MakeSpectrum(pos,&spectrum,photontype);

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

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

            if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(pos));
            if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));

           fPh_in_atmo->Add(b);
        }
	
	// for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
	if(phtype == Cerenkov) BuildLightTrack(Posi,DirTrack);
    }
    else Msg(EsafMsg::Panic)<<"<MakeTrack> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
    
    return fPh_in_atmo; 
}

//_________________________________________________________________________________________________
 Double_t TestLightSource::MakeSpectrum(const EarthVector& pos,EsafSpectrum* spectrum,string photontype) {
    //
    // 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 == "COMPLETE" && photontype == "fluo") {
       Double_t energy=80.*MeV;
       TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),energy,spectrum);
    }
    else if (spectrumtype == "COMPLETE" && photontype == "cerenkov") {
         TFormula cerenkov("cerenkov","1 /(x*x)");
         if ( spectrum !=0 ) spectrum -> Reset(&cerenkov,100,300*nm,450*nm);
	 TotalYield = 1.;
    }
    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.SpectrumType or PhotonType = " << spectrumtype<<", "<<photontype <<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(Double_t alt) {
    //
    // return Pointer to the angular distribution 
    //

   TF12 AngDist;

   Double_t EthCer = GetEnergyThreshold(alt);
   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;
    EarthVector Posf = Pos;
    /*
    Double_t L,h(0),depth(0);

    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;
	}
    }
    */
    Int_t status = Atmosphere::Get()->InvertGrammage(Pos,Dir,LgExt,Posf);
    if(status != 0) {
        Msg(EsafMsg::Warning) << "<GetLongitudinalExtension> : CANNOT FIND FINAL POSITION, InvertGrammage status = "<<status<<MsgDispatch;
	return Pos;
    }
    
    return Posf;
}

//__________________________________________________________________________________________________
 EarthVector TestLightSource::GetHmaxPos(const EarthVector& impact,const EarthVector& Direc,Double_t theta) const {
    //
    // get Hmax 3D-position -- for HmaxSPOT option only
    // Here theta MUST be LOCAL zenith angle (at impact)
    //
    
    EarthVector rtn(0,0,0);
    EarthVector Dir(Direc);
    if(Dir.Mag() != 1) Dir = Dir.Unit();
    Double_t tolerance = 0.5*m;
    if(theta > 80) Msg(EsafMsg::Warning) << "theta > 80deg NOT foreseen in <TestLightSource::GetHmaxPos>"<<MsgDispatch;
    
    // get Hmax value from LOCAL theta
    Double_t Hmax = (1.89 - 7.59*Log(Cos(theta))) * km;
    
    // get 3D-position
    EarthVector first = impact;
    EarthVector last = impact - (30*km/Cos(Dir.Theta() - Pi()))*Dir;
    EarthVector middle(1.);

    while((last.Zv() - first.Zv()) > tolerance) {
        middle = 0.5*(last+first);
        if(fabs(middle.Zv() - Hmax) < tolerance) {rtn = middle; break;}
        if(Hmax < middle.Zv()) last = middle;
        else first = middle;
    }
    if(rtn.Mag() == 0) rtn = first;
    
#ifdef DEBUG
cout<<" IN GETHMAXPOS :"<<endl;
cout<<"impact position = "<<impact<<endl;
cout<<"Hmax foreseen = "<<Hmax<<endl;
cout<<"Hmax implemented = "<<rtn.Zv()<<endl;
#endif
    
    return rtn;
}

//__________________________________________________________________________________________________
 void TestLightSource::BuildLightTrack(const EarthVector& posinit, const EarthVector& dir) {
    //
    // for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
    //
    
    // init
    ConfigFileParser* pConfig = Config::Get()->GetCF("RadiativeTransfer","BunchRadiativeTransfer");
    Double_t depthstep = pConfig->GetNum("BunchRadiativeTransfer.DepthStep")*g/cm2;
    fPh_in_atmo->ClearTrack();
    const Atmosphere* atmo = Atmosphere::Get();
    EarthVector presentpos(posinit);
    EarthVector nextpos(posinit);
    EarthVector posi(posinit);
    if(posi.IsUnderSeaLevel()) {
        Msg(EsafMsg::Warning) << "<BuildLightTrack> Starting pos is under sea level -> set to Nadir"<<MsgDispatch;
	posi.SetMag(0.);
    }
    Int_t status = -100;  // for Atmosphere::InvertGrammage() status
    
    // look at impact and total depth
    EarthVector impact = atmo->ImpactASL(posi,dir);
    if(impact.Z() == HUGE) impact = atmo->ImpactAtTOA(posi,dir);
    if(impact.Z() == HUGE) Msg(EsafMsg::Warning) << "<BuildLightTrack> : No GROUND impact nor TOA impact -> SHOULD NOT"<<MsgDispatch;
    Double_t FinalDepth = atmo->Grammage(posi,impact);
    
    // get nb of TrackLight steps
    UInt_t nb = UInt_t(floor(FinalDepth/depthstep)) + 2;  // +1 for posi,  +1 for tuning last step exit
    if(nb > 1000) {
        Msg(EsafMsg::Warning) << "<BuildLightTrack> : Tracklength is huge             = " << (impact - posi).Mag()/km<<" km"<<MsgDispatch;
        Msg(EsafMsg::Warning) << "<BuildLightTrack> : Corresponding depth is huge too = " << FinalDepth*cm2/g <<" g/cm2"<<MsgDispatch;
        Msg(EsafMsg::Warning) << "<BuildLightTrack> : Thus nb of TRACKLIGHT steps set to 1000"<<MsgDispatch;
	nb = 1000;
	depthstep = FinalDepth / (nb - 2); // -1 for posi,  -1 for tuning last step exit
    }
    fPh_in_atmo->SetNbTrackSteps(nb);
    
    // fill the track until last step entry (last step exit done by hand below
    for(UInt_t i=0; i<nb-1; i++) {
	presentpos = nextpos;
	fPh_in_atmo->FillTrack(presentpos);
	if(i == (nb-2)) break;  // to avoid line below for last step of the loop
        status = atmo->InvertGrammage(presentpos,dir,depthstep,nextpos);
	if(status != 0) Msg(EsafMsg::Warning) << "<BuildLightTrack> InvertGrammage() status is problematic -> status = " << status <<MsgDispatch;
    }
    
    // fine tune last step
    fPh_in_atmo->FillTrack(impact);


#ifdef DEBUG
cout<<"impact position = "<<impact<<endl;
    FinalDepth = atmo->Grammage(presentpos,impact);
cout<<"finaldepth w.r.t. depthstep = "<<FinalDepth/depthstep<<endl;
    if(FinalDepth/depthstep > 1) Msg(EsafMsg::Debug) << "<BuildLightTrack> : Last TRACKLIGHT step tuning failed, it counts as "<<FinalDepth/depthstep<<" times a normal step"<<MsgDispatch;    
#endif

}

//____________________________________________________________________________________________
 Double_t TestLightSource::GetEnergyThreshold(Double_t SC_alt) const {
    //
    // get the energy threshold for cerenkov emission
    //

    // Get Atmosphere Index
    const Atmosphere* atmo = Atmosphere::Get();
    Double_t delta = atmo->Index_Minus1(SC_alt);
    // When Index is not implemented in the Atmosphere use the Corsika formula
    if ( delta == 0) delta = 0.000283 * atmo->Air_Density( SC_alt )/atmo->Air_Density(0);
    
    return ElectronMassC2()/sqrt(2*delta);
}


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