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

Last change on this file since 1654 was 1610, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 15.1 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-2002
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// - blue: 2D distance of (acceptance_hadrons, acceptances_gammas)
48// to optimum (0, 1)
49// 1-sqrt(Ag*Ag + (1-Ah)*(1-Ah))
50// * Bottom Left Corner:
51// Naive quality factor: Ag/sqrt(Ah)
52// * Bottom Right Corner:
53// - black: Acceprtance Gammas vs. Acceptance Hadrons
54// - blue cross: minimum of distance to (0, 1)
55//
56////////////////////////////////////////////////////////////////////////////
57#include "MHHadronness.h"
58
59#include <TPad.h>
60#include <TGraph.h>
61#include <TStyle.h>
62#include <TCanvas.h>
63#include <TMarker.h>
64
65#include "MLog.h"
66#include "MLogManip.h"
67
68#include "MParList.h"
69#include "MBinning.h"
70#include "MHadronness.h"
71
72#include "MMcEvt.hxx"
73
74ClassImp(MHHadronness);
75
76// --------------------------------------------------------------------------
77//
78// Setup histograms, nbins is the number of bins used for the evaluation.
79// The default is 100 bins.
80//
81MHHadronness::MHHadronness(Int_t nbins, const char *name, const char *title)
82{
83 //
84 // set the name and title of this object
85 //
86 fName = name ? name : "MHHadronness";
87 fTitle = title ? title : "Gamma/Hadron Separation Quality Histograms";
88
89 fGraph = new TGraph;
90 fGraph->SetTitle("Acceptance Gammas vs. Hadrons");
91 fGraph->SetMaximum(1);
92
93 fGhness = new TH1D("Ghness", "Hadronness", nbins, 0, 1);
94 fPhness = new TH1D("Phness", "Hadronness", nbins, 0, 1);
95 fGhness->SetXTitle("Hadronness");
96 fPhness->SetXTitle("Hadronness");
97 fGhness->SetYTitle("Counts");
98 fPhness->SetYTitle("Counts");
99
100 fIntGhness = new TH1D("AccGammas", "Acceptance", nbins, 0, 1);
101 fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1);
102 fIntGhness->SetXTitle("Hadronness");
103 fIntPhness->SetXTitle("Hadronness");
104 fIntGhness->SetYTitle("Acceptance");
105 fIntPhness->SetYTitle("Acceptance");
106
107 /*
108 fQfac = new TH1D("Qfac", "Naive Quality factor", nbins, 0, 1);
109 fQfac->SetXTitle("Hadronness");
110 fQfac->SetYTitle("Quality");
111 */
112 fQfac = new TGraph;
113 fQfac->SetTitle(" Naive Quality factor ");
114
115 fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1);
116 fMinDist->SetXTitle("Hadronness");
117 fMinDist->SetYTitle("Distance");
118
119 //fQfac->SetDirectory(NULL);
120 fGhness->SetDirectory(NULL);
121 fPhness->SetDirectory(NULL);
122 fMinDist->SetDirectory(NULL);
123 fIntGhness->SetDirectory(NULL);
124 fIntPhness->SetDirectory(NULL);
125}
126
127// --------------------------------------------------------------------------
128//
129// Delete the histograms.
130//
131MHHadronness::~MHHadronness()
132{
133 delete fGhness;
134 delete fIntGhness;
135 delete fPhness;
136 delete fIntPhness;
137 delete fQfac;
138 delete fMinDist;
139 delete fGraph;
140}
141
142// --------------------------------------------------------------------------
143//
144// Setup Filling of the histograms. It needs:
145// MMcEvt and MHadronness
146//
147Bool_t MHHadronness::SetupFill(const MParList *plist)
148{
149 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
150 if (!fMcEvt)
151 {
152 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
153 return kFALSE;
154 }
155
156 fHadronness = (MHadronness*)plist->FindObject("MHadronness");
157 if (!fHadronness)
158 {
159 *fLog << err << dbginf << "MHadronness not found... aborting." << endl;
160 return kFALSE;
161 }
162
163 /*
164 MBinning* bins = (MBinning*)plist->FindObject("BinningHadronness");
165 if (!bins)
166 {
167 *fLog << err << dbginf << "BinningHadronness [MBinning] not found... aborting." << endl;
168 return kFALSE;
169 }
170
171 SetBinning(&fHist, binsalpha, binsenergy, binstheta);
172
173 fHist.Sumw2();
174 */
175
176 return kTRUE;
177}
178
179// --------------------------------------------------------------------------
180//
181// Fill the Hadronness from a MHadronness container into the corresponding
182// histogram dependant on the particle id.
183//
184// Sometimes a distance is calculated as NaN (not a number). Such events
185// are skipped at the moment.
186//
187Bool_t MHHadronness::Fill(const MParContainer *par)
188{
189 // Preliminary Workaround: FIXME!
190
191 const Double_t h = fHadronness->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 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(1, n+1);
251 const Stat_t sump = fPhness->Integral(1, n+1);
252
253 Float_t max=0;
254
255 for (Int_t i=1; i<=n; i++)
256 {
257 const Stat_t ip = fPhness->Integral(1, i)/sump;
258 const Stat_t ig = fGhness->Integral(1, i)/sumg;
259
260 fIntPhness->SetBinContent(i, ip);
261 fIntGhness->SetBinContent(i, ig);
262
263 fGraph->SetPoint(i, ip, ig);
264
265 if (ip<=0)
266 continue;
267
268 fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));
269
270 Double_t val = ig/sqrt(ip);
271 fQfac->SetPoint(i, ig, val);
272
273 if (val>max)
274 max = val;
275 }
276
277 fQfac->SetMaximum(max*1.05);
278
279 return kTRUE;
280}
281
282// --------------------------------------------------------------------------
283//
284// Search the corresponding points for the given hadron acceptance (acchad)
285// and interpolate the tow points (linear)
286//
287Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
288{
289 const Int_t n = fGraph->GetN();
290 const Double_t *x = fGraph->GetX();
291 const Double_t *y = fGraph->GetY();
292
293 Int_t i = 0;
294 while (i<n && x[i]<acchad)
295 i++;
296
297 if (i==0 || i==n)
298 return 0;
299
300 if (i==n-1)
301 i--;
302
303 const Double_t x1 = x[i-1];
304 const Double_t y1 = y[i-1];
305
306 const Double_t x2 = x[i];
307 const Double_t y2 = y[i];
308
309 return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
310}
311
312// --------------------------------------------------------------------------
313//
314// Search the corresponding points for the given gamma acceptance (accgam)
315// and interpolate the tow points (linear)
316//
317Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
318{
319 const Int_t n = fGraph->GetN();
320 const Double_t *x = fGraph->GetX();
321 const Double_t *y = fGraph->GetY();
322
323 Int_t i = 0;
324 while (i<n && y[i]<accgam)
325 i++;
326
327 if (i==0 || i==n)
328 return 0;
329
330 if (i==n-1)
331 i--;
332
333 const Double_t x1 = y[i-1];
334 const Double_t y1 = x[i-1];
335
336 const Double_t x2 = y[i];
337 const Double_t y2 = x[i];
338
339 return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
340}
341
342// --------------------------------------------------------------------------
343//
344// Print the corresponding Gammas Acceptance for a hadron acceptance of
345// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
346// acceptance and the corresponding acceptances and hadroness value is
347// printed, together with the maximum Q-factor.
348//
349void MHHadronness::Print(Option_t *) const
350{
351 *fLog << all;
352 *fLog << "Hadronness histograms:" << endl;
353 *fLog << "---------------------" << endl;
354
355 if (fGraph->GetN()==0)
356 {
357 *fLog << " <No Entries>" << endl;
358 return;
359 }
360
361 *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
362 *fLog << " acc(hadron) acc(gamma)" << endl <<endl;
363
364 *fLog << " 0.005 " << Form("%6.3f", GetGammaAcceptance(0.005)) << endl;
365 *fLog << " 0.02 " << Form("%6.3f", GetGammaAcceptance(0.02)) << endl;
366 *fLog << " 0.05 " << Form("%6.3f", GetGammaAcceptance(0.05)) << endl;
367 *fLog << " 0.1 " << Form("%6.3f", GetGammaAcceptance(0.1 )) << endl;
368 *fLog << " 0.2 " << Form("%6.3f", GetGammaAcceptance(0.2 )) << endl;
369 *fLog << " 0.3 " << Form("%6.3f", GetGammaAcceptance(0.3 )) << endl;
370 *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << " 0.1 " << endl;
371 *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << " 0.2 " << endl;
372 *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << " 0.3 " << endl;
373 *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << " 0.4 " << endl;
374 *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << " 0.5 " << endl;
375 *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << " 0.6 " << endl;
376 *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << " 0.7 " << endl;
377 *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << " 0.8 " << endl;
378 *fLog << endl;
379
380 const Int_t h = fMinDist->GetMaximumBin();
381 *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
382 *fLog << " Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
383 *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;
384
385 *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
386 *fLog << " Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
387 *fLog << endl;
388}
389
390// --------------------------------------------------------------------------
391//
392// Draw clone of all histograms. (For the Meaning see class description)
393//
394TObject *MHHadronness::DrawClone(Option_t *opt) const
395{
396 if (fGraph->GetN()==0)
397 return NULL;
398
399 TCanvas &c = *MakeDefCanvas("Hadronness", fTitle);
400 c.Divide(2, 2);
401
402 gROOT->SetSelectedPad(NULL);
403
404 c.cd(1);
405 gStyle->SetOptStat(10);
406 Getghness()->DrawCopy();
407 Getphness()->SetLineColor(kRed);
408 Getphness()->DrawCopy("same");
409 c.cd(2);
410 gStyle->SetOptStat(0);
411 Getighness()->DrawCopy();
412 Getiphness()->SetLineColor(kRed);
413 Getiphness()->DrawCopy("same");
414 fMinDist->SetLineColor(kBlue);
415 fMinDist->DrawCopy("same");
416 c.cd(3);
417 //GetQfac()->DrawCopy();
418 TGraph &g2 = (TGraph&)*fQfac->DrawClone("A*");
419 g2.SetBit(kCanDelete);
420 gPad->Modified();
421 gPad->Update();
422 if (g2.GetHistogram())
423 {
424 g2.GetXaxis()->SetRangeUser(0, 1);
425 g2.GetXaxis()->SetTitle("Acceptance Gammas");
426 g2.GetYaxis()->SetTitle("Quality");
427 g2.SetMarkerStyle(kFullDotSmall);
428 g2.Draw("P");
429
430 gPad->Modified();
431 gPad->Update();
432 }
433 c.cd(4);
434 gPad->Modified();
435 gPad->Update();
436 TGraph &g = (TGraph&)*fGraph->DrawClone("AC");
437 g.SetBit(kCanDelete);
438 gPad->Modified();
439 gPad->Update();
440 if (g.GetHistogram())
441 {
442 g.GetXaxis()->SetRangeUser(0, 1);
443 g.GetXaxis()->SetTitle("Acceptance Hadrons");
444 g.GetYaxis()->SetTitle("Acceptance Gammas");
445 g.SetMarkerStyle(kFullDotSmall);
446 g.Draw("P");
447
448 gPad->Modified();
449 gPad->Update();
450 }
451 const Int_t h = fMinDist->GetMaximumBin();
452 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
453 fIntGhness->GetBinContent(h), kStar);
454 m->SetMarkerColor(kBlue);
455 m->SetBit(kCanDelete);
456 m->Draw();
457 return &c;
458}
459
460// --------------------------------------------------------------------------
461//
462// Draw all histograms. (For the Meaning see class description)
463//
464void MHHadronness::Draw(Option_t *)
465{
466 if (fGraph->GetN()==0)
467 return;
468
469 if (!gPad)
470 MakeDefCanvas("Hadronness", fTitle);
471
472 gPad->Divide(2, 2);
473
474 gPad->cd(1);
475 gStyle->SetOptStat(10);
476 Getghness()->Draw();
477 Getphness()->SetLineColor(kRed);
478 Getphness()->Draw("same");
479 gPad->cd(2);
480 gStyle->SetOptStat(0);
481 Getighness()->Draw();
482 Getiphness()->SetLineColor(kRed);
483 Getiphness()->Draw("same");
484 fMinDist->SetLineColor(kBlue);
485 fMinDist->Draw("same");
486 gPad->cd(3);
487 //GetQfac()->Draw();
488 fQfac->Draw("A*");
489 gPad->Modified();
490 gPad->Update();
491 if (fQfac->GetHistogram())
492 {
493 fQfac->GetXaxis()->SetRangeUser(0, 1);
494 fQfac->GetXaxis()->SetTitle("Acceptance Gammas");
495 fQfac->GetYaxis()->SetTitle("Quality");
496 fQfac->SetMarkerStyle(kFullDotSmall);
497 fQfac->Draw("P");
498
499 gPad->Modified();
500 gPad->Update();
501 }
502 gPad->cd(4);
503 gPad->Modified();
504 gPad->Update();
505 fGraph->Draw("AC");
506 gPad->Modified();
507 gPad->Update();
508 if (fGraph->GetHistogram())
509 {
510 fGraph->GetXaxis()->SetRangeUser(0, 1);
511 fGraph->GetXaxis()->SetTitle("Acceptance Hadrons");
512 fGraph->GetYaxis()->SetTitle("Acceptance Gammas");
513 fGraph->SetMarkerStyle(kFullDotSmall);
514 fGraph->Draw("P");
515 gPad->Modified();
516 gPad->Update();
517 }
518 const Int_t h = fMinDist->GetMaximumBin();
519 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
520 fIntGhness->GetBinContent(h), kStar);
521 m->SetMarkerColor(kBlue);
522 m->SetBit(kCanDelete);
523 m->Draw();
524}
Note: See TracBrowser for help on using the repository browser.