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

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