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

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