source: trunk/MagicSoft/Mars/mhist/MHHadronness.cc@ 1948

Last change on this file since 1948 was 1948, checked in by wittek, 22 years ago
*** empty log message ***
File size: 15.3 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): Thomas Bretz, 5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Abelardo Moralejo <mailto:moralejo@pd.infn.it>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MHHadronness
29//
30// This is histogram is a way to evaluate the quality of a gamma/hadron
31// seperation method. It is filled from a MHadronness container, which
32// stores a hadroness for the current event. The Value must be in the
33// range [0,1]. To fill the histograms correctly the information
34// whether it is a gamma or hadron (not a gamma) must be available from
35// a MMcEvt container.
36//
37// In the constructor you can change the number of used bns for the
38// evaluation.
39//
40// The meaning of the histograms (Draw, DrawClone) are the following:
41// * Upper Left Corner:
42// - black: histogram of all hadronesses for gammas
43// - red: histogram of all hadronesses for non gammas
44// * Upper Right Corner:
45// - black: acceptance of gammas (Ag) vs. the hadroness
46// - red: acceptance of non gammas (Ah) vs. the hadroness
47// * Bottom Left Corner:
48// Naive quality factor: Ag/sqrt(Ah)
49// * Bottom Right Corner:
50// - black: Acceprtance Gammas vs. Acceptance Hadrons
51// - blue cross: minimum of distance to (0, 1)
52//
53// As a default MHHadronness searches for the container "MHadronness".
54// This can be overwritten by a different pointer specified in the
55// Fill function (No type check is done!) Use a different name in
56// MFillH.
57//
58////////////////////////////////////////////////////////////////////////////
59#include "MHHadronness.h"
60
61#include <TPad.h>
62#include <TGraph.h>
63#include <TStyle.h>
64#include <TCanvas.h>
65#include <TMarker.h>
66
67#include "MLog.h"
68#include "MLogManip.h"
69
70#include "MParList.h"
71#include "MBinning.h"
72#include "MHadronness.h"
73
74#include "MMcEvt.hxx"
75
76ClassImp(MHHadronness);
77
78// --------------------------------------------------------------------------
79//
80// Setup histograms, nbins is the number of bins used for the evaluation.
81// The default is 100 bins.
82//
83MHHadronness::MHHadronness(Int_t nbins, const char *name, const char *title)
84{
85 //
86 // set the name and title of this object
87 //
88 fName = name ? name : "MHHadronness";
89 fTitle = title ? title : "Gamma/Hadron Separation Quality Histograms";
90
91 fGraph = new TGraph;
92 fGraph->SetTitle("Acceptance Gammas vs. Hadrons");
93 fGraph->SetMaximum(1);
94
95 fGhness = new TH1D("Ghness", "Hadronness", nbins, 0, 1);
96 fPhness = new TH1D("Phness", "Hadronness", nbins, 0, 1);
97 fGhness->SetXTitle("Hadronness");
98 fPhness->SetXTitle("Hadronness");
99 fGhness->SetYTitle("Counts");
100 fPhness->SetYTitle("Counts");
101
102 fIntGhness = new TH1D("AccGammas", "Acceptance", nbins, 0, 1);
103 fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1);
104 fIntGhness->SetXTitle("Hadronness");
105 fIntPhness->SetXTitle("Hadronness");
106 fIntGhness->SetYTitle("Acceptance");
107 fIntPhness->SetYTitle("Acceptance");
108 fIntGhness->SetMaximum(1.1);
109 fIntPhness->SetMaximum(1.1);
110
111 fQfac = new TGraph;
112 fQfac->SetTitle(" Naive Quality factor ");
113
114 fGhness->SetDirectory(NULL);
115 fPhness->SetDirectory(NULL);
116 fIntGhness->SetDirectory(NULL);
117 fIntPhness->SetDirectory(NULL);
118}
119
120// --------------------------------------------------------------------------
121//
122// Delete the histograms.
123//
124MHHadronness::~MHHadronness()
125{
126 delete fGhness;
127 delete fIntGhness;
128 delete fPhness;
129 delete fIntPhness;
130 delete fQfac;
131 delete fGraph;
132}
133
134// --------------------------------------------------------------------------
135//
136// Setup Filling of the histograms. It needs:
137// MMcEvt and MHadronness
138//
139Bool_t MHHadronness::SetupFill(const MParList *plist)
140{
141 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
142 if (!fMcEvt)
143 {
144 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
145 return kFALSE;
146 }
147
148 fHadronness = (MHadronness*)plist->FindObject("MHadronness");
149
150 /*
151 MBinning* bins = (MBinning*)plist->FindObject("BinningHadronness");
152 if (!bins)
153 {
154 *fLog << err << dbginf << "BinningHadronness [MBinning] not found... aborting." << endl;
155 return kFALSE;
156 }
157
158 SetBinning(&fHist, binsalpha, binsenergy, binstheta);
159
160 fHist.Sumw2();
161 */
162
163 return kTRUE;
164}
165
166// --------------------------------------------------------------------------
167//
168// Fill the Hadronness from a MHadronness container into the corresponding
169// histogram dependant on the particle id.
170//
171// Every particle Id different than kGAMMA is considered a hadron.
172//
173// If you call Fill with a pointer (eg. from MFillH) this container is
174// used as a hadronness (Warning: No type check is done!) otherwise
175// the default is used ("MHadronness")
176//
177// Sometimes a distance is calculated as NaN (not a number). Such events
178// are skipped at the moment.
179//
180Bool_t MHHadronness::Fill(const MParContainer *par)
181{
182 // Preliminary Workaround: FIXME!
183 if (!par && !fHadronness)
184 {
185 *fLog << err << "MHHadronness::Fill: No MHadronness container specified!" << endl;
186 return kFALSE;
187 }
188
189 const MHadronness &had = par ? *(MHadronness*)par : *fHadronness;
190
191 const Double_t h = had.GetHadronness();
192
193 if (TMath::IsNaN(h))
194 return kCONTINUE;
195
196 if (fMcEvt->GetPartId()==kGAMMA)
197 fGhness->Fill(h);
198 else
199 fPhness->Fill(h);
200
201 return kTRUE;
202}
203
204Float_t MHHadronness::GetQ05() const
205{
206 Int_t n = fQfac->GetN();
207
208 Double_t val1x=0;
209 Double_t val2x=1;
210
211 Double_t val1y=0;
212 Double_t val2y=0;
213
214 for (Int_t i=1; i<=n; i++)
215 {
216 Double_t x, y;
217
218 fQfac->GetPoint(i, x, y);
219
220 if (x<0.5 && x>val1x)
221 {
222 val1x = x;
223 val1y = y;
224 }
225
226 if (x>0.5 && x<val2x)
227 {
228 val2x = x;
229 val2y = y;
230 }
231 }
232
233 return val2x-val1x == 0 ? 0 : val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
234}
235
236// --------------------------------------------------------------------------
237//
238// Finalize the histograms:
239// - integrate the hadroness histograms --> acceptance
240// - fill the Minimum Distance histogram (formular see class description)
241// - fill the Quality histogram (formular see class description)
242//
243Bool_t MHHadronness::Finalize()
244{
245 Int_t n = fGhness->GetNbinsX();
246
247 fGraph->Set(n);
248 fQfac->Set(n);
249
250 const Stat_t sumg = fGhness->Integral();
251 const Stat_t sump = fPhness->Integral();
252
253 *fLog << inf << "MHHadronness::Finalize - Sum Hadronness: gammas=" << sumg << " hadrons=" << sump << endl;
254
255 if (sumg==0 )
256 *fLog << warn << "MHHadronness::Finalize: Cannot calculate hadronness for 'gammas'." << endl;
257 if (sump==0)
258 *fLog << warn << "MHHadronness::Finalize: Cannot calculate hadronness for 'hadrons'." << endl;
259
260
261 // Normalize photon distribution
262 if (sumg>0)
263 fGhness->Scale(1./sumg);
264
265 // Normalize hadron distribution
266 if (sump>0)
267 fPhness->Scale(1./sump);
268
269 // Calculate acceptances
270 Float_t max=0;
271
272 for (Int_t i=1; i<=n; i++)
273 {
274 const Stat_t ip = fPhness->Integral(1, i);
275 const Stat_t ig = fGhness->Integral(1, i);
276
277 fIntPhness->SetBinContent(i, ip);
278 fIntGhness->SetBinContent(i, ig);
279
280 fGraph->SetPoint(i, ip, ig);
281
282 if (ip<=0)
283 continue;
284
285 const Double_t val = ig/sqrt(ip);
286 fQfac->SetPoint(i, ig, val);
287
288 if (val>max)
289 max = val;
290 }
291
292 fQfac->SetMaximum(max*1.05);
293
294 return kTRUE;
295}
296
297// --------------------------------------------------------------------------
298//
299// Search the corresponding points for the given hadron acceptance (acchad)
300// and interpolate the tow points (linear)
301//
302Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
303{
304 const Int_t n = fGraph->GetN();
305 const Double_t *x = fGraph->GetX();
306 const Double_t *y = fGraph->GetY();
307
308 Int_t i = 0;
309 while (i<n && x[i]<acchad)
310 i++;
311
312 if (i==0 || i==n)
313 return 0;
314
315 if (i==n-1)
316 i--;
317
318 const Double_t x1 = x[i-1];
319 const Double_t y1 = y[i-1];
320
321 const Double_t x2 = x[i];
322 const Double_t y2 = y[i];
323
324 return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
325}
326
327// --------------------------------------------------------------------------
328//
329// Search the corresponding points for the given gamma acceptance (accgam)
330// and interpolate the tow points (linear)
331//
332Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
333{
334 const Int_t n = fGraph->GetN();
335 const Double_t *x = fGraph->GetX();
336 const Double_t *y = fGraph->GetY();
337
338 Int_t i = 0;
339 while (i<n && y[i]<accgam)
340 i++;
341
342 if (i==0 || i==n)
343 return 0;
344
345 if (i==n-1)
346 i--;
347
348 const Double_t x1 = y[i-1];
349 const Double_t y1 = x[i-1];
350
351 const Double_t x2 = y[i];
352 const Double_t y2 = x[i];
353
354 return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
355}
356
357// --------------------------------------------------------------------------
358//
359// Print the corresponding Gammas Acceptance for a hadron acceptance of
360// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
361// acceptance and the corresponding acceptances and hadroness value is
362// printed, together with the maximum Q-factor.
363//
364void MHHadronness::Print(Option_t *) const
365{
366 *fLog << all;
367 *fLog << "Hadronness histograms:" << endl;
368 *fLog << "---------------------" << endl;
369
370 if (fGraph->GetN()==0)
371 {
372 *fLog << " <No Entries>" << endl;
373 return;
374 }
375
376 *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
377 *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h)" << endl <<endl;
378
379 *fLog << " 0.005 " << Form("%6.3f", GetGammaAcceptance(0.005)) << " " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
380 *fLog << " 0.02 " << Form("%6.3f", GetGammaAcceptance(0.02)) << " " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
381 *fLog << " 0.05 " << Form("%6.3f", GetGammaAcceptance(0.05)) << " " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
382 *fLog << " 0.1 " << Form("%6.3f", GetGammaAcceptance(0.1 )) << " " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
383 *fLog << " 0.2 " << Form("%6.3f", GetGammaAcceptance(0.2 )) << " " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
384 *fLog << " 0.3 " << Form("%6.3f", GetGammaAcceptance(0.3 )) << " " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
385 *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << " 0.1 " << endl;
386 *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << " 0.2 " << endl;
387 *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << " 0.3 " << endl;
388 *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << " 0.4 " << endl;
389 *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << " 0.5 " << endl;
390 *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << " 0.6 " << endl;
391 *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << " 0.7 " << endl;
392 *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << " 0.8 " << endl;
393 *fLog << endl;
394
395 *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
396 *fLog << " Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
397 *fLog << endl;
398}
399
400// --------------------------------------------------------------------------
401//
402// Draw clone of all histograms. (For the Meaning see class description)
403//
404TObject *MHHadronness::DrawClone(Option_t *opt) const
405{
406 if (fGraph->GetN()==0)
407 return NULL;
408
409 TCanvas &c = *MakeDefCanvas("Hadronness", fTitle);
410 c.Divide(2, 2);
411
412 gROOT->SetSelectedPad(NULL);
413
414 c.cd(1);
415 //gStyle->SetOptStat(10);
416 Getghness()->DrawCopy();
417 Getphness()->SetLineColor(kRed);
418 Getphness()->DrawCopy("same");
419
420 c.cd(2);
421 //gStyle->SetOptStat(0);
422 Getighness()->DrawCopy();
423 Getiphness()->SetLineColor(kRed);
424 Getiphness()->DrawCopy("same");
425
426 c.cd(3);
427 TGraph &g2 = (TGraph&)*fQfac->DrawClone("A*");
428 g2.SetBit(kCanDelete);
429 gPad->Modified();
430 gPad->Update();
431 if (g2.GetHistogram())
432 {
433 g2.GetXaxis()->SetRangeUser(0, 1);
434 g2.GetXaxis()->SetTitle("Acceptance Gammas");
435 g2.GetYaxis()->SetTitle("Quality");
436 g2.SetMarkerStyle(kFullDotSmall);
437 g2.Draw("P");
438
439 gPad->Modified();
440 gPad->Update();
441 }
442
443 c.cd(4);
444 gPad->Modified();
445 gPad->Update();
446 TGraph &g = (TGraph&)*fGraph->DrawClone("AC");
447 g.SetBit(kCanDelete);
448 gPad->Modified();
449 gPad->Update();
450 if (g.GetHistogram())
451 {
452 g.GetXaxis()->SetRangeUser(0, 1);
453 g.GetXaxis()->SetTitle("Acceptance Hadrons");
454 g.GetYaxis()->SetTitle("Acceptance Gammas");
455 g.SetMarkerStyle(kFullDotSmall);
456 g.Draw("P");
457
458 gPad->Modified();
459 gPad->Update();
460 }
461 /*
462 const Int_t h = fMinDist->GetMaximumBin();
463 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
464 fIntGhness->GetBinContent(h), kStar);
465 m->SetMarkerColor(kBlue);
466 m->SetBit(kCanDelete);
467 m->Draw();
468 */
469
470 gStyle->SetOptStat(1111);
471
472 return &c;
473}
474
475// --------------------------------------------------------------------------
476//
477// Draw all histograms. (For the Meaning see class description)
478//
479void MHHadronness::Draw(Option_t *)
480{
481 if (fGraph->GetN()==0)
482 return;
483
484 if (!gPad)
485 MakeDefCanvas("Hadronness", fTitle);
486
487 gPad->Divide(2, 2);
488
489 gPad->cd(1);
490 //gStyle->SetOptStat(10);
491 Getghness()->Draw();
492 Getphness()->SetLineColor(kRed);
493 Getphness()->Draw("same");
494
495 gPad->cd(2);
496 //gStyle->SetOptStat(0);
497 Getighness()->Draw();
498 Getiphness()->SetLineColor(kRed);
499 Getiphness()->Draw("same");
500
501 gPad->cd(3);
502 fQfac->Draw("A*");
503 gPad->Modified();
504 gPad->Update();
505 if (fQfac->GetHistogram())
506 {
507 fQfac->GetXaxis()->SetRangeUser(0, 1);
508 fQfac->GetXaxis()->SetTitle("Acceptance Gammas");
509 fQfac->GetYaxis()->SetTitle("Quality");
510 fQfac->SetMarkerStyle(kFullDotSmall);
511 fQfac->Draw("P");
512
513 gPad->Modified();
514 gPad->Update();
515 }
516 gPad->cd(4);
517 gPad->Modified();
518 gPad->Update();
519 fGraph->Draw("AC");
520 gPad->Modified();
521 gPad->Update();
522 if (fGraph->GetHistogram())
523 {
524 fGraph->GetXaxis()->SetRangeUser(0, 1);
525 fGraph->GetXaxis()->SetTitle("Acceptance Hadrons");
526 fGraph->GetYaxis()->SetTitle("Acceptance Gammas");
527 fGraph->SetMarkerStyle(kFullDotSmall);
528 fGraph->Draw("P");
529 gPad->Modified();
530 gPad->Update();
531 }
532}
Note: See TracBrowser for help on using the repository browser.