Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

SlastShowerGenerator - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: SlastShowerGenerator.cc,v 1.53 2005/05/12 09:15:11 stutz Exp $
// Author: Dmitry V.Naumov  03/03/2004

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: SlastShowerGenerator                                                 *
 *  Package: Shower                                                          *
 *  Coordinator: Sergio.Bottai, Dmitry Naumov                                *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// SlastShowerGenerator
//
// Slast is short of "Shower Light Attenuated to the Space Telescope".
// Originally written by Dmitry V.Naumov as a collection of F77 routines doing
// the full chain:
// Shower generation, light production (fluorescence and cherenkov), light
// attenuation in the atmosphere. 
// See SlastLightToEuso C++ / Fortran interface which is doing the full chain.
//
// SlastShowerGenerator is the simple generator written in pure C++ which is
// NOT doing the full simulation chain. Instead it returnes only the
// ShowerTrack object which is taken later by ShowerLightSource and
// RadiativeTransfer parts of the ESAF package.
//
// ESAF prodides some testing macros which can help to understand how to use
// SlastShowerGenerator.  See macros/CheckDistrs.C as an example.
//
// A simple way to use SlastShowerGenerator is directly in root:
//
// root[]   gSystem->Load("libatmosphere.so");
// root[]   gSystem->Load("libgenbase.so");
// root[]   gSystem->Load("libshowers.so");
// root[]   SlastShowerGenerator *ShowerGenerator = new SlastShowerGenerator;
// root[]   ShowerTrack *track = ShowerGenerator->Get();
// root[]   track->DrawXYZ()
//
// as an example.
//
// Now having the track object in hand one can use whatever functions avalaible
// for it. See ShowerTrack class for documentation of ShowerTrack object

#include "SlastShowerGenerator.hh"
#include "MCTruth.hh"
#include "ShowerTrack.hh"
#include "Atmosphere.hh"
#include "Config.hh"
#include "EConst.hh"
#include <TMath.h>
#include <TRandom2.h>
#include <TProfile.h>

ClassImp(SlastShowerGenerator)
using EConst::EarthRadius;
using namespace TMath;

//_________________________________________________________________________________________ 
SlastShowerGenerator::SlastShowerGenerator() : EventGenerator("slast++"), EsafMsgSource(),
                           fTrack(0), fTruth(0), fInCome(0), fOutCome(0), fQuiet(kFALSE)  {
// This is ctor
    if(!Init()) Msg(EsafMsg::Panic)<<"SlastShowerGenerator can not be initialized"<<MsgDispatch;
}


//_________________________________________________________________________________________ 
 Bool_t SlastShowerGenerator::Init()
{
    // SlastShowerGenerator is initialized reading from config/LightToEuso/GeneratorLightToEuso.cfg file
    // Edit that file to setup desired parameters of the showers such as energy, theta, phi intervals
    // or initial position in case of unique parameters.
    // Init() is called in SlastShowerGenerator ctor
    Msg(EsafMsg::Info) << "SlastShowerGenerator Enabled" << MsgDispatch;
    // Enabling the atmosphere
    Atmosphere::Get()->GetType();
    Msg(EsafMsg::Info) << "SlastShowerGenerator Enables Atmosphere " << Atmosphere::Get()->GetType() << MsgDispatch;
    // EUSO detector initializations
    
    ConfigFileParser *pConfig = Config::Get()->GetCF("General","Euso");
    fEusoAltitude= pConfig->GetNum("Euso.fAltitude")*km;
    fEusoVector.SetXYZ(0,0,fEusoAltitude);
    // General Shower Generator initializations
    pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
    fFoV       = pConfig->GetNum("GeneratorLightToEuso.FoV")*deg;
    fEnergyMin = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMin");
    fEnergyMax = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMax");
    fEnergySlope = pConfig->GetNum("GeneratorLightToEuso.EnergySlope");
    fThetaMin  = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMin")*deg;
    fThetaMax  = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMax")*deg;
    fPhiMin    = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMin")*deg;
    fPhiMax    = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMax")*deg;
    fType      = pConfig->GetStr("GeneratorLightToEuso.UhecrType");
    fDepthStep = pConfig->GetNum("GeneratorLightToEuso.DepthStep")*g/cm2;
    fFirstPointX = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorX")*km;
    fFirstPointY = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorY")*km;
    fFirstPointZ = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorZ")*km;
    fFirstType = pConfig->GetStr("GeneratorLightToEuso.InteractionType");

    if (!fQuiet) {
	MsgForm(EsafMsg::Info,"Generate %s event(s) in:",fType.c_str());
	MsgForm(EsafMsg::Info,"Energy interval [%.4E,%.4E]",fEnergyMin,fEnergyMax);
	MsgForm(EsafMsg::Info,"Theta  interval [%.4f,%.4f]",fThetaMin,fThetaMax);
	MsgForm(EsafMsg::Info,"Phi    interval [%.4f,%.4f]",fPhiMin,fPhiMax);
    }
    fAtomicMass = EConst::AtomicMass(fType);
    // Initialize SLAST specific options
    pConfig = Config::Get()->GetCF("LightToEuso","SlastLightToEuso");
    fEnergyDistribution    = pConfig->GetStr("SlastLightToEuso.EnergyDistributionParametrization");
    fShowerParametrization = pConfig->GetStr("SlastLightToEuso.ShowerParametrization");
    return kTRUE;
}

