source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/reconstruction/modules/shower/energy/src/RecoShowerTrack.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: 15.4 KB
Line 
1// $Id: RecoShowerTrack.cc 2927 2011-06-16 18:02:41Z mabl $
2// Author: Dmitry.Naumov   2005/10/11
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: RecoShowerTrack                                                      *
8 *  Package: Energy                                                          *
9 *  Coordinator: Dmitry.Naumov                                               *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// RecoShowerTrack
16//
17// <extensive class description>
18//
19//   Config file parameters
20//   ======================
21//
22//   <parameter name>: <parameter description>
23//   -Valid options: <available options>
24//
25
26#include "RecoShowerTrack.hh"
27#include "EConst.hh"
28#include "SlastShowerGenerator.hh"
29#include "ShowerTrack.hh"
30#include "FluoCalculator.hh"
31#include "Config.hh"
32#include "EsafSpectrum.hh"
33#include "LightSourceFactory.hh"
34#include "Ground.hh"
35#include "O1_ClearSkyPropagator.hh"
36#include "DetectorGeometry.hh"
37#include "RadiativeFactory.hh"
38#include "EnergyModule.hh"
39#include "RecoEvent.hh"
40#include "EEvent.hh"
41#include "ETruth.hh"
42#include "Atmosphere.hh"
43#include "EGeometry.hh"
44
45#include <TMath.h>
46#include <TCanvas.h>
47#include <TProfile.h>
48#include <TList.h>
49#include <TPolyMarker3D.h>
50
51using namespace sou;
52using namespace TMath;
53using namespace EConst;
54
55ClassImp(RecoShowerTrack)
56
57// constructors, destructor
58//_____________________________________________________________________________
59RecoShowerTrack::RecoShowerTrack(): EsafMsgSource() {
60    //
61    // Constructor
62    //
63    Constructor();
64}
65
66//_____________________________________________________________________________
67RecoShowerTrack::RecoShowerTrack(Float_t t, Float_t p, Float_t h, Float_t tfov, Float_t pfov): EsafMsgSource() {
68    //
69    // Constructor
70    //
71    Constructor();
72    SetThetaPhi(t,p);
73    SetThetaPhiFoV(tfov,pfov);
74    SetHmax(h);
75   
76}
77//_____________________________________________________________________________
78void RecoShowerTrack::Constructor() {
79    //
80    // This is default ctor
81    fHmax = -100;
82    fShowerDir.SetXYZ(0,0,0);
83    fEusoDir.SetXYZ(0,0,0);
84    string fluoname         = Config::Get()->GetCF("Reco","EnergyModule")->GetStr("EnergyModule.fFluoCalculator");
85    fFluoCalculator         = LightSourceFactory::Get()->GetFluoCalculator( fluoname );   
86    fEusoAlt                = Config::Get()->GetCF("General","Euso")->GetNum("Euso.fAltitude")*km;
87    fEnergyDistributionName = Config::Get()->GetCF("Reco","EnergyModule")->GetStr("EnergyModule.fEnergyDistributionName");
88    fEsafSpectrum    = new EsafSpectrum();
89    fShowerGenerator = new SlastShowerGenerator(kTRUE);
90    fShowerGenerator->SetQuiet();
91    fDetectorGeometry = new DetectorGeometry();
92    fGround             = RadiativeFactory::Get()->GetGround();
93    fClearSkyPropagator = new O1_ClearSkyPropagator(fGround);
94    fEnergyModule = NULL;
95    fEnergy = 1.e21;
96    fDebugList = new TList;
97
98        ConfigFileParser* pConf = Config::Get()->GetCF("Electronics","MacroCell");
99        fGTU = pConf->GetNum("MacroCell.fGtuTimeLength");
100        fGTU/=1000.;
101
102}
103
104//_____________________________________________________________________________
105RecoShowerTrack::~RecoShowerTrack() {
106    //
107    // Destructor
108    //
109    SafeDelete(fShowerGenerator);
110    SafeDelete(fDetectorGeometry);
111    SafeDelete(fFluoCalculator);
112    SafeDelete(fClearSkyPropagator);
113    SafeDelete(fDebugList);
114}
115//_____________________________________________________________________________
116void RecoShowerTrack::Clear() {
117    //
118    //
119    //
120    for (UInt_t i(0); i < fSteps.size(); i++) {
121        RecoShowerStep& rs = fSteps[i];   
122        rs.fSpectrum->Delete(); 
123    }
124    fSteps.clear();
125    for (Int_t i(0); i<fDebugList->GetEntries(); i++) {
126         TObject* h = (TObject*)fDebugList->At(i);
127         SafeDelete(h);
128    }
129    fDebugList->Clear();
130    fDebugList = new TList;
131   
132}
133
134//______________________________________________________________________________
135const RecoShowerStep& RecoShowerTrack::GetStepThetaPhi(Float_t Theta, Float_t Phi) const {
136    //
137    // Return the step of the shower corresponding to the given time
138
139    Int_t index = 0;
140    Double_t mmin = -1.;
141    TVector3 nFoV;
142
143    nFoV.SetMagThetaPhi(1.,Theta,Phi);
144
145    for (UInt_t i(0); i < fSteps.size(); i++) {
146      const RecoShowerStep& rs = fSteps[i];
147      EarthVector dir  = (rs.fDir).Unit();
148      Double_t mag = dir*nFoV;
149      if (mag > mmin) {mmin = mag; index = i;}
150    }
151
152    return (fSteps[index]);
153}
154
155//_____________________________________________________________________________
156Bool_t RecoShowerTrack::MakeTrack() {
157    //
158    //
159    //
160    EEvent  *SimuEvent  = fEnergyModule->GetRecoEvent()->GetSimuEvent();
161    Double_t radius   = SimuEvent->GetGeometry()->GetOpticsRadius();;
162    fDetectorGeometry->SetRadius(radius);
163    fDetectorGeometry->SetPos(SimuEvent->GetGeometry()->GetPos());
164    fClearSkyPropagator->CopyDetectorGeometry(fDetectorGeometry,false);
165
166    Int_t progress = 0;
167    Clear();
168    if (fHmax == -100 || !fShowerDir.Mag() || !fEusoDir.Mag()) {
169        if (fHmax == -100)     Msg(EsafMsg::Warning) << "Hmax is not setup\n" << MsgDispatch;
170        if (!fShowerDir.Mag()) Msg(EsafMsg::Warning) << "Shower Direction is not setup\n" << MsgDispatch;
171        if (!fEusoDir.Mag())   Msg(EsafMsg::Warning) << "Shower Direction is not setup\n" << MsgDispatch;
172                               Msg(EsafMsg::Panic)   << "Can not do a shower\n" << MsgDispatch;
173    }
174
175    // find 3D point of shower maximum
176
177    Float_t AB = EarthRadius() + fEusoAlt;
178    Float_t AC = EarthRadius() + fHmax;
179    Float_t cs = Cos(fEusoDir.Theta());
180    Float_t sn = Sin(fEusoDir.Theta());
181    Float_t z = fEusoAlt - AB*cs*cs + cs*Sqrt(AC*AC - AB*AB*sn*sn);
182    Float_t rho = (fEusoAlt-z)*Tan(fEusoDir.Theta());
183    Float_t x   = -rho*Cos(fEusoDir.Phi());
184    Float_t y   = -rho*Sin(fEusoDir.Phi());
185
186    Float_t tMin = 0, tMax = 0;
187
188    EarthVector MaximumPosition(x,y,z);
189    fDepthMax = Atmosphere::Get()->Grammage(MaximumPosition,fShowerDir,"dir");
190    EarthVector ISS(0,0,fEusoAlt);
191    EarthVector fInitPos;
192    // estimate for the current energy and assuming GIL parameterization the depth at the shower maximum (should be around 850 g/cm2)
193    // and jump from the point of the shower maximum along the track direction for that in the atmosphere which is far by the shower
194    // maximum depth
195    Double_t depthToMaximum = fShowerGenerator->GetXmax();
196    Int_t status = Atmosphere::Get()->InvertGrammage(MaximumPosition,fShowerDir,depthToMaximum,fInitPos,101*km);
197    Msg(EsafMsg::Debug) << "<GetFirstPoint> InvertGrammage() status is " << status << MsgDispatch;
198
199    fShowerGenerator->SetShower(fEnergy,fShowerDir.Theta(),fShowerDir.Phi(),fInitPos);
200    fShowerGenerator->SetQuiet(kTRUE);
201    fShowerGenerator->SetDepthStep(0.2);
202    ShowerTrack *fShowerTrack = (ShowerTrack*)fShowerGenerator->Get();
203
204    if (!fShowerTrack->Size()) return kFALSE; // do not process further if shower track object is empty
205
206    // Define tMin and tMax of shower evolution
207
208    {
209      const ShowerStep& s = fShowerTrack->GetStep(0);
210      Float_t     time = 0.5*(s.GetTimei() + s.GetTimef());
211      EarthVector dir  = (ISS - 0.5*(s.GetXYZi() + s.GetXYZf()));
212      tMin = (time + dir.Mag()/Clight())*1.e-3; // in microsecond
213    }
214
215    {
216     const ShowerStep& s = fShowerTrack->GetStep(fShowerTrack->Size()-1);
217      Float_t     time = 0.5*(s.GetTimei() + s.GetTimef());
218      EarthVector dir  = (ISS - 0.5*(s.GetXYZi() + s.GetXYZf()));
219      tMax = (time + dir.Mag()/Clight())*1.e-3; // in microsecond
220    }
221
222
223    Float_t GtuStep = fGTU;
224    if (fEnergyModule) 
225      GtuStep=fEnergyModule->GetRecoEvent()->GetHeader().GetGtuLength()*1.e-3;
226
227    Int_t Nbins = (Int_t)Floor((tMax - tMin)/GtuStep);
228
229    TH1F* h[10];
230    for (Int_t i(0); i < 10; i++) {
231         TString title = Form("EnergyModule_Step%d",i);
232         TH1F* hh = NULL;
233         if ((hh = (TH1F*)gDirectory->FindObject(title))) hh->Delete();
234         h[i] = new TH1F(title,"",Nbins,0,tMax - tMin);
235    }
236
237    // fill the track, find new Hmax
238    Double_t newHmax(-1), fluxMax(0);
239    for (UInt_t i(0); i < fShowerTrack->Size(); i++) {
240      Float_t age = 0, TotalYield = 0, Attenuation = 0;
241      const ShowerStep& s = fShowerTrack->GetStep(i);
242
243      // get position and compute for it: fluo yield + attenuation + time at EUSO
244      EarthVector pos  = 0.5*(s.GetXYZi() + s.GetXYZf());
245      Float_t     time = 0.5*(s.GetTimei() + s.GetTimef());
246      Float_t     X    = 0.5*(s.GetXi() + s.GetXf());
247      EarthVector dir  = ISS - pos;
248      Float_t     EusoTime = (time + dir.Mag()/Clight())*1.e-3 - tMin; // in microsecond
249
250      if (!(i % 50)) {
251        age = 0.5*(s.GetAgei() + s.GetAgef());
252        TF2* EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
253        TF12 EnergyDistributionProj("EDS_name",EnergyDistribution,age,"X");
254        Double_t temp = Atmosphere::Get()->Temperature(pos.Zv());
255        TotalYield  = 0;
256        Attenuation = 0;
257        if (temp<=0) continue;
258        TotalYield = fFluoCalculator->GetFluoYield(pos.Zv(),&EnergyDistributionProj,fEsafSpectrum);
259        Attenuation=fClearSkyPropagator->GoToDetector(*fEsafSpectrum,pos);
260      }
261
262      Float_t  Length = (s.GetXYZi() - s.GetXYZf()).Mag();
263      Double_t SolidAngle = 1./(4.*Pi()*dir.Mag2());
264      Double_t flux  = SolidAngle*s.GetNelectrons()*Attenuation*TotalYield*Length;
265
266      if (flux > fluxMax) {
267          newHmax = pos.Zv();
268          fluxMax = flux;
269      } 
270
271      h[0]->Fill(EusoTime,s.GetNelectrons());
272      h[1]->Fill(EusoTime,flux);
273      h[2]->Fill(EusoTime,pos.X()*flux);
274      h[3]->Fill(EusoTime,pos.Y()*flux);
275      h[4]->Fill(EusoTime,pos.Z()*flux);
276      h[5]->Fill(EusoTime,EusoTime*flux);
277      h[6]->Fill(EusoTime,TotalYield*flux);
278      h[7]->Fill(EusoTime,Attenuation*flux);
279      h[8]->Fill(EusoTime,age*flux);
280      h[9]->Fill(EusoTime,X*flux);
281     
282      if (100.*(i+1)/fShowerTrack->Size() >= progress) {
283        Msg(EsafMsg::Info).SetProgress(progress);
284        Msg(EsafMsg::Info) << "MakeTrack():" << MsgCount;
285        progress += 10;
286      }
287    }
288   
289    if (progress <= 100) Msg(EsafMsg::Info) << MsgDispatch;
290
291    for (Int_t i = 1; i <= Nbins; i++) {
292         Double_t flux   = h[1]->GetBinContent(i);
293         if (!flux) continue;
294         Double_t ne     = h[0]->GetBinContent(i);
295         Double_t x      = h[2]->GetBinContent(i)/flux;
296         Double_t y      = h[3]->GetBinContent(i)/flux;
297         Double_t z      = h[4]->GetBinContent(i)/flux;
298         Double_t time   = h[5]->GetBinContent(i)/flux;
299         Double_t yield  = h[6]->GetBinContent(i)/flux;
300         Double_t atten  = h[7]->GetBinContent(i)/flux;
301         Double_t age    = h[8]->GetBinContent(i)/flux;
302         Double_t X      = h[9]->GetBinContent(i)/flux;
303         EarthVector pos(x,y,z);
304
305         EarthVector dir = ISS - pos;
306         TVector3 nFoV;
307         nFoV.SetMagThetaPhi(1.,dir.Theta(),dir.Phi());
308         Float_t length = 1.*Clight()*GtuStep*1.e+3/(1.+fShowerDir*nFoV);
309
310         RecoShowerStep rs;
311         rs.fEusoTime    = time;
312         rs.fPos         = pos;
313         rs.fDir         = dir;
314         rs.fFluoYield   = yield;
315         rs.fAttenuation = atten;
316         rs.fNe          = ne;
317         rs.fFlux        = flux;
318         rs.fLength      = length;
319         rs.fX           = X*cm2/g;
320
321         TF2 *EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
322         TF12 EnergyDistributionProj("EDS_name",EnergyDistribution,age,"X");
323         Double_t temp = Atmosphere::Get()->Temperature(pos.Zv());
324         if (temp<=0) continue;
325         fFluoCalculator->GetFluoYield(pos.Zv(),&EnergyDistributionProj,
326                                       fEsafSpectrum);
327         fClearSkyPropagator->GoToDetector(*fEsafSpectrum,pos.Zv());
328
329         rs.fSpectrum = new TH1F(Form("RecoShowerTrackSpectrum%d",i),"",
330                                 100,300*nm,400*nm);
331
332         for (Int_t nt = 0; nt < 1000; nt++) 
333           rs.fSpectrum->Fill(fEsafSpectrum->GetLambda());
334
335         Add(rs);
336    }
337
338    for (Int_t i(0); i < 10; i++) h[i]->Delete();
339
340    fShowerTrack->Reset();
341    fShowerGenerator->Reset();
342
343    Msg(EsafMsg::Debug) << "Hmax set up " << fHmax/km << " Hmax found " << newHmax/km << MsgDispatch;
344    fHmax = newHmax;
345
346    return kTRUE;
347}
348
349//_____________________________________________________________________________
350void RecoShowerTrack::DoDebugHistos() {
351    //
352    //
353    //
354    Float_t tMin, tMax;
355
356    // Set tMin and tMax
357
358    {
359      RecoShowerStep& rs = fSteps[0];
360      tMin = rs.fEusoTime;
361    }
362
363    {
364      RecoShowerStep& rs = fSteps[fSteps.size()-1];
365      tMax = rs.fEusoTime;
366    }
367   
368    TH1* h[12];
369    h[0] = new TH1F("debug0","",fSteps.size(),tMin,tMax);
370    fDebugList->Add(h[0]);
371    for (Int_t i(1); i < 11; i++) {
372         TString title = Form("debug%d",i);
373         TProfile* hh = NULL;
374         if (( hh = (TProfile*)gDirectory->FindObject(title))) hh->Delete();
375         if (i >= 9)  h[i] = new TProfile(title,"",150,0,30);
376         else         h[i] = new TProfile(title,"",150,tMin,tMax);
377         fDebugList->Add(h[i]);
378    }
379    h[11] = new TH1F("debug11","",150,0,30);
380    fDebugList->Add(h[11]);
381
382    TPolyMarker3D *track = new TPolyMarker3D;
383    fDebugList->Add(track);
384
385    EEvent  *SimuEvent  = fEnergyModule->GetRecoEvent()->GetSimuEvent();
386    Double_t DetectorRadius   = SimuEvent->GetGeometry()->GetOpticsRadius();;
387    Double_t Area = TMath::Pi()*DetectorRadius*DetectorRadius;
388    Double_t EnergyScale=1;
389    if (SimuEvent)
390        EnergyScale = SimuEvent->GetTruth()->GetTrueEnergy()/(fEnergy*1.e-6);
391
392    Int_t points(0);
393    fXmin=fYmin=fZmin=1.e10; 
394    fXmax=fYmax=fZmax=-1.e10; 
395
396    for (UInt_t i(0); i < fSteps.size(); i++) {
397        RecoShowerStep& rs = fSteps[i];
398        h[0]->Fill(rs.fEusoTime-tMin,rs.fNe*EnergyScale);
399        h[1]->Fill(rs.fEusoTime-tMin,rs.fPos.X()/km);
400        h[2]->Fill(rs.fEusoTime-tMin,rs.fPos.Y()/km);
401        h[3]->Fill(rs.fEusoTime-tMin,rs.fPos.Z()/km);
402        h[4]->Fill(rs.fEusoTime-tMin,rs.fFluoYield*m);
403        h[5]->Fill(rs.fEusoTime-tMin,rs.fAttenuation);
404        h[6]->Fill(rs.fEusoTime-tMin,rs.fFlux);
405        h[7]->Fill(rs.fEusoTime-tMin,rs.fLength/km);
406        h[8]->Fill(rs.fEusoTime-tMin,rs.fX);
407
408        h[9]->Fill(rs.fPos.Zv()/km,rs.fFluoYield*m);
409        h[10]->Fill(rs.fPos.Zv()/km,rs.fAttenuation);
410        Double_t fluoEP = rs.fFlux*Area*EnergyScale*Cos(rs.fDir.Theta());
411        Int_t    nfluoEP = (Int_t)fluoEP;
412        h[11]->Fill(rs.fPos.Zv()/km,fluoEP);
413        if (nfluoEP) for (Int_t j(0); j<nfluoEP; j++) track->SetPoint(points++,rs.fPos.X()/km,rs.fPos.Y()/km,rs.fPos.Z()/km);
414        if (fXmin > rs.fPos.X()/km) fXmin = rs.fPos.X()/km;
415        if (fYmin > rs.fPos.Y()/km) fYmin = rs.fPos.Y()/km;
416        if (fZmin > rs.fPos.Z()/km) fZmin = rs.fPos.Z()/km;
417        if (fXmax < rs.fPos.X()/km) fXmax = rs.fPos.X()/km;
418        if (fYmax < rs.fPos.Y()/km) fYmax = rs.fPos.Y()/km;
419        if (fZmax < rs.fPos.Z()/km) fZmax = rs.fPos.Z()/km;
420
421    }
422
423}
424
425//_____________________________________________________________________________
426void RecoShowerTrack::Debug(TCanvas *c, Int_t id, Option_t *opt) {
427    //
428    //
429    //
430    if (!c) {
431        Warning("RecoShowerTrack::Debug","Please specify the pad where I have to draw");
432        return;
433    }
434   
435    if (id < 0 || id > 11) {
436        Warning("RecoShowerTrack::Debug","Please specify id between 0 and 11");
437        return;
438    }
439   
440    TH1* h = (TH1*)fDebugList->At(id);
441    h->Draw(opt);
442}
443
Note: See TracBrowser for help on using the repository browser.