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

Last change on this file since 1925 was 1921, checked in by wittek, 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-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
109 fQfac = new TGraph;
110 fQfac->SetTitle(" Naive Quality factor ");
111
112 fGhness->SetDirectory(NULL);
113 fPhness->SetDirectory(NULL);
114 fIntGhness->SetDirectory(NULL);
115 fIntPhness->SetDirectory(NULL);
116}
117
118// --------------------------------------------------------------------------
119//
120// Delete the histograms.
121//
122MHHadronness::~MHHadronness()
123{
124 delete fGhness;
125 delete fIntGhness;
126 delete fPhness;
127 delete fIntPhness;
128 delete fQfac;
129 delete fGraph;
130}
131
132// --------------------------------------------------------------------------
133//
134// Setup Filling of the histograms. It needs:
135// MMcEvt and MHadronness
136//
137Bool_t MHHadronness::SetupFill(const MParList *plist)
138{
139 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
140 if (!fMcEvt)
141 {
142 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
143 return kFALSE;
144 }
145
146 fHadronness = (MHadronness*)plist->FindObject("MHadronness");
147
148 /*
149 MBinning* bins = (MBinning*)plist->FindObject("BinningHadronness");
150 if (!bins)
151 {
152 *fLog << err << dbginf << "BinningHadronness [MBinning] not found... aborting." << endl;
153 return kFALSE;
154 }
155
156 SetBinning(&fHist, binsalpha, binsenergy, binstheta);
157
158 fHist.Sumw2();
159 */
160
161 return kTRUE;
162}
163
164// --------------------------------------------------------------------------
165//
166// Fill the Hadronness from a MHadronness container into the corresponding
167// histogram dependant on the particle id.
168//
169// Every particle Id different than kGAMMA is considered a hadron.
170//
171// If you call Fill with a pointer (eg. from MFillH) this container is
172// used as a hadronness (Warning: No type check is done!) otherwise
173// the default is used ("MHadronness")
174//
175// Sometimes a distance is calculated as NaN (not a number). Such events
176// are skipped at the moment.
177//
178Bool_t MHHadronness::Fill(const MParContainer *par)
179{
180 // Preliminary Workaround: FIXME!
181 if (!par && !fHadronness)
182 {
183 *fLog << err << "MHHadronness::Fill: No MHadronness container specified!" << endl;
184 return kFALSE;
185 }
186
187 const MHadronness &had = par ? *(MHadronness*)par : *fHadronness;
188
189 const Double_t h = had.GetHadronness();
190
191 if (TMath::IsNaN(h))
192 return kCONTINUE;
193
194 if (fMcEvt->GetPartId()==kGAMMA)
195 fGhness->Fill(h);
196 else
197 fPhness->Fill(h);
198
199 return kTRUE;
200}
201
202Float_t MHHadronness::GetQ05() const
203{
204 Int_t n = fQfac->GetN();
205
206 Double_t val1x=0;
207 Double_t val2x=1;
208
209 Double_t val1y=0;
210 Double_t val2y=0;
211
212 for (Int_t i=1; i<=n; i++)
213 {
214 Double_t x, y;
215
216 fQfac->GetPoint(i, x, y);
217
218 if (x<0.5 && x>val1x)
219 {
220 val1x = x;
221 val1y = y;
222 }
223
224 if (x>0.5 && x<val2x)
225 {
226 val2x = x;
227 val2y = y;
228 }
229 }
230
231 return val2x-val1x == 0 ? 0 : val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5);
232}
233
234// --------------------------------------------------------------------------
235//
236// Finalize the histograms:
237// - integrate the hadroness histograms --> acceptance
238// - fill the Minimum Distance histogram (formular see class description)
239// - fill the Quality histogram (formular see class description)
240//
241Bool_t MHHadronness::Finalize()
242{
243 Int_t n = fGhness->GetNbinsX();
244
245 fGraph->Set(n);
246 fQfac->Set(n);
247
248 const Stat_t sumg = fGhness->Integral();
249 const Stat_t sump = fPhness->Integral();
250
251 *fLog << inf << "MHHadronness::Finalize - Sum Hadronness: gammas=" << sumg << " hadrons=" << sump << endl;
252
253 if (sumg==0 || sump==0)
254 *fLog << warn << "MHHadronness::Finalize: Cannot calculate hadronness." << endl;
255
256
257 // Normalize photon distribution
258 if (sumg>0)
259 fGhness->Scale(1./sumg);
260
261 // Normalize hadron distribution
262 if (sump>0)
263 fGhness->Scale(1./sump);
264
265 // Calculate acceptances
266 Float_t max=0;
267
268
269 for (Int_t i=1; i<=n; i++)
270 {
271 const Stat_t ip = sump!=0 ? fPhness->Integral(1, i) : 0;
272 const Stat_t ig = sumg!=0 ? fGhness->Integral(1, i) : 0;
273
274 fIntPhness->SetBinContent(i, ip);
275 fIntGhness->SetBinContent(i, ig);
276
277 fGraph->SetPoint(i, ip, ig);
278
279 if (ip<=0)
280 continue;
281
282 const Double_t val = ig/sqrt(ip);
283 fQfac->SetPoint(i, ig, val);
284
285 if (val>max)
286 max = val;
287 }
288
289 fQfac->SetMaximum(max*1.05);
290
291 return kTRUE;
292}
293
294// --------------------------------------------------------------------------
295//
296// Search the corresponding points for the given hadron acceptance (acchad)
297// and interpolate the tow points (linear)
298//
299Double_t MHHadronness::GetGammaAcceptance(Double_t acchad) const
300{
301 const Int_t n = fGraph->GetN();
302 const Double_t *x = fGraph->GetX();
303 const Double_t *y = fGraph->GetY();
304
305 Int_t i = 0;
306 while (i<n && x[i]<acchad)
307 i++;
308
309 if (i==0 || i==n)
310 return 0;
311
312 if (i==n-1)
313 i--;
314
315 const Double_t x1 = x[i-1];
316 const Double_t y1 = y[i-1];
317
318 const Double_t x2 = x[i];
319 const Double_t y2 = y[i];
320
321 return (y2-y1)/(x2-x1) * (acchad-x2) + y2;
322}
323
324// --------------------------------------------------------------------------
325//
326// Search the corresponding points for the given gamma acceptance (accgam)
327// and interpolate the tow points (linear)
328//
329Double_t MHHadronness::GetHadronAcceptance(Double_t accgam) const
330{
331 const Int_t n = fGraph->GetN();
332 const Double_t *x = fGraph->GetX();
333 const Double_t *y = fGraph->GetY();
334
335 Int_t i = 0;
336 while (i<n && y[i]<accgam)
337 i++;
338
339 if (i==0 || i==n)
340 return 0;
341
342 if (i==n-1)
343 i--;
344
345 const Double_t x1 = y[i-1];
346 const Double_t y1 = x[i-1];
347
348 const Double_t x2 = y[i];
349 const Double_t y2 = x[i];
350
351 return (y2-y1)/(x2-x1) * (accgam-x2) + y2;
352}
353
354// --------------------------------------------------------------------------
355//
356// Print the corresponding Gammas Acceptance for a hadron acceptance of
357// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
358// acceptance and the corresponding acceptances and hadroness value is
359// printed, together with the maximum Q-factor.
360//
361void MHHadronness::Print(Option_t *) const
362{
363 *fLog << all;
364 *fLog << "Hadronness histograms:" << endl;
365 *fLog << "---------------------" << endl;
366
367 if (fGraph->GetN()==0)
368 {
369 *fLog << " <No Entries>" << endl;
370 return;
371 }
372
373 *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
374 *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h)" << endl <<endl;
375
376 *fLog << " 0.005 " << Form("%6.3f", GetGammaAcceptance(0.005)) << " " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
377 *fLog << " 0.02 " << Form("%6.3f", GetGammaAcceptance(0.02)) << " " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
378 *fLog << " 0.05 " << Form("%6.3f", GetGammaAcceptance(0.05)) << " " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
379 *fLog << " 0.1 " << Form("%6.3f", GetGammaAcceptance(0.1 )) << " " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
380 *fLog << " 0.2 " << Form("%6.3f", GetGammaAcceptance(0.2 )) << " " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
381 *fLog << " 0.3 " << Form("%6.3f", GetGammaAcceptance(0.3 )) << " " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
382 *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << " 0.1 " << endl;
383 *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << " 0.2 " << endl;
384 *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << " 0.3 " << endl;
385 *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << " 0.4 " << endl;
386 *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << " 0.5 " << endl;
387 *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << " 0.6 " << endl;
388 *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << " 0.7 " << endl;
389 *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << " 0.8 " << endl;
390 *fLog << endl;
391
392 *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
393 *fLog << " Acc Hadrons = " << Form("%5.1f", GetHadronAcceptance(0.5)*100) << "%" << endl;
394 *fLog << endl;
395}
396
397// --------------------------------------------------------------------------
398//
399// Draw clone of all histograms. (For the Meaning see class description)
400//
401TObject *MHHadronness::DrawClone(Option_t *opt) const
402{
403 if (fGraph->GetN()==0)
404 return NULL;
405
406 TCanvas &c = *MakeDefCanvas("Hadronness", fTitle);
407 c.Divide(2, 2);
408
409 gROOT->SetSelectedPad(NULL);
410
411 c.cd(1);
412 //gStyle->SetOptStat(10);
413 Getghness()->DrawCopy();
414 Getphness()->SetLineColor(kRed);
415 Getphness()->DrawCopy("same");
416
417 c.cd(2);
418 //gStyle->SetOptStat(0);
419 Getighness()->DrawCopy();
420 Getiphness()->SetLineColor(kRed);
421 Getiphness()->DrawCopy("same");
422
423 c.cd(3);
424 TGraph &g2 = (TGraph&)*fQfac->DrawClone("A*");
425 g2.SetBit(kCanDelete);
426 gPad->Modified();
427 gPad->Update();
428 if (g2.GetHistogram())
429 {
430 g2.GetXaxis()->SetRangeUser(0, 1);
431 g2.GetXaxis()->SetTitle("Acceptance Gammas");
432 g2.GetYaxis()->SetTitle("Quality");
433 g2.SetMarkerStyle(kFullDotSmall);
434 g2.Draw("P");
435
436 gPad->Modified();
437 gPad->Update();
438 }
439
440 c.cd(4);
441 gPad->Modified();
442 gPad->Update();
443 TGraph &g = (TGraph&)*fGraph->DrawClone("AC");
444 g.SetBit(kCanDelete);
445 gPad->Modified();
446 gPad->Update();
447 if (g.GetHistogram())
448 {
449 g.GetXaxis()->SetRangeUser(0, 1);
450 g.GetXaxis()->SetTitle("Acceptance Hadrons");
451 g.GetYaxis()->SetTitle("Acceptance Gammas");
452 g.SetMarkerStyle(kFullDotSmall);
453 g.Draw("P");
454
455 gPad->Modified();
456 gPad->Update();
457 }
458 /*
459 const Int_t h = fMinDist->GetMaximumBin();
460 TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
461 fIntGhness->GetBinContent(h), kStar);
462 m->SetMarkerColor(kBlue);
463 m->SetBit(kCanDelete);
464 m->Draw();
465 */
466
467 gStyle->SetOptStat(111100);
468
469 return &c;
470}
471
472// --------------------------------------------------------------------------
473//
474// Draw all histograms. (For the Meaning see class description)
475//
476void MHHadronness::Draw(Option_t *)
477{
478 if (fGraph->GetN()==0)
479 return;
480
481 if (!gPad)
482 MakeDefCanvas("Hadronness", fTitle);
483
484 gPad->Divide(2, 2);
485
486 gPad->cd(1);
487 //gStyle->SetOptStat(10);
488 Getghness()->Draw();
489 Getphness()->SetLineColor(kRed);
490 Getphness()->Draw("same");
491
492 gPad->cd(2);
493 //gStyle->SetOptStat(0);
494 Getighness()->Draw();
495 Getiphness()->SetLineColor(kRed);
496 Getiphness()->Draw("same");
497
498 gPad->cd(3);
499 fQfac->Draw("A*");
500 gPad->Modified();
501 gPad->Update();
502 if (fQfac->GetHistogram())
503 {
504 fQfac->GetXaxis()->SetRangeUser(0, 1);
505 fQfac->GetXaxis()->SetTitle("Acceptance Gammas");
506 fQfac->GetYaxis()->SetTitle("Quality");
507 fQfac->SetMarkerStyle(kFullDotSmall);
508 fQfac->Draw("P");
509
510 gPad->Modified();
511 gPad->Update();
512 }
513 gPad->cd(4);
514 gPad->Modified();
515 gPad->Update();
516 fGraph->Draw("AC");
517 gPad->Modified();
518 gPad->Update();
519 if (fGraph->GetHistogram())
520 {
521 fGraph->GetXaxis()->SetRangeUser(0, 1);
522 fGraph->GetXaxis()->SetTitle("Acceptance Hadrons");
523 fGraph->GetYaxis()->SetTitle("Acceptance Gammas");
524 fGraph->SetMarkerStyle(kFullDotSmall);
525 fGraph->Draw("P");
526 gPad->Modified();
527 gPad->Update();
528 }
529}
Note: See TracBrowser for help on using the repository browser.