source: trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc@ 2409

Last change on this file since 2409 was 2326, checked in by moralejo, 22 years ago
*** empty log message ***
File size: 55.6 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 expressed
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Robert Wagner 9/2002 <mailto:rw@rwagner.de>
19! Author(s): Robert Wagner 3/2003 <mailto:rw@rwagner.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27// //
28// MHOnSubtraction //
29// //
30// //
31// extracts the gamma signal from a pure ON-signal given in an //
32// ALPHA-Energy-Theta histogram. The class will take this histogram from //
33// the parameter list and will provide result histograms in the //
34// parameter list. //
35// //
36// This class still is work in progress. //
37// //
38// //
39//////////////////////////////////////////////////////////////////////////////
40
41// This part of MARS is code still evolving. Please do not change the code
42// without prior feedback by the author.
43
44#include "MHOnSubtraction.h"
45
46#include <TCanvas.h>
47#include <TF1.h>
48#include <TPaveText.h>
49#include <TLegend.h>
50#include <TStyle.h>
51
52#include "MBinning.h"
53#include "MParList.h"
54#include "MHArray.h"
55#include "MH3.h"
56
57#include "MHAlphaEnergyTheta.h"
58#include "MHAlphaEnergyTime.h"
59
60#include "MLog.h"
61#include "MLogManip.h"
62
63ClassImp(MHOnSubtraction);
64
65using namespace std;
66
67// --------------------------------------------------------------------------
68//
69// Default Constructor.
70//
71MHOnSubtraction::MHOnSubtraction(const char *name, const char *title) : fMaxSignif(0), fMaxRedChiSq(0), fSlope(20.0)
72{
73 fName = name ? name : "MHOnSubtraction";
74 fTitle = title ? title : "Extracts Gamma signal from pure ON data";
75 fHistogramType = "Theta";
76 fChiSquareHisto=0;
77 fSignificanceHisto=0;
78 fSummedAlphaPlots=0;
79 fThetaBin=0;
80
81}
82
83
84// --------------------------------------------------------------------------
85//
86// Destructor.
87//
88MHOnSubtraction::~MHOnSubtraction()
89{
90 if (fChiSquareHisto) delete fChiSquareHisto;
91 if (fSignificanceHisto) delete fSignificanceHisto;
92 if (fSummedAlphaPlots) delete fSummedAlphaPlots;
93}
94
95// -----------------------------------------------------------------------
96//
97// Calculate Significance according to Li and Ma, ApJ 272 (1983) 317, eq
98// (17). n_{On}, n_{Off} and the ratio of On-times for On and Off
99// measurements have to be provided.
100//
101// This function underestimates the true significance for it does not take
102// into account errors on the event numbers. A more accurate variation wil
103// be provided soon
104//
105Double_t MHOnSubtraction::CalcSignificance(Double_t nOn, Double_t nOff, Double_t theta)
106{
107 if (nOn<=0)
108 *fLog << warn << "Got " << nOn << " total events, " << flush;
109 if (nOff<=0)
110 *fLog << warn << "Got " << nOff << " background events, " << flush;
111 if (nOn<=0 || nOff<=0.0) {
112 *fLog << warn << "returning significance of 0.0" << endl;
113 return 0.0;
114 }
115 return sqrt(2*(nOn*log((1+theta)*nOn/(theta*(nOn+nOff)))
116 +nOff*log((1+theta)*(nOff/(nOn+nOff)))));
117}
118
119// -----------------------------------------------------------------------
120//
121// This function takes a first look at a given ALPHA distribution
122// and determines if a fit is applicable.
123//
124// Fits a given TH1 containing an ALPHA distribution with a combined
125// gaussian plus polynomial of third grade and returns the fitted
126// function. By signalRegionFactor the width of the region in which
127// signal extraction is to be performed can be specified in units of
128// sigma of the gaussian which is fitted to the distribution (default:
129// 3.5). Accordingly, FitHistogram returns the range in which it suggests
130// signal extraction (lowerBin < ALPHA < upperBin). An estimated
131// significance of the signal is also provided. draw specifies whether
132// FitHistogram should draw the distribution.
133//
134Bool_t MHOnSubtraction::FitHistogram(TH1 &alphaHisto, Double_t &sigLiMa,
135 Double_t &lowerBin, Double_t &upperBin,
136 Float_t signalRegionFactor, const Bool_t draw,
137 TString funcName)
138{
139 if (alphaHisto.GetEntries() == 0) {
140 *fLog << warn << "Histogram contains no entries. Returning. " << endl;
141 lowerBin=0;
142 upperBin=0;
143 sigLiMa=0;
144 if (draw) {
145 TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no events)");
146 lab->SetFillStyle(0);
147 lab->SetBorderSize(0);
148 lab->Draw();
149 }
150 return kFALSE;
151 }
152
153 //cosmetics
154 alphaHisto.GetXaxis()->SetTitle("ALPHA");
155 alphaHisto.GetYaxis()->SetTitle("events");
156
157 Double_t outerEvents = 0.0;
158 Double_t innerEvents = 0.0;
159 Int_t outerBins = 0;
160 Int_t innerBins = 0;
161 Int_t outerBinsZero = 0;
162 Int_t innerBinsZero = 0;
163
164 for (Int_t alphaBin = 1; alphaBin < alphaHisto.GetNbinsX(); alphaBin++) {
165 if ((alphaHisto.GetBinLowEdge(alphaBin) >= -35.)&& //inner Region
166 (alphaHisto.GetBinLowEdge(alphaBin+1) <= 35.0)) {
167 innerEvents+=alphaHisto.GetBinContent(alphaBin);
168 innerBins++;
169 if (alphaHisto.GetBinContent(alphaBin)==0.0)
170 innerBinsZero++;
171 } else {
172 if ((alphaHisto.GetBinLowEdge(alphaBin) >= -89.5)&& //outer Region
173 (alphaHisto.GetBinLowEdge(alphaBin+1) <=89.5)) {
174 outerEvents+=alphaHisto.GetBinContent(alphaBin);
175 outerBins++;
176 if (alphaHisto.GetBinContent(alphaBin)==0.0)
177 outerBinsZero++;
178 }
179 }
180 }
181
182 *fLog << dbg << "Plot contains " <<
183 outerEvents << " outer ev (" << outerEvents/outerBins << "/bin), " <<
184 innerEvents << " inner ev (" << innerEvents/innerBins << "/bin) " << endl;
185
186 if ((outerBinsZero!=0) || (innerBinsZero != 0))
187 *fLog << warn << "There are ";
188 if (outerBinsZero != 0)
189 *fLog << dbg << outerBinsZero << " empty outer bins ";
190 if (innerBinsZero != 0)
191 *fLog << dbg <<innerBinsZero << " empty inner bins ";
192 if ((outerBinsZero!=0) || (innerBinsZero != 0))
193 *fLog << endl;
194
195 if (outerEvents == 0.0) {
196 *fLog << warn << "No events in OFF region. Assuming background = 0"
197 << endl;
198 TF1 *po = new TF1("pol0"+funcName,"pol0(0)",-89.5,89.5);
199 po->SetLineWidth(2);
200 po->SetLineColor(2);
201 po->SetParNames("c");
202 po->SetParameter(0,0.0);
203 alphaHisto.GetListOfFunctions()->Add(po);
204 alphaHisto.SetLineColor(97);
205 if (draw) {
206 alphaHisto.DrawCopy(); //rwagner
207 po->Draw("SAME");
208 }
209 sigLiMa = 0;
210 lowerBin = 1;
211 upperBin = alphaHisto.GetNbinsX()-1;
212 signalRegionFactor = 0;
213
214 return kFALSE; //No gaus fit applied
215 }
216
217 if (outerEvents/outerBins < 2.5) {
218 *fLog << warn << "Only " << /*setprecision(2) <<*/ outerEvents/outerBins
219 << " events/bin in OFF region: "
220 << "Assuming this as background." << endl;
221
222 TF1 *po = new TF1("pol0"+funcName,"pol0(0)",-89.5,89.5);
223 po->SetLineWidth(2);
224 po->SetLineColor(94);
225 po->SetParNames("c");
226 po->SetParameter(0,outerEvents/outerBins);
227 alphaHisto.GetListOfFunctions()->Add(po);
228 if (draw) {
229 alphaHisto.DrawCopy(); //rwagner
230 po->Draw("SAME");
231 }
232
233 Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0);
234 Int_t maxBin = centerBin > alphaHisto.GetNbinsX() - centerBin ?
235 alphaHisto.GetNbinsX()- centerBin : centerBin;
236
237 Int_t lowerSignalEdge = centerBin;
238 for (Int_t i=3; i < maxBin; i++) {
239 if ((Double_t)alphaHisto.GetBinContent(centerBin-i)
240 < outerEvents/outerBins) {
241 lowerSignalEdge = centerBin-i+1;
242 break;
243 }
244 }
245 if (centerBin<lowerSignalEdge) lowerSignalEdge=centerBin;
246
247 Int_t upperSignalEdge = centerBin;
248 for (Int_t i=3; i < maxBin; i++) {
249 if ((Double_t)alphaHisto.GetBinContent(centerBin+i)
250 <= outerEvents/outerBins) {
251 upperSignalEdge=centerBin+i-1;
252 break;
253 }
254 }
255 if (centerBin>upperSignalEdge) upperSignalEdge=centerBin;
256
257 Double_t nOnInt = 0;
258 for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++)
259 nOnInt += alphaHisto.GetBinContent(i) - outerEvents/outerBins;
260
261 Double_t nOffInt = (upperSignalEdge - lowerSignalEdge + 1)
262 * (outerEvents/outerBins);
263
264 sigLiMa = CalcSignificance(nOnInt, nOffInt, 1);
265
266 if (sigLiMa < 3)
267 alphaHisto.SetLineColor(2);
268
269
270 *fLog << inf << "Estimated significance is " << sigLiMa
271 << " sigma " << endl;
272 *fLog << inf << "Signal region is "
273 << lowerBin << " < ALPHA < " << upperBin
274 << " (Most likely you wanna ignore this)" << endl;
275
276 return kFALSE; //No gaus fit applied
277 }
278
279 //fit combined gaus+pol3 to data
280 TF1 *gp = new TF1("gauspol3"+funcName,"gaus(0)+pol3(3)",-89.5,89.5);
281 gp->SetLineWidth(2);
282 gp->SetLineColor(2);
283 gp->SetParNames("Excess","A","sigma","a","b","c","d");
284 gp->SetParLimits(0, 200,2000);
285 gp->SetParLimits(1, -4.,4.);
286 gp->SetParLimits(2, 2.,20.);
287 // gp->SetParameter(6, 0.0000); // include for slope(0)=0 constrain
288 TString gpDrawOptions = draw ? "RIQ" : "RIQ0";
289 gpDrawOptions = "RIQ0"; // rwagner
290 alphaHisto.Fit("gauspol3"+funcName,gpDrawOptions);
291 alphaHisto.DrawCopy(); //rwagner
292 alphaHisto.SetDirectory(NULL); //rwagner
293
294 Double_t gausMean = gp->GetParameter(1);
295 Double_t gausSigma = gp->GetParameter(2);
296
297// TF1 *p3 = new TF1("p3"+funcName,"pol3(0)",-signalRegionFactor*gausSigma+gausMean,
298// signalRegionFactor*gausSigma+gausMean);
299 TF1 *p3 = new TF1("p3"+funcName,"pol3(0)",-100,
300 100);
301 p3->SetParameters(gp->GetParameter(3),gp->GetParameter(4),
302 gp->GetParameter(5),gp->GetParameter(6));
303 // p3->SetLineStyle(2);
304 p3->SetLineColor(4);
305 p3->SetLineWidth(1);
306 if (draw) p3->Draw("SAME");
307
308 TF1 *ga = new TF1("ga"+funcName,"gaus(0)",-40,40);
309 ga->SetParameters(gp->GetParameter(0),gp->GetParameter(1),gp->GetParameter(2));
310 ga->SetLineColor(93);
311 ga->SetLineWidth(2);
312 // if (draw) ga->Draw("SAME");
313
314 // identify the signal region: signalRegionFactor*gausSigma
315 // this allows us to
316 // (1) determine n_{ON}, n_{OFF}
317 Double_t scalingFactor =
318 (alphaHisto.GetXaxis()->GetBinLowEdge(alphaHisto.GetNbinsX()+1)
319 - alphaHisto.GetXaxis()->GetBinLowEdge(1)) / alphaHisto.GetNbinsX();
320
321 Double_t nOnInt = (gp->Integral(-signalRegionFactor*gausSigma+gausMean,
322 signalRegionFactor*gausSigma+gausMean) / scalingFactor);
323 Double_t nOffInt = (p3->Integral(-signalRegionFactor*gausSigma+gausMean,
324 signalRegionFactor*gausSigma+gausMean) / scalingFactor);
325
326 // (2) determine the signal region from fit in degrees
327 // we do it a bit more complicated: assuming that the binning in all
328 // histograms is the same, we want to be sure that summing up is always
329 // done over the same bins.
330
331 lowerBin = alphaHisto.GetXaxis()->FindBin(-signalRegionFactor*gausSigma+gausMean);
332 upperBin = alphaHisto.GetXaxis()->FindBin( signalRegionFactor*gausSigma+gausMean);
333
334 lowerBin = alphaHisto.GetBinLowEdge((Int_t)lowerBin);
335 upperBin = alphaHisto.GetBinLowEdge((Int_t)upperBin)+alphaHisto.GetBinWidth((Int_t)upperBin);
336
337 sigLiMa = CalcSignificance(nOnInt, nOffInt, 1);
338// if (sigLiMa < 3)
339// alphaHisto.SetLineColor(2);
340
341
342 *fLog << inf << "Fit estimates significance to be "
343 << sigLiMa << " sigma " << endl;
344
345 *fLog << inf << "Fit yields signal region to be "
346 << lowerBin << " < ALPHA < " << upperBin
347 << " (Chisquare/dof=" << gp->GetChisquare()/gp->GetNDF()
348 << ", prob=" << gp->GetProb() << ")" << endl;
349
350 return kTRUE; //returning gaus fit
351}
352
353
354
355
356
357// -----------------------------------------------------------------------
358//
359// Does the actual extraction of the gamma signal. For performance
360// reasons, fits already done by MHOnSubtraction::FitHistogram are used
361// and not redone.
362// From it, a polynomial function for the background is evaluated and the
363// gamma signal is extracted in the region given by lowerBin < ALPHA <
364// upperBin.
365// Significance of the signal is also provided. draw specifies whether
366// FitHistogram should draw the distribution.
367//
368Bool_t MHOnSubtraction::ExtractSignal(TH1 &alphaHisto, Double_t &sigLiMa,
369 Double_t &lowerBin, Double_t &upperBin,
370 Double_t &gammaSignal, Double_t &errorGammaSignal,
371 Double_t &off, Double_t &errorOff,
372 Float_t signalRegionFactor, const Bool_t draw, TString funcName,
373 TPad *drawPad, Int_t drawBase)
374{
375 TF1 *gausPol = alphaHisto.GetFunction("gauspol3"+funcName);
376 TF1 *pol = alphaHisto.GetFunction("pol0"+funcName);
377
378 if (!gausPol && !pol) {
379 *fLog << err << "Fatal: ALPHA histogram has no gauspol or pol "
380 << " fit attached to it." << endl;
381 TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no fit)");
382 lab->SetFillStyle(3000);
383 lab->SetBorderSize(0);
384 lab->DrawClone();
385 lab->SetBit(kCanDelete);
386 return kFALSE;
387 }
388
389 TF1* po;
390 po = NULL;
391
392 if (gausPol) {
393 Double_t gausMean = gausPol->GetParameter(1);
394 Double_t gausSigma = gausPol->GetParameter(2);
395 po = new TF1("po"+funcName,"pol3(0)",
396 -signalRegionFactor*gausSigma+gausMean,
397 signalRegionFactor*gausSigma+gausMean);
398 po->SetParameters(gausPol->GetParameter(3),gausPol->GetParameter(4),
399 gausPol->GetParameter(5),gausPol->GetParameter(6));
400
401 TF1 *ga = new TF1("ga"+funcName,"gaus(0)",-40,40);
402 ga->SetParameters(gausPol->GetParameter(0),gausPol->GetParameter(1),
403 gausPol->GetParameter(2));
404
405 if (draw) {
406 alphaHisto.Draw();
407 gausPol->Draw("SAME"); //Maybe not even necessary?
408
409 po->SetLineColor(4);
410 po->SetLineWidth(2);
411 po->Draw("SAME");
412
413 ga->SetLineColor(93);
414 ga->SetLineWidth(2);
415 ga->Draw("SAME");
416
417 char legendTitle[80];
418 sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f",
419 lowerBin, upperBin);
420
421 TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle);
422
423 legend->SetBorderSize(0);
424 legend->SetFillColor(10);
425 legend->SetFillStyle(0);
426 legend->AddEntry(gausPol, "combined N_{on}","l");
427 legend->AddEntry(po,"polynomial N_{bg} Signal region","l");
428 legend->AddEntry(ga, "putative gaussian N_{S}","l");
429 legend->Draw();
430 }
431 } // gausPol
432
433
434 if (pol) {
435 po = pol;
436
437 if (draw) {
438 alphaHisto.Draw();
439
440 po->SetLineColor(6);
441 po->SetLineWidth(2);
442 po->Draw("SAME");
443
444 char legendTitle[80];
445 sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f",
446 lowerBin, upperBin);
447
448 TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle);
449
450 legend->SetBorderSize(0);
451 legend->SetFillColor(10);
452 legend->SetFillStyle(0);
453 legend->AddEntry(po,"polynomial N_{bg} Signal region","l");
454 legend->Draw();
455 }
456 } // pol
457
458 Double_t nOn = 0;
459 Double_t eNOn = 0;
460
461 Int_t binNumberLow = alphaHisto.GetXaxis()->FindBin(lowerBin);
462 Int_t binNumberHigh = alphaHisto.GetXaxis()->FindBin(upperBin);
463
464 for (Int_t bin=binNumberLow; bin<binNumberHigh+1; bin++) {
465 nOn += alphaHisto.GetBinContent(bin);
466 eNOn += alphaHisto.GetBinError(bin) * alphaHisto.GetBinError(bin);
467 } //for bin
468 eNOn = sqrt(eNOn);
469
470 // Evaluate background
471
472 Double_t nOff = 0;
473 Double_t eNOff = 0;
474 for (Int_t bin=binNumberLow; bin<binNumberHigh+1; bin++) {
475 Double_t x = .5*(alphaHisto.GetBinLowEdge(bin)+alphaHisto.GetBinLowEdge(bin+1));
476 Double_t offEvts = po->Eval(x);
477 //cout << bin << ": " << offEvts << endl;
478 nOff += offEvts;
479 } //for bin
480 eNOff = sqrt(fabs(nOff));
481
482 if (nOn==0) // there should not be a negative number of signal events
483 nOff=0;
484
485 if (nOff<0) { // there should not be a negative number of off events
486 nOff=0;
487 eNOff=0;
488 }
489
490 *fLog << inf << "nEvts = " << nOn << "+-" << eNOn << ", ";
491
492 off = nOff; errorOff = eNOff;
493 gammaSignal = nOn-nOff; errorGammaSignal = sqrt(eNOn*eNOn + eNOff*eNOff);
494
495 *fLog << inf << "nBack = " << nOff << "+-" << eNOff << ", ";
496 *fLog << inf << "nSig = " << gammaSignal << "+-" << errorGammaSignal << endl;
497
498 sigLiMa = CalcSignificance(nOn, nOff, 1);
499 // Double_t sigFit=(ga->GetParameter(1))/(ga->GetParError(1)); //Mean / sigMean
500
501 *fLog << inf << "Significance: "<<sigLiMa<<" sigma "<<endl;
502
503 if (draw) {
504 TPaveText *pt = new TPaveText(0.11,.74,.57,.88,"NDC ");
505 char tx[60];
506 sprintf(tx, "Excess: %2.2f #pm %2.2f", nOn-nOff, sqrt(eNOn*eNOn + eNOff*eNOff));
507 pt->AddText(tx);
508 sprintf(tx, "Off: %2.2f #pm %2.2f", nOff, eNOff);
509 pt->AddText(tx);
510 sprintf(tx, "Significance: %2.2f #sigma", sigLiMa);
511 pt->AddText(tx);
512 pt->SetFillStyle(0);
513 pt->SetBorderSize(0);
514 pt->SetTextAlign(12);
515 pt->Draw("");
516 }
517 if (draw||drawPad) {
518
519 // Plot significance vs alpha_cut
520
521 Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0);
522 Int_t maxBin = centerBin > alphaHisto.GetNbinsX()- centerBin ?
523 alphaHisto.GetNbinsX()- centerBin : centerBin;
524
525 const Int_t totalBins = centerBin;
526 Float_t alpha[totalBins];
527 Float_t signi[totalBins];
528 Float_t maxSigni = 0;
529
530 for (Int_t i=0; i < maxBin; i++) {
531 Double_t nOn = 0; Double_t eNOn = 0;
532 Double_t nOff = 0; Double_t eNOff = 0;
533 for (Int_t bin=centerBin-i; bin<centerBin+i+1; bin++) {
534 nOn += alphaHisto.GetBinContent(bin);
535 eNOn += alphaHisto.GetBinError(bin) * alphaHisto.GetBinError(bin);
536 Double_t x = .5*(alphaHisto.GetBinLowEdge(bin)+alphaHisto.GetBinLowEdge(bin+1));
537 Double_t offEvts = po->Eval(x);
538 nOff += offEvts;
539 } //for bin
540 eNOn = sqrt(eNOn);
541 eNOff = sqrt(nOff);
542 alpha[i] =
543 (alphaHisto.GetXaxis()->GetBinLowEdge(centerBin+i+1)-
544 alphaHisto.GetXaxis()->GetBinLowEdge(centerBin-i))/2;
545 signi[i] = CalcSignificance(nOn, nOff, 1);
546 maxSigni = maxSigni > signi[i] ? maxSigni : signi[i];
547 }
548
549 if (!drawPad) {
550 TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600);
551 c3->cd();
552 }
553
554 if (drawPad)
555 drawPad->cd(drawBase-1);
556
557 TGraph* gr = new TGraph(totalBins-2,alpha,signi);
558 TString drawOpt = "L";
559
560 if (draw || (drawPad && fSigniPlotIndex == 0))
561 drawOpt += "A";
562
563 gr->Draw(drawOpt);
564 if (drawPad && fSigniPlotIndex == 0) {
565 gr->GetXaxis()->SetTitle("|#alpha_{max}|");
566 gr->GetYaxis()->SetTitle("Significance");
567 }
568 gr->SetMarkerStyle(2);
569 gr->SetMarkerSize(1);
570 gr->SetMarkerColor(4+fSigniPlotIndex);
571 gr->SetLineColor(4+fSigniPlotIndex);
572 gr->SetLineWidth(1);
573 gr->SetTitle("Significance vs ALPHA integration range");
574 gr->Draw("L");
575
576 fSigniPlotIndex++;
577
578 } //draw
579
580 return kTRUE;
581}
582
583// -----------------------------------------------------------------------
584//
585// Extraction of Gamma signal is performed on a MHAlphaEnergyTheta or
586// MHAlphaEnergyTime object expected in the ParList.
587//
588Bool_t MHOnSubtraction::Calc(MParList *parlist, const Bool_t Draw)
589{
590 //Find source histograms
591 fHistogramType = "Theta";
592 MHAlphaEnergyTheta* mhisttheta = (MHAlphaEnergyTheta*)parlist->FindObject("MHAlphaEnergyTheta");
593 if (mhisttheta) return Calc((TH3D*)mhisttheta->GetHist(), parlist, Draw);
594
595 fHistogramType = "Time";
596 MHAlphaEnergyTime* mhisttime = (MHAlphaEnergyTime*)parlist->FindObject("MHAlphaEnergyTime");
597 if (mhisttime) return Calc((TH3D*)mhisttime->GetHist(), parlist, Draw);
598
599 *fLog << err << "No MHAlphaEnergyTheta type object found in the parameter list. Aborting." << endl;
600 return kFALSE;
601}
602
603// -----------------------------------------------------------------------
604//
605// Extraction of Gamma signal is performed on a TH3D histogram in
606// ALPHA, Energy and Theta.
607//
608Bool_t MHOnSubtraction::CalcAET(TH3D *aetHisto, MParList *parlist, const Bool_t Draw)
609{
610 // Analyze aetHisto
611 // -------------------------------------
612 Int_t alphaBins = aetHisto->GetNbinsX(); // # of alpha bins
613 Int_t energyBins = aetHisto->GetNbinsY();
614 Int_t thetaBins = aetHisto->GetNbinsZ();
615
616 // We output an array of histograms to the parlist.
617 // Create a template for such a histogram, needed by MHArray
618 // -------------------------------------
619 MBinning *binsfft = new MBinning("BinningMH3X");
620 binsfft->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1),
621 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
622 parlist->AddToList(binsfft);
623 MH3 *energyHistogram = new MH3(1); // The template histogram in energy
624 // for a given theta value
625 energyHistogram->SetName("MHOnSubtractionEnergyHistogramTemplate");
626 parlist->AddToList(energyHistogram);
627
628 fThetaHistoArray = new MHArray("MHOnSubtractionEnergyHistogramTemplate", kTRUE,
629 "MHOnSubtractionGammaSignalArray",
630 "Array of histograms in energy bins");
631 fThetaHistoArray->SetupFill(parlist);
632 parlist->AddToList(fThetaHistoArray);
633
634 // Canvases---direct debug output for the time being
635 // -------------------------------------
636 TCanvas *c3 = new TCanvas("c3", "Plots by MHOnSubtraction::ExtractSignal", 800,600);
637 cout << thetaBins << " x " << energyBins << endl;
638 c3->Divide(thetaBins,energyBins);
639
640 TCanvas *c4a = new TCanvas("c4a", "Energy distributions for different ZA", 800,600);
641
642 TH1D* histalphaon[energyBins*thetaBins]; // ALPHA histograms
643
644 fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
645 fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 41, -0.5, 40.5);
646 fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha",
647 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
648 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
649
650 // -------------------------------------
651
652 fThetaLegend = new TLegend(0.83, 0.07, 0.98, 0.42, "Energy Distributions");
653 fThetaLegend->SetBorderSize(1);
654
655 Double_t overallSigLiMa = 0;
656
657 for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) {
658
659 char hname[80];
660 sprintf(hname, "Energy distribution for %s bin %d", fHistogramType.Data(),
661 thetaBin);
662 *fLog << all << "Calculating " << hname << endl;
663
664 Double_t minLowerBin = -10; // minimum integration range
665 Double_t maxUpperBin = 10; // minimum integration range
666 Double_t maxAlpha = 70; // maximum integration range
667
668 Double_t sumSigLiMa = 0;
669
670 // This loop just fixes the integration range
671 // And alerts when no significant excess is found.
672
673 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
674
675 sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin);
676 histalphaon[(thetaBin-1)*(energyBins)+energyBin] =
677 new TH1D(hname, "ALPHA histogram",
678 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
679 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
680 histalphaon[(thetaBin-1)*(energyBins)+energyBin]->Sumw2();
681 histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetTitle(hname);
682
683 sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d",
684 energyBin, thetaBin);
685 *fLog << inf << hname << endl;
686
687 // now we fill one histogram for one particular set
688 // of Energy and Theta...
689
690 if (aetHisto->InheritsFrom("TH3D")||aetHisto->InheritsFrom("TH2D")) {
691 aetHisto->GetYaxis()->SetRange(energyBin,energyBin); // E
692 if (aetHisto->InheritsFrom("TH3D"))
693 aetHisto->GetZaxis()->SetRange(thetaBin,thetaBin); // theta
694 histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto->Project3D("xe");
695 } else {
696 histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto;
697 }
698
699 *fLog << inf << "This histogram has "
700 << histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetEntries()
701 << " entries" << endl;
702
703 sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin);
704 histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetName(hname);
705
706 // Test: Find signal region and significance in histogram
707 Double_t lowerBin, upperBin, sSigLiMa;
708 FitHistogram(*histalphaon[(thetaBin-1)*(energyBins)+energyBin],
709 sSigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE);
710
711 Double_t redChiSq = histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetChisquare()/histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetNDF();
712
713 fChiSquareHisto->Fill(redChiSq);
714 fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq;
715
716 // histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetProb() <<endl;
717
718 fSummedAlphaPlots->Add(histalphaon[(thetaBin-1)*(energyBins)+energyBin], 1);
719
720 if (sSigLiMa < 3)
721 *fLog << inf << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
722 sumSigLiMa+=sSigLiMa;
723
724 // Evaluate integration range
725
726 lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
727 minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;
728
729 upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
730 maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
731
732 } //energyBin
733
734 *fLog << inf << "Integration range for this " << fHistogramType
735 << " value will be " << minLowerBin << " < ALPHA < "
736 << maxUpperBin << endl;
737
738 *fLog << inf << "Summed estd. significance for this " << fHistogramType
739 << " value is " << sumSigLiMa << endl;
740
741 // Create Histogram in Energy for this Theta value
742 cout << "STARTIGN REAL CALC1 "<< endl;
743 fThetaHistoArray->AddHistogram();
744 cout << "STARTIGN REAL CALC1a "<< endl;
745 fThetaHistoArray->SetIndex(thetaBin-1);
746
747 MH3 &mThetaHisto = *(MH3*)(fThetaHistoArray->GetH());
748 mThetaHisto.Print();
749
750 TH1D &thetaHisto = (TH1D&)(mThetaHisto.GetHist());
751 sprintf(hname,"Energy distribution for theta bin %d", thetaBin);
752 thetaHisto.SetTitle(hname);
753 thetaHisto.SetEntries(0);
754
755 // we have a rough idea of what is going on in the ALPHA plot...
756 // So now we can indeed extract the signal.
757
758// cout << "STARTIGN REAL CALC "<< endl;
759
760 sumSigLiMa = 0;
761 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
762
763 sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d",
764 energyBin, thetaBin);
765 *fLog << inf << hname << endl;
766
767 Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;
768
769 c3->cd((thetaBin-1)*(energyBins)+energyBin);
770
771 ExtractSignal(*histalphaon[(thetaBin-1)*(energyBins)+energyBin],
772 sigLiMa, minLowerBin, maxUpperBin,
773 gammaSignal, errorGammaSignal,
774 off, errorOff, (Float_t)3, kTRUE);
775
776 fSignificanceHisto->Fill(sigLiMa);
777 fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;
778
779 thetaHisto.SetBinContent(energyBin, gammaSignal*TMath::Exp(7-energyBin));
780 thetaHisto.SetBinError(energyBin, errorGammaSignal);
781 thetaHisto.SetEntries(thetaHisto.GetEntries()+gammaSignal);
782
783 sumSigLiMa += sigLiMa;
784
785 }//energyBin
786
787 *fLog << inf << "Summed significance for this " << fHistogramType
788 << " value is " << sumSigLiMa << endl;
789
790 //fit exponential function to data
791 TF1* e = new TF1("expF","expo",0,5);
792 e->SetLineWidth(3);
793 e->SetLineColor(thetaBin);
794 // e->SetParLimits(1, -0.5,3); //limits on slope
795
796 if (fSlope!=20.0) {
797 e->SetParameter(1, fSlope);
798 *fLog << inf << "Fixing slope to " << e->GetParameter(1) << endl;
799 }
800
801 TString eDrawOptions = Draw ? "RIQ" : "RIQ0";
802 cout << "doing the fit" << endl;
803 thetaHisto.Fit("expF");
804
805 Double_t expoSlope = e->GetParameter(1);
806
807 *fLog << inf << "Determined slope for energy distribution is "
808 << expoSlope << endl;
809
810 Double_t integralEvts =
811 e->Integral(aetHisto->GetYaxis()->GetBinLowEdge(1),
812 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
813
814 *fLog << inf << "Integral in E range [" <<
815 aetHisto->GetYaxis()->GetBinLowEdge(1) << ";" <<
816 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1) << "] is " <<
817 integralEvts << endl;
818
819 // Plot Energy histogram
820
821 c4a->cd(0);
822 thetaHisto.SetLineColor(thetaBin);
823 if (thetaBin==1) {
824 thetaHisto.Draw();
825 } else {
826 thetaHisto.Draw("SAME");
827 }
828
829 sprintf(hname,"Theta bin %d",thetaBin);
830 fThetaLegend->AddEntry(&thetaHisto, hname, "l");
831
832 overallSigLiMa += sumSigLiMa;
833
834 } //thetaBin
835
836 fThetaLegend->Draw();
837
838 Double_t sigLiMa, lowerBin, upperBin;
839 FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
840 fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");
841
842 *fLog << inf << "Summed overall significance is " << overallSigLiMa << endl;
843
844 *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
845 // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
846
847 // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
848 // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
849
850 return kTRUE;
851}
852
853
854
855
856
857
858
859
860// -----------------------------------------------------------------------
861//
862// Extraction of Gamma signal is performed on a TH3D histogram in
863// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms
864// in energy and Theta.
865//
866Bool_t MHOnSubtraction::Calc(TH3 *aetHisto, MParList *parlist, const Bool_t Draw)
867{
868 // Analyze aetHisto
869 // -------------------------------------
870 Int_t energyBins = aetHisto->GetNbinsY();
871 Int_t thetaBins = aetHisto->GetNbinsZ();
872
873 *fLog << "Total events: " << aetHisto->GetEntries() << endl;
874 *fLog << "Total energy bins: " << energyBins << endl;
875 *fLog << "Total theta bins: " << thetaBins << endl;
876
877 // We output results in an array of histograms to the parameter list.
878 // Create a template for such a histogram, needed by MH3
879 // -------------------------------------
880 MBinning *binsmh3x = new MBinning("BinningMH3X");
881 // dummy binning, real one follows some lines down
882 binsmh3x->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1),
883 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
884 parlist->AddToList(binsmh3x);
885
886 MBinning *binsmh3y = new MBinning("BinningMH3Y");
887 // dummy binning, real one follows some lines down
888 binsmh3y->SetEdges(thetaBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
889 aetHisto->GetZaxis()->GetBinLowEdge(thetaBins+1));
890 parlist->AddToList(binsmh3y);
891
892 MH3 *signalHisto = new MH3(2); // Signal(E,t)
893 signalHisto->SetupFill(parlist);
894 parlist->AddToList(signalHisto);
895 signalHisto->SetName("MHOnSubtractionSignal");
896 signalHisto->SetTitle("Gamma Events");
897
898 MH3 *offHisto = new MH3(2); // Off(E,t)
899 offHisto->SetupFill(parlist);
900 parlist->AddToList(offHisto);
901 offHisto->SetName("MHOnSubtractionOff");
902 offHisto->SetTitle("Background Events");
903
904 MH3 *signifHisto = new MH3(2); // Significance(E,t)
905 signifHisto->SetupFill(parlist);
906 parlist->AddToList(signifHisto);
907 signifHisto->SetName("MHOnSubtractionSignif");
908 signifHisto->SetTitle("Significance");
909
910// if (!fChiSquareHisto)
911// fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
912// if (!fSignificanceHisto)
913// fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 40.5, -0.5, 41);
914// if (!fSummedAlphaPlots)
915// fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha",
916// alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
917// aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
918
919 TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600);
920 c4->Divide(energyBins+3,thetaBins,0,0,0);
921
922 TH2D& signalTH2DHist = (TH2D&)signalHisto->GetHist();
923 TH2D& offTH2DHist = (TH2D&)offHisto->GetHist();
924 TH2D& signifTH2DHist = (TH2D&)signifHisto->GetHist();
925
926 signalTH2DHist.Reset();
927 signalTH2DHist.SetTitle("Gamma Signal");
928 signalTH2DHist.SetXTitle("E [GeV]");
929 signalTH2DHist.SetYTitle("theta [deg]");
930 signalHisto->SetBinning(&signalTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
931 signalTH2DHist.Sumw2();
932
933 offTH2DHist.Reset();
934 offTH2DHist.SetTitle("Off Signal");
935 offTH2DHist.SetXTitle("E [GeV]");
936 offTH2DHist.SetYTitle("theta [deg]");
937 offHisto->SetBinning(&offTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
938 offTH2DHist.Sumw2();
939
940 signifTH2DHist.Reset();
941 signifTH2DHist.SetTitle("Significance");
942 signifTH2DHist.SetXTitle("E [GeV]");
943 signifTH2DHist.SetYTitle("theta [deg]");
944 signifHisto->SetBinning(&signifTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
945 signifTH2DHist.Sumw2();
946
947 for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) {
948
949 *fLog << dbg << "--------------------------------------" << endl;
950 *fLog << dbg << "Processing ALPHA plots for theta bin " << thetaBin << endl;
951 *fLog << dbg << "--------------------------------------" << endl;
952 aetHisto->GetZaxis()->SetRange(thetaBin, thetaBin);
953 TH2F* aeHisto = (TH2F*)aetHisto->Project3D("yxe");
954 aeHisto->SetDirectory(NULL);
955 char hname[80];
956 sprintf(hname, "%2.1f < #theta < %2.1f",
957 aetHisto->GetZaxis()->GetBinLowEdge(thetaBin),
958 aetHisto->GetZaxis()->GetBinLowEdge(thetaBin+1));
959 *fLog << inf << "There are " << aeHisto->GetEntries()
960 << " entries in " << hname << endl;
961 aeHisto->SetTitle(hname);
962 sprintf(hname, "MHOnSubtractionThetaBin%d", thetaBin);
963 aeHisto->SetName(hname);
964
965 c4->cd((energyBins+3)*(thetaBin-1)+1);
966 aeHisto->SetMarkerColor(3);
967 aeHisto->DrawCopy();
968
969 c4->Update();
970
971 // We hand over a nonstandard parameter list, which
972 // contails lower-dimensional result histograms
973 // signalHisto, offHisto and signifHisto
974
975// cout << fLog->GetDebugLevel() << endl;
976 fLog->SetDebugLevel(1);
977
978 MParList *parlistth2 = new MParList();
979 // parlistth2->SetNullOutput();
980 parlistth2->AddToList(binsmh3x);
981
982 MH3* signalHisto2 = new MH3(1); // Signal(E)
983 signalHisto2->SetupFill(parlistth2);
984 parlistth2->AddToList(signalHisto2);
985 signalHisto2->SetName("MHOnSubtractionSignal");
986 signalHisto2->SetTitle("Gamma Events");
987
988 MH3* offHisto2 = new MH3(1); // Off(E)
989 offHisto2->SetupFill(parlistth2);
990 parlistth2->AddToList(offHisto2);
991 offHisto2->SetName("MHOnSubtractionOff");
992 offHisto2->SetTitle("Background Events");
993
994 MH3* signifHisto2 = new MH3(1); // Significance(E)
995 signifHisto2->SetupFill(parlistth2);
996 parlistth2->AddToList(signifHisto2);
997 signifHisto2->SetName("MHOnSubtractionSignif");
998 signifHisto2->SetTitle("Significance");
999
1000 fLog->SetDebugLevel(-1);
1001
1002 TH1D& signalTH1DHist = (TH1D&)signalHisto2->GetHist();
1003 TH1D& offTH1DHist = (TH1D&)offHisto2->GetHist();
1004 TH1D& signifTH1DHist = (TH1D&)signifHisto2->GetHist();
1005
1006 signalTH1DHist.Sumw2();
1007 offTH1DHist.Sumw2();
1008 signifTH1DHist.Sumw2();
1009
1010 TH2Calc(aeHisto, parlistth2, kFALSE, c4,
1011 (energyBins+3)*(thetaBin-1)+4);
1012
1013
1014 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
1015
1016 // cout << "filling " << thetaBin << " " << energyBin << ": "
1017// << signalTH1DHist.GetBinContent(energyBin) << "+-"
1018// << signalTH1DHist.GetBinError(energyBin) << " "
1019// << offTH1DHist.GetBinContent(energyBin) << "+-"
1020// << offTH1DHist.GetBinError(energyBin) << endl;
1021
1022 if (signalTH1DHist.GetBinContent(energyBin)>=0.0) {
1023 signalTH2DHist.SetBinContent(energyBin, thetaBin,
1024 signalTH1DHist.GetBinContent(energyBin));
1025 signalTH2DHist.SetBinError(energyBin, thetaBin,
1026 signalTH1DHist.GetBinError(energyBin));
1027 }
1028
1029 if (offTH1DHist.GetBinContent(energyBin)>=0.0) {
1030 offTH2DHist.SetBinContent(energyBin, thetaBin,
1031 offTH1DHist.GetBinContent(energyBin));
1032 offTH2DHist.SetBinError(energyBin, thetaBin,
1033 offTH1DHist.GetBinError(energyBin));
1034 }
1035
1036 signifTH2DHist.SetBinContent(energyBin, thetaBin,
1037 signifTH1DHist.GetBinContent(energyBin));
1038 signifTH2DHist.SetBinError(energyBin, thetaBin,
1039 signifTH1DHist.GetBinError(energyBin));
1040 } //energyBin
1041
1042
1043 c4->cd((energyBins+3)*(thetaBin-1)+2);
1044 sprintf(hname,"c4_%d",(energyBins+3)*(thetaBin-1)+2);
1045 TPad* tp = (TPad*)(gROOT->FindObject(hname));
1046 tp->SetLogx();
1047 signalTH1DHist.SetLineColor(2);
1048 signalTH1DHist.DrawCopy();
1049 offTH1DHist.SetLineColor(4);
1050 offTH1DHist.DrawCopy("SAME");
1051 c4->Update();
1052
1053 signalTH2DHist.SetEntries(signalTH2DHist.GetEntries()+signalTH1DHist.GetEntries());
1054 offTH2DHist.SetEntries(offTH2DHist.GetEntries()+offTH1DHist.GetEntries());
1055 signifTH2DHist.SetEntries(signifTH2DHist.GetEntries()+signifTH1DHist.GetEntries());
1056
1057 delete signifHisto2;
1058 delete offHisto2;
1059 delete signalHisto2;
1060 delete parlistth2;
1061
1062 }
1063
1064
1065 c4->Update();
1066
1067 return kTRUE;
1068}
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120/*
1121
1122// -----------------------------------------------------------------------
1123//
1124// Extraction of Gamma signal is performed on a TH3D histogram in
1125// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms
1126// in energy and Theta.
1127//
1128Bool_t MHOnSubtraction::CalcLightCurve(TH3 *aetHisto, MParList *parlist, const Bool_t draw)
1129{
1130 // Analyze aetHisto
1131 // -------------------------------------
1132 Int_t alphaBins = aetHisto->GetNbinsX();
1133 Int_t energyBins = aetHisto->GetNbinsY();
1134 Int_t timeBins = aetHisto->GetNbinsZ();
1135
1136 *fLog << "Total events: " << aetHisto->GetEntries() << endl;
1137 *fLog << "Total energy bins: " << energyBins << endl;
1138 *fLog << "Total time bins: " << timeBins << endl;
1139
1140 // We output results in an array of histograms to the parameter list.
1141 // Create a template for such a histogram, needed by MH3
1142 // -------------------------------------
1143 MBinning *binsmh3x = new MBinning("BinningMH3X");
1144 // dummy binning, real one follows some lines down
1145 binsmh3x->SetEdges(timeBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
1146 aetHisto->GetZaxis()->GetBinLowEdge(timeBins+1));
1147 parlist->AddToList(binsmh3x);
1148
1149 MH3 *signalHisto = new MH3(1); // Signal(t)
1150 signalHisto->SetupFill(parlist);
1151 parlist->AddToList(signalHisto);
1152 signalHisto->SetName("MHOnSubtractionSignal");
1153 signalHisto->SetTitle("Gamma Events");
1154
1155 MH3 *offHisto = new MH3(1); // Off(t)
1156 offHisto->SetupFill(parlist);
1157 parlist->AddToList(offHisto);
1158 offHisto->SetName("MHOnSubtractionOff");
1159 offHisto->SetTitle("Background Events");
1160
1161 MH3 *signifHisto = new MH3(1); // Significance(t)
1162 signifHisto->SetupFill(parlist);
1163 parlist->AddToList(signifHisto);
1164 signifHisto->SetName("MHOnSubtractionSignif");
1165 signifHisto->SetTitle("Significance");
1166
1167 TH1D& signalTH1DHist = (TH1D&)signalHisto->GetHist();
1168 TH1D& offTH1DHist = (TH1D&)offHisto->GetHist();
1169 TH1D& signifTH1DHist = (TH1D&)signifHisto->GetHist();
1170
1171 signalTH1DHist.Reset();
1172 signalTH1DHist.SetTitle("Gamma Signal");
1173 signalTH1DHist.SetXTitle("Time [MJD]");
1174 signalTH1DHist.SetYTitle("Events");
1175 signalHisto->SetBinning(&signalTH1DHist, aetHisto->GetZaxis());
1176 signalTH1DHist.Sumw2();
1177
1178 offTH1DHist.Reset();
1179 offTH1DHist.SetTitle("Background Signal");
1180 offTH1DHist.SetXTitle("Time [MJD]");
1181 offTH1DHist.SetYTitle("Events");
1182 offHisto->SetBinning(&offTH1DHist, aetHisto->GetZaxis());
1183 offTH1DHist.Sumw2();
1184
1185 signifTH1DHist.Reset();
1186 signifTH1DHist.SetTitle("Significance");
1187 signifTH1DHist.SetXTitle("Time [MJD]");
1188 signifTH1DHist.SetYTitle("Significance #sigma");
1189 signifHisto->SetBinning(&signifTH1DHist, aetHisto->GetZaxis());
1190 signifTH1DHist.Sumw2();
1191
1192 //The following histogram is an additional histogram needed for
1193 //the light curve
1194
1195 TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600);
1196 c4->Divide(8,7,0,0,0);
1197
1198
1199
1200 // The following loop fixes a common integration region for each
1201 // energy bin and alerts when no significant excess is found.
1202
1203 Double_t minLowerBin = -10; // minimum integration range
1204 Double_t maxUpperBin = 10; // minimum integration range
1205 Double_t maxAlpha = 70; // maximum integration range
1206 Double_t sumSigLiMa = 0;
1207
1208 TH1F* alphaHisto[timeBins];
1209
1210 for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {
1211
1212 *fLog << dbg << "--------------------------------------" << endl;
1213 *fLog << dbg << "Processing ALPHA plots for time bin " << timeBin << endl;
1214 *fLog << dbg << "--------------------------------------" << endl;
1215
1216 aetHisto->GetZaxis()->SetRange(timeBin, timeBin);
1217
1218 char hname[80];
1219 sprintf(hname, "MHOnSubtractionTimeBin%d", timeBin);
1220 alphaHisto[timeBin-1] =
1221 new TH1F(hname, "ALPHA histogram",
1222 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
1223 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
1224 alphaHisto[timeBin-1]->SetDirectory(NULL);
1225 alphaHisto[timeBin-1]->Sumw2();
1226 alphaHisto[timeBin-1] = (TH1F*)aetHisto->Project3D("xe");
1227 alphaHisto[timeBin-1]->SetName(hname);
1228 alphaHisto[timeBin-1]->SetDirectory(NULL);
1229
1230 sprintf(hname, "%6.0f < t < %6.0f",
1231 aetHisto->GetZaxis()->GetBinLowEdge(timeBin),
1232 aetHisto->GetZaxis()->GetBinLowEdge(timeBin+1));
1233 *fLog << inf << "There are " << alphaHisto[timeBin-1]->GetEntries()
1234 << " entries in " << hname << endl;
1235 alphaHisto[timeBin-1]->SetTitle(hname);
1236
1237 // Find signal region and significance in histogram
1238 Double_t lowerBin, upperBin, sSigLiMa;
1239 char funcName[20];
1240 sprintf(funcName, "%2d", timeBin);
1241
1242 Bool_t drawFit = kTRUE;
1243
1244 c4->cd(timeBin);
1245// alphaHisto[timeBin-1]->SetMarkerColor(3);
1246 alphaHisto[timeBin-1]->DrawCopy();
1247
1248 c4->Update();
1249
1250 fSigniPlotColor = 0;
1251 ;
1252 Bool_t fit = FitHistogram(*alphaHisto[timeBin-1], sSigLiMa,
1253 lowerBin, upperBin, (Float_t)3.0, drawFit,
1254 funcName);
1255
1256 if (fit) {
1257 if (sSigLiMa < 3)
1258 *fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
1259 sumSigLiMa+=sSigLiMa;
1260
1261 // Evaluate integration range
1262 lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
1263 minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;
1264 upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
1265 maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
1266 }
1267
1268
1269
1270
1271
1272
1273 } //timeBin
1274
1275 *fLog << inf << "=> Integration range will be " << minLowerBin << " < ALPHA < "
1276 << maxUpperBin << "," << endl;
1277 *fLog << inf << " Summed estimated significance is " << sumSigLiMa << endl;
1278
1279 // we have an idea of what is going on in the ALPHA plot...
1280 // So now we can indeed extract the signal.
1281
1282 sumSigLiMa = 0;
1283 for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {
1284 *fLog << inf << "Processing ALPHA distribution for time bin " << timeBin << endl;
1285 if (alphaHisto[timeBin-1]->GetEntries() == 0) {
1286 *fLog << warn << "No attempt to extract a signal from 0 events." << endl;
1287 if (draw) {
1288 TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-");
1289 lab->SetBorderSize(0);
1290 lab->Draw();
1291 }
1292 } else {
1293 char funcName[20];
1294 sprintf(funcName, "%2d", timeBin);
1295
1296 Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;
1297
1298 Bool_t drawFit = kFALSE;
1299
1300 ExtractSignal(*alphaHisto[timeBin-1],
1301 sigLiMa, minLowerBin, maxUpperBin,
1302 gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
1303 funcName, NULL, 1);
1304
1305 sumSigLiMa += sigLiMa;
1306
1307 fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;
1308
1309 // Fill Signal
1310 TH1D &h = (TH1D&)(signalHisto->GetHist());
1311 h.SetBinContent(timeBin, gammaSignal);
1312 h.SetBinError(timeBin, errorGammaSignal);
1313 h.SetEntries(h.GetEntries()+gammaSignal);
1314
1315 // Fill Off Events
1316 TH1D &ho = (TH1D&)(offHisto->GetHist());
1317 ho.SetBinContent(timeBin, off);
1318 ho.SetBinError(timeBin, errorOff);
1319 ho.SetEntries(ho.GetEntries()+off);
1320
1321 // Fill Significance
1322 TH1D &hg = (TH1D&)(signifHisto->GetHist());
1323 hg.SetBinContent(timeBin, sigLiMa);
1324 hg.SetEntries(hg.GetEntries()+sigLiMa);
1325 }
1326
1327 }//timeBin
1328
1329 *fLog << inf << "Summed significance is " << sumSigLiMa << endl;
1330
1331
1332 if (draw) {
1333 Double_t sigLiMa, lowerBin, upperBin;
1334 FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
1335 fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");
1336
1337 *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
1338 // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
1339
1340 // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
1341 // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
1342
1343 fChiSquareHisto->DrawClone();
1344 fSignificanceHisto->DrawClone();
1345 fSummedAlphaPlots->DrawClone();
1346 }
1347
1348 return kTRUE;
1349}
1350
1351
1352
1353
1354
1355*/
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432// -----------------------------------------------------------------------
1433//
1434// Extraction of Gamma signal is performed on a TH2 histogram in
1435// ALPHA and Energy. Output to parameter list is TH1 histogram in
1436// energy
1437//
1438Bool_t MHOnSubtraction::TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t draw, TPad *drawPad, Int_t drawBase)
1439{
1440
1441 fSigniPlotColor = 0;
1442
1443 // Analyze aeHisto
1444 // -------------------------------------
1445 const Int_t alphaBins = aeHisto->GetNbinsX(); // # of alpha bins
1446 const Int_t energyBins = aeHisto->GetNbinsY();
1447
1448 // Check availability of result histograms
1449 // -------------------------------------
1450
1451 MH3* signalHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignal","MH3");
1452 MH3* offHisto = (MH3*)parlist->FindObject("MHOnSubtractionOff","MH3");
1453 MH3* signifHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignif","MH3");
1454
1455 if (signalHisto && offHisto && signifHisto) {
1456// *fLog << dbg << "Result histograms are present in parameter list " <<
1457// "and will be used further." << endl;
1458 } else {
1459
1460 *fLog << warn << "Result histograms are not present in parameter " <<
1461 "list and thus are going to be created now." << endl;
1462
1463 MBinning *binsmh3x = new MBinning("BinningMH3X");
1464 binsmh3x->SetEdgesLog(energyBins,aeHisto->GetYaxis()->GetBinLowEdge(1),
1465 aeHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
1466 parlist->AddToList(binsmh3x);
1467
1468 signalHisto = new MH3(1); // Signal(E)
1469 signalHisto->SetName("MHOnSubtractionSignal");
1470 signalHisto->SetTitle("Extracted Signal");
1471 parlist->AddToList(signalHisto);
1472 signalHisto->SetBinning(&((TH1D&)signalHisto->GetHist()),
1473 aeHisto->GetYaxis());
1474
1475 offHisto = new MH3(1); // Off(E)
1476 offHisto->SetName("MHOnSubtractionOff");
1477 offHisto->SetTitle("Off Signal");
1478 parlist->AddToList(offHisto);
1479 offHisto->SetBinning(&((TH1D&)offHisto->GetHist()),
1480 aeHisto->GetYaxis());
1481
1482 signifHisto = new MH3(1); // Significance(E)
1483 signifHisto->SetName("MHOnSubtractionSignif");
1484 signifHisto->SetTitle("Significance");
1485 parlist->AddToList(signifHisto);
1486 signifHisto->SetBinning(&((TH1D&)signifHisto->GetHist()),
1487 aeHisto->GetYaxis());
1488
1489 }
1490 if (fChiSquareHisto==0x0)
1491 fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
1492 if (fSignificanceHisto==0x0)
1493 fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 41, -0.5, 40.5);
1494 if (fSummedAlphaPlots==0x0)
1495 fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha",
1496 alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1),
1497 aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
1498
1499 // The following loop fixes a common integration region for each
1500 // energy bin and alerts when no significant excess is found.
1501
1502 Double_t minLowerBin = -10; // minimum integration range
1503 Double_t maxUpperBin = 10; // minimum integration range
1504 Double_t maxAlpha = 70; // maximum integration range
1505 Double_t sumSigLiMa = 0;
1506
1507 TH1F* alphaHisto[energyBins];
1508
1509 fSigniPlotIndex = 0; // cf. ExtractSignal
1510
1511 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
1512 *fLog << inf << "Preprocessing ALPHA distribution for E bin " << energyBin << endl;
1513 char hname[80];
1514 sprintf(hname, "MHOnSubtractionAlpha%d", energyBin);
1515 alphaHisto[energyBin-1] =
1516 new TH1F(hname, "ALPHA histogram",
1517 alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1),
1518 aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
1519 alphaHisto[energyBin-1]->SetDirectory(NULL);
1520 alphaHisto[energyBin-1]->Sumw2();
1521 alphaHisto[energyBin-1] = (TH1F*)aeHisto->ProjectionX("xe", energyBin, energyBin);
1522 alphaHisto[energyBin-1]->SetName(hname);
1523 alphaHisto[energyBin-1]->SetDirectory(NULL);
1524
1525 sprintf(hname, "%2.1f < E < %2.1f",
1526 aeHisto->GetYaxis()->GetBinLowEdge(energyBin),
1527 aeHisto->GetYaxis()->GetBinLowEdge(energyBin+1));
1528 alphaHisto[energyBin-1]->SetTitle(hname);
1529 // alphaHisto[energyBin-1]->DrawCopy();
1530 alphaHisto[energyBin-1]->SetDirectory(NULL);
1531
1532 // Find signal region and significance in histogram
1533 Double_t lowerBin, upperBin, sSigLiMa;
1534 char funcName[20];
1535 sprintf(funcName, "%2d", energyBin);
1536
1537 Bool_t drawFit = kFALSE;
1538
1539 if (drawPad) {
1540 drawPad->cd(drawBase+energyBin-1);
1541 drawFit = kTRUE;
1542 }
1543
1544 Bool_t fit = FitHistogram(*alphaHisto[energyBin-1], sSigLiMa,
1545 lowerBin, upperBin, (Float_t)3.0, drawFit,
1546 funcName);
1547
1548 if (fit) {
1549 Double_t redChiSq = alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetChisquare()/
1550 alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetNDF();
1551 fChiSquareHisto->Fill(redChiSq);
1552 fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq;
1553 fSummedAlphaPlots->Add(alphaHisto[energyBin-1], 1);
1554 if (sSigLiMa < 3)
1555 *fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
1556 sumSigLiMa+=sSigLiMa;
1557
1558 // Evaluate integration range
1559 lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
1560 minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;
1561 upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
1562 maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
1563 } else {
1564 fChiSquareHisto->Fill(0);
1565 }
1566
1567// debugOut->Update();
1568
1569 } //energyBin
1570
1571 *fLog << inf << "=> Integration range for this theta bin will be " << minLowerBin << " < ALPHA < "
1572 << maxUpperBin << "," << endl;
1573 *fLog << inf << " Summed estimated significance is " << sumSigLiMa << endl;
1574
1575 // we have an idea of what is going on in the ALPHA plot...
1576 // So now we can indeed extract the signal.
1577
1578 sumSigLiMa = 0;
1579 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
1580 *fLog << inf << "Processing ALPHA distribution for E bin " << energyBin << endl;
1581 if (alphaHisto[energyBin-1]->GetEntries() == 0) {
1582 *fLog << warn << "No attempt to extract a signal from 0 events." << endl;
1583 if (draw) {
1584 TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-");
1585 lab->SetBorderSize(0);
1586 lab->Draw();
1587 }
1588 } else {
1589 char funcName[20];
1590 sprintf(funcName, "%2d", energyBin);
1591
1592 Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;
1593
1594 Bool_t drawFit = kFALSE;
1595
1596// if (drawPad) {
1597// cout << "Going to use pad " <<drawBase+energyBin-1 << endl;
1598// drawPad->cd(drawBase+energyBin-1);
1599// drawFit = kTRUE;
1600// }
1601
1602
1603 ExtractSignal(*alphaHisto[energyBin-1],
1604 sigLiMa, minLowerBin, maxUpperBin,
1605 gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
1606 funcName, drawPad, drawBase);
1607
1608 sumSigLiMa += sigLiMa;
1609 fSignificanceHisto->Fill(sigLiMa);
1610 fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;
1611
1612 // Fill Signal
1613 TH1D &h = (TH1D&)(signalHisto->GetHist());
1614 h.SetBinContent(energyBin, gammaSignal);
1615 h.SetBinError(energyBin, errorGammaSignal);
1616 h.SetEntries(h.GetEntries()+gammaSignal);
1617
1618 // Fill Off Events
1619 TH1D &ho = (TH1D&)(offHisto->GetHist());
1620 ho.SetBinContent(energyBin, off);
1621 ho.SetBinError(energyBin, errorOff);
1622 ho.SetEntries(ho.GetEntries()+off);
1623
1624 // Fill Significance
1625 TH1D &hg = (TH1D&)(signifHisto->GetHist());
1626 hg.SetBinContent(energyBin, sigLiMa);
1627 hg.SetEntries(hg.GetEntries()+sigLiMa);
1628 }
1629
1630 }//energyBin
1631
1632 *fLog << inf << "Summed significance is " << sumSigLiMa << endl;
1633
1634
1635 if (draw) {
1636 Double_t sigLiMa, lowerBin, upperBin;
1637 FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
1638 fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");
1639
1640 *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
1641 // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
1642
1643 // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
1644 // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
1645
1646 fChiSquareHisto->DrawClone();
1647 fSignificanceHisto->DrawClone();
1648 fSummedAlphaPlots->DrawClone();
1649 }
1650
1651 return kTRUE;
1652}
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669// -----------------------------------------------------------------------
1670//
1671// Extraction of Gamma signal is performed on a TH1 histogram in ALPHA
1672//
1673//
1674Bool_t MHOnSubtraction::Calc(TH1 *aHisto, MParList *parlist,
1675 const Bool_t draw, Float_t signalRegion)
1676{
1677 // Find signal region and estimate significance
1678 Double_t lowerBin, upperBin, sigLiMa;
1679 FitHistogram(*aHisto, sigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE);
1680
1681 //if (sigLiMa < 3)
1682 // *fLog << err << "No significant excess (sigma=" << sigLiMa <<")!"<<endl;
1683
1684 if (signalRegion!=0) {
1685 lowerBin = -signalRegion;
1686 upperBin = signalRegion;
1687 }
1688
1689 Double_t gammaSignal, errorGammaSignal, off, errorOff;
1690
1691 ExtractSignal(*aHisto, sigLiMa, lowerBin, upperBin,
1692 gammaSignal, errorGammaSignal, off, errorOff,
1693 (Float_t)3.0, draw);
1694
1695 fSignificance = sigLiMa;
1696
1697 *fLog << inf << "Signal is "
1698 << gammaSignal << "+-" << errorGammaSignal << endl
1699 << inf << "Significance is " << sigLiMa << endl;
1700
1701 return kTRUE;
1702}
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721// -------------------------------------------------------------------------
1722//
1723// Set the binnings and prepare the filling of the histograms
1724//
1725Bool_t MHOnSubtraction::SetupFill(const MParList *pliston)
1726{
1727 return kTRUE;
1728}
1729
1730// -------------------------------------------------------------------------
1731//
1732// Dummy Fill
1733//
1734Bool_t MHOnSubtraction::Fill(const MParContainer *pcont, const Stat_t w)
1735{
1736 return kTRUE;
1737}
1738
1739// -------------------------------------------------------------------------
1740//
1741// Draw a copy of the histogram
1742//
1743TObject *MHOnSubtraction::DrawClone(Option_t *opt) const
1744{
1745 TCanvas &c = *MakeDefCanvas("OnSubtraction", "Results from OnSubtraction");
1746 c.Divide(2,2);
1747
1748 gROOT->SetSelectedPad(NULL);
1749 gStyle->SetOptStat(0);
1750
1751 c.cd(1);
1752 fSummedAlphaPlots->DrawCopy(opt);
1753
1754 c.cd(2);
1755 fSignificanceHisto->DrawCopy(opt);
1756
1757 c.cd(3);
1758 TString drawopt="";
1759 for (fThetaHistoArray->SetIndex(0);
1760 MH3* h=(MH3*)(fThetaHistoArray->GetH());
1761 fThetaHistoArray->IncIndex())
1762 {
1763 // Get the current histogram
1764 TH1D& hist = (TH1D&)h->GetHist();
1765 fSummedAlphaPlots->SetTitle("Energy distributions for different theta bins");
1766 hist.Draw(drawopt);
1767 drawopt="SAME";
1768 }
1769 fThetaLegend->Draw();
1770
1771 c.cd(4);
1772 fChiSquareHisto->DrawCopy(opt);
1773
1774 c.Modified();
1775 c.Update();
1776
1777 return &c;
1778}
1779
1780// -------------------------------------------------------------------------
1781//
1782// Draw the histogram
1783//
1784void MHOnSubtraction::Draw(Option_t *opt)
1785{
1786 if (!gPad)
1787 MakeDefCanvas("OnSubtraction", "Results from OnSubtraction");
1788
1789 gPad->Divide(2,3);
1790
1791// Ok, at some point this will be implemented
1792
1793 gPad->Modified();
1794 gPad->Update();
1795}
1796
1797
1798
1799
1800
Note: See TracBrowser for help on using the repository browser.