source: trunk/Mars/mhflux/MHEnergyEst.cc@ 9627

Last change on this file since 9627 was 9565, checked in by tbretz, 15 years ago
*** empty log message ***
File size: 16.7 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 *fLog << " Bias of Delta E/E: " << Form("%+4.2f%%", 100*f->GetParameter(1)) << endl;
257 *fLog << " Sigma of Delta E/E: " << Form("%4.2f%%", 100*f->GetParameter(2)) << endl;
258 *fLog << " Res of Delta E/E: " << Form("%4.2f%%", 100*TMath::Hypot(f->GetParameter(1), f->GetParameter(2))) << endl;
259
260 delete h;
261}
262
263// --------------------------------------------------------------------------
264//
265// Return Correction Coefficients (weights)
266// hist = E_mc/E_est
267//
268void MHEnergyEst::GetWeights(TH1D &hist) const
269{
270 // Project into EnergyEst_ey
271 // the "e" ensures that errors are calculated
272 TH1D *h1 = (TH1D*)fHEnergy.Project3D("rtlprmft_ex"); // E_Est
273 TH1D *h2 = (TH1D*)fHEnergy.Project3D("rtlprmft_ey"); // E_mc
274
275 h2->Copy(hist);
276 hist.Divide(h1);
277
278 delete h1;
279 delete h2;
280
281 hist.SetNameTitle("EnergyRatio", "Ratio between estimated and monte carlo energy");
282 hist.SetXTitle("E [GeV]");
283 hist.SetYTitle("N_{mc}/N_{est} [1]");
284 hist.SetDirectory(0);
285}
286
287void MHEnergyEst::Paint(Option_t *opt)
288{
289 TVirtualPad *pad = gPad;
290
291 pad->cd(1);
292
293 TH1 *hx=0;
294 TH1 *hy=0;
295
296 if (pad->GetPad(1))
297 {
298 pad->GetPad(1)->cd(1);
299
300 if ((hx = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ex"))))
301 {
302 hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
303 hx = fHEnergy.Project3D("ex");
304 }
305
306 if ((hy = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ey"))))
307 {
308 hy->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
309 hy = fHEnergy.Project3D("ey");
310 }
311
312 if (hx && hy)
313 {
314 hx->SetLineColor(kBlue);
315 hx->SetMarkerColor(kBlue);
316 hy->SetMaximum();
317 hy->SetMaximum(TMath::Max(hx->GetMaximum(), hy->GetMaximum())*1.2);
318 if (hy->GetMaximum()>0)
319 gPad->SetLogy();
320 }
321
322 if (pad->GetPad(1)->GetPad(2))
323 {
324 pad->GetPad(1)->GetPad(2)->cd(1);
325 if ((hx = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ez"))))
326 {
327 hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
328 fHEnergy.Project3D("ez");
329 }
330
331 pad->GetPad(1)->GetPad(2)->cd(2);
332 if (gPad->FindObject("ResEnergyEst_py"))
333 {
334 hx = (TH1D*)fHResolutionEst.ProjectionY("_py", -1, -1, "e");
335 TPaveStats *stats = dynamic_cast<TPaveStats*>(hx->FindObject("stats"));
336 if (stats)
337 {
338 stats->SetBit(BIT(17)); // TStyle.cxx: kTakeStyle=BIT(17)
339 stats->SetX1NDC(0.63);
340 stats->SetY1NDC(0.68);
341 }
342
343 hx->Fit("gaus", "Q", "", -0.25, 1.0);
344 hx->GetFunction("gaus")->SetLineColor(kBlue);
345 hx->GetFunction("gaus")->SetLineWidth(2);
346 gPad=NULL;
347 gStyle->SetOptFit(101);
348 }
349 }
350 }
351
352 if (pad->GetPad(2))
353 {
354 pad->GetPad(2)->cd(1);
355 if (gPad->FindObject("EnergyEst_yx"))
356 {
357 TH2D *hyx = static_cast<TH2D*>(fHEnergy.Project3D("yx"));
358 UpdateProf(*hyx, kTRUE);
359 }
360
361 TLine *l = (TLine*)gPad->FindObject("TLine");
362 if (l)
363 {
364 const Float_t min = TMath::Max(fHEnergy.GetXaxis()->GetXmin(), fHEnergy.GetYaxis()->GetXmin());
365 const Float_t max = TMath::Min(fHEnergy.GetXaxis()->GetXmax(), fHEnergy.GetYaxis()->GetXmax());
366
367 l->SetX1(min);
368 l->SetX2(max);
369 l->SetY1(min);
370 l->SetY2(max);
371 }
372
373 pad->GetPad(2)->cd(2);
374 UpdateProf(fHResolutionMC, kFALSE);
375
376 pad->GetPad(2)->cd(3);
377 UpdateProf(fHResolutionEst, kFALSE);
378 }
379}
380
381void MHEnergyEst::UpdateProf(TH2 &h, Bool_t logy)
382{
383 const TString pname = MString::Format("Prof%s", h.GetName());
384
385 if (!gPad->FindObject(pname))
386 return;
387
388 TH1D *hx = h.ProfileX(pname, -1, -1, "s");
389 hx->SetLineColor(kBlue);
390 hx->SetMarkerColor(kBlue);
391
392 if (logy && hx->GetMaximum()>0)
393 gPad->SetLogy();
394}
395
396TH1 *MHEnergyEst::MakeProj(const char *how)
397{
398 TH1 *p = fHEnergy.Project3D(how);
399 p->SetDirectory(NULL);
400 p->SetBit(kCanDelete);
401 p->SetBit(TH1::kNoStats);
402 p->SetMarkerStyle(kFullDotMedium);
403 p->SetLineColor(kBlue);
404
405 return p;
406}
407
408TH1 *MHEnergyEst::MakeProf(TH2 &h)
409{
410 TH1 *p = h.ProfileX(MString::Format("Prof%s", h.GetName()), -1, -1, "s");
411 p->SetDirectory(NULL);
412 p->SetBit(kCanDelete);
413 p->SetLineWidth(2);
414 p->SetLineColor(kBlue);
415 p->SetFillStyle(4000);
416 p->SetStats(kFALSE);
417
418 return p;
419}
420
421// --------------------------------------------------------------------------
422//
423// Draw the histogram
424//
425void MHEnergyEst::Draw(Option_t *opt)
426{
427 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
428
429 // Do the projection before painting the histograms into
430 // the individual pads
431 AppendPad("");
432
433 pad->SetBorderMode(0);
434 pad->Divide(2, 1, 1e-10, 1e-10);
435
436 TH1 *h;
437
438 // ----------------------------------------
439
440 pad->cd(1);
441 gPad->SetBorderMode(0);
442
443 gPad->Divide(1, 2, 1e-10, 1e-10);
444
445 TVirtualPad *pad2 = gPad;
446
447 // ----------------------------------------
448
449 pad2->cd(1);
450 gPad->SetBorderMode(0);
451
452 gPad->SetGridx();
453 gPad->SetGridy();
454 gPad->SetLogx();
455
456 h = MakeProj("ey");
457 h->SetTitle("Energy disribution: Monte Carlo E_{mc} (black), Estimated E_{est} (blue)");
458 h->SetXTitle("E [GeV]"); // E_mc
459 h->SetYTitle("Counts");
460 h->Draw();
461
462 h->GetXaxis()->SetMoreLogLabels();
463 h->GetXaxis()->SetNoExponent();
464
465 h = MakeProj("ex");
466 h->SetLineColor(kBlue);
467 h->SetMarkerColor(kBlue);
468 h->Draw("same");
469
470 // ----------------------------------------
471
472 pad2->cd(2);
473 gPad->SetBorderMode(0);
474
475 TVirtualPad *pad3 = gPad;
476 pad3->Divide(2, 1, 1e-10, 1e-10);
477 pad3->cd(1);
478 gPad->SetBorderMode(0);
479 gPad->SetGridx();
480 gPad->SetGridy();
481
482 h = MakeProj("ez");
483 h->SetTitle("Zenith Angle Distribution");
484 h->GetXaxis()->SetMoreLogLabels();
485 h->GetXaxis()->SetNoExponent();
486 h->Draw();
487
488 // ----------------------------------------
489
490 pad3->cd(2);
491 gPad->SetBorderMode(0);
492 gPad->SetGridx();
493 gPad->SetGridy();
494
495 h = fHResolutionEst.ProjectionY("_py");
496 h->SetTitle("Distribution of \\Delta E/E_{est}");
497 h->SetDirectory(NULL);
498 h->SetBit(kCanDelete);
499 h->GetXaxis()->SetRangeUser(-1.3, 1.3);
500 h->Draw();
501 // ----------------------------------------
502
503 pad->cd(2);
504 gPad->SetBorderMode(0);
505
506 gPad->Divide(1, 3, 1e-10, 1e-10);
507 pad2 = gPad;
508
509 // ----------------------------------------
510
511 pad2->cd(1);
512 gPad->SetBorderMode(0);
513 gPad->SetLogy();
514 gPad->SetLogx();
515 gPad->SetGridx();
516 gPad->SetGridy();
517
518 // Results in crashes....
519 //gROOT->GetListOfCleanups()->Add(gPad); // WHY?
520
521 TH2D *h2 = (TH2D*)fHEnergy.Project3D("yx");
522 h2->SetDirectory(NULL);
523 h2->SetBit(kCanDelete);
524 h2->SetFillColor(kBlue);
525 h2->SetTitle("Estimated energy E_{mc} vs Monte Carlo energy E_{est}");
526
527 h2->Draw("");
528 MakeProf(*h2)->Draw("E0 same");
529
530 h2->GetXaxis()->SetMoreLogLabels();
531 h2->GetXaxis()->SetNoExponent();
532
533 TLine line;
534 line.DrawLine(0,0,1,1);
535
536 line.SetLineColor(kBlue);
537 line.SetLineWidth(2);
538 line.SetLineStyle(kDashed);
539
540 // ----------------------------------------
541
542 pad2->cd(2);
543 gPad->SetBorderMode(0);
544 gPad->SetLogx();
545 gPad->SetGridx();
546 gPad->SetGridy();
547 fHResolutionMC.Draw();
548 MakeProf(fHResolutionMC)->Draw("E0 same");
549
550 fHResolutionMC.GetXaxis()->SetMoreLogLabels();
551 fHResolutionMC.GetXaxis()->SetNoExponent();
552
553 line.DrawLine(fHResolutionMC.GetXaxis()->GetXmin(), 0, fHResolutionMC.GetXaxis()->GetXmax(), 0);
554
555 // ----------------------------------------
556
557 pad2->cd(3);
558 gPad->SetBorderMode(0);
559 gPad->SetLogx();
560 gPad->SetGridx();
561 gPad->SetGridy();
562 fHResolutionEst.Draw();
563 MakeProf(fHResolutionEst)->Draw("E0 same");
564
565 fHResolutionEst.GetXaxis()->SetMoreLogLabels();
566 fHResolutionEst.GetXaxis()->SetNoExponent();
567
568 line.DrawLine(fHResolutionEst.GetXaxis()->GetXmin(), 0, fHResolutionEst.GetXaxis()->GetXmax(), 0);
569}
570
571// --------------------------------------------------------------------------
572//
573// Returns the mapped value from the Matrix
574//
575Double_t MHEnergyEst::GetVal(Int_t i) const
576{
577 return (*fMatrix)[fMap[i]];
578}
579
580// --------------------------------------------------------------------------
581//
582// You can use this function if you want to use a MHMatrix instead of the
583// given containers. This function adds all necessary columns to the
584// given matrix. Afterward you should fill the matrix with the corresponding
585// data (eg from a file by using MHMatrix::Fill). If you now loop
586// through the matrix (eg using MMatrixLoop) MEnergyEstParam2::Process
587// will take the values from the matrix instead of the containers.
588//
589void MHEnergyEst::InitMapping(MHMatrix *mat)
590{
591 if (fMatrix)
592 return;
593
594 fMatrix = mat;
595
596 fMap[0] = fMatrix->AddColumn("MMcEvt.fEnergy");
597 fMap[1] = fMatrix->AddColumn("MMcEvt.fImpact/100");
598 fMap[2] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta*kRad2Deg");
599}
600
601void MHEnergyEst::StopMapping()
602{
603 fMatrix = NULL;
604}
Note: See TracBrowser for help on using the repository browser.