//_________________________________________________________________________________________ 
SlastShowerGenerator::~SlastShowerGenerator() {
//  This is dtor
    SafeDelete(fTrack);
    SafeDelete(fTruth);
}

//_________________________________________________________________________________________ 
 void SlastShowerGenerator::SetShower(Double_t e,Double_t t,Double_t p,EarthVector pos) {
    // This is a special method to setup desired showertrack parameters from the code (not from the config file)
    // This can be a usefull function in the reconstruction code
    
    SetEnergy(fEnergyMin = fEnergyMax = e); 
    SetTheta(fThetaMin = fThetaMax = t);
    SetPhi(fPhiMin = fPhiMax = p);
    fFirstPointX = pos.X();
    fFirstPointY = pos.Y();
    fFirstPointZ = pos.Z();
}

//_________________________________________________________________________________________ 
 void SlastShowerGenerator::GetX1() {
    // Returns atmosphere depth before the first interaction. 
    // Attention: units are ESAF internal. In order to get in human way one has to multiple by cm2/g 
    if ( fFirstType == "POS" ) {
        fX1 = Atmosphere::Get()->Grammage(fInitPoint,-fOmega,"dir");
        // cout << " fInitPoint " << fInitPoint.X() <<" "<<fInitPoint.Y()<<" "<<fInitPoint.Z()<<endl;
        // cout << " fX1 POS " << fX1*cm2/g << endl;;
    }
    else if (fFirstType == "X1" ) {
        ConfigFileParser *pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
        fX1      = pConfig->GetNum("GeneratorLightToEuso.InteractionX1")*g/cm2;
        // cout << " fX1 X1 " << fX1*cm2/g << endl;;
    }
    else
        fX1 =  -Log(gRandom->Rndm())*HadronInteractionLength(fAtomicMass,fEnergy)*g/cm2;
        // cout <<" fX1 std " << fX1*cm2/g<<endl;
}


//_________________________________________________________________________________________ 
 Double_t SlastShowerGenerator::HadronInteractionLength(Double_t A, Double_t E) {
    // SlastShowerGenerator internal function which computes the hadronic interactio  length.
    // Used by GetX1() method
    Float_t R0 = 1.287e-13, Aair = 14.483,  q = 0.93, AN = 14.00674, AO = 15.9994;
    Float_t AA = 39.948, AC = 12.011, CN = 0.7847, CO = 0.2105, CA = 0.0047, CC = 0.0001;
    Float_t S_A = Pi()*Power(R0,2);
    Float_t F_A1 = Power(S_A*Na(),-1);
    Float_t F_A2 = Aair*F_A1;
    Float_t X = Power(A,1./3);
    Float_t Y = 1./X;
    Float_t B = Power(AN,1./3);
    Float_t S = CN*Power(B+X-q*(1./B+Y),2); // Nitrogen (N2);
    B  = Power(AO,1./3);
    S += CO*Power(B+X-q*(1./B+Y),2);     // Oxigen (O2)
    B  = Power(AA,1./3);
    S += CA*Power(B+X-q*(1./B+Y),2);     // Argon (Ar)
    B  = Power(AC,1./3);
    S += CC*Power(B+X-q*(1./B+Y),2);     // Carbon (in CO2)
    return F_A2/S*NucleonAirCrossSection(1.e11/A)/NucleonAirCrossSection(E/A);
}


