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

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 23.3 KB
Line 
1// $Id: ShowerTrack.cc 2808 2009-01-29 10:01:14Z naumov $
2// Author: Alessandro Thea, Sergio Bottai, Dmitry Naumov  23/11/2003
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: ShowerStep                                                           *
8 *  Package: Shower                                                          *
9 *  Coordinator: Sergio.Bottai, Dmitry.Naumov                                *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// ShowerTrack class serves as an universal container of Atmospheric Air Shower
16// produced by Ultra High Energy Cosmic Ray entering the atmosphere.
17// The ShowerTrack object should be produced by interfaces to specific generators
18// like CORSIKA, AIRES, UNISIM, SLAST, etc
19// In the ShowerTrack object the relevant information is saved:
20// vector of ShowerStep's, UHECR energy, theta, phi, first interaction depth,
21// Energy Threshold of electrons, Unit vector of the track direction
22// first interaction point in MES
23// Impact point on the Earth if any
24
25// In addition ShowerTrack provides usefull functional utilities like: Histogram with energy, lateral etc distributions
26// For the debugging purposes ShowerTrack object can draw itself as XY or XYZ view directly in root:
27// root[]   ShowerTrack *track = ShowerGenerator->Get();
28// root[]   track->DrawXYZ() as an example.
29// root[]   track->DrawXY() as an example.
30// Check a usefull macro: macros/selftest/CheckDistrs.C
31#include "PhysicsData.hh"
32#include "ShowerTrack.hh"
33#include "EarthVector.hh"
34
35#include <stdexcept>
36#include "EVector.hh"
37#include <TF12.h>
38#include <TMath.h>
39#include <TROOT.h>
40using namespace sou;
41using namespace TMath;
42
43ClassImp(ShowerTrack)
44
45//______________________________________________________________________________   
46Double_t ShowerTrackHillas(Double_t *x, Double_t *par) {
47    //
48    // A.M. Hillas, J.Phys.G:Nucl.Phys..8(1982) 1461-1473
49    // This function is the derivate of energy of the original distribution done by D.V.Naumov
50    //
51    Double_t E = x[0], age = x[1];
52    Double_t E0 = 26;
53    if (age>=0.4)
54        E0 = 44 - 17*Power((age-1.46),2);
55    return 1.e8*Power((0.89*E0-1.2)/(E+E0),age)*age*(1.e4+2*E0+E*(2+age))/
56        (E+E0)/Power(1.e4+E*age,3);
57}
58
59
60//______________________________________________________________________________
61Double_t ShowerTrackGiller(Double_t *x, Double_t *par) {
62    //
63    // M.Giller, A.Kacerzyk et al (ICRC 2003 + J.Phys.G Nucl. Part. Phys. 30 (2004))
64    // Normalization : spectrum MUST be integrated using dE  (not dln(E))
65    //
66   
67    Double_t E = x[0], age = x[1];
68    Double_t a = 1.005, b = 0.06, c = 189, d = 7.06*age + 12.48, C = 0.111*age + 0.134, Ecr = 80;
69    return C / E *(1-a*Exp(-d*E/Ecr))*Power(1+E/Ecr,-(age+b*Log(E/Ecr/c)));
70}
71
72//______________________________________________________________________________
73Double_t ShowerTrackNerling(Double_t *x, Double_t *par) {
74    //
75    // Nerling et al. Astropart. 24 (2006) 421-437, and Nerling PhD report
76    //
77    // even if a parametrization in age of normalization coeff. is given in the article, it depends on threshold energy
78    // so here spectrum IS NOT NORMALIZED !! Normalization SHOULD BE DONE numerically in the rest of the code
79    // Normalization : spectrum MUST be integrated using dE  (not dln(E))
80    //
81
82    // This function should be yet worked around due to the wrong normalization
83    Double_t E = x[0]/MeV, age = x[1];
84    Double_t a1 = 6.42522 - 1.53183*age, a2 = 168.168 - 42.1368*age;
85    return 1. / (E+a1) / Power(E+a2,age);
86}
87
88//______________________________________________________________________________
89Double_t ShowerTrackMelot(Double_t *x, Double_t *par) {
90    //
91    //
92    //
93
94    return 0;
95}
96
97//______________________________________________________________________________
98Double_t ShowerTrackNKG1(Double_t *x, Double_t *par) {
99    //
100    // dNe/dx derived from 'standard' NGK formula for Lateral distribution. x
101    // is the distance to shower axis divided by Ro=moliere radius : x=D/Ro .
102    // The integration of NGK1 over x from 0. to infinity is normalized to 1 .
103    // The integral from radius r1 to radius r2 give the fraction of total
104    // number of electrons which lie inside such interval.
105    //
106   
107    Double_t D = x[0], age = x[1];
108    Double_t e1=2.0;
109    Double_t e2=4.5;
110    return Gamma(e2-age)/Gamma(age)/Gamma(e2-2.*age)*Power(D,age-(e1-1.0))
111        *Power((1.0+D),age-e2);
112}
113
114//______________________________________________________________________________
115Double_t ShowerTrackNKG2(Double_t *x, Double_t *par) {
116    //
117    // dNe/dr derived from 'standard' NGK formula for Lateral distribution. r is
118    // the distance to shower axis Rm=moliere radius is a parameter . The
119    // integration of NGK2 over r from 0. to infinity is normalized to 1 . The
120    // integral from radius r1 to radius r2 give the fraction of total number of
121    // electrons which lie inside such interval.*/
122    //
123   
124    Double_t R = x[0], age = x[1], Rm=par[0];
125    Double_t e1=2.0;
126    Double_t e2=4.5;
127    Double_t D=R/Rm;
128    return Gamma(e2-age)/Gamma(age)/Gamma(e2-2.*age)*Power(D,age-(e1-1.0))
129        *Power((1.0+D),age-e2)/Rm;
130}
131
132//______________________________________________________________________________
133Double_t ShowerTrackNKGhadron(Double_t *x, Double_t *par) {
134    //
135    // idem ShowerTrackNKG2, but Rm is rescale by a factor 0.5577
136    // Scaling factor comes from Gora et al. Astropart. 24 (2006)
137    // this new NKG agrees well with JNC function, see Capedevielle et al. ICRC2001
138    //
139   
140    Double_t R = x[0], age = x[1], Rm=par[0]*0.5577;
141    Double_t e1=2.0;
142    Double_t e2=4.5;
143    Double_t D=R/Rm;
144    Double_t result = Gamma(e2-age)/Gamma(age)/Gamma(e2-2.*age)*Power(D,age-(e1-1.0))
145        *Power((1.0+D),age-e2)/Rm;
146    if (IsNaN(result)) {
147        cout << "ShowerTrackNKGhadron returnes NAN. I return ZERO "  << endl;
148        return 0;
149    }
150    return result;
151}
152
153//______________________________________________________________________________
154Double_t ShowerTrackBaltru(Double_t *x, Double_t *par) {
155    //
156    // dNe/dtheta(theta,Et) from Baltrusaitis et al. J.Phys.G:Nucl. Phys. 13 (1987)
157    // where theta is the angle between the electrons and the shower axis and Et
158    // (MeV) is the energy thr for electrons considered. The integration of
159    // dNe/dtheta(theta,Et) over dtheta from 0 to pi is normalized to 1.  The
160    // integral from theta1 to theta2 give the fraction of total number of
161    // electrons which lie inside such angular interval (distribution in phi is
162    // supposed uniform).
163    //
164   
165    Double_t theta = x[0], Et = x[1];
166    Double_t a = 0.85;
167    Double_t b = 0.66;
168    Double_t theta0 = a*Power(Et,-b);
169    Double_t pigreco = Pi();
170    return Exp(-theta/theta0)/theta0/(1.-Exp(-pigreco/theta0) );
171}
172
173#ifdef __APPLE__
174//G.Barrand : on Mac the below global variables induces a crash at startup of
175//            Simu. Then we create the objects "when needed" in the
176//            ShowerTrack::Get methods. Sight, using global variables,
177//            as long as the singleton pattern, is NEVER a good idea.       
178TF2* ShowerTrack::fEnergyAgeDistributionDefault = 0;
179TF2* ShowerTrack::fEnergyAgeDistributionHillas  = 0;
180TF2* ShowerTrack::fEnergyAgeDistributionGiller  = 0;
181TF2* ShowerTrack::fEnergyAgeDistributionNerling = 0;
182TF2* ShowerTrack::fEnergyAgeDistributionMelot   = 0;
183TF2* ShowerTrack::fLateralDistributionDefalt = 0;
184TF2* ShowerTrack::fLateralDistribution1 = 0;
185TF2* ShowerTrack::fLateralDistribution2 = 0;
186TF2* ShowerTrack::fLateralDistribution3 = 0;
187TF2* ShowerTrack::fAngularDistribution1 = 0;
188#else
189// Set of static distributions. Note! Default does not mean that this is the best rather it is
190// a convinient choice for the debugging study
191TF2* ShowerTrack::fEnergyAgeDistributionDefault = new TF2("EneAgeDefault",ShowerTrackGiller,0.1,1000,0,2);
192TF2* ShowerTrack::fEnergyAgeDistributionHillas  = new TF2("hillas",ShowerTrackHillas,0.1,1000,0,2);
193TF2* ShowerTrack::fEnergyAgeDistributionGiller  = new TF2("giller",ShowerTrackGiller,0.1,1000,0,2);
194TF2* ShowerTrack::fEnergyAgeDistributionNerling = new TF2("nerling",ShowerTrackNerling,0.1,1000,0,2);
195TF2* ShowerTrack::fEnergyAgeDistributionMelot   = new TF2("melot",ShowerTrackMelot,0.1,1000,0,2);
196
197TF2* ShowerTrack::fLateralDistributionDefalt = new TF2("LateralDistributionDefalt",ShowerTrackNKG1,0.001,50,0,2);
198TF2* ShowerTrack::fLateralDistribution1 = new TF2("NKG1",ShowerTrackNKG1,0.001,50,0,2);
199TF2* ShowerTrack::fLateralDistribution2 = new TF2("NKG",ShowerTrackNKG2,0.001,5000,0,2,1);
200TF2* ShowerTrack::fLateralDistribution3 = new TF2("NKGhadron",ShowerTrackNKGhadron,0.001,5000,0,2,1);
201
202TF2* ShowerTrack::fAngularDistribution1 = new TF2("baltru",ShowerTrackBaltru,0.,Pi(),.5,1000.);
203#endif
204
205//______________________________________________________________________________
206ShowerTrack::ShowerTrack(): PhysicsData("shower"), EsafMsgSource() {
207    //
208    // Constructor
209    //
210    SetNumEnergyHistos(20);
211    SetNumLateralHistos(20);
212}
213//______________________________________________________________________________
214ShowerTrack::~ShowerTrack() {
215    //
216    // Destructor
217    //
218   Reset();
219
220}
221//______________________________________________________________________________
222TF2* ShowerTrack::GetEnergyDistribution(TString  name) {
223    //
224    // Return Pointer to the energy distribution of electrons in the shower
225    //
226
227#ifdef __APPLE__
228  //G.Barrand :
229  if(!fEnergyAgeDistributionDefault)
230    fEnergyAgeDistributionDefault = new TF2("EneAgeDefault",ShowerTrackGiller,0.1,1000,m,2);
231  if(!fEnergyAgeDistributionHillas)
232    fEnergyAgeDistributionHillas  = new TF2("hillas",ShowerTrackHillas,0.1,1000,0,2);
233  if(!fEnergyAgeDistributionGiller)
234    fEnergyAgeDistributionGiller  = new TF2("giller",ShowerTrackGiller,0.1,1000,0,2);
235  if(!fEnergyAgeDistributionNerling)
236    fEnergyAgeDistributionNerling = new TF2("nerling",ShowerTrackNerling,0.1,1000,0,2);
237  if(!fEnergyAgeDistributionMelot)
238    fEnergyAgeDistributionMelot   = new TF2("melot",ShowerTrackMelot,0.1,1000,0,2);
239#endif
240
241    if( name == "hillas" )
242        return fEnergyAgeDistributionHillas;
243    else if ( name == "giller")
244        return fEnergyAgeDistributionGiller;
245    else if ( name == "nerling")
246        return fEnergyAgeDistributionNerling;
247    else if (name == "default" || name == "")
248        return fEnergyAgeDistributionDefault;
249    else
250        throw runtime_error(("ShowerTrack does not know this parametrization: "+name).Data());
251}
252//______________________________________________________________________________
253TF2* ShowerTrack::GetLateralDistribution(TString  name) {
254    //
255    // Return Pointer to the lateral distribution of electrons in the shower
256    //
257#ifdef __APPLE__
258  //G.Barrand :
259  if(!fLateralDistributionDefalt)
260    fLateralDistributionDefalt = new TF2("LateralDistributionDefalt",ShowerTrackNKG1,0.001,50,0,2);
261  if(!fLateralDistribution1)
262    fLateralDistribution1 = new TF2("NKG1",ShowerTrackNKG1,0.001,50,0,2);
263  if(!fLateralDistribution2)
264    fLateralDistribution2 = new TF2("NKG",ShowerTrackNKG2,0.001,5000,0,2,1);
265  if(!fLateralDistribution3)
266    fLateralDistribution3 = new TF2("NKGhadron",ShowerTrackNKGhadron,0.001,5000,0,2,1);
267#endif
268
269    if( name == "NKG1" )
270        return fLateralDistribution1;
271    else if (name == "default" || name == "")
272        return fLateralDistribution1;
273    else if (name == "NKG")
274        return fLateralDistribution2;
275    else if (name == "NKGhadron")
276        return fLateralDistribution3;
277    else
278        throw runtime_error(("ShowerTrack does not know this parametrization: "+name).Data());
279}
280//______________________________________________________________________________
281TF2* ShowerTrack::GetAngularDistribution(TString  name) {
282    //
283    // Return Pointer to the angular distribution of electrons in the shower
284    //
285#ifdef __APPLE__
286  //G.Barrand :
287  if(!fAngularDistribution1)
288    fAngularDistribution1 = new TF2("baltru",ShowerTrackBaltru,0.,Pi(),.5,1000.);
289#endif
290
291    if( name == "baltru" )
292        return fAngularDistribution1;
293    else if (name == "default" || name == "")
294        return fAngularDistribution1;
295    else
296        throw runtime_error(("ShowerTrack does not know this parametrization: "+name).Data());
297}
298
299//______________________________________________________________________________
300const ShowerStep& ShowerTrack::GetStep( UInt_t i ) const {
301    //
302    // Return the i-th step of the shower
303    //
304    if ( i >= fSteps.size() )
305        throw runtime_error("Index out of range in ShowerTrack::GetStep");
306    return (fSteps[i]);
307}
308//______________________________________________________________________________
309const ShowerStep& ShowerTrack::GetLastStep( ) const {
310    //
311    // Return last step of the shower
312    //
313    return (fSteps[Size()-1]);
314}
315
316//______________________________________________________________________________
317const ShowerStep& ShowerTrack::operator[]( UInt_t i ) const {
318    //
319    // Return the i-th step using an operator
320    //
321        if ( i >= fSteps.size() )
322                throw runtime_error("Index out of range in ShowerTrack::[]");
323        return (fSteps[i]);
324}
325//______________________________________________________________________________
326const vector<ShowerStep>& ShowerTrack::GetSteps() const {
327    //
328    // Returns all steps in the track
329    //
330        return fSteps; 
331}
332
333//______________________________________________________________________________
334void ShowerTrack::AtmUpdateTrack() {
335    // Update the Track to the density profile of the atmosphere
336    // currently setted in ESAF. Leave the track direction in space unchanged,
337    // leave all the shower physics in g/cm^2 unchanged. Change the
338    // points in physical space. This is done in order to be able to handle
339    // shower track created by shower generators outside ESAF which are using
340    // a different atmospheric density profile respect to the one currently
341    // setted in ESAF.
342
343}
344//______________________________________________________________________________
345void ShowerTrack::Reset() {
346  // delete each ShowerStep object
347  if ( Clear() ) fSteps.clear();
348  // reset hit ground information
349  fHitGround = kFALSE;
350}
351//______________________________________________________________________________
352Bool_t ShowerTrack::Clear() {
353    vector<ShowerStep>::iterator step;
354    vector<ShowerStep>::iterator step2;
355
356    for (step = fSteps.begin(); step != fSteps.end(); step++) {
357        TH1F* Energy_tmp      = step->fHistoEnergy;
358        TH1F* Lateral_tmp     = step->fHistoLateral;
359        TH2F* EneAng_tmp      = step->fHistoEneAng;
360        TH2F* RadPhiEle_tmp   = step->fHistoRadPhiEle;
361        TH2F* RadDTimeEle_tmp = step->fHistoRadDTimeEle;
362        TH2F* RadPhiEloss_tmp = step->fHistoRadPhiEloss;
363        TH1F* AngCher_tmp     = step->fHistoAngCher;
364        if (Energy_tmp||Lateral_tmp||EneAng_tmp||RadPhiEle_tmp||
365            RadDTimeEle_tmp||RadPhiEloss_tmp ||AngCher_tmp ) {
366            // scan all other steps to find the same pointer (if any)
367            for (step2 = step+1; step2 != fSteps.end(); step2++) {
368                if(Energy_tmp == step2->fHistoEnergy)           step2->fHistoEnergy      = NULL;
369                if(Lateral_tmp == step2->fHistoLateral)         step2->fHistoLateral     = NULL;
370                if(EneAng_tmp == step2->fHistoEneAng)           step2->fHistoEneAng      = NULL;
371                if(RadPhiEle_tmp == step2->fHistoRadPhiEle)     step2->fHistoRadPhiEle   = NULL;
372                if(RadDTimeEle_tmp == step2->fHistoRadDTimeEle) step2->fHistoRadDTimeEle = NULL;
373                if(RadPhiEloss_tmp == step2->fHistoRadPhiEloss) step2->fHistoRadPhiEloss = NULL;
374                if(AngCher_tmp == step2->fHistoAngCher)         step2->fHistoAngCher     = NULL;
375            }
376        }
377
378        // save pointers to histos
379        SafeDelete(step->fHistoEnergy);
380        SafeDelete(step->fHistoLateral);
381        SafeDelete(step->fHistoEneAng);
382        SafeDelete(step->fHistoRadPhiEle);
383        SafeDelete(step->fHistoRadDTimeEle);
384        SafeDelete(step->fHistoRadPhiEloss);
385        SafeDelete(step->fHistoAngCher);
386    }
387    Int_t i(0);
388    for (step = fSteps.begin(); step != fSteps.end(); step++) {
389        if (step->fHistoEnergy)      Msg(EsafMsg::Warning) << "Energy histo not deleted " << i << MsgDispatch;
390        if (step->fHistoLateral)     Msg(EsafMsg::Warning) << "Lateral histo not deleted " << i << MsgDispatch;
391        if (step->fHistoEneAng)      Msg(EsafMsg::Warning) << "EneAng histo not deleted " << i << MsgDispatch;
392        if (step->fHistoRadDTimeEle) Msg(EsafMsg::Warning) << "RadTimeEle histo not deleted " << i << MsgDispatch;
393        if (step->fHistoRadPhiEloss) Msg(EsafMsg::Warning) << "RadPhiEloss histo not deleted " << i << MsgDispatch;
394        if (step->fHistoAngCher)     Msg(EsafMsg::Warning) << "AngCher histo not deleted " << i << MsgDispatch;
395        i++;
396    }
397    return kTRUE;
398}
399
400//______________________________________________________________________________
401void ShowerTrack::Add( ShowerStep step ) {
402    step.SetStepID(Size() + 1);
403    step.SetParentTrack(this); 
404    fSteps.push_back(step);
405}
406//______________________________________________________________________________
407void ShowerTrack::DrawXYZ(Option_t *opt) {
408    TString option(opt);
409    Int_t n = GetNumStep();
410
411    if (n == 0) {
412        Msg(EsafMsg::Debug) << "DrawXYZ() informs: The track is empty. Check why" << MsgDispatch;
413        return;
414    }
415   
416    Double_t kHuge = 1.e20;
417    Double_t Xmin(kHuge), Xmax(-kHuge), Ymin(kHuge), Ymax(-kHuge), Zmin(kHuge), Zmax(-kHuge), Nmax(1);
418    Double_t Threshold(1.e-2);
419   
420    for (Int_t i=0; i<n; i++) {
421        Double_t xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
422        Double_t ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
423        Double_t zmean = (GetStep(i).GetXYZi().Z() + GetStep(i).GetXYZf().Z())/2/km;
424        Double_t ne    =  GetStep(i).GetNelectrons();
425        if (ne/Nmax > Threshold) {
426            Xmin = (Xmin < xmean ? Xmin : xmean); 
427            Xmax = (Xmax > xmean ? Xmax : xmean); 
428            Ymin = (Ymin < ymean ? Ymin : ymean); 
429            Ymax = (Ymax > ymean ? Ymax : ymean); 
430            Zmin = (Zmin < zmean ? Zmin : zmean); 
431            Zmax = (Zmax > zmean ? Zmax : zmean); 
432        }
433        if (Nmax < GetStep(i).GetNelectrons() )
434            Nmax = GetStep(i).GetNelectrons();
435    }
436
437    // build the histogram
438    TString name = "ShowerTrackXYZ";
439    TH3F* th = (TH3F*)gROOT->FindObject(name);
440    SafeDelete(th);   
441   
442    if (!th) {
443        Double_t OffSet = 1.1;
444        Double_t x1(-5),x2(5),y1(-5),y2(5),z1(0),z2(30); // make 10x10x30 kms box for default
445        if (Xmin < x1) x1 = Xmin*OffSet;
446        if (Xmax > x2) x2 = Xmax*OffSet;
447        if (Ymin < y1) y1 = Ymin*OffSet;
448        if (Ymax > y2) y2 = Ymax*OffSet;       
449        Int_t xbins = Int_t((x2-x1)/1);  // make roughly 1 km bins size
450        Int_t ybins = Int_t((y2-y1)/1);  // make roughly 1 km bins size
451        th = new TH3F( name, "XYZ development", xbins,x1,x2,ybins,y1,y2,30,z1,z2); 
452        th->SetStats(0);
453        th->SetXTitle("X [km]");
454        th->SetYTitle("Y [km]");
455        th->SetZTitle("Z [km]");
456    }
457    for (Int_t i=0; i<n; i++) {
458        double xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
459        double ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
460        double zmean = (GetStep(i).GetXYZi().Z() + GetStep(i).GetXYZf().Z())/2/km;
461        double Ne    = GetStep(i).GetNelectrons();
462        th->Fill(xmean,ymean,zmean,1.+Ne/Nmax);
463    }
464    th->Draw(option);
465}
466
467//______________________________________________________________________________
468void ShowerTrack::DrawXY(Option_t *opt) {
469    TString option(opt);
470
471    Int_t n = GetNumStep();
472
473    if (n == 0) {
474        Msg(EsafMsg::Debug) << "DrawXY() informs: The track is empty. Check why" << MsgDispatch;
475        return;
476    }
477   
478    Double_t kHuge = 1.e20;
479    Double_t Xmin(kHuge), Xmax(-kHuge), Ymin(kHuge), Ymax(-kHuge),Nmax(0);
480    Double_t Threshold(1.e-2);
481   
482    for (Int_t i=0; i<n; i++) {
483        Double_t xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
484        Double_t ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
485        Double_t ne    =  GetStep(i).GetNelectrons();
486        if (ne/Nmax > Threshold) {
487            Xmin = (Xmin < xmean ? Xmin : xmean);
488            Xmax = (Xmax > xmean ? Xmax : xmean); 
489            Ymin = (Ymin < ymean ? Ymin : ymean); 
490            Ymax = (Ymax > ymean ? Ymax : ymean); 
491        }
492        if (Nmax < GetStep(i).GetNelectrons() )
493            Nmax = GetStep(i).GetNelectrons();
494    }
495    // build the histogram
496    TString name = "ShowerTrackXY";
497    TH2F* th = (TH2F*)gROOT->FindObject(name);
498    SafeDelete(th);   
499   
500    if (!th) {
501        Double_t OffSet = 1.1;
502        Double_t x1(-5),x2(5),y1(-5),y2(5); // make 10x10 kms box for default
503        if (Xmin < x1) x1 = Xmin*OffSet;
504        if (Xmax > x2) x2 = Xmax*OffSet;
505        if (Ymin < y1) y1 = Ymin*OffSet;
506        if (Ymax > y2) y2 = Ymax*OffSet; 
507        Int_t xbins = Int_t((x2-x1)/1);  // make roughly 1 km bins size
508        Int_t ybins = Int_t((y2-y1)/1);  // make roughly 1 km bins size
509             
510        th = new TH2F( name, "XY development", xbins,x1,x2,ybins,y1,y2); 
511        th->SetStats(0);
512        th->SetXTitle("X [km]");
513        th->SetYTitle("Y [km]");
514    }
515    for (Int_t i=0; i<n; i++) {
516        double xmean = (GetStep(i).GetXYZi().X() + GetStep(i).GetXYZf().X())/2/km;
517        double ymean = (GetStep(i).GetXYZi().Y() + GetStep(i).GetXYZf().Y())/2/km;
518        double Ne    = GetStep(i).GetNelectrons();
519        th->Fill(xmean,ymean,Ne);
520    }
521    th->Draw(option);
522}
523
524//______________________________________________________________________________
525EarthVector ShowerTrack::FirstPos() const {
526    //
527    // return the entry position of the first ShowerStep
528    //
529    EarthVector rtn(0,0,HUGE);
530    if(fSteps.size()) {
531        const EVector& entry = fSteps[0].GetXYZi();
532        rtn.SetXYZ(entry.X(),entry.Y(),entry.Z());
533    }
534   
535    return rtn;
536}
537
538//______________________________________________________________________________
539EarthVector ShowerTrack::LastPos() const {
540    //
541    // return the exit position of the last ShowerStep
542    //
543    EarthVector rtn(0,0,HUGE);
544    size_t last = fSteps.size();
545    if(last) {
546        const EVector& exit = fSteps[last-1].GetXYZf();
547        rtn.SetXYZ(exit.X(),exit.Y(),exit.Z());
548    }
549   
550    return rtn;
551}
552
553//______________________________________________________________________________
554Bool_t ShowerTrack::CheckTrack() {
555    //
556    // Check that the track does not contain Nans
557    //
558    vector <ShowerStep>::iterator it;
559    for (it = fSteps.begin(); it != fSteps.end(); ) {
560         // begin to check the content of the step
561         Bool_t problem = kFALSE;
562         EVector posi = (*it).GetXYZi();
563         EVector posf = (*it).GetXYZf();
564         Double_t nel = (*it).GetNelectrons();
565         Double_t ti  = (*it).GetTimei();
566         Double_t tf  = (*it).GetTimef();
567         if (TMath::IsNaN(posi.X()) || TMath::IsNaN(posi.Y()) || TMath::IsNaN(posi.Z())) problem = kTRUE;
568         if (TMath::IsNaN(posf.X()) || TMath::IsNaN(posf.Y()) || TMath::IsNaN(posf.Z())) problem = kTRUE;
569         if (TMath::IsNaN(nel) || TMath::IsNaN(ti) || TMath::IsNaN(tf))                  problem = kTRUE;
570         if (problem) {
571             fSteps.erase(it);
572             Msg(EsafMsg::Warning) << "A step with NAN is found. Erased. "  << MsgDispatch;
573         }
574         else it++;
575    }
576    return kTRUE;
577}
578
Note: See TracBrowser for help on using the repository browser.