Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

ShowerTrack - source file

// $Id: ShowerTrack.cc,v 1.40 2005/10/02 14:17:55 thea Exp $
// Author: Alessandro Thea, Sergio Bottai, Dmitry Naumov  23/11/2003

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

//_____________________________________________________________________________
//
// ShowerTrack class serves as an universal container of Atmospheric Air Shower 
// produced by Ultra High Energy Cosmic Ray entering the atmosphere.
// The ShowerTrack object should be produced by interfaces to specific generators
// like CORSIKA, AIRES, UNISIM, SLAST, etc
// In the ShowerTrack object the relevant information is saved:
// vector of ShowerStep's, UHECR energy, theta, phi, first interaction depth, 
// Energy Threshold of electrons, Unit vector of the track direction
// first interaction point in MES
// Impact point on the Earth if any

// In addition ShowerTrack provides usefull functional utilities like: Histogram with energy, lateral etc distributions
// For the debugging purposes ShowerTrack object can draw itself as XY or XYZ view directly in root:
// root[]   ShowerTrack *track = ShowerGenerator->Get();
// root[]   track->DrawXYZ() as an example. 
// root[]   track->DrawXY() as an example.
// Check a usefull macro: macros/selftest/CheckDistrs.C
#include "PhysicsData.hh"
#include "ShowerTrack.hh"
#include "EarthVector.hh"

#include <stdexcept>
#include "EVector.hh"
#include <TF12.h>

using namespace sou;

ClassImp(ShowerTrack)

//______________________________________________________________________________    
Double_t ShowerTrackHillas(Double_t *x, Double_t *par) {
    //
    // A.M. Hillas, J.Phys.G:Nucl.Phys..8(1982) 1461-1473
    // This function is the derivate of energy of the original distribution done by D.V.Naumov
    //
    Double_t E = x[0], age = x[1];
    Double_t E0 = 26;
    if (age>=0.4)
        E0 = 44 - 17*TMath::Power((age-1.46),2);
    return 1.e8*TMath::Power((0.89*E0-1.2)/(E+E0),age)*age*(1.e4+2*E0+E*(2+age))/
        (E+E0)/TMath::Power(1.e4+E*age,3);
}


//______________________________________________________________________________
Double_t ShowerTrackGiller(Double_t *x, Double_t *par) {
    // 
    // M.Giller, A.Kacerzyk et al (ICRC 2003)
    //
    Double_t E = x[0], age = x[1];
    Double_t a = 1.005, b = 0.06, c = 189, d = 7.06*age + 12.48, C = 0.111*age + 0.134, Ecr = 80;
    return C/E*(1-a*TMath::Exp(-d*E/Ecr))*TMath::Power(1+E/Ecr,-(age+b*TMath::Log(E/Ecr/c)));
}

//______________________________________________________________________________
Double_t ShowerTrackNerling(Double_t *x, Double_t *par) {
    //
    // F.Nerling, R.Engel, C.Guerard et al (ICRC 2003). Also M.Risse (ICRC 2003, ICRC 2001)
    // Note that the numbers are different at the poster and in the paper by M.Risse
    // We adopt the paper values
    // 

    // This function should be yet worked around due to the wrong normalization
    Double_t E = x[0], age = x[1];
    Double_t a1 = 6.879 - 2.092*age, a2 = 122.0;
    return 1/(E+a1)/TMath::Power(E+a2,age);
}

//______________________________________________________________________________
Double_t ShowerTrackMelot(Double_t *x, Double_t *par) {
    //
    //
    //

    return 0;
}

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

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


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

// Set of static distributions. Note! Default does not mean that this is the best rather it is
// a convinient choice for the debugging study
TF2* ShowerTrack::fEnergyAgeDistributionDefault = new TF2("EneAgeDefault",ShowerTrackGiller,0.1,1000,0,2);
TF2* ShowerTrack::fEnergyAgeDistributionHillas  = new TF2("hillas",ShowerTrackHillas,0.1,1000,0,2);
TF2* ShowerTrack::fEnergyAgeDistributionGiller  = new TF2("giller",ShowerTrackGiller,0.1,1000,0,2);
TF2* ShowerTrack::fEnergyAgeDistributionNerling = new TF2("nerling",ShowerTrackNerling,0.1,1000,0,2);
TF2* ShowerTrack::fEnergyAgeDistributionMelot   = new TF2("melot",ShowerTrackMelot,0.1,1000,0,2);