//_________________________________________________________________________________________ 
 Double_t SlastShowerGenerator::NucleonAirCrossSection(Double_t E) {
    /*
    Computes NUCLEN + AIR CROSS SECTION:
        ... USES AN EMPIRICAL PARAMETRIZATION OF TOTAL INELASTIC 
        ... NUCLEON AIR CROSS SECTION (in mbarn)
        ... INPUT:  NUCLEON ENERGY in eV
        ... OUTPUT: NUCLEON AIR CROSS SECTION (in mbarn)
        ... REFERENCES: 
        ... 1. Mielke H.H., Foller, M, Knapp J. 
        ...    // J.Phys.G: Nucl.Part.Phys. 1994. V 20, P637
        ... 2. V.A. Naumov, T.S. Sinegovskaya (eq (17))
    */
    return 290 - 8.7*Log(E/1.e9) + 1.14*Power(Log(E/1.e9),2);
}


//_________________________________________________________________________________________ 
 const ShowerStep& SlastShowerGenerator::GetShowerStep() {
    // Filling the ShowerStep object with relevant information
    fStep.Clear();
    fStep.fXi = fXcurrent;
    fStep.fXf = fXnext;
    fStep.fXYZi = fCurrentPoint;
    fStep.fXYZf = fNextPoint;
    fStep.fTimei = fTimecurrent;
    fStep.fTimef = fTimenext;
    GetShowerParametrization(fXcurrent);
    fStep.fAgei = fAge;
    GetShowerParametrization(fXnext);
    fStep.fAgef = fAge;
    GetShowerParametrization((fXnext+fXcurrent)/2);
    fStep.fNelectrons = fNe;
    return fStep;
}


//_________________________________________________________________________________________ 
 void SlastShowerGenerator::GetShowerParametrization(Double_t x) {
    // Assigning of fNe - number of electrons in the current step according to the Shower Parameterization.
    //  GIL stands for Greizen-Ilina-Linsley which is the QGSJET model 
    // parameterization.
    // GFA stands for Gaussian Function in Age
    // GHF stands for Gaisser Hillas Function
    if (fShowerParametrization == "GIL"){
        GIL(x);}
    else if (fShowerParametrization == "GFA"){
        GFA(x);}  
    else if (fShowerParametrization == "GHF"){
        GHF(x);}  
    else 
        throw invalid_argument("SlastShowerGenerator does not know about ShowerParametrization "+
               fShowerParametrization);
}


//_________________________________________________________________________________________ 
 Bool_t SlastShowerGenerator::GetFarFirstPoint() {
    // This method generates the first ``far'' point which is ensured to be visible by EUSO FoV.
    // 1. Generate point uniformly on the Earth surface underneath of EUSO
    Double_t x,y,z, Zmin, H = fEusoVector.Z() + EarthRadius();
    Zmin = H*Power(Sin(fFoV),2) + 
        Cos(fFoV)*Sqrt( Power(EarthRadius(),2) - 
                Power(H,2)*Power(sin(fFoV),2) );
    while (kTRUE) {    
           gRandom->Sphere(x,y,z,EConst::EarthRadius());
        if (z > Zmin) break;
    }

    fInitPoint.SetXYZ(x,y,z);
    fEarthImpact = fInitPoint;
    fEarthImpact.SetZ(fInitPoint.Z()-EarthRadius());

    // 2. Generate random Theta and Phi 
    fPhi = (fPhiMax - fPhiMin)*gRandom->Rndm() + fPhiMin;

    Double_t r1 = Cos(2*fThetaMin);
    Double_t r2 = Cos(2*fThetaMax);    
    fTheta = 0.5*ACos(r1 - gRandom->Rndm()*(r1-r2));

    //    fOmega.SetXYZ(-Sin(fTheta)*Cos(fPhi),-Sin(fTheta)*Sin(fPhi),-Cos(fTheta));
    EarthVector v1 = fInitPoint.Unit();
    EarthVector vrot;
    fOmega.SetXYZ(1,0,0);
    EarthVector Uz(0,0,1);
    vrot = Uz.Cross(v1);
    fOmega.Rotate(v1.Theta(),vrot);
    fOmega.Rotate(fPhi,v1);
    vrot = v1.Cross(fOmega);
    fOmega.Rotate(TMath::Pi()/2-fTheta,vrot);

    cout << "fOmega1 " <<-Sin(fTheta)*Cos(fPhi)<<" "<<-Sin(fTheta)*Sin(fPhi)<<" "<<-Cos(fTheta)<<endl;
    cout << "fOmega2 " <<fOmega.X()<<" "<<fOmega.Y()<<" "<<fOmega.Z()<<endl;

    // 3. Jump to the higher altitude along the generated direction to the sphere of Radius rEarth + h
    Float_t h = 80*km;
    // 4. Find first need track length
    Double_t L = fOmega*fInitPoint + Sqrt(h*(2*fInitPoint.Mag() + h) + Power(-fOmega*fInitPoint,2));
    fInitPoint -= fOmega*L;    
    fInitPoint.SetZ(fInitPoint.Z()-EarthRadius());
    fHitGround = 1;
    return kTRUE;
}


