source: JEM-EUSO/esaf_cc_at_lal/packages/common/eventviewer/src/EShowerHistoPainter.cc @ 114

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

actual version of ESAF at CCin2p3
File size: 9.8 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: EShowerHistoPainter.cc 2819 2009-03-13 11:48:16Z naumov $
3// Author: Naumov Dmitry   Jul, 26 2004
4
5#include "EShowerHistoPainter.hh"
6#include "EShower.hh"
7#include "EShowerStep.hh"
8#include "ESystemOfUnits.hh"
9#include "EGeometry.hh"
10#include "EDetector.hh"
11#include "EConst.hh"
12
13#include <TTimer.h>
14#include <TPad.h>
15#include <TH1F.h>
16#include <TPaveText.h>
17#include <TMath.h>
18#include <TList.h>
19
20ClassImp(EShowerHistoPainter)
21
22using namespace sou;
23using namespace EConst;
24
25// TOOL function : give local altitude
26Double_t Zv(const TVector3& v) {
27    return sqrt( pow(6.370949e6*m + v.Z(),2) + v.Perp2() ) - 6.370949e6*m;
28}
29
30
31//_____________________________________________________________________________
32EShowerHistoPainter::EShowerHistoPainter( EShower* sh, const EGeometry* geom, EDetector* det) {
33    //
34    // Constructor
35    //
36
37    fShower = sh;
38    fGeometry = geom;
39    fDetector = det;
40    if (!fShower || !fGeometry || !fDetector) 
41       {Error("EShowerHistoPainter","must be defined all three objects: shower, geometry and detector");}
42    fHistos = new TList();
43
44    fTimeBins = 0;
45    fAltBins = 0;
46    fNumFrames = 100;
47    fFrame = -1;
48
49    Build();
50}
51
52//_____________________________________________________________________________
53EShowerHistoPainter::~EShowerHistoPainter() {
54    //
55    // Destructor
56    //
57   
58    if ( fTimeBins ) delete [] fTimeBins;
59    if ( fAltBins ) delete [] fAltBins;
60    if ( fHistos ) fHistos->Delete();
61
62}
63
64//_____________________________________________________________________________
65void EShowerHistoPainter::Build() {
66    //
67    // Build
68    //
69   
70    if ( !fShower ) {
71        Info("Build()","EShower object is NULL. Painter made zombie.");
72        MakeZombie();
73        return;
74    }
75
76    fNeMax = 0;
77    fAgeMax = 0;
78    fDepthMax = 0;
79   
80    Int_t n = fShower->GetNumSteps();
81    if ( n <= 0 ) {
82        Info("Build()","Empty Shower ( NumSteps = 0 ). Painter made zombie.");
83        MakeZombie();
84        return;
85    }
86
87    fTimeBins = new Double_t[n+1];
88    fAltBins = new Double_t[n+1];
89
90    for (Int_t i=0; i<n; i++) {
91
92        fTimeBins[i] = fShower->GetStep(i)->GetTimei()/microsecond;
93        fTimeBins[i+1] = fShower->GetStep(i)->GetTimef()/microsecond;
94
95//        if ( i > 0 &&  fTimeBins[i-1] > fTimeBins[i] ) {
96//            Printf("Porcapupazza i = %d, fTimeBins[i-1] = %f, fTimeBins[i] = %f",fTimeBins[i-1], fTimeBins[i]);
97//        }
98
99        if( fNeMax < fShower->GetStep(i)->GetNumElectrons() )
100            fNeMax = fShower->GetStep(i)->GetNumElectrons();
101        if ( fAgeMax <  0.5*(fShower->GetStep(i)->GetAgei()
102                    + fShower->GetStep(i)->GetAgef()))
103            fAgeMax = 0.5*(fShower->GetStep(i)->GetAgei()
104                    + fShower->GetStep(i)->GetAgef());
105        if ( fDepthMax < (fShower->GetStep(i)->GetXi()
106                    + fShower->GetStep(i)->GetXf())/2/g*cm2)
107            fDepthMax = (fShower->GetStep(i)->GetXi()
108                    + fShower->GetStep(i)->GetXf())/2/g*cm2;
109       
110        fAltBins[n - i] = Zv(fShower->GetStep(i)->GetPosi())/km;
111        fAltBins[n-1 - i] = Zv(fShower->GetStep(i)->GetPosf())/km;
112    }
113   
114   
115    {
116      EShowerStep *s = fShower->GetStep(0);
117      MinGtu = s->GetTimei()+(fGeometry->GetPos()-s->GetPosi()).Mag()/Clight();
118      MinGtu+= s->GetTimef()+(fGeometry->GetPos()-s->GetPosf()).Mag()/Clight();
119      MinGtu = 0.5*MinGtu/microsecond; // in microsecond
120
121    }
122
123    {
124      EShowerStep *s = fShower->GetStep(fShower->GetNumSteps()-1);
125      MaxGtu = s->GetTimei()+(fGeometry->GetPos()-s->GetPosi()).Mag()/Clight();
126      MaxGtu+= s->GetTimef()+(fGeometry->GetPos()-s->GetPosf()).Mag()/Clight();
127      MaxGtu = 0.5*MaxGtu/microsecond; // in microsecond
128    }
129
130    Float_t GtuStep = fDetector->GetElectronics()->GetGtuLength()*1.e-3;
131    GtuBins = (Int_t)Floor((MaxGtu - MinGtu)/GtuStep);
132
133    // build histograms
134    TH1F *th0 = GetTimeHisto("Shower_Ne", "Ne vs time");
135    TH1F *th1 = GetTimeHisto("Shower_Age", "Age vs time");
136    TH1F *th2 = GetTimeHisto("Shower_Depth", "Depth vs time");
137    th0->SetYTitle("Ne");
138    th1->SetYTitle("Age");
139    th2->SetYTitle("Slant depth");
140   
141    TH1F *alt0 = GetAltHisto("Shower_Ne_alt", "Ne vs altitude");
142    alt0->SetYTitle("Ne");
143   
144    // fill same information but versus gtu
145    TH1F* gtu0 = GetGtuHisto("Shower_Ne_GTU","Ne vs GTU");
146    gtu0->SetYTitle("Ne");
147
148    for (Int_t i=0; i<n; i++) {
149        EShowerStep* s = fShower->GetStep(i);
150        Double_t Ne    = s->GetNumElectrons();
151        Double_t Age   = (s->GetAgei()+s->GetAgef())/2;
152        Double_t Depth = (s->GetXi()+s->GetXf())/2/g*cm2;
153        Double_t alt = (Zv(s->GetPosi())+ Zv(s->GetPosf()))/2/km; 
154
155        TVector3 pos  = 0.5*(s->GetPosi() + s->GetPosf());
156        Double_t time  = (s->GetTimei()+s->GetTimef())/2;
157        TVector3 dir  = fGeometry->GetPos() - pos;
158        Float_t  EusoTimei = (s->GetTimei() + (fGeometry->GetPos() - s->GetPosi()).Mag()/Clight())/microsecond - MinGtu; // in microsecond
159        Float_t  EusoTimef = (s->GetTimef() + (fGeometry->GetPos() - s->GetPosi()).Mag()/Clight())/microsecond - MinGtu; // in microsecond
160        Float_t  EusoTime = 0.5*(EusoTimei+EusoTimef);
161
162        th0->Fill(time/microsecond,Ne);
163        th1->Fill(time/microsecond,Age);
164        th2->Fill(time/microsecond,Depth);
165       
166        alt0->Fill(alt,Ne);
167
168        gtu0->Fill(EusoTime,Ne);
169    }
170
171
172
173   
174
175}
176
177//_____________________________________________________________________________
178TH1F *EShowerHistoPainter::GetTimeHisto(const char *name, const char *title) {
179    //
180    //
181    //
182
183    if ( title == 0 ) title = name;
184
185    TH1F* th = (TH1F*)fHistos->FindObject(name);
186   
187    if (!th) {
188        th = new TH1F( name, title, fShower->GetNumSteps(), fTimeBins);
189        fHistos->Add(th);
190        th->SetStats(0);
191        th->SetXTitle("Development Time, #{mu}s");
192        th->SetFillColor(9);
193        th->SetDirectory(0);
194    }
195
196    return th;
197}
198
199//_____________________________________________________________________________
200TH1F *EShowerHistoPainter::GetGtuHisto(const char *name, const char *title) {
201    //
202    //
203    //
204
205    if ( title == 0 ) title = name;
206
207    TH1F* th = (TH1F*)fHistos->FindObject(name);
208   
209    if (!th) {
210        th = new TH1F( name, title, GtuBins, 0, MaxGtu-MinGtu);
211        fHistos->Add(th);
212        th->SetStats(0);
213        th->SetXTitle("GTU Time, #{mu}s");
214        th->SetFillColor(9);
215        th->SetDirectory(0);
216    }
217
218    return th;
219}
220
221//_____________________________________________________________________________
222TH1F *EShowerHistoPainter::GetAltHisto(const char *name, const char *title) {
223    //
224    //
225    //
226
227    if ( title == 0 ) title = name;
228
229    TH1F* th = (TH1F*)fHistos->FindObject(name);
230   
231    if (!th) {
232        th = new TH1F( name, title, fShower->GetNumSteps(), fAltBins);
233        fHistos->Add(th);
234        th->SetStats(0);
235        th->SetXTitle("Altitude (km)");
236        th->SetFillColor(9);
237        th->SetDirectory(0);
238    }
239
240    return th;
241}
242
243//______________________________________________________________________________
244void EShowerHistoPainter::Draw(Option_t *opt) {
245    //
246    //
247    //
248   
249    if ( IsZombie() ) return;
250
251    TString option(opt);
252    TH1F* th(0);
253
254    if ( option == "" ) option = "Ne";
255   
256    if(option == "Ne_alt") {
257        th = GetAltHisto("Shower_Ne_alt");
258        if(th) th->Draw();
259    }
260
261    if ( option == "Ne" ) {
262        th = GetTimeHisto("Shower_Ne");
263        if ( th ) th->Draw();
264    } 
265    else if ( option == "Age" ) {
266        th = GetTimeHisto("Shower_Age");
267        if ( th ) th->Draw();
268           
269    } 
270    else if ( option == "Depth" ) {
271        th = GetTimeHisto("Shower_Depth");
272        if ( th ) th->Draw();           
273    }
274    else if (option == "Truth") {
275        TPaveText* hdemo = new TPaveText(.05,.05,.95,.95);
276        hdemo->SetTextAlign(12);
277        hdemo->SetTextFont(52);
278        hdemo->SetTextColor(4);
279        hdemo->AddText("");
280        hdemo->AddText("This event was simulated with:");
281        hdemo->AddText("");
282        hdemo->AddText(Form("Energy %.4e MeV",fShower->GetEnergy()));
283        hdemo->AddText(Form("Theta  %.2f    deg",fShower->GetTheta()*TMath::RadToDeg()));
284        hdemo->AddText(Form("Phi    %.2f    deg",fShower->GetPhi()*TMath::RadToDeg()));
285        hdemo->AddText(Form("X1     %.2f    g/cm^{2}",fShower->GetX1()*cm2/g));
286        hdemo->AddText(Form("Xinit  %.2f    km",fShower->GetInitPos().X()/km));
287        hdemo->AddText(Form("Yinit  %.2f    km",fShower->GetInitPos().Y()/km));
288        hdemo->AddText(Form("Zinit  %.2f    km",fShower->GetInitPos().Z()/km));
289        hdemo->AddText("");
290        hdemo->Draw();
291    }
292}
293
294//_____________________________________________________________________________
295void EShowerHistoPainter::NextFrame() {
296    //
297    // Switch to the next frame and update histograms
298    //
299
300    fFrame++;
301    UpdateHistos();
302}
303
304//______________________________________________________________________________
305void EShowerHistoPainter::UpdateHistos() {
306    //
307    // Update hitograms at the current frame.
308    //
309
310    if ( IsZombie() ) return;
311
312    Int_t n;
313    if ( fFrame == -1 )
314        n = fShower->GetNumSteps();
315    else
316        n = (Int_t)(((Float_t)fShower->GetNumSteps()/fNumFrames)*fFrame);
317
318
319    TH1F* thNe    = GetTimeHisto("Shower_Ne");
320    TH1F* thAge   = GetTimeHisto("Shower_Age");
321    TH1F* thDepth = GetTimeHisto("Shower_Depth");
322
323    for (Int_t i=0; i<fShower->GetNumSteps() ; i++) {
324        Double_t Ne = fShower->GetStep(i)->GetNumElectrons();
325        Double_t Age  = (fShower->GetStep(i)->GetAgei()
326                       + fShower->GetStep(i)->GetAgef())/2;
327        Double_t Depth  = (fShower->GetStep(i)->GetXi()
328                         + fShower->GetStep(i)->GetXf())/2/g*cm2;
329
330
331        thNe->SetBinContent(i+1, i < n ? Ne : 0 );
332        thAge->SetBinContent(i+1, i < n ? Age : 0 );
333        thDepth->SetBinContent(i+1, i < n ? Depth : 0 );
334    }
335
336    thNe->SetMaximum( fNeMax*1.1 );
337    thAge->SetMaximum( fAgeMax*1.1 );
338    thDepth->SetMaximum( fDepthMax*1.1 );
339
340}
Note: See TracBrowser for help on using the repository browser.