1 | // ESAF : Euso Simulation and Analysis Framework |
---|
2 | // $Id: HmaxByShapeMethod.cc 2946 2011-06-25 13:37:26Z mabl $ |
---|
3 | // Dmitry V. Naumov created Jun, 6 2004 |
---|
4 | |
---|
5 | #include "HmaxByShapeMethod.hh" |
---|
6 | #include "Atmosphere.hh" |
---|
7 | #include "utils.hh" |
---|
8 | #include "RecoEvent.hh" |
---|
9 | #include "EDetector.hh" |
---|
10 | #include "EConst.hh" |
---|
11 | #include "EnergyTypes.hh" |
---|
12 | #include "EnergyViewer.hh" |
---|
13 | |
---|
14 | #include <TH1F.h> |
---|
15 | #include <TF1.h> |
---|
16 | // #include <TMinuit.h> |
---|
17 | #include <TVirtualFitter.h> |
---|
18 | |
---|
19 | #include <iostream> |
---|
20 | |
---|
21 | using namespace sou; |
---|
22 | using namespace TMath; |
---|
23 | using namespace EConst; |
---|
24 | using namespace std; |
---|
25 | |
---|
26 | ClassImp(HmaxByShapeMethod) |
---|
27 | // fFluorescence function |
---|
28 | //_____________________________________________________________________________ |
---|
29 | Double_t HSM_fFluorescence(Double_t *x, Double_t *par) { |
---|
30 | return par[0]*exp(-0.5*pow(x[0]-par[1],2)/pow(par[2],2)); |
---|
31 | } |
---|
32 | |
---|
33 | |
---|
34 | // Sum of background and peak functions |
---|
35 | //_____________________________________________________________________________ |
---|
36 | Double_t HSM_fitFunction(Double_t *x, Double_t *par) { |
---|
37 | return HSM_fFluorescence(x,&par[0]); |
---|
38 | } |
---|
39 | |
---|
40 | //_____________________________________________________________________________ |
---|
41 | float HSM_DensityToAltitude(float h) { |
---|
42 | // HmaxByShapeMethod* hm = (HmaxByShapeMethod*)gMinuit->GetObjectFit(); |
---|
43 | HmaxByShapeMethod* hm = dynamic_cast<HmaxByShapeMethod*>( (TVirtualFitter::GetFitter())->GetUserFunc() ); |
---|
44 | return Atmosphere::Get()->Air_Density(h) - hm->GetDensity(); |
---|
45 | } |
---|
46 | |
---|
47 | //_____________________________________________________________________________ |
---|
48 | HmaxByShapeMethod::HmaxByShapeMethod(EnergyModule* em): EsafMsgSource() { |
---|
49 | // |
---|
50 | // ctor |
---|
51 | // |
---|
52 | if (!em) Msg(EsafMsg::Panic)<< "Pointer to EnergyModule is NULL. I do not want to work like this" <<MsgDispatch; |
---|
53 | |
---|
54 | fEnergyModule = em; |
---|
55 | fEnergy = 1.e14; |
---|
56 | fFitFcn = NULL; |
---|
57 | fSigma = -HUGE; |
---|
58 | fAltitude = -HUGE; |
---|
59 | fDensity = -HUGE; |
---|
60 | } |
---|
61 | |
---|
62 | //_____________________________________________________________________________ |
---|
63 | HmaxByShapeMethod::HmaxByShapeMethod(): EsafMsgSource() { |
---|
64 | // |
---|
65 | // ctor |
---|
66 | // |
---|
67 | fEnergyModule = NULL; |
---|
68 | fEnergy = 1.e14; |
---|
69 | fFitFcn = NULL; |
---|
70 | fSigma = -HUGE; |
---|
71 | fAltitude = -HUGE; |
---|
72 | fDensity = -HUGE; |
---|
73 | } |
---|
74 | |
---|
75 | //_____________________________________________________________________________ |
---|
76 | HmaxByShapeMethod::~HmaxByShapeMethod() { |
---|
77 | // |
---|
78 | // dtor |
---|
79 | // |
---|
80 | Clear(); |
---|
81 | } |
---|
82 | |
---|
83 | //_____________________________________________________________________________ |
---|
84 | void HmaxByShapeMethod::Clear() { |
---|
85 | // |
---|
86 | // clean |
---|
87 | // |
---|
88 | SafeDelete(fFitFcn); |
---|
89 | } |
---|
90 | |
---|
91 | //_____________________________________________________________________________ |
---|
92 | Bool_t HmaxByShapeMethod::Fit(TH1F* data) { |
---|
93 | // |
---|
94 | // fit entry histogramm and estimate altitude of the shower maximum. The entry histo is assumed to be unfolded of inefficients |
---|
95 | // |
---|
96 | Clear(); |
---|
97 | RecoEvent *ev = fEnergyModule->GetRecoEvent(); |
---|
98 | Int_t fFirstHit = ev->GetStartGtu(); |
---|
99 | Int_t fLastHit = ev->GetEndGtu(); |
---|
100 | Float_t gtulength = fEnergyModule->GetRecoEvent()->GetHeader().GetGtuLength(); |
---|
101 | FitHisto(data,fFirstHit,fLastHit,gtulength); |
---|
102 | Int_t debug = fEnergyModule->GetDebugLevel(); |
---|
103 | #ifdef USE_VISUALIZATION |
---|
104 | if (debug>=1) { |
---|
105 | EnergyViewer *fViewer = fEnergyModule->GetEnergyViewer(); |
---|
106 | fViewer->DisplayAltitude(data,fFitFcn); |
---|
107 | } |
---|
108 | #endif |
---|
109 | return kTRUE; |
---|
110 | } |
---|
111 | |
---|
112 | //_____________________________________________________________________________ |
---|
113 | void HmaxByShapeMethod::FitHisto(TH1F* data, Int_t first, Int_t last, Float_t gtulength) { |
---|
114 | // |
---|
115 | // fit entry histogramm and estimate altitude of the shower maximum. The entry histo is assumed to be unfolded of inefficients |
---|
116 | // |
---|
117 | //gMinuit->SetObjectFit(this); |
---|
118 | |
---|
119 | Double_t par[3]={1,1,1}; |
---|
120 | fFitFcn = new TF1("fFitFcn",HSM_fitFunction,first,last,3); |
---|
121 | fFitFcn->SetNpx(500); |
---|
122 | fFitFcn->SetLineWidth(2); |
---|
123 | fFitFcn->SetLineColor(kMagenta); |
---|
124 | |
---|
125 | fFitFcn->SetParLimits(0,0,data->GetMaximum()); |
---|
126 | fFitFcn->SetParLimits(1,0,last); |
---|
127 | fFitFcn->SetParLimits(2,0,0.5*last); |
---|
128 | fFitFcn->SetParameter(1,fFluorMaxTime); |
---|
129 | |
---|
130 | data->Fit("fFitFcn","RQ"); |
---|
131 | fChi2 = data->GetFunction("fFitFcn")->GetChisquare()/data->GetFunction("fFitFcn")->GetNDF(); |
---|
132 | data->GetFunction("fFitFcn")->GetParameters(&par[0]); |
---|
133 | fFitFcn->SetParameters(par); |
---|
134 | fFluorMaxTime = par[1]; |
---|
135 | TVirtualFitter *min = TVirtualFitter::Fitter(0,4); |
---|
136 | min->SetUserFunc(this); |
---|
137 | |
---|
138 | // Compute altitude of the shower maximum |
---|
139 | Float_t x0 = 37.15*g/cm2; |
---|
140 | Float_t t0 = 1.7+0.76*Log(fEnergy/81.); // energy in MeV !!! |
---|
141 | fSigma = par[2]*gtulength; |
---|
142 | fDensity = Sqrt(2*t0)*x0*(1.+fShowerDir*fShowerMaxDir)/Clight()/fSigma; |
---|
143 | DensityToAltitude(); |
---|
144 | Msg(EsafMsg::Debug) << "Reconstructed altitude of the shower maximum "<< fAltitude/km <<MsgDispatch; |
---|
145 | Msg(EsafMsg::Debug) << "Reconstructed density at the shower maximum " << fDensity/g*cm3 << " sigma " << fSigma/microsecond << MsgDispatch; |
---|
146 | } |
---|
147 | |
---|
148 | //_____________________________________________________________________________ |
---|
149 | void HmaxByShapeMethod::DensityToAltitude() { |
---|
150 | // |
---|
151 | // Finds the altitude which corresponds to the given density |
---|
152 | Float_t h1 = 0*km, h2 = 100*km; |
---|
153 | Float_t (*f)(float)= &HSM_DensityToAltitude; |
---|
154 | zbrac(f,&h1,&h2); |
---|
155 | fAltitude = rtbis(f,h1,h2,0.01*km); |
---|
156 | } |
---|
157 | |
---|
158 | |
---|
159 | //_____________________________________________________________________________ |
---|
160 | void HmaxByShapeMethod::Set(Float_t energy,Float_t timemax, TVector3 showerdir, TVector3 showermaxdir){ |
---|
161 | // |
---|
162 | // Set up the module |
---|
163 | // |
---|
164 | SetEnergy(energy); |
---|
165 | SetFluoMaxTime(timemax); |
---|
166 | SetShowerDir(showerdir); |
---|
167 | SetShowerMaxDir(showermaxdir); |
---|
168 | } |
---|