source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/generators/showers/src/SlastShowerGenerator.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: SlastShowerGenerator.cc 2991 2011-09-26 10:08:50Z mabl $
// 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 "EsafRandom.hh"
#include "MCTruth.hh"
#include "ShowerTrack.hh"
#include "Atmosphere.hh"
#include "Config.hh"
#include "EConst.hh"
#include <TMath.h>
#include <TRandom.h>
#include <TProfile.h>
#include <TH1F.h>
#include <TGraph.h>
#include "NumbersFileParser.hh"

ClassImp(SlastShowerGenerator)
using namespace TMath;
using namespace sou;
using namespace EConst;

//_________________________________________________________________________________________ 
SlastShowerGenerator::SlastShowerGenerator(Bool_t quiet) : EventGenerator("slast++"), EsafMsgSource(),
                           fTrack(0), fTruth(0), fInCome(0), fOutCome(0), fQuiet(kFALSE)  {
    // This is ctor
    fQuiet = quiet; 
    fSpectrumRdnm = 0;
    fEnergy = -1;
    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
    //
    
    // General Shower Generator initializations
    ConfigFileParser *pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
    fFoV          = pConfig->GetNum("GeneratorLightToEuso.FoV")*deg;
    fSpectrumType = pConfig->GetStr("GeneratorLightToEuso.SpectrumType");
    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");
    fEusoAltitude = pConfig->GetNum("GeneratorLightToEuso.altitude")*km; //TOFIX : to be changed when more general FoV handling achieved
    fEusoVector.SetXYZ(0,0,fEusoAltitude);
    fRejectFakeEvents = pConfig->GetBool("GeneratorLightToEuso.RejectFakeEvents");
    fRejectNoXmax = pConfig->GetBool("GeneratorLightToEuso.RejectNoXmax");
    if(!fRejectFakeEvents && fRejectNoXmax)
         Msg(EsafMsg::Panic) << "<Init> if RejectFakeEvents not required, RejectNoXmax cannot be required" << MsgDispatch;
    
    // check if bounds are > 1deg -- EXCEPTED in "POS" option of fFirstType
    if((fThetaMax - fThetaMin) < kTolerance && fFirstType != "POS") {
        Msg(EsafMsg::Warning) << "<Init> theta bounds width cannot be zero : theta max increased by 1 degree" << MsgDispatch;
        fThetaMax += 1*DegToRad();
    }
    if((fPhiMax - fPhiMin) < kTolerance && fFirstType != "POS") {
        Msg(EsafMsg::Warning) << "<Init> phi bounds width cannot be zero : phi max increased by 1 degree" << MsgDispatch;
        fPhiMax += 1*DegToRad();
    }

    if (!fQuiet) {
        Msg(EsafMsg::Info) << "SlastShowerGenerator Enabled" << MsgDispatch;
        Msg(EsafMsg::Info) << "SlastShowerGenerator Enables Atmosphere " << Atmosphere::Get()->GetType() << MsgDispatch;
        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);
    }
    else {
        MsgForm(EsafMsg::Debug,"Generate %s event(s) in:",fType.c_str());
        MsgForm(EsafMsg::Debug,"Energy interval [%.4E,%.4E]",fEnergyMin,fEnergyMax);
        MsgForm(EsafMsg::Debug,"Theta  interval [%.4f,%.4f]",fThetaMin,fThetaMax);
        MsgForm(EsafMsg::Debug,"Phi    interval [%.4f,%.4f]",fPhiMin,fPhiMax);

    }
    fAtomicMass = AtomicMass(fType);
    
    // if non-analytical RC spectrum
    if(fSpectrumType == "GZKHiRes2005") {
        // GZK shape from best fit of HiRes data ( Phys. Lett. B 619 (2005) )
        // graph is recopied "by hand" using g3data software
        // between 1e19-1e21 eV only
        if(fEnergyMin < 1e19 || fEnergyMax > 1e21)
            Msg(EsafMsg::Panic) << "<Init> GZKHiRes2005 spectrum stands within [1e19 - 1e21] only !" << MsgDispatch;
        NumbersFileParser nf("config/Tools/GZKHiRes2005.dat",2);
        vector<Double_t> log10E_v   = nf.GetCol(0);  // expressed in eV
        vector<Double_t> fluxE3_v   = nf.GetCol(1);  // expressed in 10^24 eV2.m-2.s-1.sr-1
        TArrayD log10E(log10E_v.size()), fluxE3(log10E_v.size());
        for(size_t i=0; i < log10E_v.size(); i++) {
            log10E[i] = log10E_v[i];
            fluxE3[i] = fluxE3_v[i];
        }
        TGraph* gr = new TGraph(log10E_v.size(),log10E.GetArray(),fluxE3.GetArray());
        fSpectrumRdnm = new TH1F("GZK","GZK",10000,Log10(fEnergyMin),Log10(fEnergyMax));
        for(Int_t i=0; i < 10000; i++) {
            fSpectrumRdnm->SetBinContent(i+1,gr->Eval(fSpectrumRdnm->GetBinCenter(i+1),0,"") * log(10.) / pow(10,fSpectrumRdnm->GetBinCenter(i+1)*2. - 24.));
        }
    }
    
    
    // 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);
    SafeDelete(fSpectrumRdnm);
}

//_________________________________________________________________________________________ 
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();
    fFirstType = "POS";
}

