1 | /* ======================================================================== *\
|
---|
2 | !
|
---|
3 | ! *
|
---|
4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
8 | ! *
|
---|
9 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
10 | ! * documentation for any purpose is hereby granted without fee,
|
---|
11 | ! * provided that the above copyright notice appear in all copies and
|
---|
12 | ! * that both that copyright notice and this permission notice appear
|
---|
13 | ! * in supporting documentation. It is provided "as is" without express
|
---|
14 | ! * or implied warranty.
|
---|
15 | ! *
|
---|
16 | !
|
---|
17 | !
|
---|
18 | ! Author(s): Wolfgang Wittek 5/2002 <mailto:wittek@mppmu.mpg.de>
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2002
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | //////////////////////////////////////////////////////////////////////////////
|
---|
26 | // //
|
---|
27 | // MHFlux //
|
---|
28 | // //
|
---|
29 | // calculates absolute photon fluxes //
|
---|
30 | // from the distributions of the estimated energy //
|
---|
31 | // for the different bins in some variable 'Var' //
|
---|
32 | // (Var = Theta or time) //
|
---|
33 | // //
|
---|
34 | //////////////////////////////////////////////////////////////////////////////
|
---|
35 |
|
---|
36 | #include "MHFlux.h"
|
---|
37 |
|
---|
38 | #include <TStyle.h>
|
---|
39 |
|
---|
40 | #include <TF1.h>
|
---|
41 | #include <TH2.h>
|
---|
42 | #include <TProfile.h>
|
---|
43 |
|
---|
44 |
|
---|
45 | #include <TCanvas.h>
|
---|
46 |
|
---|
47 | #include "MTime.h"
|
---|
48 |
|
---|
49 | #include "MBinning.h"
|
---|
50 | #include "MParList.h"
|
---|
51 |
|
---|
52 | #include "MLog.h"
|
---|
53 | #include "MLogManip.h"
|
---|
54 |
|
---|
55 | ClassImp(MHFlux);
|
---|
56 |
|
---|
57 |
|
---|
58 | // --------------------------------------------------------------------------
|
---|
59 | //
|
---|
60 | // Default Constructor. It sets the variable name (Theta or time)
|
---|
61 | // and the units for the variable
|
---|
62 | //
|
---|
63 | MHFlux::MHFlux(const TH2D &h2d, const Bool_t Draw,
|
---|
64 | const char *varname, const char *unit)
|
---|
65 | : fHOrig(), fHUnfold(), fHFlux()
|
---|
66 | {
|
---|
67 | if (varname == NULL || unit == NULL)
|
---|
68 | {
|
---|
69 | *fLog << "MHFlux : varname or unit not defined" << endl;
|
---|
70 | }
|
---|
71 |
|
---|
72 | fVarname = varname;
|
---|
73 | fUnit = unit;
|
---|
74 |
|
---|
75 | TString strg(varname);
|
---|
76 | strg += unit;
|
---|
77 |
|
---|
78 | // char txt[100];
|
---|
79 |
|
---|
80 | // original distribution of E-est for different bins
|
---|
81 | // of the variable (Theta or time)
|
---|
82 | // sprintf(txt, "gammas vs. E-est and %s",varname);
|
---|
83 |
|
---|
84 | TString strg1 = "no.of gammas vs. E-est and ";
|
---|
85 | strg1 += varname;
|
---|
86 |
|
---|
87 | ((TH2D&)h2d).Copy(fHOrig);
|
---|
88 |
|
---|
89 | fHOrig.SetName("E-est");
|
---|
90 | fHOrig.SetTitle(strg1);
|
---|
91 |
|
---|
92 | fHOrig.SetDirectory(NULL);
|
---|
93 | fHOrig.SetXTitle("E-est [GeV] ");
|
---|
94 | fHOrig.SetYTitle(strg);
|
---|
95 | fHOrig.Sumw2();
|
---|
96 |
|
---|
97 | SetBinning((TH2*)&fHOrig, (TH2*)&h2d);
|
---|
98 |
|
---|
99 |
|
---|
100 | // unfolded distribution of E-unfold for different bins
|
---|
101 | // of the variable (Theta or time)
|
---|
102 | // sprintf(txt, "gammas vs. E-unfold and %s",varname);
|
---|
103 | TString strg2 = "no.of gammas vs. E-unfold and ";
|
---|
104 | strg2 += varname;
|
---|
105 |
|
---|
106 | fHUnfold.SetName("E-unfolded");
|
---|
107 | fHUnfold.SetTitle(strg2);
|
---|
108 |
|
---|
109 | fHUnfold.SetDirectory(NULL);
|
---|
110 | fHUnfold.SetXTitle("E-unfold [GeV] ");
|
---|
111 | fHUnfold.SetYTitle(strg);
|
---|
112 | fHUnfold.Sumw2();
|
---|
113 |
|
---|
114 | SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig);
|
---|
115 |
|
---|
116 |
|
---|
117 | // absolute photon flux vs. E-unfold
|
---|
118 | // for different bins of the variable (Theta or time)
|
---|
119 | //
|
---|
120 | // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname);
|
---|
121 | TString strg3 = "gamma flux [1/(s m2 GeV) vs. E-unfold and ";
|
---|
122 | strg3 += varname;
|
---|
123 |
|
---|
124 | fHFlux.SetName("photon flux");
|
---|
125 | fHFlux.SetTitle(strg3);
|
---|
126 |
|
---|
127 | fHFlux.SetDirectory(NULL);
|
---|
128 | fHFlux.SetXTitle("E-unfold [GeV] ");
|
---|
129 | fHFlux.SetYTitle(strg);
|
---|
130 | fHFlux.Sumw2();
|
---|
131 |
|
---|
132 | SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold);
|
---|
133 |
|
---|
134 |
|
---|
135 | // copy fHOrig into fHUnfold in case no unfolding is done
|
---|
136 | const Int_t nEunf = fHUnfold.GetNbinsX();
|
---|
137 | const Int_t nVar = fHUnfold.GetNbinsY();
|
---|
138 | for (int m=1; m<=nEunf; m++)
|
---|
139 | {
|
---|
140 | for (int n=1; n<=nVar; n++)
|
---|
141 | {
|
---|
142 | Double_t cont = fHOrig.GetBinContent(m,n);
|
---|
143 | Double_t dcont = fHOrig.GetBinError(m,n);
|
---|
144 | fHUnfold.SetBinContent(m,n,cont);
|
---|
145 | fHUnfold.SetBinError(m,n,dcont);
|
---|
146 | }
|
---|
147 | }
|
---|
148 | //..............................................
|
---|
149 | // draw the No.of photons vs. E-est
|
---|
150 | // for the individual bins of the variable Var
|
---|
151 |
|
---|
152 | if (Draw == kTRUE)
|
---|
153 | {
|
---|
154 | const Int_t nVar = fHOrig.GetNbinsY();
|
---|
155 |
|
---|
156 | for (int n=1; n<=nVar; n++)
|
---|
157 | {
|
---|
158 | TString strg0("Orig-");
|
---|
159 | strg0 += fVarname;
|
---|
160 | TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E");
|
---|
161 |
|
---|
162 | char txt0[100];
|
---|
163 | strg0 = fVarname;
|
---|
164 | strg0 += "-bin %d";
|
---|
165 | sprintf(txt0, strg0, n);
|
---|
166 |
|
---|
167 | TString strg1("No.of photons vs. E-est for ");
|
---|
168 | strg1 += txt0;
|
---|
169 | new TCanvas(txt0,strg1);
|
---|
170 | // TCanvas &c = *MakeDefCanvas(txt0, strg1);
|
---|
171 | // gROOT->SetSelectedPad(NULL);
|
---|
172 |
|
---|
173 | gPad->SetLogx();
|
---|
174 |
|
---|
175 | h.SetName(txt0);
|
---|
176 | h.SetTitle(strg1);
|
---|
177 | h.SetXTitle("E-est [GeV] ");
|
---|
178 | h.SetYTitle("No.of photons");
|
---|
179 | h.DrawCopy();
|
---|
180 |
|
---|
181 | // c.Modified();
|
---|
182 | // c.Update();
|
---|
183 | }
|
---|
184 | }
|
---|
185 | //........................
|
---|
186 | }
|
---|
187 |
|
---|
188 | // -------------------------------------------------------------------------
|
---|
189 | //
|
---|
190 | // Dummy Fill (has to be included because in base class MH Fill is set to 0);
|
---|
191 | // without the dummy Fill one gets the error message :
|
---|
192 | //
|
---|
193 | // Error: Can't call MHFlux::MHFlux(evttime,"time","[s]") in current scope
|
---|
194 | // FILE:macros/flux.C LINE:465
|
---|
195 | // Possible candidates are...
|
---|
196 | // filename line:size busy function type and name (in MHFlux)
|
---|
197 | // filename line:size busy function type and name (in MH)
|
---|
198 | // filename line:size busy function type and name (in MParContainer)
|
---|
199 | // filename line:size busy function type and name (in TObject)
|
---|
200 | //
|
---|
201 | Bool_t MHFlux::Fill(const MParContainer *par)
|
---|
202 | {
|
---|
203 | return kTRUE;
|
---|
204 | }
|
---|
205 |
|
---|
206 |
|
---|
207 | // -------------------------------------------------------------------------
|
---|
208 | //
|
---|
209 | // Unfold the distribution in E-est
|
---|
210 | //
|
---|
211 | void MHFlux::Unfold(const Bool_t Draw)
|
---|
212 | {
|
---|
213 | //..............................................
|
---|
214 | // draw the No.of photons vs. E-unfold
|
---|
215 | // for the individual bins of the variable Var
|
---|
216 |
|
---|
217 | if (Draw == kTRUE)
|
---|
218 | {
|
---|
219 | const Int_t nVar = fHUnfold.GetNbinsY();
|
---|
220 |
|
---|
221 | for (int n=1; n<=nVar; n++)
|
---|
222 | {
|
---|
223 | TString strg0("Unfold-");
|
---|
224 | strg0 += fVarname;
|
---|
225 | TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
|
---|
226 |
|
---|
227 | char txt0[100];
|
---|
228 | strg0 = fVarname;
|
---|
229 | strg0 += "-bin %d";
|
---|
230 | sprintf(txt0, strg0, n);
|
---|
231 |
|
---|
232 | TString strg1("No.of photons vs. E-unfold for ");
|
---|
233 | strg1 += txt0;
|
---|
234 | new TCanvas(txt0,strg1);
|
---|
235 |
|
---|
236 | // TCanvas &c = *MakeDefCanvas(txt0, strg1);
|
---|
237 | // gROOT->SetSelectedPad(NULL);
|
---|
238 |
|
---|
239 | gPad->SetLogx();
|
---|
240 |
|
---|
241 | h.SetName(txt0);
|
---|
242 | h.SetTitle(strg1);
|
---|
243 | h.SetXTitle("E-unfold [GeV] ");
|
---|
244 | h.SetYTitle("No.of photons");
|
---|
245 | h.DrawCopy();
|
---|
246 |
|
---|
247 | // c.Modified();
|
---|
248 | // c.Update();
|
---|
249 | }
|
---|
250 | }
|
---|
251 | //........................
|
---|
252 | }
|
---|
253 |
|
---|
254 |
|
---|
255 | // -------------------------------------------------------------------------
|
---|
256 | //
|
---|
257 | // Calculate photon flux by dividing the distribution in Eunf (fHUnfold) by
|
---|
258 | // the width of the energy interval (deltaE)
|
---|
259 | // the effective ontime (*teff)
|
---|
260 | // and the effective collection area (*aeff)
|
---|
261 | //
|
---|
262 | void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar,
|
---|
263 | const TH2D *aeff, const Bool_t Draw)
|
---|
264 | {
|
---|
265 | // Note that fHUnfold has bins in Eunf and Var
|
---|
266 | // *teff has bins in Var (the same bins in Var as fHUnfold)
|
---|
267 | // *thetabar has bins in Var (the same bins in Var as fHUnfold)
|
---|
268 | // *aeff has bins in Etru and Theta
|
---|
269 | // (where in general the binning in Etru is different
|
---|
270 | // from the binning in Eunf)
|
---|
271 | // The variable Var may be 'time' or 'Theta'
|
---|
272 |
|
---|
273 | // Draw = kTRUE means the differential flux vs E-unf should be drawn
|
---|
274 | // for the individual bins of the variable Var
|
---|
275 |
|
---|
276 | //....................................
|
---|
277 | // define dummy histogram *aeff
|
---|
278 | ((TH1*)aeff)->Sumw2();
|
---|
279 | MBinning binsetru("BinningEtru");
|
---|
280 | binsetru.SetEdgesLog(10, 10, 1e3);
|
---|
281 |
|
---|
282 | MBinning binsthetatru("BinningThetatru");
|
---|
283 | binsthetatru.SetEdges(7, -2.5, 32.5);
|
---|
284 | //SetBinning((TH1*)aeff, &binsetru, &binsthetatru);
|
---|
285 | SetBinning((TH2*)aeff, &binsetru, &binsthetatru);
|
---|
286 |
|
---|
287 | const Int_t netru = aeff->GetNbinsX();
|
---|
288 | const Int_t ntheta = aeff->GetNbinsY();
|
---|
289 |
|
---|
290 | for (int j=1; j<=netru; j++)
|
---|
291 | {
|
---|
292 | for (int k=1; k<=ntheta; k++)
|
---|
293 | {
|
---|
294 | Double_t cont = 10000.0;;
|
---|
295 | ((TH1*)aeff)->SetBinContent(j,k,cont);
|
---|
296 |
|
---|
297 | Double_t dcont = 100.0;
|
---|
298 | ((TH1*)aeff)->SetBinError(j,k,dcont);
|
---|
299 | }
|
---|
300 | }
|
---|
301 | // *fLog << "Dummy aeff : netru =" << netru << ", ntheta = " << ntheta << endl;
|
---|
302 | //....................................
|
---|
303 |
|
---|
304 | // number of Eunf and Var bins (histograms : fHUnfold, fHFlux)
|
---|
305 | const Int_t nEunf = fHFlux.GetNbinsX();
|
---|
306 | const Int_t nVar = fHFlux.GetNbinsY();
|
---|
307 |
|
---|
308 | // number of Etru and Theta bins (histogram *aeff of collection area)
|
---|
309 |
|
---|
310 | const Int_t nEtru = aeff->GetNbinsX();
|
---|
311 | const Int_t nTheta = aeff->GetNbinsY();
|
---|
312 |
|
---|
313 | //*fLog << "nEunf =" << nEunf << ", nVar = " << nVar << endl;
|
---|
314 | //*fLog << "nEtru =" << nEtru << ", nTheta = " << nTheta << endl;
|
---|
315 |
|
---|
316 | //...................................................................
|
---|
317 | // calculate effective collection area
|
---|
318 | // for the Eunf and Var bins of the histogram fHUnfold
|
---|
319 | // from the histogram *aeff, which has bins in Etru and Theta
|
---|
320 | // the result is the histogram fHAeff
|
---|
321 | //
|
---|
322 |
|
---|
323 | TH2D fHAeff;
|
---|
324 | fHAeff.Sumw2();
|
---|
325 | SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
|
---|
326 |
|
---|
327 | Int_t errflag;
|
---|
328 | Double_t c0, c1, c2;
|
---|
329 | Double_t t1, t2, t3;
|
---|
330 | Double_t a1, a2, a3;
|
---|
331 |
|
---|
332 | Double_t *aeffbar;
|
---|
333 | aeffbar = new Double_t[nEtru];
|
---|
334 | Double_t *daeffbar;
|
---|
335 | daeffbar = new Double_t[nEtru];
|
---|
336 |
|
---|
337 | Double_t aeffEunfVar;
|
---|
338 | Double_t daeffEunfVar;
|
---|
339 |
|
---|
340 |
|
---|
341 | //------ start n loop ------
|
---|
342 | for (int n=1; n<=nVar; n++)
|
---|
343 | {
|
---|
344 | Double_t Thetabar = thetabar->GetBinContent(n);
|
---|
345 | Double_t cosThetabar = cos(Thetabar);
|
---|
346 |
|
---|
347 | // determine Theta bins (k1, k2, k3) for interpolation in Theta
|
---|
348 | // k0 denotes the Theta bin from whicvh the error is copied
|
---|
349 | Int_t k0=0, k1=0, k2=0, k3=0;
|
---|
350 | for (int k=3; k<=nTheta; k++)
|
---|
351 | {
|
---|
352 | Double_t Thetalow = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(k);
|
---|
353 | if (Thetabar < Thetalow)
|
---|
354 | {
|
---|
355 | k1 = k-2;
|
---|
356 | k2 = k-1;
|
---|
357 | k3 = k;
|
---|
358 | k0 = k2;
|
---|
359 | break;
|
---|
360 | }
|
---|
361 | }
|
---|
362 |
|
---|
363 | if (k3 == 0)
|
---|
364 | {
|
---|
365 | k1 = nTheta-2;
|
---|
366 | k2 = nTheta-1;
|
---|
367 | k3 = nTheta;
|
---|
368 | k0 = k2;
|
---|
369 | }
|
---|
370 |
|
---|
371 | if (Thetabar < ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2))
|
---|
372 | k0 = 1;
|
---|
373 | else if (Thetabar > ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta))
|
---|
374 | k0 = nTheta;
|
---|
375 |
|
---|
376 | Double_t Thetamin = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(1);
|
---|
377 | Double_t Thetamax = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta+1);
|
---|
378 | if (Thetabar < Thetamin || Thetabar > Thetamax)
|
---|
379 | {
|
---|
380 | *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar
|
---|
381 | << ", Thetamin =" << Thetamin
|
---|
382 | << ", Thetamax =" << Thetamax << endl;
|
---|
383 | }
|
---|
384 |
|
---|
385 | //*fLog << "Var bin " << n << ":" << endl;
|
---|
386 | //*fLog << "Thetabar= " << Thetabar
|
---|
387 | // << ", k0= " << k0
|
---|
388 | // << ", k1= " << k1
|
---|
389 | // << ", k2= " << k2
|
---|
390 | // << ", k3= " << k3 << endl;
|
---|
391 |
|
---|
392 |
|
---|
393 | // calculate effective collection area at Theta = Thetabar
|
---|
394 | // by quadratic interpolation in cos(Theta);
|
---|
395 | // do this for each bin of Etru
|
---|
396 | //
|
---|
397 |
|
---|
398 | for (int j=1; j<=nEtru; j++)
|
---|
399 | {
|
---|
400 | c0 = 0.0;
|
---|
401 | c1 = 0.0;
|
---|
402 | c2 = 0.0;
|
---|
403 |
|
---|
404 | t1 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k1) );
|
---|
405 | t2 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k2) );
|
---|
406 | t3 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k3) );
|
---|
407 |
|
---|
408 | a1 = aeff->GetBinContent(j,k1);
|
---|
409 | a2 = aeff->GetBinContent(j,k2);
|
---|
410 | a3 = aeff->GetBinContent(j,k3);
|
---|
411 |
|
---|
412 | Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
|
---|
413 | aeffbar[j] = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
|
---|
414 | daeffbar[j] = aeff->GetBinError(j,k0);
|
---|
415 |
|
---|
416 | //*fLog << "Etru bin " << j << ": tbar= " << Thetabar
|
---|
417 | // << ", abar= " << aeffbar[j]
|
---|
418 | // << ", dabar= " << daeffbar[j] << endl;
|
---|
419 | }
|
---|
420 |
|
---|
421 | //--- start m loop ---
|
---|
422 | // calculate effective collection area at (E = Ebar, Theta = Thetabar)
|
---|
423 | // by quadratic interpolation in log10(Etru)
|
---|
424 | // do this for each bin of Eunf
|
---|
425 | //
|
---|
426 | for (int m=1; m<=nEunf; m++)
|
---|
427 | {
|
---|
428 | Double_t log10Ebar = 0.5 *
|
---|
429 | ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m) )
|
---|
430 | +log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
|
---|
431 | Double_t Ebar = pow(10.0, log10Ebar);
|
---|
432 |
|
---|
433 | // determine Etru bins (j1, j2, j3) for interpolation in E
|
---|
434 | // j0 denotes the Etru bin from which the error is copied
|
---|
435 | Int_t j0=0, j1=0, j2=0, j3=0;
|
---|
436 |
|
---|
437 | for (int j=3; j<=nEtru; j++)
|
---|
438 | {
|
---|
439 | Double_t Elow = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(j);
|
---|
440 | if (Ebar < Elow)
|
---|
441 | {
|
---|
442 | j1 = j-2;
|
---|
443 | j2 = j-1;
|
---|
444 | j3 = j;
|
---|
445 | j0 = j2;
|
---|
446 | break;
|
---|
447 | }
|
---|
448 | }
|
---|
449 |
|
---|
450 | if (j3 == 0)
|
---|
451 | {
|
---|
452 | j1 = nEtru-2;
|
---|
453 | j2 = nEtru-1;
|
---|
454 | j3 = nEtru;
|
---|
455 | j0 = j2;
|
---|
456 | }
|
---|
457 |
|
---|
458 | if (Ebar < ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(2))
|
---|
459 | j0 = 1;
|
---|
460 | else if (Ebar > ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru))
|
---|
461 | j0 = nEtru;
|
---|
462 |
|
---|
463 | Double_t Etrumin = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(1);
|
---|
464 | Double_t Etrumax = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru+1);
|
---|
465 | if (Ebar < Etrumin || Ebar > Etrumax)
|
---|
466 | {
|
---|
467 | *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
|
---|
468 | << ", Etrumin =" << Etrumin
|
---|
469 | << ", Etrumax =" << Etrumax << endl;
|
---|
470 | }
|
---|
471 |
|
---|
472 | //*fLog << "Var bin " << n << ":" << endl;
|
---|
473 | //*fLog << "Ebar= " << Ebar
|
---|
474 | // << ", j1= " << j1
|
---|
475 | // << ", j2= " << j2
|
---|
476 | // << ", j3= " << j3 << endl;
|
---|
477 |
|
---|
478 |
|
---|
479 | c0=0.0;
|
---|
480 | c1=0.0;
|
---|
481 | c2=0.0;
|
---|
482 |
|
---|
483 | t1 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1) )
|
---|
484 | +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1+1)) );
|
---|
485 |
|
---|
486 | t2 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2) )
|
---|
487 | +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2+1)) );
|
---|
488 |
|
---|
489 | t3 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3) )
|
---|
490 | +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3+1)) );
|
---|
491 |
|
---|
492 |
|
---|
493 | a1 = aeffbar[j1];
|
---|
494 | a2 = aeffbar[j2];
|
---|
495 | a3 = aeffbar[j3];
|
---|
496 |
|
---|
497 | Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
|
---|
498 | aeffEunfVar = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
|
---|
499 | daeffEunfVar = daeffbar[j0];
|
---|
500 |
|
---|
501 | //*fLog << "Eunf bin " << m << ": Ebar= " << Ebar
|
---|
502 | // << ", aeffEunfVar = " << aeffEunfVar
|
---|
503 | // << ", daeffEunfVar = " << daeffEunfVar << endl;
|
---|
504 |
|
---|
505 | fHAeff.SetBinContent(m,n,aeffEunfVar);
|
---|
506 | fHAeff.SetBinError(m,n,daeffEunfVar);
|
---|
507 | }
|
---|
508 | //--- end m loop ---
|
---|
509 | }
|
---|
510 | //------ end n loop ------
|
---|
511 | delete aeffbar;
|
---|
512 |
|
---|
513 | //...................................................................
|
---|
514 | // now calculate the absolute gamma flux
|
---|
515 | //
|
---|
516 | for (int m=1; m<=nEunf; m++)
|
---|
517 | {
|
---|
518 | Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
|
---|
519 |
|
---|
520 | for (int n=1; n<=nVar; n++)
|
---|
521 | {
|
---|
522 | Double_t Ngam = fHUnfold.GetBinContent(m,n);
|
---|
523 | Double_t dNgam = fHUnfold.GetBinError(m,n);
|
---|
524 |
|
---|
525 | Double_t Aeff = fHAeff.GetBinContent(m,n);
|
---|
526 | Double_t dAeff = fHAeff.GetBinError(m,n);
|
---|
527 |
|
---|
528 | Double_t Effon = teff->GetBinContent(n);
|
---|
529 | Double_t dEffon = teff->GetBinError(n);
|
---|
530 |
|
---|
531 | Double_t Cont, dCont;
|
---|
532 | if (Ngam > 0.0 && DeltaE > 0.0 && Effon > 0.0 && Aeff > 0.0)
|
---|
533 | {
|
---|
534 | Cont = Ngam / (DeltaE * Effon * Aeff);
|
---|
535 | dCont = Cont * sqrt( dNgam*dNgam / (Ngam*Ngam)
|
---|
536 | + dEffon*dEffon / (Effon*Effon)
|
---|
537 | + dAeff*dAeff / (Aeff*Aeff) );
|
---|
538 | }
|
---|
539 | else
|
---|
540 | {
|
---|
541 | Cont = 1.e-20;
|
---|
542 | dCont = 1.e-20;
|
---|
543 | }
|
---|
544 |
|
---|
545 | fHFlux.SetBinContent(m,n,Cont);
|
---|
546 | fHFlux.SetBinError(m,n,dCont);
|
---|
547 |
|
---|
548 | //*fLog << "Eunf bin " << m << ", Var bin " << n
|
---|
549 | // << ": Ngam = " << Ngam << ", Flux = "
|
---|
550 | // << Cont << ", dFlux = " << dCont << endl;
|
---|
551 | //*fLog << ", DeltaE = " << DeltaE << ", Effon = " << Effon
|
---|
552 | // << ", Aeff = " << Aeff << endl;
|
---|
553 | }
|
---|
554 | }
|
---|
555 |
|
---|
556 | //..............................................
|
---|
557 | // draw the differential photon flux vs. E-unf
|
---|
558 | // for the individual bins of the variable Var
|
---|
559 |
|
---|
560 | if (Draw == kTRUE)
|
---|
561 | {
|
---|
562 | for (int n=1; n<=nVar; n++)
|
---|
563 | {
|
---|
564 | TString strg0("Flux-");
|
---|
565 | strg0 += fVarname;
|
---|
566 | TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
|
---|
567 |
|
---|
568 | char txt[100];
|
---|
569 | TString strg1("Photon flux vs. E-unfold for ");
|
---|
570 | TString strg2 = fVarname;
|
---|
571 | strg2 += "-bin %d";
|
---|
572 | sprintf(txt, strg2, n);
|
---|
573 | TString strg3 = strg1 + txt;
|
---|
574 |
|
---|
575 | new TCanvas(txt, strg3);
|
---|
576 | // TCanvas &c = *MakeDefCanvas(txt, txt);
|
---|
577 | // gROOT->SetSelectedPad(NULL);
|
---|
578 |
|
---|
579 | gPad->SetLogx();
|
---|
580 |
|
---|
581 | h.SetName(txt);
|
---|
582 | h.SetTitle(strg3);
|
---|
583 | h.SetXTitle("E-unfold [GeV] ");
|
---|
584 | h.SetYTitle("photons / (s m2 GeV)");
|
---|
585 | h.DrawCopy();
|
---|
586 |
|
---|
587 | // c.Modified();
|
---|
588 | // c.Update();
|
---|
589 | }
|
---|
590 | }
|
---|
591 | //........................
|
---|
592 | }
|
---|
593 |
|
---|
594 | // -------------------------------------------------------------------------
|
---|
595 | //
|
---|
596 | // Draw copies of the histograms
|
---|
597 | //
|
---|
598 | TObject *MHFlux::DrawClone(Option_t *opt) const
|
---|
599 | {
|
---|
600 | TCanvas &c = *MakeDefCanvas("flux", "Orig - Unfold - Flux plots");
|
---|
601 | c.Divide(2, 2);
|
---|
602 |
|
---|
603 | gROOT->SetSelectedPad(NULL);
|
---|
604 |
|
---|
605 | c.cd(1);
|
---|
606 | ((TH2*)&fHOrig)->DrawCopy("");
|
---|
607 | gPad->SetLogx();
|
---|
608 |
|
---|
609 | c.cd(2);
|
---|
610 | ((TH2*)&fHUnfold)->DrawCopy("");
|
---|
611 | gPad->SetLogx();
|
---|
612 |
|
---|
613 | c.cd(3);
|
---|
614 | ((TH2*)&fHFlux)->DrawCopy("");
|
---|
615 | gPad->SetLogx();
|
---|
616 |
|
---|
617 | c.Modified();
|
---|
618 | c.Update();
|
---|
619 |
|
---|
620 | return &c;
|
---|
621 | }
|
---|
622 |
|
---|
623 | // -------------------------------------------------------------------------
|
---|
624 | //
|
---|
625 | // Draw the histograms
|
---|
626 | //
|
---|
627 | void MHFlux::Draw(Option_t *opt)
|
---|
628 | {
|
---|
629 | if (!gPad)
|
---|
630 | MakeDefCanvas("flux", "orig-unfold-flux plots");
|
---|
631 |
|
---|
632 | gPad->Divide(2,2);
|
---|
633 |
|
---|
634 | gPad->cd(1);
|
---|
635 | fHOrig.Draw(opt);
|
---|
636 |
|
---|
637 | gPad->cd(2);
|
---|
638 | fHUnfold.Draw(opt);
|
---|
639 |
|
---|
640 | gPad->cd(3);
|
---|
641 | fHFlux.Draw(opt);
|
---|
642 |
|
---|
643 | gPad->Modified();
|
---|
644 | gPad->Update();
|
---|
645 | }
|
---|
646 |
|
---|
647 | // -------------------------------------------------------------------------
|
---|
648 | //
|
---|
649 | // Quadratic interpolation
|
---|
650 | //
|
---|
651 | // *** calculate the parameters of a parabula
|
---|
652 | // y = a + b*x + c*x**2 = F(x)
|
---|
653 | // such that yi = F(xi) for (i=1,3)
|
---|
654 | //
|
---|
655 | void MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
|
---|
656 | Double_t y1, Double_t y2, Double_t y3,
|
---|
657 | Double_t *a, Double_t *b, Double_t *c, Int_t *errflag)
|
---|
658 | {
|
---|
659 | double ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33;
|
---|
660 | double det,det1;
|
---|
661 | //double yt1,yt2,yt3;
|
---|
662 |
|
---|
663 | det = x2*x3*x3 + x1*x2*x2 + x3*x1*x1
|
---|
664 | - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
|
---|
665 |
|
---|
666 | if (det != 0.0)
|
---|
667 | {
|
---|
668 | *errflag = 0;
|
---|
669 | det1 = 1.0/det;
|
---|
670 |
|
---|
671 | ai11 = x2*x3*x3 - x3*x2*x2;
|
---|
672 | ai12 = x3*x1*x1 - x1*x3*x3;
|
---|
673 | ai13 = x1*x2*x2 - x2*x1*x1;
|
---|
674 |
|
---|
675 | ai21 = x2*x2 - x3*x3;
|
---|
676 | ai22 = x3*x3 - x1*x1;
|
---|
677 | ai23 = x1*x1 - x2*x2;
|
---|
678 |
|
---|
679 | ai31 = x3 - x2;
|
---|
680 | ai32 = x1 - x3;
|
---|
681 | ai33 = x2 - x1;
|
---|
682 |
|
---|
683 | *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1;
|
---|
684 | *b = (ai21*y1 + ai22*y2 + ai23*y3) * det1;
|
---|
685 | *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1;
|
---|
686 |
|
---|
687 | //yt1 = *a + *b * x1 + *c * x1*x1;
|
---|
688 | //yt2 = *a + *b * x2 + *c * x2*x2;
|
---|
689 | //yt3 = *a + *b * x3 + *c * x3*x3;
|
---|
690 |
|
---|
691 | //*fLog << "x1 = " << x1 << ", x2 = " << x2 << ", x3 = " << x3 << endl;
|
---|
692 | //*fLog << "y1 = " << y1 << ", y2 = " << y2 << ", y3 = " << y3 << endl;
|
---|
693 | //*fLog << "yt1 = " << yt1 << ", yt2 = " << yt2
|
---|
694 | // << ", yt3 = " << yt3 << endl;
|
---|
695 | //*fLog << "*a = " << *a << ", *b = " << *b << ", *c= " << *c
|
---|
696 | // << ", *errflag = " << *errflag << endl;
|
---|
697 |
|
---|
698 | return;
|
---|
699 | }
|
---|
700 |
|
---|
701 | *errflag = 1;
|
---|
702 | *a = 0.0;
|
---|
703 | *b = 0.0;
|
---|
704 | *c = 0.0;
|
---|
705 | return;
|
---|
706 | }
|
---|