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

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