// $Id: ShowerTrack.cc,v 1.38 2005/04/19 14:27:32 naumov 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 "SystemOfUnits.hh"
#include "EarthVector.hh"
#include <stdexcept>
#include "EVector.hh"
#include <TF12.h>
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((s-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;
}