//_________________________________________________________________________________________ 
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");
    }
    else if (fFirstType == "X1" ) {
        ConfigFileParser *pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
        fX1      = pConfig->GetNum("GeneratorLightToEuso.InteractionX1")*g/cm2;
    }
    else {
        TRandom *rndm = EsafRandom::Get();
        fX1 =  -Log(rndm->Rndm())*HadronInteractionLength(fAtomicMass,fEnergy)*g/cm2;
    }
}

//_________________________________________________________________________________________ 
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.e-9) + 1.14*Power(Log(E*1.e-9),2);
}

//_________________________________________________________________________________________ 
const ShowerStep& SlastShowerGenerator::GetShowerStep() {
    // Filling the ShowerStep object with relevant information
    fStep.Clear();
    fStep.fXi = fXcurrent + fX1;
    fStep.fXf = fXnext + fX1;
    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 Msg(EsafMsg::Panic) << "<GetShowerParametrization> Only GIL formula available so far" <<MsgDispatch;
    //TOFIX : GFA and GHF not coded correctly, they do not included X1 fluctuations
    /*
    else if (fShowerParametrization == "GFA") GFA(x + fX1);  
    else if (fShowerParametrization == "GHF") GHF(x + fX1); 
    else Msg(EsafMsg::Panic) << "<GetShowerParametrization> Wrong argument for ShowerParametrization :"
         << fShowerParametrization <<MsgDispatch;
    */
}