TF2* ShowerTrack::fLateralDistributionDefalt = new TF2("LateralDistributionDefalt",ShowerTrackNKG1,0.001,50,0,2);
TF2* ShowerTrack::fLateralDistribution1 = new TF2("NKG1",ShowerTrackNKG1,0.001,50,0,2);
TF2* ShowerTrack::fLateralDistribution2 = new TF2("NKG2",ShowerTrackNKG2,0.001,5000,0,2,1);

TF2* ShowerTrack::fAngularDistribution1 = new TF2("baltru",ShowerTrackBaltru,0.,TMath::Pi(),.5,1000.);

//______________________________________________________________________________
 ShowerTrack::ShowerTrack(): PhysicsData("shower"), EsafMsgSource() {
    // 
    // Constructor
    // 
    SetNumEnergyHistos(20);
    SetNumLateralHistos(20);
}
//______________________________________________________________________________
 ShowerTrack::~ShowerTrack() {
    // 
    // Destructor
    // 
   Reset();

}
//______________________________________________________________________________
 TF2* ShowerTrack::GetEnergyDistribution(TString  name) {
    // 
    // Return Pointer to the energy distribution of electrons in the shower
    // 
    if( name == "hillas" )
        return fEnergyAgeDistributionHillas;
    else if ( name == "giller")
        return fEnergyAgeDistributionGiller;
    else if ( name == "nerling")
        return fEnergyAgeDistributionNerling;
    else if (name == "default" || name == "")
        return fEnergyAgeDistributionDefault;
    else
        throw runtime_error("ShowerTrack does not know this parametrization: ");
}
//______________________________________________________________________________
 TF2* ShowerTrack::GetLateralDistribution(TString  name) {
    // 
    // Return Pointer to the lateral distribution of electrons in the shower
    // 
    if( name == "NKG1" )
        return fLateralDistribution1;
    else if (name == "default" || name == "")
        return fLateralDistribution1;
    else if (name == "NKG2")
        return fLateralDistribution2;
    else
        throw runtime_error("ShowerTrack does not know this parametrization: ");
}
//______________________________________________________________________________
 TF2* ShowerTrack::GetAngularDistribution(TString  name) {
    // 
    // Return Pointer to the angular distribution of electrons in the shower
    // 
    if( name == "baltru" )
        return fAngularDistribution1;
    else if (name == "default" || name == "")
        return fAngularDistribution1;
    else
        throw runtime_error("ShowerTrack does not know this parametrization: ");
}

//______________________________________________________________________________
 const ShowerStep& ShowerTrack::GetStep( UInt_t i ) const {
    // 
    // Return the i-th step of the shower
    // 
    if ( i >= fSteps.size() )
	throw runtime_error("Index out of range in ShowerTrack::GetStep");
    return (fSteps[i]);
}
//______________________________________________________________________________
 const ShowerStep& ShowerTrack::GetLastStep( ) const {
    // 
    // Return last step of the shower
    // 
    return (fSteps[Size()-1]);
}

//______________________________________________________________________________
const ShowerStep& ShowerTrack::operator[]( UInt_t i ) const {
    // 
    // Return the i-th step using an operator 
    // 
	if ( i >= fSteps.size() )
		throw runtime_error("Index out of range in ShowerTrack::[]");
	return (fSteps[i]);
}
//______________________________________________________________________________
const vector<ShowerStep>& ShowerTrack::GetSteps() const {
    // 
    // Returns all steps in the track
    // 
	return fSteps;	
}