//_________________________________________________________________________________________ 
 void SlastShowerGenerator::GetFirstPoint(){
    // Once ``far'' first point is found this method find the real first point according to the generated X1
    Double_t L = 0.1*km, depth(0), h(0);
    depth = Atmosphere::Get()->Grammage(fInitPoint,-fOmega,"dir");

    while(true) {
        if ( fFirstType == "POS" ) break;

        fInitPoint += fOmega*L;
	if (!fInitPoint.IsUnderSeaLevel())
             h      = fInitPoint.Zv();
	else 
	     break;     
        depth      += Atmosphere::Get()->Air_Density(h)*L;
        if(depth>=fX1) break;
    }
    fXcurrent = depth-fX1;
    fXnext    = fXcurrent;
    fCurrentPoint = fInitPoint;
    fNextPoint    = fCurrentPoint;
    fTimecurrent  = 0;
    fTimenext     = 0;
}


//_________________________________________________________________________________________ 
 Bool_t SlastShowerGenerator::GetNextStep() {
    // This method checks if next point is visible, not outside of FoV and makes the next step if 
    // it is OK.
    if ((fEusoVector-fCurrentPoint).Theta()<=fFoV) 
        fInCome = kTRUE;

    if (fInCome && ((fEusoVector-fCurrentPoint).Theta()>fFoV || fCurrentPoint.IsUnderSeaLevel()) ) 
        fOutCome = kTRUE;

    if (fOutCome) {
        Msg(EsafMsg::Debug) << " fOutCome true " << MsgDispatch;
        return kFALSE;
    }
    if (fCurrentPoint.IsUnderSeaLevel()) return kFALSE;
    fCurrentPoint = fNextPoint;
    fNextPoint    = fCurrentPoint;
    Double_t L, h(0), depth(0);
    Int_t cycling(0);
    while(1) {
        cycling++;
        L              = fDepthStep/Atmosphere::Get()->Air_Density(h)/10;
        fNextPoint    += fOmega*L;
	if (!fNextPoint.IsUnderSeaLevel())
             h         = fNextPoint.Zv();
	else 
	     return kFALSE;     
        depth         += Atmosphere::Get()->Air_Density(h)*L;
	
        if (depth>=fDepthStep)  break;
        if (cycling>100000) {
            MsgForm(EsafMsg::Debug,"NextPoint not found : cycling > 10000 with depth %f ",depth/g*cm2);
            return kFALSE;
        }
    }
    if (h<0){
       Msg(EsafMsg::Debug) << " Shower reaches the ground " << MsgDispatch;
       return kFALSE;
    }
    fXcurrent = fXnext;
    fXnext    = fXcurrent + depth;
    fTimecurrent = (fCurrentPoint-fInitPoint).Mag()/EConst::C();
    fTimenext    = (fNextPoint-fInitPoint).Mag()/EConst::C();
    return kTRUE;
}


