source: trunk/MagicSoft/Mars/mhist/MHFlux.cc@ 1618

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