source: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc@ 9292

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