//_________________________________________________________________________________________ 
 void SlastShowerGenerator::GIL(Double_t x)
{
    // Assigning of fNe - number of electrons in the current step according to GIL parameterization.
    // (GIL stands for Greizen-Ilina-Linsley which is the QGSJET model parameterization.
    
    Float_t t0 = 37.15;
    Float_t t  = x*cm2/g/t0;
    Float_t t_max = 1.7 + 0.76*(Log(fEnergy/0.81e8) - Log(fAtomicMass));
    fAge  = 2*t/(t+t_max);
    fNe   = 0;
    if(fAge<=0) return;
    fNe = fEnergy/1.45e9*Exp(t-t_max-2*t*Log(fAge));
}


//_________________________________________________________________________________________ 
 void SlastShowerGenerator::GFA(Double_t x)
{
    // Assigning of fNe - number of electrons in the current step according to Gaussian Function in Age parameterization.
    // Formulae from C. Song, Astropart. Physics 22(2004)151
    // Xmax, Nmax, sigma values fitted and extrapolated 

    Float_t Xmax,Nmax;
    Float_t sigma,p0,p1;
    if (fAtomicMass==1) {
        sigma =  0.3206 - 0.006597*Log10(fEnergy);
        p0    = -8.9475;
        p1    =  0.9874;
    }
    else {
        sigma =  0.470-0.0135*Log10(fEnergy);
        p0    = -9.043;
        p1    =  0.9923;
    }
    //  Xmax=(-3426.*sigma+1410.)*g/cm2;
    Xmax = ((200*Log10(fEnergy)-7026)*sigma+1410.)*g/cm2;
    Nmax = pow(10,p1*Log10(fEnergy)+p0);
    fAge = 3*x/(x+2.*Xmax);
    if(fAge<=0) return;
    fNe = Nmax*Exp((-1./(2.*sigma*sigma))*Power(fAge-1.,2));
    return;
}

//_________________________________________________________________________________________ 
 void SlastShowerGenerator::GHF(Double_t x)
{
    // Assigning of fNe - number of electrons in the current step according to
    // Gaisser Hillas    parametrisation
    // Formulae from C. Song, Astropart. Physics 22(2004)151
    // Xmax, Nmax, lambda, X1, X0 values fitted and extrapolated 

    Float_t Xmax,Nmax;
    Float_t X,X1,X0,lambda;
    if (fAtomicMass==1) {
        X1   =  106.2 -  3.15*Log10(fEnergy);
        X0   =  401.9 - 25.79*Log10(fEnergy);
        Nmax =  Exp(-20.95 + 2.292*Log10(fEnergy));
        Xmax = -258.1 + 54.81*Log10(fEnergy);
        lambda= 98.6 -1.922*Log10(fEnergy);
    }
    else {
        X1   =   29.72 - 1.031*Log10(fEnergy);
        X0   =  471.1 - 29.49*Log10(fEnergy);
        Nmax =  Exp(-21.68 + 2.329*Log10(fEnergy));
        Xmax = -430.7 + 59.08*Log10(fEnergy);
        lambda= 182.4 -6.1*Log10(fEnergy);
    }
    Xmax = Xmax*g/cm2;
    X1   = X1*g/cm2;
    X0   = X0*g/cm2;
    //  X=x+fX1;
    X=x;
    lambda = lambda*g/cm2;
    fNe = Nmax * Power((X-X0)/(Xmax-X0),(Xmax-X0)/lambda) * Exp((Xmax-X)/lambda);
    fAge = 3*x/(x+2.*Xmax); //should be defined ?

    return;
}


