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

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