source: releases/Mars.2014.05.26/mhflux/MHEnergyEst.cc@ 18558

Last change on this file since 18558 was 10281, checked in by tbretz, 14 years ago
File size: 16.9 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 4/2002 <mailto:wittek@mppmu.mpg.de>
19! Author(s): Abelardo Moralejo 5/2003 <mailto:moralejo@pd.infn.it>
20! Author(s): Thomas Bretz 1/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2008
23!
24!
25\* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28//
29// MHEnergyEst
30//
31// calculates the migration matrix E-est vs. E-true
32// for different bins in Theta
33//
34// Class Version 2:
35// - fHResolution
36// + fHResolutionEst
37// + fHResolutionMC
38//
39//////////////////////////////////////////////////////////////////////////////
40#include "MHEnergyEst.h"
41
42#include <TF1.h>
43#include <TLine.h>
44#include <TMath.h>
45#include <TCanvas.h>
46#include <TStyle.h>
47#include <TProfile.h>
48#include <TPaveStats.h>
49
50#include "MString.h"
51
52#include "MMcEvt.hxx"
53
54#include "MBinning.h"
55#include "MParList.h"
56#include "MParameters.h"
57#include "MHMatrix.h"
58
59#include "MLog.h"
60#include "MLogManip.h"
61
62ClassImp(MHEnergyEst);
63
64using namespace std;
65
66// --------------------------------------------------------------------------
67//
68// Default Constructor. It sets name and title of the histogram.
69//
70MHEnergyEst::MHEnergyEst(const char *name, const char *title)
71 : fMcEvt(0), fEnergy(0), fResult(0), fMatrix(0)
72{
73 //
74 // set the name and title of this object
75 //
76 fName = name ? name : "MHEnergyEst";
77 fTitle = title ? title : "Histogram for the result of the energy reconstruction";
78
79 //fNameEnergy = "MEnergyEst";
80 //fNameResult = "MinimizationValue";
81
82 fHEnergy.SetDirectory(NULL);
83 fHEnergy.SetName("EnergyEst");
84 fHEnergy.SetTitle("Histogram in E_{est}, E_{mc} and \\Theta");
85 fHEnergy.SetXTitle("E_{est} [GeV]");
86 fHEnergy.SetYTitle("E_{mc} [GeV]");
87 fHEnergy.SetZTitle("\\Theta [\\circ]");
88
89 fHResolutionEst.SetDirectory(NULL);
90 fHResolutionEst.SetName("ResEnergyEst");
91 fHResolutionEst.SetTitle("Histogram in \\Delta E/E_{est} vs E_{est}");
92 fHResolutionEst.SetXTitle("E_{est} [GeV]");
93 fHResolutionEst.SetYTitle("1-E_{mc}/E_{est}");
94
95 fHResolutionMC.SetDirectory(NULL);
96 fHResolutionMC.SetName("ResEnergyMC");
97 fHResolutionMC.SetTitle("Histogram in \\Delta E/E_{mc} vs E_{mc}");
98 fHResolutionMC.SetXTitle("E_{mc} [GeV]");
99 fHResolutionMC.SetYTitle("E_{est}/E_{mc}-1");
100
101 fHImpact.SetDirectory(NULL);
102 fHImpact.SetName("Impact");
103 fHImpact.SetTitle("\\Delta E/E vs Impact parameter");
104 fHImpact.SetXTitle("Impact parameter [m]");
105 fHImpact.SetYTitle("E_{est}/E_{mc}-1");
106
107 MBinning binsi, binse, binst, binsr;
108 binse.SetEdgesLog(21, 6.3, 100000);
109 binst.SetEdgesASin(51, -0.005, 0.505);
110 binsr.SetEdges(75, -1.75, 1.75);
111
112 // Use the binning in impact to do efficiency studies
113 binsi.SetEdges(1, 0, 1000);
114
115 SetBinning(fHEnergy, binse, binse, binst);
116 SetBinning(fHImpact, binsi, binsr);
117
118 SetBinning(fHResolutionEst, binse, binsr);
119 SetBinning(fHResolutionMC, binse, binsr);
120
121 // For some unknown reasons this must be called after
122 // the binning has been initialized at least once
123 fHEnergy.Sumw2();
124 fHImpact.Sumw2();
125 fHResolutionEst.Sumw2();
126 fHResolutionMC.Sumw2();
127}
128
129// --------------------------------------------------------------------------
130//
131// Set the binnings and prepare the filling of the histograms
132//
133Bool_t MHEnergyEst::SetupFill(const MParList *plist)
134{
135 if (!fMatrix)
136 {
137 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
138 if (!fMcEvt)
139 {
140 *fLog << err << "MMcEvt not found... aborting." << endl;
141 return kFALSE;
142 }
143 }
144
145 fEnergy = (MParameterD*)plist->FindObject("MEnergyEst", "MParameterD");
146 if (!fEnergy)
147 {
148 *fLog << err << "MEnergyEst [MParameterD] not found... aborting." << endl;
149 return kFALSE;
150 }
151
152 fResult = (MParameterD*)const_cast<MParList*>(plist)->FindCreateObj("MParameterD", "MinimizationValue");
153 if (!fResult)
154 return kFALSE;
155
156 MBinning binst, binse, binsi, binsr;
157 binse.SetEdges(fHEnergy, 'x');
158 binst.SetEdges(fHEnergy, 'z');
159 binsi.SetEdges(fHImpact, 'x');
160 binsr.SetEdges(fHImpact, 'y');
161
162 binst.SetEdges(*plist, "BinningTheta");
163 binse.SetEdges(*plist, "BinningEnergyEst");
164 binsi.SetEdges(*plist, "BinningImpact");
165 binsr.SetEdges(*plist, "BinningEnergyRes");
166
167 SetBinning(fHEnergy, binse, binse, binst);
168 SetBinning(fHImpact, binsi, binsr);
169
170 SetBinning(fHResolutionEst, binse, binsr);
171 SetBinning(fHResolutionMC, binse, binsr);
172
173 fChisq = 0;
174 fBias = 0;
175
176 fHEnergy.Reset();
177 fHImpact.Reset();
178
179 fHResolutionEst.Reset();
180 fHResolutionMC.Reset();
181
182 return kTRUE;
183}
184
185// --------------------------------------------------------------------------
186//
187// Fill the histogram
188//
189Int_t MHEnergyEst::Fill(const MParContainer *par, const Stat_t w)
190{
191 const Double_t eest = fEnergy->GetVal();
192 const Double_t etru = fMatrix ? GetVal(0) : fMcEvt->GetEnergy();
193 const Double_t imp = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
194 const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
195
196 const Double_t resEst = (eest-etru)/eest;
197 const Double_t resMC = (eest-etru)/etru;
198
199 fHEnergy.Fill(eest, etru, theta, w);
200 fHImpact.Fill(imp, resEst, w);
201
202 fHResolutionEst.Fill(eest, resEst, w);
203 fHResolutionMC.Fill(etru, resMC, w);
204
205 // For the fit we use a different quantity
206 //const Double_t res = TMath::Log10(eest/etru);
207 const Double_t res = eest-etru;
208
209 fChisq += res*res;
210 fBias += res;
211
212 return kTRUE;
213}
214
215// --------------------------------------------------------------------------
216//
217// Divide chisq and bias by number of executions
218// Print result
219//
220Bool_t MHEnergyEst::Finalize()
221{
222 fChisq /= GetNumExecutions();
223 fBias /= GetNumExecutions();
224
225 fResult->SetVal(fChisq);
226
227 *fLog << all << endl;
228 Print();
229
230 return kTRUE;
231}
232
233// --------------------------------------------------------------------------
234//
235// Print result
236//
237void MHEnergyEst::Print(Option_t *o) const
238{
239 const Double_t res = TMath::Sqrt(fChisq-fBias*fBias);
240 if (!TMath::Finite(res))
241 {
242 *fLog << all << "MHEnergyEst::Print: ERROR - Resolution is not finite (eg. NaN)." << endl;
243 return;
244 }
245
246 TH1D *h = (TH1D*)fHResolutionEst.ProjectionY("Dummy", -1, -1, "s");
247 h->Fit("gaus", "Q0", "", -1.0, 0.25);
248
249 TF1 *f = h->GetFunction("gaus");
250
251 *fLog << all << "F=" << fChisq << endl;
252 *fLog << "Results from Histogram:" << endl;
253 *fLog << " Mean of Delta E/E: " << Form("%+4.2f%%", 100*h->GetMean()) << endl;
254 *fLog << " RMS of Delta E/E: " << Form("%4.2f%%", 100*h->GetRMS()) << endl;
255 *fLog << "Results from Histogram-Fit:" << endl;
256 if (!f)
257 *fLog << " <fit did not succeed>" << endl;
258 else
259 {
260 *fLog << " Bias of Delta E/E: " << Form("%+4.2f%%", 100*f->GetParameter(1)) << endl;
261 *fLog << " Sigma of Delta E/E: " << Form("%4.2f%%", 100*f->GetParameter(2)) << endl;
262 *fLog << " Res of Delta E/E: " << Form("%4.2f%%", 100*TMath::Hypot(f->GetParameter(1), f->GetParameter(2))) << endl;
263 }
264
265 delete h;
266}
267
268// --------------------------------------------------------------------------
269//
270// Return Correction Coefficients (weights)
271// hist = E_mc/E_est
272//
273void MHEnergyEst::GetWeights(TH1D &hist) const
274{
275 // Project into EnergyEst_ey
276 // the "e" ensures that errors are calculated
277 TH1D *h1 = (TH1D*)fHEnergy.Project3D("rtlprmft_ex"); // E_Est
278 TH1D *h2 = (TH1D*)fHEnergy.Project3D("rtlprmft_ey"); // E_mc
279
280 h2->Copy(hist);
281 hist.Divide(h1);
282
283 delete h1;
284 delete h2;
285
286 hist.SetNameTitle("EnergyRatio", "Ratio between estimated and monte carlo energy");
287 hist.SetXTitle("E [GeV]");
288 hist.SetYTitle("N_{mc}/N_{est} [1]");
289 hist.SetDirectory(0);
290}
291
292void MHEnergyEst::Paint(Option_t *opt)
293{
294 TVirtualPad *pad = gPad;
295
296 pad->cd(1);
297
298 TH1 *hx=0;
299 TH1 *hy=0;
300
301 if (pad->GetPad(1))
302 {
303 pad->GetPad(1)->cd(1);
304
305 if ((hx = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ex"))))
306 {
307 hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
308 hx = fHEnergy.Project3D("ex");
309 }
310
311 if ((hy = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ey"))))
312 {
313 hy->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
314 hy = fHEnergy.Project3D("ey");
315 }
316
317 if (hx && hy)
318 {
319 hx->SetLineColor(kBlue);
320 hx->SetMarkerColor(kBlue);
321 hy->SetMaximum();
322 hy->SetMaximum(TMath::Max(hx->GetMaximum(), hy->GetMaximum())*1.2);
323 if (hy->GetMaximum()>0)
324 gPad->SetLogy();
325 }
326
327 if (pad->GetPad(1)->GetPad(2))
328 {
329 pad->GetPad(1)->GetPad(2)->cd(1);
330 if ((hx = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ez"))))
331 {
332 hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
333 fHEnergy.Project3D("ez");
334 }
335
336 pad->GetPad(1)->GetPad(2)->cd(2);
337 if (gPad->FindObject("ResEnergyEst_py"))
338 {
339 hx = (TH1D*)fHResolutionEst.ProjectionY("_py", -1, -1, "e");
340 TPaveStats *stats = dynamic_cast<TPaveStats*>(hx->FindObject("stats"));
341 if (stats)
342 {
343 stats->SetBit(BIT(17)); // TStyle.cxx: kTakeStyle=BIT(17)
344 stats->SetX1NDC(0.63);
345 stats->SetY1NDC(0.68);
346 }
347
348 if (hx->GetEntries()>0)
349 {
350 hx->Fit("gaus", "Q", "", -0.25, 1.0);
351 if (hx->GetFunction("gaus"))
352 {
353 hx->GetFunction("gaus")->SetLineColor(kBlue);
354 hx->GetFunction("gaus")->SetLineWidth(2);
355 gPad=NULL;
356 gStyle->SetOptFit(101);
357 }
358 }
359 }
360 }
361 }
362
363 if (pad->GetPad(2))
364 {
365 pad->GetPad(2)->cd(1);
366 if (gPad->FindObject("EnergyEst_yx"))
367 {
368 TH2D *hyx = static_cast<TH2D*>(fHEnergy.Project3D("yx"));
369 UpdateProf(*hyx, kTRUE);
370 }
371
372 TLine *l = (TLine*)gPad->FindObject("TLine");
373 if (l)
374 {
375 const Float_t min = TMath::Max(fHEnergy.GetXaxis()->GetXmin(), fHEnergy.GetYaxis()->GetXmin());
376 const Float_t max = TMath::Min(fHEnergy.GetXaxis()->GetXmax(), fHEnergy.GetYaxis()->GetXmax());
377
378 l->SetX1(min);
379 l->SetX2(max);
380 l->SetY1(min);
381 l->SetY2(max);
382 }
383
384 pad->GetPad(2)->cd(2);
385 UpdateProf(fHResolutionMC, kFALSE);
386
387 pad->GetPad(2)->cd(3);
388 UpdateProf(fHResolutionEst, kFALSE);
389 }
390}
391
392void MHEnergyEst::UpdateProf(TH2 &h, Bool_t logy)
393{
394 const TString pname = MString::Format("Prof%s", h.GetName());
395
396 if (!gPad->FindObject(pname))
397 return;
398
399 TH1D *hx = h.ProfileX(pname, -1, -1, "s");
400 hx->SetLineColor(kBlue);
401 hx->SetMarkerColor(kBlue);
402
403 if (logy && hx->GetMaximum()>0)
404 gPad->SetLogy();
405}
406
407TH1 *MHEnergyEst::MakeProj(const char *how)
408{
409 TH1 *p = fHEnergy.Project3D(how);
410 p->SetDirectory(NULL);
411 p->SetBit(kCanDelete);
412 p->SetBit(TH1::kNoStats);
413 p->SetMarkerStyle(kFullDotMedium);
414 p->SetLineColor(kBlue);
415
416 return p;
417}
418
419TH1 *MHEnergyEst::MakeProf(TH2 &h)
420{
421 TH1 *p = h.ProfileX(MString::Format("Prof%s", h.GetName()), -1, -1, "s");
422 p->SetDirectory(NULL);
423 p->SetBit(kCanDelete);
424 p->SetLineWidth(2);
425 p->SetLineColor(kBlue);
426 p->SetFillStyle(4000);
427 p->SetStats(kFALSE);
428
429 return p;
430}
431
432// --------------------------------------------------------------------------
433//
434// Draw the histogram
435//
436void MHEnergyEst::Draw(Option_t *opt)
437{
438 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
439
440 // Do the projection before painting the histograms into
441 // the individual pads
442 AppendPad("");
443
444 pad->SetBorderMode(0);
445 pad->Divide(2, 1, 1e-10, 1e-10);
446
447 TH1 *h;
448
449 // ----------------------------------------
450
451 pad->cd(1);
452 gPad->SetBorderMode(0);
453
454 gPad->Divide(1, 2, 1e-10, 1e-10);
455
456 TVirtualPad *pad2 = gPad;
457
458 // ----------------------------------------
459
460 pad2->cd(1);
461 gPad->SetBorderMode(0);
462
463 gPad->SetGridx();
464 gPad->SetGridy();
465 gPad->SetLogx();
466
467 h = MakeProj("ey");
468 h->SetTitle("Energy disribution: Monte Carlo E_{mc} (black), Estimated E_{est} (blue)");
469 h->SetXTitle("E [GeV]"); // E_mc
470 h->SetYTitle("Counts");
471 h->Draw();
472
473 h->GetXaxis()->SetMoreLogLabels();
474 h->GetXaxis()->SetNoExponent();
475
476 h = MakeProj("ex");
477 h->SetLineColor(kBlue);
478 h->SetMarkerColor(kBlue);
479 h->Draw("same");
480
481 // ----------------------------------------
482
483 pad2->cd(2);
484 gPad->SetBorderMode(0);
485
486 TVirtualPad *pad3 = gPad;
487 pad3->Divide(2, 1, 1e-10, 1e-10);
488 pad3->cd(1);
489 gPad->SetBorderMode(0);
490 gPad->SetGridx();
491 gPad->SetGridy();
492
493 h = MakeProj("ez");
494 h->SetTitle("Zenith Angle Distribution");
495 h->GetXaxis()->SetMoreLogLabels();
496 h->GetXaxis()->SetNoExponent();
497 h->Draw();
498
499 // ----------------------------------------
500
501 pad3->cd(2);
502 gPad->SetBorderMode(0);
503 gPad->SetGridx();
504 gPad->SetGridy();
505
506 h = fHResolutionEst.ProjectionY("_py");
507 h->SetTitle("Distribution of \\Delta E/E_{est}");
508 h->SetDirectory(NULL);
509 h->SetBit(kCanDelete);
510 h->GetXaxis()->SetRangeUser(-1.3, 1.3);
511 h->Draw();
512 // ----------------------------------------
513
514 pad->cd(2);
515 gPad->SetBorderMode(0);
516
517 gPad->Divide(1, 3, 1e-10, 1e-10);
518 pad2 = gPad;
519
520 // ----------------------------------------
521
522 pad2->cd(1);
523 gPad->SetBorderMode(0);
524 gPad->SetLogy();
525 gPad->SetLogx();
526 gPad->SetGridx();
527 gPad->SetGridy();
528
529 // Results in crashes....
530 //gROOT->GetListOfCleanups()->Add(gPad); // WHY?
531
532 TH2D *h2 = (TH2D*)fHEnergy.Project3D("yx");
533 h2->SetDirectory(NULL);
534 h2->SetBit(kCanDelete);
535 h2->SetFillColor(kBlue);
536 h2->SetTitle("Estimated energy E_{mc} vs Monte Carlo energy E_{est}");
537
538 h2->Draw("");
539 MakeProf(*h2)->Draw("E0 same");
540
541 h2->GetXaxis()->SetMoreLogLabels();
542 h2->GetXaxis()->SetNoExponent();
543
544 TLine line;
545 line.DrawLine(0,0,1,1);
546
547 line.SetLineColor(kBlue);
548 line.SetLineWidth(2);
549 line.SetLineStyle(kDashed);
550
551 // ----------------------------------------
552
553 pad2->cd(2);
554 gPad->SetBorderMode(0);
555 gPad->SetLogx();
556 gPad->SetGridx();
557 gPad->SetGridy();
558 fHResolutionMC.Draw();
559 MakeProf(fHResolutionMC)->Draw("E0 same");
560
561 fHResolutionMC.GetXaxis()->SetMoreLogLabels();
562 fHResolutionMC.GetXaxis()->SetNoExponent();
563
564 line.DrawLine(fHResolutionMC.GetXaxis()->GetXmin(), 0, fHResolutionMC.GetXaxis()->GetXmax(), 0);
565
566 // ----------------------------------------
567
568 pad2->cd(3);
569 gPad->SetBorderMode(0);
570 gPad->SetLogx();
571 gPad->SetGridx();
572 gPad->SetGridy();
573 fHResolutionEst.Draw();
574 MakeProf(fHResolutionEst)->Draw("E0 same");
575
576 fHResolutionEst.GetXaxis()->SetMoreLogLabels();
577 fHResolutionEst.GetXaxis()->SetNoExponent();
578
579 line.DrawLine(fHResolutionEst.GetXaxis()->GetXmin(), 0, fHResolutionEst.GetXaxis()->GetXmax(), 0);
580}
581
582// --------------------------------------------------------------------------
583//
584// Returns the mapped value from the Matrix
585//
586Double_t MHEnergyEst::GetVal(Int_t i) const
587{
588 return (*fMatrix)[fMap[i]];
589}
590
591// --------------------------------------------------------------------------
592//
593// You can use this function if you want to use a MHMatrix instead of the
594// given containers. This function adds all necessary columns to the
595// given matrix. Afterward you should fill the matrix with the corresponding
596// data (eg from a file by using MHMatrix::Fill). If you now loop
597// through the matrix (eg using MMatrixLoop) MEnergyEstParam2::Process
598// will take the values from the matrix instead of the containers.
599//
600void MHEnergyEst::InitMapping(MHMatrix *mat)
601{
602 if (fMatrix)
603 return;
604
605 fMatrix = mat;
606
607 fMap[0] = fMatrix->AddColumn("MMcEvt.fEnergy");
608 fMap[1] = fMatrix->AddColumn("MMcEvt.fImpact/100");
609 fMap[2] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta*kRad2Deg");
610}
611
612void MHEnergyEst::StopMapping()
613{
614 fMatrix = NULL;
615}
Note: See TracBrowser for help on using the repository browser.