Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

SlastShowerGenerator - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: SlastShowerGenerator.cc,v 1.74 2005/11/14 15:48:20 moreggia 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 "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>

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; 
    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
    if (!fQuiet) {
    Msg(EsafMsg::Info) << "SlastShowerGenerator Enabled" << MsgDispatch;
    Msg(EsafMsg::Info) << "SlastShowerGenerator Enables Atmosphere " << Atmosphere::Get()->GetType() << MsgDispatch;
    }   
    // General Shower Generator initializations
    ConfigFileParser *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");
    fEusoAltitude = pConfig->GetNum("GeneratorLightToEuso.altitude")*km; //TOFIX : to be changed when more general FoV handling achieved
    fEusoVector.SetXYZ(0,0,fEusoAltitude);

    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 = 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();
    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");
        // 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 {
        TRandom *rndm = EsafRandom::Get();
        fX1 =  -Log(rndm->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.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 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
    //

    TRandom *rndm = EsafRandom::Get();
    //Double_t x,y,z, Zmin, H = fEusoVector.Z() + EarthRadius();
    Double_t x,y,z, Zmin, H = fEusoVector.Z() + EarthRadius();
    ConfigFileParser* pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
    string ImpactMode = pConfig->GetStr("GeneratorLightToEuso.ImpactMode");
    string thetaMode  = pConfig->GetStr("GeneratorLightToEuso.ThetaMode");
    
    // if impact random on a limited surface
    if(ImpactMode == "random") {
            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 = EarthRadius() - (x*x + y*y)/(2*EarthRadius()); // needs -EarthRadius() before setting z !!  // approched value
	    //if(sqrt(x*x + y*y) > 20*km) break; //TOFIX: needs an infinite loop //GRNDetec
    }
	
    // within all Euso FoV
    else if(ImpactMode == "none") {
            Zmin         = H*Power(Sin(fFoV),2) + Cos(fFoV)*Sqrt( Power(EarthRadius(),2) - Power(H,2)*Power(sin(fFoV),2) );
            z            = (EarthRadius() - Zmin)*rndm->Rndm() + Zmin;
	    Double_t phi = rndm->Rndm()*2*TMath::Pi();
	    x            = EarthRadius()*Sin(ACos(z/EarthRadius()))*Cos(phi);
	    y            = EarthRadius()*Sin(ACos(z/EarthRadius()))*Sin(phi);	    
    }
     
    else Msg(EsafMsg::Panic) << "Wrong argument for GeneratorLightToEuso.ImpactMode config parameter" <<MsgDispatch;

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

    // 2. Generate random Theta and Phi with respect to the MSS
    fPhi = (fPhiMax - fPhiMin)*rndm->Rndm() + fPhiMin;
    
    // sin(2x)
    if(thetaMode == "sinus") {
        Double_t r1 = Cos(2*fThetaMin);
        Double_t r2 = Cos(2*fThetaMax);    
        fTheta = 0.5*ACos(r1 - rndm->Rndm()*(r1-r2));
    }
    // linear
    else if(thetaMode == "linear") fTheta = (fThetaMax - fThetaMin)*rndm->Rndm() + fThetaMin;
    // Hmax almost linear at E=1e20eV with USStd  ->  theta in [0,70]
    // comes from Hmax = 1.89 - 7.59*log(cos(th))
    else if(thetaMode == "FlatHmax") fTheta = ACos( Exp((1.89 - 10)*rndm->Rndm() / 7.59));
    else Msg(EsafMsg::Panic) << "Wrong argument for theta Mode  ->  " <<thetaMode << MsgDispatch;

    EarthVector v1 = fInitPoint.Unit();  // here init point = impact in earth-centered frame
    EarthVector vrot;
    fOmega.SetXYZ(1,0,0);
    EarthVector Uz(0,0,1);
    vrot = Uz.Cross(v1);
    if(vrot.Mag() > 0) {
        fOmega.Rotate(v1.Theta(),vrot);
        fOmega.Rotate(fPhi,v1);
        vrot = v1.Cross(fOmega);
        fOmega.Rotate(TMath::Pi()/2+fTheta,vrot);
    }
    else fOmega.SetMagThetaPhi(1.,TMath::Pi() - fTheta,fPhi + TMath::Pi());

    // 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) //TOFIX : shower dvpt should not be stopped outside FoV
        fInCome = kTRUE;

    //if (fInCome && ((fEusoVector-fCurrentPoint).Theta()>fFoV || fCurrentPoint.IsUnderSeaLevel()) ) 
        //fOutCome = kTRUE;                         //TOFIX : shower dvpt should not be stopped outside FoV

    //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)*0.1;
        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()/Clight();
    fTimenext    = (fNextPoint-fInitPoint).Mag()/Clight();
    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

    TRandom *rndm = EsafRandom::Get();
    if ( fEnergySlope == -1.) {   // flat differential spectrum in log scale
        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) );
    }

    // Check if only a specific point should be generated
    if( fFirstType == "POS") {
        fInitPoint.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);
        fPhi = (fPhiMax - fPhiMin)*rndm->Rndm() + fPhiMin;
        Double_t r1 = Cos(2*fThetaMin);
        Double_t r2 = Cos(2*fThetaMax);    
        fTheta = 0.5*ACos(r1 - rndm->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);
    Int_t particode = Int_t(floor(fAtomicMass+0.5));
    fTruth->SetParticle(particode);
    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 Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005