//______________________________________________________________________________
 void ShowerTrack::AtmUpdateTrack() {
    // Update the Track to the density profile of the atmosphere 
    // currently setted in ESAF. Leave the track direction in space unchanged,
    // leave all the shower physics in g/cm^2 unchanged. Change the
    // points in physical space. This is done in order to be able to handle
    // shower track created by shower generators outside ESAF which are using
    // a different atmospheric density profile respect to the one currently
    // setted in ESAF. 

}
//______________________________________________________________________________
 void ShowerTrack::Reset() {
  // delete each ShowerStep object
  if ( Clear() ) fSteps.clear();
  // reset hit ground information
  fHitGround = kFALSE;
}
//______________________________________________________________________________
 Bool_t ShowerTrack::Clear() {
    vector<ShowerStep>::iterator step;
    vector<ShowerStep>::iterator step2;
    
    for (step = fSteps.begin(); step != fSteps.end(); step++) {
        TH1F* Energy_tmp      = step->fHistoEnergy;
        TH1F* Lateral_tmp     = step->fHistoLateral;
        TH2F* EneAng_tmp      = step->fHistoEneAng;
        TH2F* RadPhiEle_tmp   = step->fHistoRadPhiEle;
        TH2F* RadDTimeEle_tmp = step->fHistoRadDTimeEle;
        TH2F* RadPhiEloss_tmp = step->fHistoRadPhiEloss;
        TH1F* AngCher_tmp     = step->fHistoAngCher;
        if (Energy_tmp||Lateral_tmp||EneAng_tmp||RadPhiEle_tmp||
            RadDTimeEle_tmp||RadPhiEloss_tmp ||AngCher_tmp ) {
            // scan all other steps to find the same pointer (if any)
            for (step2 = step+1; step2 != fSteps.end(); step2++) {
                if(Energy_tmp == step2->fHistoEnergy)           step2->fHistoEnergy      = NULL;
                if(Lateral_tmp == step2->fHistoLateral)         step2->fHistoLateral     = NULL;
                if(EneAng_tmp == step2->fHistoEneAng)           step2->fHistoEneAng      = NULL;
                if(RadPhiEle_tmp == step2->fHistoRadPhiEle)     step2->fHistoRadPhiEle   = NULL;
                if(RadDTimeEle_tmp == step2->fHistoRadDTimeEle) step2->fHistoRadDTimeEle = NULL;
                if(RadPhiEloss_tmp == step2->fHistoRadPhiEloss) step2->fHistoRadPhiEloss = NULL;
                if(AngCher_tmp == step2->fHistoAngCher)         step2->fHistoAngCher     = NULL;
            }
        }

        // save pointers to histos
        SafeDelete(step->fHistoEnergy);
        SafeDelete(step->fHistoLateral);
        SafeDelete(step->fHistoEneAng);
        SafeDelete(step->fHistoRadPhiEle);
        SafeDelete(step->fHistoRadDTimeEle);
        SafeDelete(step->fHistoRadPhiEloss);
        SafeDelete(step->fHistoAngCher);
    }
    Int_t i(0);
    for (step = fSteps.begin(); step != fSteps.end(); step++) {
        if (step->fHistoEnergy)      Msg(EsafMsg::Warning) << "Energy histo not deleted " << i << MsgDispatch;
        if (step->fHistoLateral)     Msg(EsafMsg::Warning) << "Lateral histo not deleted " << i << MsgDispatch;
        if (step->fHistoEneAng)      Msg(EsafMsg::Warning) << "EneAng histo not deleted " << i << MsgDispatch;
        if (step->fHistoRadDTimeEle) Msg(EsafMsg::Warning) << "RadTimeEle histo not deleted " << i << MsgDispatch;
        if (step->fHistoRadPhiEloss) Msg(EsafMsg::Warning) << "RadPhiEloss histo not deleted " << i << MsgDispatch;
        if (step->fHistoAngCher)     Msg(EsafMsg::Warning) << "AngCher histo not deleted " << i << MsgDispatch;
        i++;
    }
    return kTRUE;
}

//______________________________________________________________________________
 void ShowerTrack::Add( ShowerStep step ) { 
    step.SetStepID(Size() + 1);
    step.SetParentTrack(this); 
    fSteps.push_back(step);
}
//______________________________________________________________________________
 void ShowerTrack::DrawXYZ(Option_t *opt) {
    TString option(opt);
    Int_t n = GetNumStep();

    if (n == 0) {
        Msg(EsafMsg::Debug) << "DrawXYZ() informs: The track is empty. Check why" << MsgDispatch;
        return;
    }
    
    Double_t kHuge = 1.e20;
    Double_t Xmin(kHuge), Xmax(-kHuge), Ymin(kHuge), Ymax(-kHuge), Zmin(kHuge), Zmax(-kHuge), Nmax(1);
    Double_t Threshold(1.e-2);
    
    for (Int_t i=0; i<n; i++) {
        Double_t xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
        Double_t ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
        Double_t zmean = (GetStep(i).GetXYZi().Z() + GetStep(i).GetXYZf().Z())/2/km;
        Double_t ne    =  GetStep(i).GetNelectrons();
        if (ne/Nmax > Threshold) {
            Xmin = (Xmin < xmean ? Xmin : xmean); 
            Xmax = (Xmax > xmean ? Xmax : xmean); 
            Ymin = (Ymin < ymean ? Ymin : ymean); 
            Ymax = (Ymax > ymean ? Ymax : ymean); 
            Zmin = (Zmin < zmean ? Zmin : zmean); 
            Zmax = (Zmax > zmean ? Zmax : zmean); 
        }
        if (Nmax < GetStep(i).GetNelectrons() )
            Nmax = GetStep(i).GetNelectrons();
    }

    // build the histogram 
    TString name = "ShowerTrackXYZ";
    TH3F* th = (TH3F*)gROOT->FindObject(name);
    SafeDelete(th);    
    
    if (!th) {
        Double_t OffSet = 1.1;
        Double_t x1(-5),x2(5),y1(-5),y2(5),z1(0),z2(30); // make 10x10x30 kms box for default
        if (Xmin < x1) x1 = Xmin*OffSet;
        if (Xmax > x2) x2 = Xmax*OffSet;
        if (Ymin < y1) y1 = Ymin*OffSet;
        if (Ymax > y2) y2 = Ymax*OffSet;        
	Int_t xbins = Int_t((x2-x1)/1);  // make roughly 1 km bins size
	Int_t ybins = Int_t((y2-y1)/1);  // make roughly 1 km bins size
        th = new TH3F( name, "XYZ development", xbins,x1,x2,ybins,y1,y2,30,z1,z2); 
        th->SetStats(0);
        th->SetXTitle("X [km]");
        th->SetYTitle("Y [km]");
        th->SetZTitle("Z [km]");
    }
    for (Int_t i=0; i<n; i++) {
        double xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
        double ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
        double zmean = (GetStep(i).GetXYZi().Z() + GetStep(i).GetXYZf().Z())/2/km;
        double Ne    = GetStep(i).GetNelectrons();
        th->Fill(xmean,ymean,zmean,1.+Ne/Nmax);
    }
    th->Draw(option);
}