//_________________________________________________________________________________________ 
Bool_t SlastShowerGenerator::FindImpactOnTOA() {
    //
    // Generates point uniformly on top of atmosphere on the upper part of the earth hemisphere
    // theta, phi angles are taken at TOA impact, then translated into MES frame for further steps of event simulation
    //

    TRandom *rndm = EsafRandom::Get();
    ConfigFileParser* pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
    string ImpactMode = pConfig->GetStr("GeneratorLightToEuso.ImpactMode");
    string thetaMode  = pConfig->GetStr("GeneratorLightToEuso.ThetaMode");
    const Atmosphere* atmo = Atmosphere::Get();
    
    
#ifdef DEBUG
Int_t counterr = 0;
#endif
    
    
    // find impact on TOA, then find impact at sea level if any. Otherwise, set earthimpact to ******* TOA outgoing point **********
    // principle : - uniform distribution on a sphere of radius EarthRadius + TOA_altitude
    //             - sin(2theta) distribution, phi uniform distribution (LOCAL angles, at TOA impact centered frame)
    //             - check if defined track cut the EUSO FoV cone, if NOT  --> re-try
    if(ImpactMode == "TOA") {
        Double_t TOA_alt  = atmo->GetTOAAltitude();
        Double_t Euso_alt = fEusoVector.Z();
        Double_t FoV_radius = Euso_alt * Tan(fFoV + 1*deg);
        Double_t L1 = Sqrt( (EarthRadius() + TOA_alt)*(EarthRadius() + TOA_alt) - EarthRadius()*EarthRadius() );
        Double_t Zmin,x,y,z;
        Zmin = EarthRadius()*Sqrt(1 + 2*TOA_alt/EarthRadius() + (TOA_alt*TOA_alt - FoV_radius*FoV_radius)/(EarthRadius()*EarthRadius()));
        Zmin -= 2*L1 * Sqrt(1 - EarthRadius()*EarthRadius() / (EarthRadius() + TOA_alt) / (EarthRadius() + TOA_alt));

#ifdef DEBUG
//cout << "********** Zmin = "<<Zmin <<endl;

// *********** if EUSO tilted mode, Zmin is different. The rest of the code (FoV intersection etc..) should be adapted too *********
/*
else if(ImpactMode == "TOA_tiltedEuso") {
Double_t L2 = Sqrt( (EarthRadius() + Euso_alt)*(EarthRadius() + Euso_alt) - EarthRadius()*EarthRadius() );
Zmin = EarthRadius() + Euso_alt;
Zmin -= (L1 + L2) * Sqrt(1 - EarthRadius()*EarthRadius() / (EarthRadius() + Euso_alt) / (EarthRadius() + Euso_alt));
}
*/
#endif
        
        // find a good candidate shower, which has an opportunity to be seen by EUSO
        Bool_t eventisfound = kFALSE;
        Int_t infloop = 0;
        while(!eventisfound) {
#ifdef DEBUG
counterr++;
#endif
            // prevent infinite loop
            if(infloop++ > 100000) {
                Msg(EsafMsg::Warning) << "<FindImpactOnTOA> TOA impact mode : Inifinite loop broken by hand" << MsgDispatch;
                return kFALSE;
            }
            
            // 1. get an impact on TOA, uniformly distributed on the upper part of a sphere
            //    cut in Z (i.e. Zmin, expressed in earth-centered frame) is determined from geometrical considerations
            z            = (EarthRadius()+TOA_alt - Zmin)*rndm->Rndm() + Zmin;
            Double_t phi = rndm->Rndm()*TwoPi();
            x            = EarthRadius()*Sin(ACos(z / (EarthRadius()+TOA_alt) ))*Cos(phi);
            y            = EarthRadius()*Sin(ACos(z / (EarthRadius()+TOA_alt) ))*Sin(phi);          


            // 2. Generate random Theta and Phi w.r.t. local frame
            // sin(2x) OR linear OR Hmax modes
            // last one achieves flat Hmax distribution, at E=1e20eV with USStd  ->  theta in [0,70], comes from Hmax = 1.89 - 7.59*log(cos(th))
            fPhi_local = (TwoPi() - 0)*rndm->Rndm() + 0;
            if(thetaMode == "sinus") {
                Double_t r1 = 1.;  // cos(2*0)
                Double_t r2 = -1.; // cos(2*halfpi)
                fTheta_local = 0.5*ACos(r1 - rndm->Rndm()*(r1-r2));
            }
            else if(thetaMode == "linear") fTheta_local = (PiOver2() - 0)*rndm->Rndm() + 0;
            else if(thetaMode == "FlatHmax") fTheta_local = ACos( Exp((1.89 - 10)*rndm->Rndm() / 7.59));
            else Msg(EsafMsg::Panic) << "Wrong argument for theta Mode  ->  " <<thetaMode << MsgDispatch;
        
        
            // 3. Change fOmega from local frame to MES frame, using impact on TOA as the reference vector
            EarthVector TOA_impact_try(x,y,z);  // in earth-centered frame
            EarthVector v1 = TOA_impact_try.Unit();
            EarthVector vrot;
            fOmega.SetXYZ(1,0,0);
            EarthVector Uz(0,0,1);
            vrot = Uz.Cross(v1);
            if(vrot.Mag() > kTolerance) {
                fOmega.Rotate(v1.Theta(),vrot);
                fOmega.Rotate(fPhi_local,v1);
                vrot = v1.Cross(fOmega);
                fOmega.Rotate(PiOver2()+fTheta_local,vrot);
            }
            else fOmega.SetMagThetaPhi(1.,Pi() - fTheta_local,fPhi_local + Pi());
            
            // set angles in MES frame + check if they are within required bounds
            fTheta_mes = Pi() - fOmega.Theta();
            if(fTheta_mes < fThetaMin || fTheta_mes > fThetaMax) continue;
            
            // to get phi in [0,2pi]  (root returns it within [-pi,pi])
            if(fOmega.Phi() < 0) fPhi_mes = (fOmega.Phi() + TwoPi()) - Pi();
            else                 fPhi_mes = fOmega.Phi() + Pi();
            if(fPhi_mes < fPhiMin || fPhi_mes > fPhiMax) continue;
            
            
            // returns in MES frame
            TOA_impact_try.SetZ(TOA_impact_try.Z() - EarthRadius()); 


            fOutOfFoV = kFALSE;
            EarthVector intersection, impactASL;
            // 4. FoV handling
            if(IsInFoV(TOA_impact_try))  eventisfound = kTRUE;  // look if TOA impact is in FoV
            else if(!fRejectFakeEvents) {
                // event always kept in this case
                EarthVector intersection = FoVintersection(TOA_impact_try,fOmega);
                EarthVector impactASL    = atmo->ImpactASL(TOA_impact_try,fOmega);
                // check
                //    - if track cuts the EUSO FoV
                //    - if intersection is above ground
                //    - if intersection is below TOA
                //    - if track has reached ground before intersection
                intersection = FoVintersection(TOA_impact_try,fOmega);
                impactASL    = atmo->ImpactASL(TOA_impact_try,fOmega);
                if(intersection.Z() == -HUGE || intersection.Z() == HUGE || intersection.IsUnderSeaLevel()) fOutOfFoV = kTRUE;
                // check if intersection is below TOA altitude
                else if(intersection.Zv() > TOA_alt) fOutOfFoV = kTRUE;
                // check if track reaches ground before its intersection with euso FoV cone
                else if((impactASL - TOA_impact_try).Mag() < (intersection - TOA_impact_try).Mag()) fOutOfFoV = kTRUE;
                eventisfound = kTRUE;
            }
            else {
                // track is kept
                //    - if track cuts the EUSO FoV
                //    - if intersection is above ground
                //    - if intersection is below TOA
                //    - if track has reached ground before intersection
                intersection = FoVintersection(TOA_impact_try,fOmega);
                impactASL    = atmo->ImpactASL(TOA_impact_try,fOmega);
                if(intersection.Z() == -HUGE || intersection.Z() == HUGE || intersection.IsUnderSeaLevel()) continue;
                // check if intersection is below TOA altitude
                else if(intersection.Zv() > TOA_alt) continue;
                // check if track reaches ground before its intersection with euso FoV cone
                else if((impactASL - TOA_impact_try).Mag() < (intersection - TOA_impact_try).Mag()) continue;
                else eventisfound = kTRUE;
                fOutOfFoV = kFALSE;  // flag is disabled in 'fRejectFakeEvents == kTRUE' mode
            }
        }
        
        
        // 5. a valid TOA impact has been found
        fTOAImpact.SetXYZ(x,y,z - EarthRadius());

#ifdef DEBUG
Msg(EsafMsg::Debug) <<"********** Nb of trials = "<<counterr<<MsgDispatch;
/*
Msg(EsafMsg::Debug) <<"********** TOA found = "<<fTOAImpact<<MsgDispatch;
Msg(EsafMsg::Debug) <<"********** TOA found = "<<fTOAImpact.Mag()<<", "<<fTOAImpact.Theta()*RadToDeg()<<", "<<fTOAImpact.Phi()*RadToDeg()<<MsgDispatch;

Msg(EsafMsg::Debug) <<"********** LOCAL theta,phi = "<<fTheta_local*RadToDeg()<<", "<<fPhi_local*RadToDeg()<<MsgDispatch;
Msg(EsafMsg::Debug) <<"********** direction found = "<<fOmega<<MsgDispatch;
Msg(EsafMsg::Debug) <<"********** direction found = "<<fOmega.Mag()<<", "<<fOmega.Theta()*RadToDeg()<<", "<<fOmega.Phi()*RadToDeg()<<MsgDispatch;
*/
#endif


        // 4. Find impact at sea level if exists
        //    If it does not, set impact to outgoing point at TOA
        fEarthImpact = atmo->ImpactASL(fTOAImpact,fOmega);
        fHitGround = 0;
        if(fEarthImpact.Z() == -HUGE)      Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactASL says fTOAImpact is under sea level" << MsgDispatch;
        else if(fEarthImpact.Z() == HUGE)  Msg(EsafMsg::Debug) << "<FindImpactOnTOA> No ImpactASL found" << MsgDispatch;
        else fHitGround = 1;
        
        if(fHitGround == 0) {
            fEarthImpact = atmo->ImpactAtTOA(fTOAImpact,fOmega);
            if(fEarthImpact.Z() == -HUGE)      Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactAtTOA says fTOAImpact is under sea level" << MsgDispatch;
            else if(fEarthImpact.Z() == HUGE)  Msg(EsafMsg::Warning) << "<FindImpactOnTOA> No ImpactAtTOA found. IT SHOULD NOT OCCUR HERE" << MsgDispatch;
        }
    }
        
        
        
        
        
        
    // find impact at sea level, then find impact on TOA
    else if(ImpactMode == "ASL") {
        Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ASL impact mode : Cosmic Ray LOCAL direction is taken at sea level" << MsgDispatch;
        Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ASL impact mode : DOES NOT CHECK IF TRACK CUTS THE FOV" << MsgDispatch;

        Int_t infloop = 0;
        while(kTRUE) {
            if(infloop++ > 100000) {
                Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ASL impact mode : Inifinite loop broken by hand" << MsgDispatch;
                return kFALSE;
            }
            
            // 1. find impact at sea level
            Double_t x,y,z;
            Double_t Xmin = pConfig->GetNum("GeneratorLightToEuso.ImpactXmin")*km;
            Double_t Ymin = pConfig->GetNum("GeneratorLightToEuso.ImpactYmin")*km;
            Double_t Xmax = pConfig->GetNum("GeneratorLightToEuso.ImpactXmax")*km;
            Double_t Ymax = pConfig->GetNum("GeneratorLightToEuso.ImpactYmax")*km;
            x = rndm->Rndm()*(Xmax - Xmin) + Xmin;
            y = rndm->Rndm()*(Ymax - Ymin) + Ymin;
            z = sqrt(EarthRadius()*EarthRadius() - (x*x + y*y)); // needs -EarthRadius() before setting z !!

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


            // 2. Generate random Theta and Phi with respect to the MSS (local frame)
            // sin(2x) OR linear OR Hmax modes
            // last one achieving flat Hmax distribution, at E=1e20eV with USStd  ->  theta in [0,70], comes from Hmax = 1.89 - 7.59*log(cos(th))
            fPhi_local = (TwoPi() - 0)*rndm->Rndm() + 0;
            if(thetaMode == "sinus") {
                Double_t r1 = 1.;  // cos(2*0)
                Double_t r2 = -1.; // cos(2*halfpi)
                fTheta_local = 0.5*ACos(r1 - rndm->Rndm()*(r1-r2));
            }
            else if(thetaMode == "linear") fTheta_local = (PiOver2() - 0)*rndm->Rndm() + 0;
            else if(thetaMode == "FlatHmax") fTheta_local = ACos( Exp((1.89 - 10)*rndm->Rndm() / 7.59));
            else Msg(EsafMsg::Panic) << "Wrong argument for theta Mode  ->  " <<thetaMode << MsgDispatch;


            // 3. Change fOmega from local frame (MSS) to MES frame, using impact on ground as the reference vector
            EarthVector v1(x,y,z);
            v1 = v1.Unit();
            EarthVector vrot;
            fOmega.SetXYZ(1,0,0);
            EarthVector Uz(0,0,1);
            vrot = Uz.Cross(v1);
            if(vrot.Mag() > kTolerance) {
                fOmega.Rotate(v1.Theta(),vrot);
                fOmega.Rotate(fPhi_local,v1);
                vrot = v1.Cross(fOmega);
                fOmega.Rotate(Pi()/2+fTheta_local,vrot);
            }
            else fOmega.SetMagThetaPhi(1.,Pi() - fTheta_local,fPhi_local + Pi());


            // set angles in MES frame + check if they are within required bounds
            fTheta_mes = Pi() - fOmega.Theta();
            fOutOfFoV = kFALSE;
            // to get phi in [0,2pi]  (root returns it within [-pi,pi]
            if(fOmega.Phi() < 0) fPhi_mes = (fOmega.Phi() + TwoPi()) - Pi();
            else                 fPhi_mes = fOmega.Phi() + Pi();
            if(fPhi_mes < fPhiMin || fPhi_mes > fPhiMax || fTheta_mes < fThetaMin || fTheta_mes > fThetaMax) continue;
            else break;
        }
            

        // 4. Find TOA impact
        fTOAImpact = atmo->ImpactAtTOA(fEarthImpact,EarthVector(-fOmega));
        if(fTOAImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactAtTOA says fEarthImpact is under sea level" << MsgDispatch;
        if(fTOAImpact.Z() == HUGE)  Msg(EsafMsg::Warning) << "<FindImpactOnTOA> No ImpactAtTOA found : it Should not occur here" << MsgDispatch;
        
        
        // 4. refine value of impact at sea level (should exist)
        //    If it does not, set impact to outgoing point at TOA
        fEarthImpact = atmo->ImpactASL(fTOAImpact,fOmega);
        fHitGround = 0;
        if(fEarthImpact.Z() == -HUGE)      Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactASL says fTOAImpact is under sea level" << MsgDispatch;
        else if(fEarthImpact.Z() == HUGE)  Msg(EsafMsg::Debug) << "<FindImpactOnTOA> No ImpactASL found" << MsgDispatch;
        else fHitGround = 1;
        if(fHitGround == 0) {
            Msg(EsafMsg::Warning) << "<FindImpactASL> No earth impact found. IT SHOULD NOT OCCUR IN ASL mode" << MsgDispatch;
            fEarthImpact = atmo->ImpactAtTOA(fTOAImpact,fOmega);
            if(fEarthImpact.Z() == -HUGE)      Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactAtTOA says fTOAImpact is under sea level" << MsgDispatch;
            else if(fEarthImpact.Z() == HUGE)  Msg(EsafMsg::Warning) << "<FindImpactOnTOA> No ImpactAtTOA found. IT SHOULD NOT OCCUR HERE" << MsgDispatch;
        }
    }
     
    else Msg(EsafMsg::Panic) << "Wrong argument for GeneratorLightToEuso.ImpactMode config parameter" <<MsgDispatch;

    
    return kTRUE;
}

//_____________________________________________________________________________________________________
Bool_t SlastShowerGenerator::IsInFoV(const EarthVector& pos) const {
    //
    // returns true if pos is within Euso FoV cone
    //

    // 1. change of frame : from MES to NES (local frame with origin on EUSO + Z-axis directed toward the Earth)
    //    translation --> EUSO vector translation + rotation around x-axis to change Z direction to the opposite
    //    thus MES/NES Y-axes are in the opposite sens
    EarthVector Ux(1,0,0);
    EarthVector pos_e = pos; 
    pos_e -= fEusoVector; 
    pos_e.Rotate(Pi(),Ux);
    
    if(pos_e.Theta() < fFoV) return kTRUE;
    else return kFALSE;
}

//_____________________________________________________________________________________________________
EarthVector SlastShowerGenerator::FoVintersection(const EarthVector& pos, const EarthVector& dir) const {
    //
    // finds the intersection between the given track and the Euso FoV cone, if exists
    // if does not exist, returns (0,0,HUGE)
    //
    // intersection considered only above ground
    //
    
    Double_t w = Cos(fFoV)*Cos(fFoV);
    EarthVector direc = dir.Unit();

    // 1. if pos is underground
    EarthVector rtn;
    if( (pos.Zv() + kAltitudeTolerance) < 0 ) {
        Msg(EsafMsg::Warning) << "<FoVintersection> Starting position is underground, SHOULD NOT OCCUR" << MsgDispatch;
        rtn.SetXYZ(0,0,-HUGE);
        return rtn;
    }
    
    // 2. change of frame : from MES to NES (local frame with origin on EUSO + Z-axis directed toward the Earth)
    //    translation --> EUSO vector translation + rotation around x-axis to change Z direction to the opposite
    //    thus MES/NES Y-axes are in the opposite sens
    EarthVector Ux(1,0,0);
    EarthVector pos_e = pos; 
    EarthVector direc_e = direc;

    pos_e -= fEusoVector;  
    direc_e.Rotate(Pi(),Ux);
    pos_e.Rotate(Pi(),Ux);
    

    // 3. now intersection calculation
    Double_t mag(0);
    // trinome defining solutions
    Double_t a = direc_e.Z()*direc_e.Z() - w;
    Double_t b = 2*pos_e.Z()*direc_e.Z() - 2*w*pos_e*direc_e;
    Double_t c = pos_e.Z()*pos_e.Z() - pos_e.Mag2()*w;
    pair<Int_t,Double_t*>& p = findRoots(a,b,c);  // special cases with a=0 and b=c=0 handled
    if(p.first == 0) {
        rtn.SetXYZ(0,0,HUGE);
        return rtn;
    }
    if(p.first == 1) mag = p.second[0];
    else if(p.first == 2) {
        if((p.second[0]+kAltitudeTolerance) * (p.second[1]+kAltitudeTolerance) < 0) mag = Max(p.second[0],p.second[1]);
        else mag = Min(p.second[0],p.second[1]);
    }
    // if intersection stands before pos, it is rejected
    if(mag < -kAltitudeTolerance) {
        rtn.SetXYZ(0,0,HUGE);
        return rtn;
    }
    
    rtn = pos_e + mag*direc_e;    
    
    
    // 4. change of frame : from NES to MES frame
    rtn.Rotate(-Pi(),Ux);
    rtn += fEusoVector;

    return rtn;
}

//_________________________________________________________________________________________ 
Bool_t SlastShowerGenerator::GetFirstPoint(){
    //
    // X1 is known --> find the first point fInitPoint from TOA impact
    // if event goes out atmosphere / OR / reaches ground without interacting, returns kFALSE
    //
    
    Bool_t rtn = kTRUE;
    
    if(fFirstType == "POS") {
        fXcurrent     = 0;
        fXnext        = fXcurrent;
        fCurrentPoint = fInitPoint;
        fNextPoint    = fCurrentPoint;
        fTimecurrent  = 0;
        fTimenext     = 0;
        return rtn;
    }
    
    Int_t status = Atmosphere::Get()->InvertGrammage(fTOAImpact,fOmega,fX1,fInitPoint);
    if(status == -1) Msg(EsafMsg::Info) << "<GetFirstPoint> InvertGrammage() status is -1, it Should not" << MsgDispatch;
    
    else if(status == 0) {
        fXcurrent     = 0;
        fXnext        = fXcurrent;
        fCurrentPoint = fInitPoint;
        fNextPoint    = fCurrentPoint;
        fTimecurrent  = 0;
        fTimenext     = 0;
    }
    
    else if(status == 1)
             Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status is 1, it SHOULD NOT" << MsgDispatch;
             
    else if(status == 2) {
        if(Atmosphere::Get()->Grammage(fTOAImpact,Atmosphere::Get()->ImpactAtTOA(fTOAImpact,fOmega)) < fX1+kTolerance)
             Msg(EsafMsg::Info)    << "<GetFirstPoint> TOA reached : Primary has not interacted in the atmosphere" << MsgDispatch;
        else Msg(EsafMsg::Warning) << "<GetFirstPoint> PROBLEM : TOA reached, BUT primary SHOULD have interacted in the atmosphere" << MsgDispatch;
        rtn = kFALSE;
    }
    
    else if(status == 3) {
        if(Atmosphere::Get()->Grammage(fTOAImpact,Atmosphere::Get()->ImpactASL(fTOAImpact,fOmega)) < fX1+kTolerance)
             Msg(EsafMsg::Info)    << "<GetFirstPoint> Primary has hit the ground without interacting in the atmosphere " << MsgDispatch;
        else Msg(EsafMsg::Warning) << "<GetFirstPoint> PROBLEM : ground reached, BUT primary SHOULD have interacted before" << MsgDispatch;
        rtn = kFALSE;
    }
    
    else     Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status unknown : "<< status << MsgDispatch;
    
    return rtn;
}

//_________________________________________________________________________________________ 
Bool_t SlastShowerGenerator::IsXmaxInFoV(){
    //
    // Assess Xmax position and tell if it is in FoV
    // used to reject event with Xmax out of FoV, without simulating the dvpt
    //
    
    // inits
    Bool_t rtn = kTRUE;
    
    // get Xmax according to shower dvpt model
    if (fShowerParametrization != "GIL")
        Msg(EsafMsg::Panic) << "<GetShowerParametrization> Only GIL formula available so far" <<MsgDispatch;
    //TOFIX : GFA and GHF not coded correctly, they do not included X1 fluctuations
    // when these options are reactivated, change the code below accordingly
    Float_t x0 = 37.15*g/cm2;
    Float_t Xmax = (1.7 + 0.76*(Log(fEnergy/0.81e8) - Log(fAtomicMass))) * x0 + fX1;
    EarthVector Maxpos(1);

    
    Int_t status = Atmosphere::Get()->InvertGrammage(fTOAImpact,fOmega,Xmax,Maxpos);
    if(status == -1) Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status is -1, it Should not" << MsgDispatch;
    else if(status == 1) Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status is 1, it SHOULD NOT" << MsgDispatch;
    // TOA or sea level reached before
    else if(status == 2 || status == 3) rtn = kFALSE;
    else if(status == 0) rtn = IsInFoV(Maxpos);
    else Msg(EsafMsg::Panic) << "<GetFirstPoint> InvertGrammage() status unknown : "<< status << MsgDispatch;
    
    return rtn;
}

//_________________________________________________________________________________________ 
Bool_t SlastShowerGenerator::GetNextStep() {
    // This method checks if next point is visible and makes the next step if it is OK.
    
    //TOFIX : old flag
    fInCome = kTRUE;

    if(fCurrentPoint.IsUnderSeaLevel()) return kFALSE;
    fCurrentPoint = fNextPoint;
    
    // search next shower step position
    Int_t status = Atmosphere::Get()->InvertGrammage(fCurrentPoint,fOmega,fDepthStep,fNextPoint);
    if(status == -1) {
       Msg(EsafMsg::Warning) << "<GetNextStep> InvertGrammage() status is -1, it Should not" << MsgDispatch;
       return kFALSE;
    }
    else if(status == 0) {
        fXcurrent = fXnext;
        fXnext    = fXcurrent + fDepthStep;
        fTimecurrent = (fCurrentPoint-fInitPoint).Mag()/Clight();
        fTimenext    = (fNextPoint-fInitPoint).Mag()/Clight();
    }
    else if(status == 1) {
       Msg(EsafMsg::Warning) << "<GetNextStep> InvertGrammage() status is 1, it Should not" << MsgDispatch;
       return kFALSE;
    }
    else if(status == 2) {
       Msg(EsafMsg::Info) << "<GetNextStep> Shower goes out atmosphere : TOA reached " << MsgDispatch;
       return kFALSE;
    }
    else if(status == 3) {
       Msg(EsafMsg::Info) << "<GetNextStep> Shower reaches the ground " << MsgDispatch;
       return kFALSE;
    }
    else Msg(EsafMsg::Warning) << "<GetNextStep> InvertGrammage() status unknown : "<< status << MsgDispatch;
    
    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 x0 = 37.15;
    Float_t t  = x*cm2/g/x0;
    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) {
        Msg(EsafMsg::Warning)<<"<GIL> Asking for depth <= 0. Zero returned"<<MsgDispatch;
        fNe = 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
    // sigma values fitted and extrapolated, Xmax, Nmax from GHF used
    // Factor 1.1 added to Ne because of threshold energy used for corsika simulation (10% particles missing)
    //
    // WARNING !!!     not usable yet, cause cannot currently handle X1 fluctuations
    //

    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;
        Nmax =  Exp(-20.95 + 2.292*Log10(fEnergy));
        Xmax = -258.1 + 54.81*Log10(fEnergy);
    }
    else if(fAtomicMass == 56) {
        sigma =  0.470-0.0135*Log10(fEnergy);
        p0    = -9.043;
        p1    =  0.9923;
        Nmax =  Exp(-21.68 + 2.329*Log10(fEnergy));
        Xmax = -430.7 + 59.08*Log10(fEnergy);
    }
    else throw invalid_argument(Form("Gaussian parameterization in age not available for primaries A = %f",fAtomicMass));

    Xmax = Xmax*g/cm2;
    fAge = 3*x/(x+2.*Xmax);
    if(fAge<=0) {
        Msg(EsafMsg::Warning)<<"<GFA> Asking for depth <= 0. Zero returned"<<MsgDispatch;
        fNe = 0.;
        return;
    }
    fNe = 1.11*Nmax * Exp( -(fAge-1)*(fAge-1) / (2.*sigma*sigma) );
    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 
    // Factor 1.1 added to Ne because of threshold energy used for corsika simulation (10% particles missing)
    //
    // WARNING !!!     not usable yet, cause cannot currently handle X1 fluctuations
    //

    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 if(fAtomicMass == 56) {
        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);
    }
    else throw invalid_argument(Form("Gaussian parameterization in age not available for primaries A = %f",fAtomicMass));
    
    Xmax = Xmax*g/cm2;
    X1   = X1*g/cm2;
    X0   = X0*g/cm2;
    //  X=x+fX1;
    X=x;
    lambda = lambda*g/cm2;
    fNe = 1.11*Nmax * Power((X-X0)/(Xmax-X0),(Xmax-X0)/lambda) * Exp((Xmax-X)/lambda);
    fAge = 3*x/(x+2.*Xmax); //should be defined ?
    
    if(fAge<=0) {
        Msg(EsafMsg::Warning)<<"<GHF> Asking for depth <= 0. Zero returned"<<MsgDispatch;
        fNe = 0.;
        return;
    }

    return;
}

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

    // generate random Energy, Theta, Phi and the interaction point
    // energy slope is for differential spectrum
    TRandom *rndm = EsafRandom::Get();
    
    // if power law
    if(fSpectrumType == "powerlaw") {
        // flat differential spectrum in log scale
        if ( fEnergySlope == -1.) {
            fEnergy = fEnergyMin * Exp(rndm->Rndm() * Log(fEnergyMax/fEnergyMin) );
        }
        else {
            Double_t e1 = Power(fEnergyMin,fEnergySlope + 1);
            Double_t e2 = Power(fEnergyMax,fEnergySlope + 1);
            fEnergy = Power( e1 + (e2-e1)*rndm->Rndm(),1/(fEnergySlope + 1) );
        }
    }
    else if(fSpectrumType == "GZKHiRes2005") {
        fEnergy = -HUGE;
        while((fEnergy < fEnergyMin) || (fEnergy > fEnergyMax)) {
            Double_t log10E = fSpectrumRdnm->GetRandom();
            fEnergy = Power(10.,log10E);
        }
    }
    else Msg(EsafMsg::Panic)<<"<Get> Wrong argument for spectrum type == "<< fSpectrumType <<MsgDispatch;



    // Check if only a specific point should be generated
    if( fFirstType == "POS") {
        fInitPoint.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);
        fPhi_mes = (fPhiMax - fPhiMin)*rndm->Rndm() + fPhiMin;
        fTheta_mes = (fThetaMax - fThetaMin)*rndm->Rndm() + fThetaMin;
        fOmega.SetXYZ(-Sin(fTheta_mes)*Cos(fPhi_mes),-Sin(fTheta_mes)*Sin(fPhi_mes),-Cos(fTheta_mes));
        Msg(EsafMsg::Info) << "<Get()> Pos option, LOCAL theta and phi set equal to MES values" << MsgDispatch;
        fTheta_local = fTheta_mes;
        fPhi_local = fPhi_mes;
        
        // find TOA impact and earth impact
        fTOAImpact = Atmosphere::Get()->ImpactAtTOA(fInitPoint,EarthVector(-fOmega));
        fEarthImpact = Atmosphere::Get()->ImpactASL(fInitPoint,fOmega);
        fHitGround = 0;
        if(fEarthImpact.Z() == -HUGE)      Msg(EsafMsg::Warning) << "<Get(), POS option> ImpactASL says fInitPoint is under sea level" << MsgDispatch;
        else if(fEarthImpact.Z() == HUGE)  Msg(EsafMsg::Debug) << "<Get(), POS option> Shower does not reach the ground" << MsgDispatch;
        else fHitGround = 1;
        if(fHitGround == 0) {
            fEarthImpact = Atmosphere::Get()->ImpactAtTOA(fInitPoint,fOmega);
            if(fEarthImpact.Z() == -HUGE)      Msg(EsafMsg::Warning) << "<Get(), POS option> ImpactAtTOA says fInitPoint is under sea level" << MsgDispatch;
            else if(fEarthImpact.Z() == HUGE)  Msg(EsafMsg::Warning) << "<Get(), POS option> No ImpactAtTOA found. IT SHOULD NOT OCCUR HERE" << MsgDispatch;
        }
        
        // disable standard geometrical flags
        fRejectFakeEvents = kFALSE;
        fOutOfFoV = kFALSE;
        fRejectNoXmax = kFALSE;

        DevelopShower();
    }
    
    
    // STANDARD WAY
    else {
        Int_t cycle(0);
        while (kTRUE) {
            if ( cycle++>1000000 )                       break;
            if ( FindImpactOnTOA() && DevelopShower())   break;
        }
    }
    if (!fQuiet) {
        MsgForm(EsafMsg::Info,"ShowerTrack produced : Energy = %.2E MeV, Theta_local = %.2f deg, Phi_local = %.2f deg",
                fTrack->GetEnergy(), fTheta_local*RadToDeg(),fPhi_local*RadToDeg());
        MsgForm(EsafMsg::Info,"         Impact on TOA : (%.2f,%.2f,%.2f) km with MES angles : Theta_mes = %.2f, Phi_mes = %.2f",
                fTOAImpact.X()/km,fTOAImpact.Y()/km,fTOAImpact.Z()/km,fTheta_mes*RadToDeg(),fPhi_mes*RadToDeg());
        MsgForm(EsafMsg::Info,"         First interaction X1 = %.3f at : (%.2f,%.2f,%.2f) km gives a shower with %d steps",
                fTrack->fX1*cm2/g,fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
    }
    else {
        MsgForm(EsafMsg::Debug,"SlastShowerGenerator produced a ShowerTrack: Energy = %.2E MeV Theta = %.2f deg Phi = %.2f deg",
                fTrack->GetEnergy(), fTrack->GetTheta()*RadToDeg(),fTrack->GetPhi()*RadToDeg());
        MsgForm(EsafMsg::Debug,"                              originated with X1 = %.3f at: (%.2f,%.2f,%.2f) km with %d steps",
                fTrack->fX1*cm2/g,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,HUGE);
    fOmega.SetXYZ(0,0,HUGE);
    fCurrentPoint.SetXYZ(0,0,HUGE);
    fNextPoint.SetXYZ(0,0,HUGE);
    fEarthImpact.SetXYZ(0,0,HUGE);
    fTOAImpact.SetXYZ(0,0,HUGE);
    fTheta_local  = -HUGE;
    fPhi_local    = -HUGE;
    fTheta_mes  = -HUGE;
    fPhi_mes    = -HUGE;
    fX1     = -HUGE;
    fEnergy = -HUGE;
    
}

//_________________________________________________________________________________________ 
void SlastShowerGenerator::GetShowerInfo() {
    // Fills ShowerTrack parameters
    fTrack->fEnergy    = fEnergy*eV;
    fTrack->fTheta     = fTheta_mes;
    fTrack->fPhi       = fPhi_mes;
    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();
    Bool_t has_interacted = GetFirstPoint();
    
    // init
    if(!fTrack) fTrack = new ShowerTrack();
    else fTrack->Reset();

    // 1. if fake events not allowed, re-try for another event
    //    NB : fOutOfFoV flag is not relevant here, cause no showers should be in this case. But kept for safety
    if(fRejectFakeEvents && (!has_interacted || fOutOfFoV)) return kFALSE;
    
    // 1.1 if events without Xmax in the FoV are not allowed, re-try for another event
    if(fRejectNoXmax) {
        if(!IsXmaxInFoV()) {
            Msg(EsafMsg::Info)<<"Shower Xmax is not in FoV --> retry for another event"<<MsgDispatch;
            return kFALSE;
        }
    }

    // 2. if fake events allowed + CR has not interacted
    // --> X1 related data fixed to HUGE values
    if(!has_interacted) {
        fX1 = -HUGE;
        fInitPoint.SetXYZ(0,0,HUGE);
        GetShowerInfo();
        Msg(EsafMsg::Info)<<"Cosmic ray has not interacted in the atmosphere"<<MsgDispatch;
        return kTRUE;
    }

    // 3. if fake events allowed + events is out of FoV, shower dvpt avoided
    if(!fRejectFakeEvents && fOutOfFoV) {
        GetShowerInfo();
        Msg(EsafMsg::Info)<<"ShowerTrack does NOT cut the FoV"<<MsgDispatch;
        return kTRUE;
    }

    // 4. "standard" candidate
    //    Develop a shower of shower steps
    Int_t cycle(0);
    while(kTRUE) {
        if(cycle++>100000) return kFALSE;
        if(!GetNextStep()) break;
        if( fInCome) {
            ShowerStep s = GetShowerStep();
            fTrack->Add(s);
        }
    }
    GetShowerInfo();
    Msg(EsafMsg::Info)<<"ShowerTrack cut the FoV"<<MsgDispatch;
    
    return kTRUE;
}

//_________________________________________________________________________________________ 
MCTruth* SlastShowerGenerator::GetTruth() {
    //
    // Method returning the Truth object 
    //
    
    // init
    if(!fTruth) fTruth = new MCTruth();
    
    // can be filled even if showertrack has no steps
    fTruth->SetEnergy(fTrack->GetEnergy());
    fTruth->SetThetaPhi(fTrack->GetTheta(),fTrack->GetPhi());
    fTruth->SetFirstInt(fInitPoint,fX1);
    Int_t particode = Int_t(floor(fAtomicMass+0.5));
    fTruth->SetParticle(particode);
    fTruth->SetTOAImpact(fTOAImpact);
    fTruth->SetEarthImpact(fEarthImpact);
    
    // filled only if steps exist    
    if (fTrack->Size() > 0) {
       const ShowerStep& s = fTrack->GetLastStep();
       fTruth->SetEarthAge(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;
}

//_________________________________________________________________________________________ 
Double_t SlastShowerGenerator::GetXmax() {
    //
    // returns Xmax of the shower according to GIL parametrization (used by Reco)
    //
    Float_t x0 = 37.15*g/cm2;
    Float_t fEnergyScale;
    if (fEnergy<0) fEnergyScale = 1.e20;
    else           fEnergyScale = fEnergy;
    Float_t t_max = 1.7 + 0.76*(Log(fEnergyScale/0.81e8) - Log(fAtomicMass));
    return t_max*x0;
}