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

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