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

Last change on this file since 2313 was 2282, checked in by wittek, 21 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(lowerBin);
335 upperBin = alphaHisto.GetBinLowEdge(upperBin)+alphaHisto.GetBinWidth(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
391 if (gausPol) {
392 Double_t gausMean = gausPol->GetParameter(1);
393 Double_t gausSigma = gausPol->GetParameter(2);
394 po = new TF1("po"+funcName,"pol3(0)",
395 -signalRegionFactor*gausSigma+gausMean,
396 signalRegionFactor*gausSigma+gausMean);
397 po->SetParameters(gausPol->GetParameter(3),gausPol->GetParameter(4),
398 gausPol->GetParameter(5),gausPol->GetParameter(6));
399
400 TF1 *ga = new TF1("ga"+funcName,"gaus(0)",-40,40);
401 ga->SetParameters(gausPol->GetParameter(0),gausPol->GetParameter(1),
402 gausPol->GetParameter(2));
403
404 if (draw) {
405 alphaHisto.Draw();
406 gausPol->Draw("SAME"); //Maybe not even necessary?
407
408 po->SetLineColor(4);
409 po->SetLineWidth(2);
410 po->Draw("SAME");
411
412 ga->SetLineColor(93);
413 ga->SetLineWidth(2);
414 ga->Draw("SAME");
415
416 char legendTitle[80];
417 sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f",
418 lowerBin, upperBin);
419
420 TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle);
421
422 legend->SetBorderSize(0);
423 legend->SetFillColor(10);
424 legend->SetFillStyle(0);
425 legend->AddEntry(gausPol, "combined N_{on}","l");
426 legend->AddEntry(po,"polynomial N_{bg} Signal region","l");
427 legend->AddEntry(ga, "putative gaussian N_{S}","l");
428 legend->Draw();
429 }
430 } // gausPol
431
432
433 if (pol) {
434 po = pol;
435
436 if (draw) {
437 alphaHisto.Draw();
438
439 po->SetLineColor(6);
440 po->SetLineWidth(2);
441 po->Draw("SAME");
442
443 char legendTitle[80];
444 sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f",
445 lowerBin, upperBin);
446
447 TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle);
448
449 legend->SetBorderSize(0);
450 legend->SetFillColor(10);
451 legend->SetFillStyle(0);
452 legend->AddEntry(po,"polynomial N_{bg} Signal region","l");
453 legend->Draw();
454 }
455 } // pol
456
457 Double_t nOn = 0;
458 Double_t eNOn = 0;
459
460 Int_t binNumberLow = alphaHisto.GetXaxis()->FindBin(lowerBin);
461 Int_t binNumberHigh = alphaHisto.GetXaxis()->FindBin(upperBin);
462
463 for (Int_t bin=binNumberLow; bin<binNumberHigh+1; bin++) {
464 nOn += alphaHisto.GetBinContent(bin);
465 eNOn += alphaHisto.GetBinError(bin) * alphaHisto.GetBinError(bin);
466 } //for bin
467 eNOn = sqrt(eNOn);
468
469 // Evaluate background
470
471 Double_t nOff = 0;
472 Double_t eNOff = 0;
473 for (Int_t bin=binNumberLow; bin<binNumberHigh+1; bin++) {
474 Double_t x = .5*(alphaHisto.GetBinLowEdge(bin)+alphaHisto.GetBinLowEdge(bin+1));
475 Double_t offEvts = po->Eval(x);
476 //cout << bin << ": " << offEvts << endl;
477 nOff += offEvts;
478 } //for bin
479 eNOff = sqrt(fabs(nOff));
480
481 if (nOn==0) // there should not be a negative number of signal events
482 nOff=0;
483
484 if (nOff<0) { // there should not be a negative number of off events
485 nOff=0;
486 eNOff=0;
487 }
488
489 *fLog << inf << "nEvts = " << nOn << "+-" << eNOn << ", ";
490
491 off = nOff; errorOff = eNOff;
492 gammaSignal = nOn-nOff; errorGammaSignal = sqrt(eNOn*eNOn + eNOff*eNOff);
493
494 *fLog << inf << "nBack = " << nOff << "+-" << eNOff << ", ";
495 *fLog << inf << "nSig = " << gammaSignal << "+-" << errorGammaSignal << endl;
496
497 sigLiMa = CalcSignificance(nOn, nOff, 1);
498 // Double_t sigFit=(ga->GetParameter(1))/(ga->GetParError(1)); //Mean / sigMean
499
500 *fLog << inf << "Significance: "<<sigLiMa<<" sigma "<<endl;
501
502 if (draw) {
503 TPaveText *pt = new TPaveText(0.11,.74,.57,.88,"NDC ");
504 char tx[60];
505 sprintf(tx, "Excess: %2.2f #pm %2.2f", nOn-nOff, sqrt(eNOn*eNOn + eNOff*eNOff));
506 pt->AddText(tx);
507 sprintf(tx, "Off: %2.2f #pm %2.2f", nOff, eNOff);
508 pt->AddText(tx);
509 sprintf(tx, "Significance: %2.2f #sigma", sigLiMa);
510 pt->AddText(tx);
511 pt->SetFillStyle(0);
512 pt->SetBorderSize(0);
513 pt->SetTextAlign(12);
514 pt->Draw("");
515 }
516 if (draw||drawPad) {
517
518 // Plot significance vs alpha_cut
519
520 Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0);
521 Int_t maxBin = centerBin > alphaHisto.GetNbinsX()- centerBin ?
522 alphaHisto.GetNbinsX()- centerBin : centerBin;
523
524 const Int_t totalBins = centerBin;
525 Float_t alpha[totalBins];
526 Float_t signi[totalBins];
527 Float_t maxSigni = 0;
528
529 for (Int_t i=0; i < maxBin; i++) {
530 Double_t nOn = 0; Double_t eNOn = 0;
531 Double_t nOff = 0; Double_t eNOff = 0;
532 for (Int_t bin=centerBin-i; bin<centerBin+i+1; bin++) {
533 nOn += alphaHisto.GetBinContent(bin);
534 eNOn += alphaHisto.GetBinError(bin) * alphaHisto.GetBinError(bin);
535 Double_t x = .5*(alphaHisto.GetBinLowEdge(bin)+alphaHisto.GetBinLowEdge(bin+1));
536 Double_t offEvts = po->Eval(x);
537 nOff += offEvts;
538 } //for bin
539 eNOn = sqrt(eNOn);
540 eNOff = sqrt(nOff);
541 alpha[i] =
542 (alphaHisto.GetXaxis()->GetBinLowEdge(centerBin+i+1)-
543 alphaHisto.GetXaxis()->GetBinLowEdge(centerBin-i))/2;
544 signi[i] = CalcSignificance(nOn, nOff, 1);
545 maxSigni = maxSigni > signi[i] ? maxSigni : signi[i];
546 }
547
548 if (!drawPad) {
549 TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600);
550 c3->cd();
551 }
552
553 if (drawPad)
554 drawPad->cd(drawBase-1);
555
556 TGraph* gr = new TGraph(totalBins-2,alpha,signi);
557 TString drawOpt = "L";
558
559 if (draw || (drawPad && fSigniPlotIndex == 0))
560 drawOpt += "A";
561
562 gr->Draw(drawOpt);
563 if (drawPad && fSigniPlotIndex == 0) {
564 gr->GetXaxis()->SetTitle("|#alpha_{max}|");
565 gr->GetYaxis()->SetTitle("Significance");
566 }
567 gr->SetMarkerStyle(2);
568 gr->SetMarkerSize(1);
569 gr->SetMarkerColor(4+fSigniPlotIndex);
570 gr->SetLineColor(4+fSigniPlotIndex);
571 gr->SetLineWidth(1);
572 gr->SetTitle("Significance vs ALPHA integration range");
573 gr->Draw("L");
574
575 fSigniPlotIndex++;
576
577 } //draw
578
579 return kTRUE;
580}
581
582// -----------------------------------------------------------------------
583//
584// Extraction of Gamma signal is performed on a MHAlphaEnergyTheta or
585// MHAlphaEnergyTime object expected in the ParList.
586//
587Bool_t MHOnSubtraction::Calc(MParList *parlist, const Bool_t Draw)
588{
589 //Find source histograms
590 fHistogramType = "Theta";
591 MHAlphaEnergyTheta* mhisttheta = (MHAlphaEnergyTheta*)parlist->FindObject("MHAlphaEnergyTheta");
592 if (mhisttheta) return Calc((TH3D*)mhisttheta->GetHist(), parlist, Draw);
593
594 fHistogramType = "Time";
595 MHAlphaEnergyTime* mhisttime = (MHAlphaEnergyTime*)parlist->FindObject("MHAlphaEnergyTime");
596 if (mhisttime) return Calc((TH3D*)mhisttime->GetHist(), parlist, Draw);
597
598 *fLog << err << "No MHAlphaEnergyTheta type object found in the parameter list. Aborting." << endl;
599 return kFALSE;
600}
601
602// -----------------------------------------------------------------------
603//
604// Extraction of Gamma signal is performed on a TH3D histogram in
605// ALPHA, Energy and Theta.
606//
607Bool_t MHOnSubtraction::CalcAET(TH3D *aetHisto, MParList *parlist, const Bool_t Draw)
608{
609 // Analyze aetHisto
610 // -------------------------------------
611 Int_t alphaBins = aetHisto->GetNbinsX(); // # of alpha bins
612 Int_t energyBins = aetHisto->GetNbinsY();
613 Int_t thetaBins = aetHisto->GetNbinsZ();
614
615 // We output an array of histograms to the parlist.
616 // Create a template for such a histogram, needed by MHArray
617 // -------------------------------------
618 MBinning *binsfft = new MBinning("BinningMH3X");
619 binsfft->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1),
620 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
621 parlist->AddToList(binsfft);
622 MH3 *energyHistogram = new MH3(1); // The template histogram in energy
623 // for a given theta value
624 energyHistogram->SetName("MHOnSubtractionEnergyHistogramTemplate");
625 parlist->AddToList(energyHistogram);
626
627 fThetaHistoArray = new MHArray("MHOnSubtractionEnergyHistogramTemplate", kTRUE,
628 "MHOnSubtractionGammaSignalArray",
629 "Array of histograms in energy bins");
630 fThetaHistoArray->SetupFill(parlist);
631 parlist->AddToList(fThetaHistoArray);
632
633 // Canvases---direct debug output for the time being
634 // -------------------------------------
635 TCanvas *c3 = new TCanvas("c3", "Plots by MHOnSubtraction::ExtractSignal", 800,600);
636 cout << thetaBins << " x " << energyBins << endl;
637 c3->Divide(thetaBins,energyBins);
638
639 TCanvas *c4a = new TCanvas("c4a", "Energy distributions for different ZA", 800,600);
640
641 TH1D* histalphaon[energyBins*thetaBins]; // ALPHA histograms
642
643 fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
644 fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 40.5, -0.5, 41);
645 fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha",
646 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
647 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
648
649 // -------------------------------------
650
651 fThetaLegend = new TLegend(0.83, 0.07, 0.98, 0.42, "Energy Distributions");
652 fThetaLegend->SetBorderSize(1);
653
654 Double_t overallSigLiMa = 0;
655
656 for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) {
657
658 char hname[80];
659 sprintf(hname, "Energy distribution for %s bin %d", fHistogramType.Data(),
660 thetaBin);
661 *fLog << all << "Calculating " << hname << endl;
662
663 Double_t minLowerBin = -10; // minimum integration range
664 Double_t maxUpperBin = 10; // minimum integration range
665 Double_t maxAlpha = 70; // maximum integration range
666
667 Double_t sumSigLiMa = 0;
668
669 // This loop just fixes the integration range
670 // And alerts when no significant excess is found.
671
672 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
673
674 sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin);
675 histalphaon[(thetaBin-1)*(energyBins)+energyBin] =
676 new TH1D(hname, "ALPHA histogram",
677 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
678 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
679 histalphaon[(thetaBin-1)*(energyBins)+energyBin]->Sumw2();
680 histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetTitle(hname);
681
682 sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d",
683 energyBin, thetaBin);
684 *fLog << inf << hname << endl;
685
686 // now we fill one histogram for one particular set
687 // of Energy and Theta...
688
689 if (aetHisto->InheritsFrom("TH3D")||aetHisto->InheritsFrom("TH2D")) {
690 aetHisto->GetYaxis()->SetRange(energyBin,energyBin); // E
691 if (aetHisto->InheritsFrom("TH3D"))
692 aetHisto->GetZaxis()->SetRange(thetaBin,thetaBin); // theta
693 histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto->Project3D("xe");
694 } else {
695 histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto;
696 }
697
698 *fLog << inf << "This histogram has "
699 << histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetEntries()
700 << " entries" << endl;
701
702 sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin);
703 histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetName(hname);
704
705 // Test: Find signal region and significance in histogram
706 Double_t lowerBin, upperBin, sSigLiMa;
707 FitHistogram(*histalphaon[(thetaBin-1)*(energyBins)+energyBin],
708 sSigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE);
709
710 Double_t redChiSq = histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetChisquare()/histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetNDF();
711
712 fChiSquareHisto->Fill(redChiSq);
713 fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq;
714
715 // histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetProb() <<endl;
716
717 fSummedAlphaPlots->Add(histalphaon[(thetaBin-1)*(energyBins)+energyBin], 1);
718
719 if (sSigLiMa < 3)
720 *fLog << inf << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
721 sumSigLiMa+=sSigLiMa;
722
723 // Evaluate integration range
724
725 lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
726 minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;
727
728 upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
729 maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
730
731 } //energyBin
732
733 *fLog << inf << "Integration range for this " << fHistogramType
734 << " value will be " << minLowerBin << " < ALPHA < "
735 << maxUpperBin << endl;
736
737 *fLog << inf << "Summed estd. significance for this " << fHistogramType
738 << " value is " << sumSigLiMa << endl;
739
740 // Create Histogram in Energy for this Theta value
741 cout << "STARTIGN REAL CALC1 "<< endl;
742 fThetaHistoArray->AddHistogram();
743 cout << "STARTIGN REAL CALC1a "<< endl;
744 fThetaHistoArray->SetIndex(thetaBin-1);
745
746 MH3 &mThetaHisto = *(MH3*)(fThetaHistoArray->GetH());
747 mThetaHisto.Print();
748
749 TH1D &thetaHisto = (TH1D&)(mThetaHisto.GetHist());
750 sprintf(hname,"Energy distribution for theta bin %d", thetaBin);
751 thetaHisto.SetTitle(hname);
752 thetaHisto.SetEntries(0);
753
754 // we have a rough idea of what is going on in the ALPHA plot...
755 // So now we can indeed extract the signal.
756
757// cout << "STARTIGN REAL CALC "<< endl;
758
759 sumSigLiMa = 0;
760 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
761
762 sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d",
763 energyBin, thetaBin);
764 *fLog << inf << hname << endl;
765
766 Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;
767
768 c3->cd((thetaBin-1)*(energyBins)+energyBin);
769
770 ExtractSignal(*histalphaon[(thetaBin-1)*(energyBins)+energyBin],
771 sigLiMa, minLowerBin, maxUpperBin,
772 gammaSignal, errorGammaSignal,
773 off, errorOff, (Float_t)3, kTRUE);
774
775 fSignificanceHisto->Fill(sigLiMa);
776 fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;
777
778 thetaHisto.SetBinContent(energyBin, gammaSignal*TMath::Exp(7-energyBin));
779 thetaHisto.SetBinError(energyBin, errorGammaSignal);
780 thetaHisto.SetEntries(thetaHisto.GetEntries()+gammaSignal);
781
782 sumSigLiMa += sigLiMa;
783
784 }//energyBin
785
786 *fLog << inf << "Summed significance for this " << fHistogramType
787 << " value is " << sumSigLiMa << endl;
788
789 //fit exponential function to data
790 TF1* e = new TF1("expF","expo",0,5);
791 e->SetLineWidth(3);
792 e->SetLineColor(thetaBin);
793 // e->SetParLimits(1, -0.5,3); //limits on slope
794
795 if (fSlope!=20.0) {
796 e->SetParameter(1, fSlope);
797 *fLog << inf << "Fixing slope to " << e->GetParameter(1) << endl;
798 }
799
800 TString eDrawOptions = Draw ? "RIQ" : "RIQ0";
801 cout << "doing the fit" << endl;
802 thetaHisto.Fit("expF");
803
804 Double_t expoSlope = e->GetParameter(1);
805
806 *fLog << inf << "Determined slope for energy distribution is "
807 << expoSlope << endl;
808
809 Double_t integralEvts =
810 e->Integral(aetHisto->GetYaxis()->GetBinLowEdge(1),
811 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
812
813 *fLog << inf << "Integral in E range [" <<
814 aetHisto->GetYaxis()->GetBinLowEdge(1) << ";" <<
815 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1) << "] is " <<
816 integralEvts << endl;
817
818 // Plot Energy histogram
819
820 c4a->cd(0);
821 thetaHisto.SetLineColor(thetaBin);
822 if (thetaBin==1) {
823 thetaHisto.Draw();
824 } else {
825 thetaHisto.Draw("SAME");
826 }
827
828 sprintf(hname,"Theta bin %d",thetaBin);
829 fThetaLegend->AddEntry(&thetaHisto, hname, "l");
830
831 overallSigLiMa += sumSigLiMa;
832
833 } //thetaBin
834
835 fThetaLegend->Draw();
836
837 Double_t sigLiMa, lowerBin, upperBin;
838 FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
839 fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");
840
841 *fLog << inf << "Summed overall significance is " << overallSigLiMa << endl;
842
843 *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
844 // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
845
846 // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
847 // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
848
849 return kTRUE;
850}
851
852
853
854
855
856
857
858
859// -----------------------------------------------------------------------
860//
861// Extraction of Gamma signal is performed on a TH3D histogram in
862// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms
863// in energy and Theta.
864//
865Bool_t MHOnSubtraction::Calc(TH3 *aetHisto, MParList *parlist, const Bool_t Draw)
866{
867 // Analyze aetHisto
868 // -------------------------------------
869 Int_t energyBins = aetHisto->GetNbinsY();
870 Int_t thetaBins = aetHisto->GetNbinsZ();
871
872 *fLog << "Total events: " << aetHisto->GetEntries() << endl;
873 *fLog << "Total energy bins: " << energyBins << endl;
874 *fLog << "Total theta bins: " << thetaBins << endl;
875
876 // We output results in an array of histograms to the parameter list.
877 // Create a template for such a histogram, needed by MH3
878 // -------------------------------------
879 MBinning *binsmh3x = new MBinning("BinningMH3X");
880 // dummy binning, real one follows some lines down
881 binsmh3x->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1),
882 aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
883 parlist->AddToList(binsmh3x);
884
885 MBinning *binsmh3y = new MBinning("BinningMH3Y");
886 // dummy binning, real one follows some lines down
887 binsmh3y->SetEdges(thetaBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
888 aetHisto->GetZaxis()->GetBinLowEdge(thetaBins+1));
889 parlist->AddToList(binsmh3y);
890
891 MH3 *signalHisto = new MH3(2); // Signal(E,t)
892 signalHisto->SetupFill(parlist);
893 parlist->AddToList(signalHisto);
894 signalHisto->SetName("MHOnSubtractionSignal");
895 signalHisto->SetTitle("Gamma Events");
896
897 MH3 *offHisto = new MH3(2); // Off(E,t)
898 offHisto->SetupFill(parlist);
899 parlist->AddToList(offHisto);
900 offHisto->SetName("MHOnSubtractionOff");
901 offHisto->SetTitle("Background Events");
902
903 MH3 *signifHisto = new MH3(2); // Significance(E,t)
904 signifHisto->SetupFill(parlist);
905 parlist->AddToList(signifHisto);
906 signifHisto->SetName("MHOnSubtractionSignif");
907 signifHisto->SetTitle("Significance");
908
909// if (!fChiSquareHisto)
910// fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
911// if (!fSignificanceHisto)
912// fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 40.5, -0.5, 41);
913// if (!fSummedAlphaPlots)
914// fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha",
915// alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
916// aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
917
918 TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600);
919 c4->Divide(energyBins+3,thetaBins,0,0,0);
920
921 TH2D& signalTH2DHist = (TH2D&)signalHisto->GetHist();
922 TH2D& offTH2DHist = (TH2D&)offHisto->GetHist();
923 TH2D& signifTH2DHist = (TH2D&)signifHisto->GetHist();
924
925 signalTH2DHist.Reset();
926 signalTH2DHist.SetTitle("Gamma Signal");
927 signalTH2DHist.SetXTitle("E [GeV]");
928 signalTH2DHist.SetYTitle("theta [deg]");
929 signalHisto->SetBinning(&signalTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
930 signalTH2DHist.Sumw2();
931
932 offTH2DHist.Reset();
933 offTH2DHist.SetTitle("Off Signal");
934 offTH2DHist.SetXTitle("E [GeV]");
935 offTH2DHist.SetYTitle("theta [deg]");
936 offHisto->SetBinning(&offTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
937 offTH2DHist.Sumw2();
938
939 signifTH2DHist.Reset();
940 signifTH2DHist.SetTitle("Significance");
941 signifTH2DHist.SetXTitle("E [GeV]");
942 signifTH2DHist.SetYTitle("theta [deg]");
943 signifHisto->SetBinning(&signifTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
944 signifTH2DHist.Sumw2();
945
946 for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) {
947
948 *fLog << dbg << "--------------------------------------" << endl;
949 *fLog << dbg << "Processing ALPHA plots for theta bin " << thetaBin << endl;
950 *fLog << dbg << "--------------------------------------" << endl;
951 aetHisto->GetZaxis()->SetRange(thetaBin, thetaBin);
952 TH2F* aeHisto = (TH2F*)aetHisto->Project3D("yxe");
953 aeHisto->SetDirectory(NULL);
954 char hname[80];
955 sprintf(hname, "%2.1f < #theta < %2.1f",
956 aetHisto->GetZaxis()->GetBinLowEdge(thetaBin),
957 aetHisto->GetZaxis()->GetBinLowEdge(thetaBin+1));
958 *fLog << inf << "There are " << aeHisto->GetEntries()
959 << " entries in " << hname << endl;
960 aeHisto->SetTitle(hname);
961 sprintf(hname, "MHOnSubtractionThetaBin%d", thetaBin);
962 aeHisto->SetName(hname);
963
964 c4->cd((energyBins+3)*(thetaBin-1)+1);
965 aeHisto->SetMarkerColor(3);
966 aeHisto->DrawCopy();
967
968 c4->Update();
969
970 // We hand over a nonstandard parameter list, which
971 // contails lower-dimensional result histograms
972 // signalHisto, offHisto and signifHisto
973
974// cout << fLog->GetDebugLevel() << endl;
975 fLog->SetDebugLevel(1);
976
977 MParList *parlistth2 = new MParList();
978 // parlistth2->SetNullOutput();
979 parlistth2->AddToList(binsmh3x);
980
981 MH3* signalHisto2 = new MH3(1); // Signal(E)
982 signalHisto2->SetupFill(parlistth2);
983 parlistth2->AddToList(signalHisto2);
984 signalHisto2->SetName("MHOnSubtractionSignal");
985 signalHisto2->SetTitle("Gamma Events");
986
987 MH3* offHisto2 = new MH3(1); // Off(E)
988 offHisto2->SetupFill(parlistth2);
989 parlistth2->AddToList(offHisto2);
990 offHisto2->SetName("MHOnSubtractionOff");
991 offHisto2->SetTitle("Background Events");
992
993 MH3* signifHisto2 = new MH3(1); // Significance(E)
994 signifHisto2->SetupFill(parlistth2);
995 parlistth2->AddToList(signifHisto2);
996 signifHisto2->SetName("MHOnSubtractionSignif");
997 signifHisto2->SetTitle("Significance");
998
999 fLog->SetDebugLevel(-1);
1000
1001 TH1D& signalTH1DHist = (TH1D&)signalHisto2->GetHist();
1002 TH1D& offTH1DHist = (TH1D&)offHisto2->GetHist();
1003 TH1D& signifTH1DHist = (TH1D&)signifHisto2->GetHist();
1004
1005 signalTH1DHist.Sumw2();
1006 offTH1DHist.Sumw2();
1007 signifTH1DHist.Sumw2();
1008
1009 TH2Calc(aeHisto, parlistth2, kFALSE, c4,
1010 (energyBins+3)*(thetaBin-1)+4);
1011
1012
1013 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
1014
1015 // cout << "filling " << thetaBin << " " << energyBin << ": "
1016// << signalTH1DHist.GetBinContent(energyBin) << "+-"
1017// << signalTH1DHist.GetBinError(energyBin) << " "
1018// << offTH1DHist.GetBinContent(energyBin) << "+-"
1019// << offTH1DHist.GetBinError(energyBin) << endl;
1020
1021 if (signalTH1DHist.GetBinContent(energyBin)>=0.0) {
1022 signalTH2DHist.SetBinContent(energyBin, thetaBin,
1023 signalTH1DHist.GetBinContent(energyBin));
1024 signalTH2DHist.SetBinError(energyBin, thetaBin,
1025 signalTH1DHist.GetBinError(energyBin));
1026 }
1027
1028 if (offTH1DHist.GetBinContent(energyBin)>=0.0) {
1029 offTH2DHist.SetBinContent(energyBin, thetaBin,
1030 offTH1DHist.GetBinContent(energyBin));
1031 offTH2DHist.SetBinError(energyBin, thetaBin,
1032 offTH1DHist.GetBinError(energyBin));
1033 }
1034
1035 signifTH2DHist.SetBinContent(energyBin, thetaBin,
1036 signifTH1DHist.GetBinContent(energyBin));
1037 signifTH2DHist.SetBinError(energyBin, thetaBin,
1038 signifTH1DHist.GetBinError(energyBin));
1039 } //energyBin
1040
1041
1042 c4->cd((energyBins+3)*(thetaBin-1)+2);
1043 sprintf(hname,"c4_%d",(energyBins+3)*(thetaBin-1)+2);
1044 TPad* tp = (TPad*)(gROOT->FindObject(hname));
1045 tp->SetLogx();
1046 signalTH1DHist.SetLineColor(2);
1047 signalTH1DHist.DrawCopy();
1048 offTH1DHist.SetLineColor(4);
1049 offTH1DHist.DrawCopy("SAME");
1050 c4->Update();
1051
1052 signalTH2DHist.SetEntries(signalTH2DHist.GetEntries()+signalTH1DHist.GetEntries());
1053 offTH2DHist.SetEntries(offTH2DHist.GetEntries()+offTH1DHist.GetEntries());
1054 signifTH2DHist.SetEntries(signifTH2DHist.GetEntries()+signifTH1DHist.GetEntries());
1055
1056 delete signifHisto2;
1057 delete offHisto2;
1058 delete signalHisto2;
1059 delete parlistth2;
1060
1061 }
1062
1063
1064 c4->Update();
1065
1066 return kTRUE;
1067}
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// Extraction of Gamma signal is performed on a TH3D histogram in
1124// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms
1125// in energy and Theta.
1126//
1127Bool_t MHOnSubtraction::CalcLightCurve(TH3 *aetHisto, MParList *parlist, const Bool_t draw)
1128{
1129 // Analyze aetHisto
1130 // -------------------------------------
1131 Int_t alphaBins = aetHisto->GetNbinsX();
1132 Int_t energyBins = aetHisto->GetNbinsY();
1133 Int_t timeBins = aetHisto->GetNbinsZ();
1134
1135 *fLog << "Total events: " << aetHisto->GetEntries() << endl;
1136 *fLog << "Total energy bins: " << energyBins << endl;
1137 *fLog << "Total time bins: " << timeBins << endl;
1138
1139 // We output results in an array of histograms to the parameter list.
1140 // Create a template for such a histogram, needed by MH3
1141 // -------------------------------------
1142 MBinning *binsmh3x = new MBinning("BinningMH3X");
1143 // dummy binning, real one follows some lines down
1144 binsmh3x->SetEdges(timeBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
1145 aetHisto->GetZaxis()->GetBinLowEdge(timeBins+1));
1146 parlist->AddToList(binsmh3x);
1147
1148 MH3 *signalHisto = new MH3(1); // Signal(t)
1149 signalHisto->SetupFill(parlist);
1150 parlist->AddToList(signalHisto);
1151 signalHisto->SetName("MHOnSubtractionSignal");
1152 signalHisto->SetTitle("Gamma Events");
1153
1154 MH3 *offHisto = new MH3(1); // Off(t)
1155 offHisto->SetupFill(parlist);
1156 parlist->AddToList(offHisto);
1157 offHisto->SetName("MHOnSubtractionOff");
1158 offHisto->SetTitle("Background Events");
1159
1160 MH3 *signifHisto = new MH3(1); // Significance(t)
1161 signifHisto->SetupFill(parlist);
1162 parlist->AddToList(signifHisto);
1163 signifHisto->SetName("MHOnSubtractionSignif");
1164 signifHisto->SetTitle("Significance");
1165
1166 TH1D& signalTH1DHist = (TH1D&)signalHisto->GetHist();
1167 TH1D& offTH1DHist = (TH1D&)offHisto->GetHist();
1168 TH1D& signifTH1DHist = (TH1D&)signifHisto->GetHist();
1169
1170 signalTH1DHist.Reset();
1171 signalTH1DHist.SetTitle("Gamma Signal");
1172 signalTH1DHist.SetXTitle("Time [MJD]");
1173 signalTH1DHist.SetYTitle("Events");
1174 signalHisto->SetBinning(&signalTH1DHist, aetHisto->GetZaxis());
1175 signalTH1DHist.Sumw2();
1176
1177 offTH1DHist.Reset();
1178 offTH1DHist.SetTitle("Background Signal");
1179 offTH1DHist.SetXTitle("Time [MJD]");
1180 offTH1DHist.SetYTitle("Events");
1181 offHisto->SetBinning(&offTH1DHist, aetHisto->GetZaxis());
1182 offTH1DHist.Sumw2();
1183
1184 signifTH1DHist.Reset();
1185 signifTH1DHist.SetTitle("Significance");
1186 signifTH1DHist.SetXTitle("Time [MJD]");
1187 signifTH1DHist.SetYTitle("Significance #sigma");
1188 signifHisto->SetBinning(&signifTH1DHist, aetHisto->GetZaxis());
1189 signifTH1DHist.Sumw2();
1190
1191 //The following histogram is an additional histogram needed for
1192 //the light curve
1193
1194 TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600);
1195 c4->Divide(8,7,0,0,0);
1196
1197
1198
1199 // The following loop fixes a common integration region for each
1200 // energy bin and alerts when no significant excess is found.
1201
1202 Double_t minLowerBin = -10; // minimum integration range
1203 Double_t maxUpperBin = 10; // minimum integration range
1204 Double_t maxAlpha = 70; // maximum integration range
1205 Double_t sumSigLiMa = 0;
1206
1207 TH1F* alphaHisto[timeBins];
1208
1209 for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {
1210
1211 *fLog << dbg << "--------------------------------------" << endl;
1212 *fLog << dbg << "Processing ALPHA plots for time bin " << timeBin << endl;
1213 *fLog << dbg << "--------------------------------------" << endl;
1214
1215 aetHisto->GetZaxis()->SetRange(timeBin, timeBin);
1216
1217 char hname[80];
1218 sprintf(hname, "MHOnSubtractionTimeBin%d", timeBin);
1219 alphaHisto[timeBin-1] =
1220 new TH1F(hname, "ALPHA histogram",
1221 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
1222 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
1223 alphaHisto[timeBin-1]->SetDirectory(NULL);
1224 alphaHisto[timeBin-1]->Sumw2();
1225 alphaHisto[timeBin-1] = (TH1F*)aetHisto->Project3D("xe");
1226 alphaHisto[timeBin-1]->SetName(hname);
1227 alphaHisto[timeBin-1]->SetDirectory(NULL);
1228
1229 sprintf(hname, "%6.0f < t < %6.0f",
1230 aetHisto->GetZaxis()->GetBinLowEdge(timeBin),
1231 aetHisto->GetZaxis()->GetBinLowEdge(timeBin+1));
1232 *fLog << inf << "There are " << alphaHisto[timeBin-1]->GetEntries()
1233 << " entries in " << hname << endl;
1234 alphaHisto[timeBin-1]->SetTitle(hname);
1235
1236 // Find signal region and significance in histogram
1237 Double_t lowerBin, upperBin, sSigLiMa;
1238 char funcName[20];
1239 sprintf(funcName, "%2d", timeBin);
1240
1241 Bool_t drawFit = kTRUE;
1242
1243 c4->cd(timeBin);
1244// alphaHisto[timeBin-1]->SetMarkerColor(3);
1245 alphaHisto[timeBin-1]->DrawCopy();
1246
1247 c4->Update();
1248
1249 fSigniPlotColor = 0;
1250 ;
1251 Bool_t fit = FitHistogram(*alphaHisto[timeBin-1], sSigLiMa,
1252 lowerBin, upperBin, (Float_t)3.0, drawFit,
1253 funcName);
1254
1255 if (fit) {
1256 if (sSigLiMa < 3)
1257 *fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
1258 sumSigLiMa+=sSigLiMa;
1259
1260 // Evaluate integration range
1261 lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
1262 minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;
1263 upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
1264 maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
1265 }
1266
1267
1268
1269
1270
1271
1272 } //timeBin
1273
1274 *fLog << inf << "=> Integration range will be " << minLowerBin << " < ALPHA < "
1275 << maxUpperBin << "," << endl;
1276 *fLog << inf << " Summed estimated significance is " << sumSigLiMa << endl;
1277
1278 // we have an idea of what is going on in the ALPHA plot...
1279 // So now we can indeed extract the signal.
1280
1281 sumSigLiMa = 0;
1282 for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {
1283 *fLog << inf << "Processing ALPHA distribution for time bin " << timeBin << endl;
1284 if (alphaHisto[timeBin-1]->GetEntries() == 0) {
1285 *fLog << warn << "No attempt to extract a signal from 0 events." << endl;
1286 if (draw) {
1287 TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-");
1288 lab->SetBorderSize(0);
1289 lab->Draw();
1290 }
1291 } else {
1292 char funcName[20];
1293 sprintf(funcName, "%2d", timeBin);
1294
1295 Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;
1296
1297 Bool_t drawFit = kFALSE;
1298
1299 ExtractSignal(*alphaHisto[timeBin-1],
1300 sigLiMa, minLowerBin, maxUpperBin,
1301 gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
1302 funcName, NULL, 1);
1303
1304 sumSigLiMa += sigLiMa;
1305
1306 fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;
1307
1308 // Fill Signal
1309 TH1D &h = (TH1D&)(signalHisto->GetHist());
1310 h.SetBinContent(timeBin, gammaSignal);
1311 h.SetBinError(timeBin, errorGammaSignal);
1312 h.SetEntries(h.GetEntries()+gammaSignal);
1313
1314 // Fill Off Events
1315 TH1D &ho = (TH1D&)(offHisto->GetHist());
1316 ho.SetBinContent(timeBin, off);
1317 ho.SetBinError(timeBin, errorOff);
1318 ho.SetEntries(ho.GetEntries()+off);
1319
1320 // Fill Significance
1321 TH1D &hg = (TH1D&)(signifHisto->GetHist());
1322 hg.SetBinContent(timeBin, sigLiMa);
1323 hg.SetEntries(hg.GetEntries()+sigLiMa);
1324 }
1325
1326 }//timeBin
1327
1328 *fLog << inf << "Summed significance is " << sumSigLiMa << endl;
1329
1330
1331 if (draw) {
1332 Double_t sigLiMa, lowerBin, upperBin;
1333 FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
1334 fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");
1335
1336 *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
1337 // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
1338
1339 // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
1340 // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
1341
1342 fChiSquareHisto->DrawClone();
1343 fSignificanceHisto->DrawClone();
1344 fSummedAlphaPlots->DrawClone();
1345 }
1346
1347 return kTRUE;
1348}
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// Extraction of Gamma signal is performed on a TH2 histogram in
1434// ALPHA and Energy. Output to parameter list is TH1 histogram in
1435// energy
1436//
1437Bool_t MHOnSubtraction::TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t draw, TPad *drawPad, Int_t drawBase)
1438{
1439
1440 fSigniPlotColor = 0;
1441
1442 // Analyze aeHisto
1443 // -------------------------------------
1444 const Int_t alphaBins = aeHisto->GetNbinsX(); // # of alpha bins
1445 const Int_t energyBins = aeHisto->GetNbinsY();
1446
1447 // Check availability of result histograms
1448 // -------------------------------------
1449
1450 MH3* signalHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignal","MH3");
1451 MH3* offHisto = (MH3*)parlist->FindObject("MHOnSubtractionOff","MH3");
1452 MH3* signifHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignif","MH3");
1453
1454 if (signalHisto && offHisto && signifHisto) {
1455// *fLog << dbg << "Result histograms are present in parameter list " <<
1456// "and will be used further." << endl;
1457 } else {
1458
1459 *fLog << warn << "Result histograms are not present in parameter " <<
1460 "list and thus are going to be created now." << endl;
1461
1462 MBinning *binsmh3x = new MBinning("BinningMH3X");
1463 binsmh3x->SetEdgesLog(energyBins,aeHisto->GetYaxis()->GetBinLowEdge(1),
1464 aeHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
1465 parlist->AddToList(binsmh3x);
1466
1467 signalHisto = new MH3(1); // Signal(E)
1468 signalHisto->SetName("MHOnSubtractionSignal");
1469 signalHisto->SetTitle("Extracted Signal");
1470 parlist->AddToList(signalHisto);
1471 signalHisto->SetBinning(&((TH1D&)signalHisto->GetHist()),
1472 aeHisto->GetYaxis());
1473
1474 offHisto = new MH3(1); // Off(E)
1475 offHisto->SetName("MHOnSubtractionOff");
1476 offHisto->SetTitle("Off Signal");
1477 parlist->AddToList(offHisto);
1478 offHisto->SetBinning(&((TH1D&)offHisto->GetHist()),
1479 aeHisto->GetYaxis());
1480
1481 signifHisto = new MH3(1); // Significance(E)
1482 signifHisto->SetName("MHOnSubtractionSignif");
1483 signifHisto->SetTitle("Significance");
1484 parlist->AddToList(signifHisto);
1485 signifHisto->SetBinning(&((TH1D&)signifHisto->GetHist()),
1486 aeHisto->GetYaxis());
1487
1488 }
1489 if (fChiSquareHisto==0x0)
1490 fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
1491 if (fSignificanceHisto==0x0)
1492 fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 40.5, -0.5, 41);
1493 if (fSummedAlphaPlots==0x0)
1494 fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha",
1495 alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1),
1496 aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
1497
1498 // The following loop fixes a common integration region for each
1499 // energy bin and alerts when no significant excess is found.
1500
1501 Double_t minLowerBin = -10; // minimum integration range
1502 Double_t maxUpperBin = 10; // minimum integration range
1503 Double_t maxAlpha = 70; // maximum integration range
1504 Double_t sumSigLiMa = 0;
1505
1506 TH1F* alphaHisto[energyBins];
1507
1508 fSigniPlotIndex = 0; // cf. ExtractSignal
1509
1510 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
1511 *fLog << inf << "Preprocessing ALPHA distribution for E bin " << energyBin << endl;
1512 char hname[80];
1513 sprintf(hname, "MHOnSubtractionAlpha%d", energyBin);
1514 alphaHisto[energyBin-1] =
1515 new TH1F(hname, "ALPHA histogram",
1516 alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1),
1517 aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
1518 alphaHisto[energyBin-1]->SetDirectory(NULL);
1519 alphaHisto[energyBin-1]->Sumw2();
1520 alphaHisto[energyBin-1] = (TH1F*)aeHisto->ProjectionX("xe", energyBin, energyBin);
1521 alphaHisto[energyBin-1]->SetName(hname);
1522 alphaHisto[energyBin-1]->SetDirectory(NULL);
1523
1524 sprintf(hname, "%2.1f < E < %2.1f",
1525 aeHisto->GetYaxis()->GetBinLowEdge(energyBin),
1526 aeHisto->GetYaxis()->GetBinLowEdge(energyBin+1));
1527 alphaHisto[energyBin-1]->SetTitle(hname);
1528 // alphaHisto[energyBin-1]->DrawCopy();
1529 alphaHisto[energyBin-1]->SetDirectory(NULL);
1530
1531 // Find signal region and significance in histogram
1532 Double_t lowerBin, upperBin, sSigLiMa;
1533 char funcName[20];
1534 sprintf(funcName, "%2d", energyBin);
1535
1536 Bool_t drawFit = kFALSE;
1537
1538 if (drawPad) {
1539 drawPad->cd(drawBase+energyBin-1);
1540 drawFit = kTRUE;
1541 }
1542
1543 Bool_t fit = FitHistogram(*alphaHisto[energyBin-1], sSigLiMa,
1544 lowerBin, upperBin, (Float_t)3.0, drawFit,
1545 funcName);
1546
1547 if (fit) {
1548 Double_t redChiSq = alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetChisquare()/
1549 alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetNDF();
1550 fChiSquareHisto->Fill(redChiSq);
1551 fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq;
1552 fSummedAlphaPlots->Add(alphaHisto[energyBin-1], 1);
1553 if (sSigLiMa < 3)
1554 *fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
1555 sumSigLiMa+=sSigLiMa;
1556
1557 // Evaluate integration range
1558 lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
1559 minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;
1560 upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
1561 maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
1562 } else {
1563 fChiSquareHisto->Fill(0);
1564 }
1565
1566// debugOut->Update();
1567
1568 } //energyBin
1569
1570 *fLog << inf << "=> Integration range for this theta bin will be " << minLowerBin << " < ALPHA < "
1571 << maxUpperBin << "," << endl;
1572 *fLog << inf << " Summed estimated significance is " << sumSigLiMa << endl;
1573
1574 // we have an idea of what is going on in the ALPHA plot...
1575 // So now we can indeed extract the signal.
1576
1577 sumSigLiMa = 0;
1578 for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
1579 *fLog << inf << "Processing ALPHA distribution for E bin " << energyBin << endl;
1580 if (alphaHisto[energyBin-1]->GetEntries() == 0) {
1581 *fLog << warn << "No attempt to extract a signal from 0 events." << endl;
1582 if (draw) {
1583 TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-");
1584 lab->SetBorderSize(0);
1585 lab->Draw();
1586 }
1587 } else {
1588 char funcName[20];
1589 sprintf(funcName, "%2d", energyBin);
1590
1591 Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;
1592
1593 Bool_t drawFit = kFALSE;
1594
1595// if (drawPad) {
1596// cout << "Going to use pad " <<drawBase+energyBin-1 << endl;
1597// drawPad->cd(drawBase+energyBin-1);
1598// drawFit = kTRUE;
1599// }
1600
1601
1602 ExtractSignal(*alphaHisto[energyBin-1],
1603 sigLiMa, minLowerBin, maxUpperBin,
1604 gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
1605 funcName, drawPad, drawBase);
1606
1607 sumSigLiMa += sigLiMa;
1608 fSignificanceHisto->Fill(sigLiMa);
1609 fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;
1610
1611 // Fill Signal
1612 TH1D &h = (TH1D&)(signalHisto->GetHist());
1613 h.SetBinContent(energyBin, gammaSignal);
1614 h.SetBinError(energyBin, errorGammaSignal);
1615 h.SetEntries(h.GetEntries()+gammaSignal);
1616
1617 // Fill Off Events
1618 TH1D &ho = (TH1D&)(offHisto->GetHist());
1619 ho.SetBinContent(energyBin, off);
1620 ho.SetBinError(energyBin, errorOff);
1621 ho.SetEntries(ho.GetEntries()+off);
1622
1623 // Fill Significance
1624 TH1D &hg = (TH1D&)(signifHisto->GetHist());
1625 hg.SetBinContent(energyBin, sigLiMa);
1626 hg.SetEntries(hg.GetEntries()+sigLiMa);
1627 }
1628
1629 }//energyBin
1630
1631 *fLog << inf << "Summed significance is " << sumSigLiMa << endl;
1632
1633
1634 if (draw) {
1635 Double_t sigLiMa, lowerBin, upperBin;
1636 FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
1637 fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");
1638
1639 *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
1640 // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
1641
1642 // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
1643 // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
1644
1645 fChiSquareHisto->DrawClone();
1646 fSignificanceHisto->DrawClone();
1647 fSummedAlphaPlots->DrawClone();
1648 }
1649
1650 return kTRUE;
1651}
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668// -----------------------------------------------------------------------
1669//
1670// Extraction of Gamma signal is performed on a TH1 histogram in ALPHA
1671//
1672//
1673Bool_t MHOnSubtraction::Calc(TH1 *aHisto, MParList *parlist,
1674 const Bool_t draw, Float_t signalRegion)
1675{
1676 // Find signal region and estimate significance
1677 Double_t lowerBin, upperBin, sigLiMa;
1678 FitHistogram(*aHisto, sigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE);
1679
1680 //if (sigLiMa < 3)
1681 // *fLog << err << "No significant excess (sigma=" << sigLiMa <<")!"<<endl;
1682
1683 if (signalRegion!=0) {
1684 lowerBin = -signalRegion;
1685 upperBin = signalRegion;
1686 }
1687
1688 Double_t gammaSignal, errorGammaSignal, off, errorOff;
1689
1690 ExtractSignal(*aHisto, sigLiMa, lowerBin, upperBin,
1691 gammaSignal, errorGammaSignal, off, errorOff,
1692 (Float_t)3.0, draw);
1693
1694 fSignificance = sigLiMa;
1695
1696 *fLog << inf << "Signal is "
1697 << gammaSignal << "+-" << errorGammaSignal << endl
1698 << inf << "Significance is " << sigLiMa << endl;
1699
1700 return kTRUE;
1701}
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720// -------------------------------------------------------------------------
1721//
1722// Set the binnings and prepare the filling of the histograms
1723//
1724Bool_t MHOnSubtraction::SetupFill(const MParList *pliston)
1725{
1726 return kTRUE;
1727}
1728
1729// -------------------------------------------------------------------------
1730//
1731// Dummy Fill
1732//
1733Bool_t MHOnSubtraction::Fill(const MParContainer *pcont, const Stat_t w)
1734{
1735 return kTRUE;
1736}
1737
1738// -------------------------------------------------------------------------
1739//
1740// Draw a copy of the histogram
1741//
1742TObject *MHOnSubtraction::DrawClone(Option_t *opt) const
1743{
1744 TCanvas &c = *MakeDefCanvas("OnSubtraction", "Results from OnSubtraction");
1745 c.Divide(2,2);
1746
1747 gROOT->SetSelectedPad(NULL);
1748 gStyle->SetOptStat(0);
1749
1750 c.cd(1);
1751 fSummedAlphaPlots->DrawCopy(opt);
1752
1753 c.cd(2);
1754 fSignificanceHisto->DrawCopy(opt);
1755
1756 c.cd(3);
1757 TString drawopt="";
1758 for (fThetaHistoArray->SetIndex(0);
1759 MH3* h=(MH3*)(fThetaHistoArray->GetH());
1760 fThetaHistoArray->IncIndex())
1761 {
1762 // Get the current histogram
1763 TH1D& hist = (TH1D&)h->GetHist();
1764 fSummedAlphaPlots->SetTitle("Energy distributions for different theta bins");
1765 hist.Draw(drawopt);
1766 drawopt="SAME";
1767 }
1768 fThetaLegend->Draw();
1769
1770 c.cd(4);
1771 fChiSquareHisto->DrawCopy(opt);
1772
1773 c.Modified();
1774 c.Update();
1775
1776 return &c;
1777}
1778
1779// -------------------------------------------------------------------------
1780//
1781// Draw the histogram
1782//
1783void MHOnSubtraction::Draw(Option_t *opt)
1784{
1785 if (!gPad)
1786 MakeDefCanvas("OnSubtraction", "Results from OnSubtraction");
1787
1788 gPad->Divide(2,3);
1789
1790// Ok, at some point this will be implemented
1791
1792 gPad->Modified();
1793 gPad->Update();
1794}
1795
1796
1797
1798
1799
Note: See TracBrowser for help on using the repository browser.