//______________________________________________________________________________
 void ShowerTrack::DrawXY(Option_t *opt) {
    TString option(opt);

    Int_t n = GetNumStep();

    if (n == 0) {
        Msg(EsafMsg::Debug) << "DrawXY() informs: The track is empty. Check why" << MsgDispatch;
        return;
    }
    
    Double_t kHuge = 1.e20;
    Double_t Xmin(kHuge), Xmax(-kHuge), Ymin(kHuge), Ymax(-kHuge),Nmax(0);
    Double_t Threshold(1.e-2);
    
    for (Int_t i=0; i<n; i++) {
        Double_t xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
        Double_t ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
        Double_t ne    =  GetStep(i).GetNelectrons();
        if (ne/Nmax > Threshold) {
            Xmin = (Xmin < xmean ? Xmin : xmean); 
            Xmax = (Xmax > xmean ? Xmax : xmean); 
            Ymin = (Ymin < ymean ? Ymin : ymean); 
            Ymax = (Ymax > ymean ? Ymax : ymean); 
        }
        if (Nmax < GetStep(i).GetNelectrons() )
            Nmax = GetStep(i).GetNelectrons();
    }
    // build the histogram 
    TString name = "ShowerTrackXY";
    TH2F* th = (TH2F*)gROOT->FindObject(name);
    SafeDelete(th);    
    
    if (!th) {
        Double_t OffSet = 1.1;
        Double_t x1(-5),x2(5),y1(-5),y2(5); // make 10x10 kms box for default
        if (Xmin < x1) x1 = Xmin*OffSet;
        if (Xmax > x2) x2 = Xmax*OffSet;
        if (Ymin < y1) y1 = Ymin*OffSet;
        if (Ymax > y2) y2 = Ymax*OffSet;  
	Int_t xbins = Int_t((x2-x1)/1);  // make roughly 1 km bins size
	Int_t ybins = Int_t((y2-y1)/1);  // make roughly 1 km bins size
	      
        th = new TH2F( name, "XY development", xbins,x1,x2,ybins,y1,y2); 
        th->SetStats(0);
        th->SetXTitle("X [km]");
        th->SetYTitle("Y [km]");
    }
    for (Int_t i=0; i<n; i++) {
        double xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
        double ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
        double Ne    = GetStep(i).GetNelectrons();
        th->Fill(xmean,ymean,Ne);
    }
    th->Draw(option);
}

//______________________________________________________________________________
 EarthVector ShowerTrack::FirstPos() const {
    //
    // return the entry position of the first ShowerStep
    //
    EarthVector rtn(0,0,HUGE);
    if(fSteps.size()) {
        const EVector& entry = fSteps[0].GetXYZi();
	rtn.SetXYZ(entry.X(),entry.Y(),entry.Z());
    }
    
    return rtn;
}

//______________________________________________________________________________
 EarthVector ShowerTrack::LastPos() const {
    //
    // return the exit position of the last ShowerStep
    //
    EarthVector rtn(0,0,HUGE);
    size_t last = fSteps.size();
    if(last) {
        const EVector& exit = fSteps[last-1].GetXYZf();
	rtn.SetXYZ(exit.X(),exit.Y(),exit.Z());
    }
    
    return rtn;
}


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