//_________________________________________________________________________________________ 
 PhysicsData* SlastShowerGenerator::Get() {
    // method returning the ShowerTrack object

    // generate random Energy, Theta, Phi and the interaction point
    // energy slope is for differential spectrum
    if ( fEnergySlope == -1.) {   // flat differential spectrum in log scale
        fEnergy = fEnergyMin * Exp(gRandom->Rndm() * Log(fEnergyMax/fEnergyMin) );
    }
    else {
        Double_t e1 = Power(fEnergyMin,fEnergySlope + 1);
        Double_t e2 = Power(fEnergyMax,fEnergySlope + 1);
        fEnergy = Power( e1 + (e2-e1)*gRandom->Rndm(),1/(fEnergySlope + 1) );
    }
    // Check if only a specific point should be generated
    if( fFirstType == "POS") {
        fInitPoint.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);
        fPhi = (fPhiMax - fPhiMin)*gRandom->Rndm() + fPhiMin;
        Double_t r1 = Cos(2*fThetaMin);
        Double_t r2 = Cos(2*fThetaMax);    
        fTheta = 0.5*ACos(r1 - gRandom->Rndm()*(r1-r2));

        fOmega.SetXYZ(-Sin(fTheta)*Cos(fPhi),-Sin(fTheta)*Sin(fPhi),-Cos(fTheta));
        DevelopShower();
        //if (DevelopShower())
        //            fTrack->fInitPos.SetXYZ(fInitPoint.X(),fInitPoint.Y(),fInitPoint.Z()); // Reset Initial Position to required
    }
    else {
        Int_t cycle(0);
        while (kTRUE) {
            if ( cycle++>1000 )                          break;
            if ( GetFarFirstPoint() && DevelopShower())  break;
        }
    }
    if (!fQuiet) {
        MsgForm(EsafMsg::Info,"SlastShowerGenerator produced a ShowerTrack: Energy = %.2E MeV Theta = %.2f deg Phi = %.2f deg",
                fTrack->GetEnergy(), fTrack->GetTheta()*RadToDeg(),fTrack->GetPhi()*RadToDeg());
        MsgForm(EsafMsg::Info,"                              originated at: (%.2f,%.2f,%.2f) km with %d steps",
                fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
    }
    return fTrack;
}

//_________________________________________________________________________________________ 
 void SlastShowerGenerator::Reset() {
    // Reset ShowerTrack and SlastShowerGenerator internal parameters
    fHitGround = 0;
    fInCome = kFALSE;
    fOutCome = kFALSE;
    fInitPoint.SetXYZ(0,0,0);
    fOmega.SetXYZ(0,0,0);
    fCurrentPoint.SetXYZ(0,0,0);
    fNextPoint.SetXYZ(0,0,0);
    fEarthImpact.SetXYZ(0,0,0);
}

//_________________________________________________________________________________________ 
 void SlastShowerGenerator::GetShowerInfo() {
    // Fills ShowerTrack parameters
    fTrack->fEnergy    = fEnergy*eV;
    fTrack->fTheta     = fTheta;
    fTrack->fPhi       = fPhi;
    fTrack->fX1        = fX1;
    fTrack->fEthrEl    = 0;
    fTrack->fDirVers   = fOmega;
    fTrack->fInitPos   = fInitPoint;
    fTrack->fHitGround = fHitGround;
    fTrack->fEarthImpact = fEarthImpact;
}

//_________________________________________________________________________________________ 
 Bool_t SlastShowerGenerator::DevelopShower() {
    // Method developing the shower
    GetX1();
    GetFirstPoint();
    // Develop a shower of shower steps
    if( !fTrack)
        fTrack = new ShowerTrack();
    else
        fTrack->Reset();
    Int_t cycle(0);
    while (kTRUE) {
        if(cycle++>100000) 
            return kFALSE;
        if(!GetNextStep())
            break;
        if( fInCome) {
            ShowerStep s = GetShowerStep();
            fTrack->Add(s);
        }
    }
    GetShowerInfo();
    return kTRUE;
}

//_________________________________________________________________________________________ 
 MCTruth* SlastShowerGenerator::GetTruth() {
    // Method returning the Truth object 
    if (!fTrack->Size()) return NULL;
    if (!fTruth) 
        fTruth = new MCTruth();
    fTruth->SetEnergy(fTrack->GetEnergy());
    fTruth->SetThetaPhi(fTrack->GetTheta(),fTrack->GetPhi());
    fTruth->SetFirstInt(fInitPoint,fX1);
    const ShowerStep& s = fTrack->GetLastStep();
    fTruth->SetEarthImpact(fEarthImpact,s.GetAgef());
    // get the shower maximum
    Float_t MaxElectrons(0);
    Int_t MaxIndex(-1);
    for (UInt_t i(0); i<fTrack->Size(); i++) {
        const ShowerStep& s = fTrack->GetStep(i);
        if (MaxElectrons < s.GetNelectrons()) {
            MaxIndex = i;
            MaxElectrons = s.GetNelectrons();
        }
    }
    if (MaxIndex != -1) {
        const ShowerStep& s = fTrack->GetStep(MaxIndex);
        fTruth->SetShowerMax(0.5*(s.GetXYZi() + s.GetXYZf()),0.5*(s.GetXi() + s.GetXf()));
    }
    return fTruth;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31