source: trunk/Mars/mhist/MHHadronness.cc@ 15369

Last change on this file since 15369 was 9856, checked in by tbretz, 14 years ago
Fixed scaling of the intgeral histogram in MHHadroness.
File size: 17.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-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: Acceptance 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// If you are using filtercuts which gives you only two discrete values
59// of the hadronness (0.25 and 0.75) you may want to use MHHadronness(2).
60// If you request Q05() from such a MHHadronness instance you will get
61// Acc_g/sqrt(Acc_h)
62//
63////////////////////////////////////////////////////////////////////////////
64#include "MHHadronness.h"
65
66#include <TH1.h>
67#include <TPad.h>
68#include <TMath.h>
69#include <TGraph.h>
70#include <TStyle.h>
71#include <TCanvas.h>
72#include <TMarker.h>
73
74#include "MParList.h"
75#include "MBinning.h"
76#include "MHMatrix.h"
77#include "MParameters.h"
78
79#include "MLog.h"
80#include "MLogManip.h"
81
82#include "MMcEvt.hxx"
83
84ClassImp(MHHadronness);
85
86using namespace std;
87
88// --------------------------------------------------------------------------
89//
90// Setup histograms, nbins is the number of bins used for the evaluation.
91// The default is 100 bins.
92//
93MHHadronness::MHHadronness(Int_t nbins, const char *name, const char *title)
94 : fMatrix(NULL)
95{
96 //
97 // set the name and title of this object
98 //
99 fName = name ? name : "MHHadronness";
100 fTitle = title ? title : "Gamma/Hadron Separation Quality Histograms";
101
102 fGraph = new TGraph;
103 fGraph->SetTitle("Acceptance Gammas vs. Hadrons");
104 fGraph->SetMarkerStyle(kFullDotSmall);
105
106 fGhness = new TH1D("Ghness", "Acceptance vs. Hadronness (Gammas)", nbins, 0, 1.0001);
107 fPhness = new TH1D("Phness", "Acceptance vs. Hadronness (Hadrons)", nbins, 0, 1.0001);
108 fGhness->SetXTitle("Hadronness");
109 fPhness->SetXTitle("Hadronness");
110 fGhness->SetYTitle("Acceptance");
111 fPhness->SetYTitle("Acceptance");
112 fPhness->SetLineColor(kRed);
113 fGhness->SetDirectory(NULL);
114 fPhness->SetDirectory(NULL);
115
116 fIntGhness = new TH1D("AccGammas", "Integral Acceptance vs. Hadronness (Gammas)", nbins, 0, 1);
117 fIntPhness = new TH1D("AccHadrons", "Integral Acceptance vs. Hadronness (Hadrons)", nbins, 0, 1);
118 fIntGhness->SetXTitle("Hadronness");
119 fIntPhness->SetXTitle("Hadronness");
120 fIntGhness->SetYTitle("Acceptance");
121 fIntPhness->SetYTitle("Acceptance");
122 fIntGhness->SetMaximum(1);
123 fIntPhness->SetMaximum(1);
124 fIntGhness->SetMinimum(0);
125 fIntPhness->SetMinimum(0);
126 fIntGhness->SetDirectory(NULL);
127 fIntPhness->SetDirectory(NULL);
128 fIntPhness->SetLineColor(kRed);
129 fIntGhness->SetBit(TH1::kNoStats);
130 fIntPhness->SetBit(TH1::kNoStats);
131
132 fQfac = new TGraph;
133 fQfac->SetTitle(" Naive Quality factor ");
134 fQfac->SetMarkerStyle(kFullDotSmall);
135}
136
137// --------------------------------------------------------------------------
138//
139// Delete the histograms.
140//
141MHHadronness::~MHHadronness()
142{
143 delete fGhness;
144 delete fIntGhness;
145 delete fPhness;
146 delete fIntPhness;
147 delete fQfac;
148 delete fGraph;
149}
150
151// --------------------------------------------------------------------------
152//
153// Setup Filling of the histograms. It needs:
154// MMcEvt and MHadronness
155//
156Bool_t MHHadronness::SetupFill(const MParList *plist)
157{
158 if (!fMatrix)
159 {
160 fMcEvt = (MMcEvt*)plist->FindObject(AddSerialNumber("MMcEvt"));
161 if (!fMcEvt)
162 {
163 *fLog << err << dbginf << AddSerialNumber("MMcEvt");
164 *fLog << " not found... aborting." << endl;
165 return kFALSE;
166 }
167 }
168
169 fHadronness = (MParameterD*)plist->FindObject("MHadronness");
170
171// fGhness->Reset();
172// fPhness->Reset();
173
174 /*
175 MBinning* bins = (MBinning*)plist->FindObject("BinningHadronness");
176 if (!bins)
177 {
178 *fLog << err << dbginf << "BinningHadronness [MBinning] not found... aborting." << endl;
179 return kFALSE;
180 }
181
182 SetBinning(&fHist, binsalpha, binsenergy, binstheta);
183
184 fHist.Sumw2();
185 */
186
187 return kTRUE;
188}
189
190// --------------------------------------------------------------------------
191//
192// Fill the Hadronness from a MHadronness container into the corresponding
193// histogram dependant on the particle id.
194//
195// Every particle Id different than kGAMMA is considered a hadron.
196//
197// If you call Fill with a pointer (eg. from MFillH) this container is
198// used as a hadronness (Warning: No type check is done!) otherwise
199// the default is used ("MHadronness")
200//
201// Sometimes a distance is calculated as NaN (not a number). Such events
202// are skipped at the moment.
203//
204Int_t MHHadronness::Fill(const MParContainer *par, const Stat_t w)
205{
206 // Preliminary Workaround: FIXME!
207 if (!par && !fHadronness)
208 {
209 *fLog << err << "MHHadronness::Fill: No MHadronness container specified!" << endl;
210 return kERROR;
211 }
212
213 const MParameterD &had = par ? *(MParameterD*)par : *fHadronness;
214
215 const Double_t h = TMath::Min(TMath::Max(had.GetVal(), 0.), 1.);
216
217 if (!TMath::Finite(h))
218 return kTRUE; // Use kCONTINUE with extreme care!
219
220 const Int_t particleid = fMatrix ? (Int_t)(*fMatrix)[fMap] : fMcEvt->GetPartId();
221
222 if (particleid==MMcEvt::kGAMMA)
223 fGhness->Fill(h, w);
224 else
225 fPhness->Fill(h, w);
226
227 return kTRUE;
228}
229
230// --------------------------------------------------------------------------
231//
232// Returns the quality factor at gamma acceptance 0.5.
233//
234// If the histogram containes only two bins we return the
235// naive quality: ag/sqrt(ah)
236// with ag the acceptance of gammas and ah the acceptance of hadrons.
237//
238// You can use this (nbins=2) in case of some kind of filter cuts giving
239// only a result: gamma yes/no (means discrete values of hadronness 0.25
240// or 0.75)
241//
242// FIXME: In the later case weights cannot be used!
243//
244Float_t MHHadronness::GetQ05() const
245{
246 if (fGhness->GetNbinsX()==2)
247 {
248 // acceptance of all gamma-like gammas (h<0.5)
249 const Double_t ig = fGhness->GetBinContent(1);
250
251 // acceptance of all gamma-like hadrons (h<0.5)
252 const Double_t ip = fPhness->GetBinContent(1);
253
254 if (ip==0)
255 return 0; // FIXME!
256
257 // naive quality factor
258 const Double_t q = ig / sqrt(ip);
259
260 *fLog << all << ip << "/" << ig << ": " << q << endl;
261
262 return q;
263 }
264
265 const Int_t n = fQfac->GetN();
266
267 Double_t val1x=0;
268 Double_t val2x=1;
269
270 Double_t val1y=0;
271 Double_t val2y=0;
272
273 for (Int_t i=1; i<=n; i++)
274 {
275 Double_t x, y;
276
277 fQfac->GetPoint(i, x, y);
278
279 if (x<0.5 && x>val1x)
280 {
281 val1x = x;
282 val1y = y;
283 }
284
285 if (x>0.5 && x<val2x)
286 {
287 val2x = x;
288 val2y = y;
289 }
290 }
291
292 //*fLog << dbg << val1x << "/" << val1y << " " << val2x << "/" << val2y << endl;
293
294 return val2x-val1x == 0 ? 0 : val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
295}
296
297// --------------------------------------------------------------------------
298//
299// Finalize the histograms:
300// - integrate the hadroness histograms --> acceptance
301// - fill the Minimum Distance histogram (formular see class description)
302// - fill the Quality histogram (formular see class description)
303//
304void MHHadronness::CalcGraph(Double_t sumg, Double_t sump)
305{
306 Int_t n = fGhness->GetNbinsX();
307
308 fGraph->Set(n);
309 fQfac->Set(n);
310
311 // Calculate acceptances
312 Float_t max=0;
313
314 for (Int_t i=1; i<=n; i++)
315 {
316 const Stat_t ip = fPhness->Integral(1, i)/sump;
317 const Stat_t ig = fGhness->Integral(1, i)/sumg;
318
319 fIntPhness->SetBinContent(i, ip);
320 fIntGhness->SetBinContent(i, ig);
321
322 fGraph->SetPoint(i, ip, ig);
323
324 if (ip<=0)
325 continue;
326
327 const Double_t val = ig/sqrt(ip);
328 fQfac->SetPoint(i, ig, val);
329
330 if (val>max)
331 max = val;
332 }
333
334 fQfac->SetMaximum(max*1.05);
335}
336
337Bool_t MHHadronness::Finalize()
338{
339 const Stat_t sumg = fGhness->Integral();
340 const Stat_t sump = fPhness->Integral();
341
342 *fLog << inf << "Sum Hadronness: gammas=" << sumg << " hadrons=" << sump << endl;
343
344 // Normalize photon distribution
345 if (sumg>0)
346 fGhness->Scale(1./sumg);
347 else
348 *fLog << warn << "Cannot calculate hadronness for 'gammas'." << endl;
349
350 // Normalize hadron distribution
351 if (sump>0)
352 fPhness->Scale(1./sump);
353 else
354 *fLog << warn << "Cannot calculate hadronness for 'hadrons'." << endl;
355
356 CalcGraph(1, 1);
357
358 return kTRUE;
359}
360
361void MHHadronness::Paint(Option_t *opt)
362{
363 Stat_t sumg = fGhness->Integral();
364 Stat_t sump = fPhness->Integral();
365
366 // Normalize photon distribution
367 if (sumg<=0)
368 sumg=1;
369
370 // Normalize hadron distribution
371 if (sump<=0)
372 sump=1;
373
374 CalcGraph(sumg, sump);
375}
376
377// --------------------------------------------------------------------------
378//
379// Search the corresponding points for the given hadron acceptance (acchad)
380// and interpolate the tow points (linear)
381//
382Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
383{
384 const Int_t n = fGraph->GetN();
385 const Double_t *x = fGraph->GetX();
386 const Double_t *y = fGraph->GetY();
387
388 Int_t i = 0;
389 while (i<n && x[i]<acchad)
390 i++;
391
392 if (i==0 || i==n)
393 return 0;
394
395 if (i==n-1)
396 i--;
397
398 const Double_t x1 = x[i-1];
399 const Double_t y1 = y[i-1];
400
401 const Double_t x2 = x[i];
402 const Double_t y2 = y[i];
403
404 return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
405}
406
407// --------------------------------------------------------------------------
408//
409// Search the corresponding points for the given gamma acceptance (accgam)
410// and interpolate the tow points (linear)
411//
412Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
413{
414 const Int_t n = fGraph->GetN();
415 const Double_t *x = fGraph->GetX();
416 const Double_t *y = fGraph->GetY();
417
418 Int_t i = 0;
419 while (i<n && y[i]<accgam)
420 i++;
421
422 if (i==0 || i==n)
423 return 0;
424
425 if (i==n-1)
426 i--;
427
428 const Double_t x1 = y[i-1];
429 const Double_t y1 = x[i-1];
430
431 const Double_t x2 = y[i];
432 const Double_t y2 = x[i];
433
434 return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
435}
436
437// --------------------------------------------------------------------------
438//
439// Search the hadronness corresponding to a given hadron acceptance.
440//
441Double_t MHHadronness::GetHadronness(Double_t acchad) const
442{
443 for (int i=1; i<fIntPhness->GetNbinsX()+1; i++)
444 if (fIntPhness->GetBinContent(i)>acchad)
445 return fIntPhness->GetBinLowEdge(i);
446
447 return -1;
448}
449
450// --------------------------------------------------------------------------
451//
452// Print the corresponding Gammas Acceptance for a hadron acceptance of
453// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
454// acceptance and the corresponding acceptances and hadroness value is
455// printed, together with the maximum Q-factor.
456//
457void MHHadronness::Print(Option_t *) const
458{
459 *fLog << all;
460 *fLog << underline << GetDescriptor() << endl;
461
462 if (fGraph->GetN()==0)
463 {
464 *fLog << " <No Entries>" << endl;
465 return;
466 }
467
468 *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
469 *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h) h" << endl <<endl;
470
471 *fLog << " 0.005 " << Form("%6.3f", GetGammaAcceptance(0.005)) << " " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << " " << GetHadronness(0.005) << endl;
472 *fLog << " 0.02 " << Form("%6.3f", GetGammaAcceptance(0.02)) << " " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << " " << GetHadronness(0.02) << endl;
473 *fLog << " 0.05 " << Form("%6.3f", GetGammaAcceptance(0.05)) << " " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << " " << GetHadronness(0.05) << endl;
474 *fLog << " 0.1 " << Form("%6.3f", GetGammaAcceptance(0.1 )) << " " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << " " << GetHadronness(0.1) << endl;
475 *fLog << " 0.2 " << Form("%6.3f", GetGammaAcceptance(0.2 )) << " " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << " " << GetHadronness(0.2) << endl;
476 *fLog << " 0.3 " << Form("%6.3f", GetGammaAcceptance(0.3 )) << " " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << " " << GetHadronness(0.3) << endl;
477 *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << " 0.1 " << endl;
478 *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << " 0.2 " << endl;
479 *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << " 0.3 " << endl;
480 *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << " 0.4 " << endl;
481 *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << " 0.5 " << endl;
482 *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << " 0.6 " << endl;
483 *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << " 0.7 " << endl;
484 *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << " 0.8 " << endl;
485 *fLog << endl;
486
487 *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
488 *fLog << " Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
489 *fLog << endl;
490}
491
492// --------------------------------------------------------------------------
493//
494// Draw all histograms. (For the Meaning see class description)
495//
496void MHHadronness::Draw(Option_t *)
497{
498 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas("Hadronness", fTitle);
499 pad->SetBorderMode(0);
500
501 AppendPad("");
502
503 pad->Divide(2, 2);
504
505 TH1 *h;
506
507 pad->cd(1);
508 gPad->SetBorderMode(0);
509 //gStyle->SetOptStat(10);
510 MH::DrawSame(*fGhness, *fPhness, "Hadronness"); // Displ both stat boxes
511
512 pad->cd(2);
513 gPad->SetBorderMode(0);
514 gPad->SetGridx();
515 gPad->SetGridy();
516 fIntGhness->Draw();
517 fIntPhness->Draw("same");
518
519 pad->cd(3);
520 gPad->SetBorderMode(0);
521 gPad->SetGridx();
522 gPad->SetGridy();
523 fQfac->Draw("A*");
524 gPad->Modified();
525 gPad->Update();
526 if ((h=fQfac->GetHistogram()))
527 {
528 h->GetXaxis()->SetRangeUser(0, 1);
529 h->SetXTitle("Acceptance Gammas");
530 h->SetYTitle("Quality");
531 fQfac->Draw("P");
532 gPad->Modified();
533 gPad->Update();
534 }
535
536 pad->cd(4);
537 gPad->SetBorderMode(0);
538 gPad->SetGridx();
539 gPad->SetGridy();
540 fGraph->Draw("AC");
541 gPad->Modified();
542 gPad->Update();
543 if ((h=fGraph->GetHistogram()))
544 {
545 h->GetXaxis()->SetRangeUser(0, 1);
546 h->SetXTitle("Acceptance Hadrons");
547 h->SetYTitle("Acceptance Gammas");
548 fGraph->SetMaximum(1);
549 fGraph->Draw("P");
550 gPad->Modified();
551 gPad->Update();
552 }
553}
554
555// --------------------------------------------------------------------------
556//
557// You can use this function if you want to use a MHMatrix instead of
558// MMcEvt. This function adds all necessary columns to the
559// given matrix. Afterward you should fill the matrix with the corresponding
560// data (eg from a file by using MHMatrix::Fill). If you now loop
561// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
562// will take the values from the matrix instead of the containers.
563//
564void MHHadronness::InitMapping(MHMatrix *mat)
565{
566 if (fMatrix)
567 return;
568
569 fMatrix = mat;
570
571 TString str = AddSerialNumber("MMcEvt");
572 str += ".fPartId";
573
574 fMap = fMatrix->AddColumn(str);
575}
576
577void MHHadronness::StopMapping()
578{
579 fMatrix = NULL;
580}
Note: See TracBrowser for help on using the repository browser.