source: trunk/Mars/mhflux/MAlphaFitter.cc@ 18716

Last change on this file since 18716 was 9576, checked in by tbretz, 15 years ago
*** empty log message ***
File size: 29.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2008
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MAlphaFitter
28//
29// Create a single Alpha-Plot. The alpha-plot is fitted online. You can
30// check the result when it is filles in the MStatusDisplay
31// For more information see MHFalseSource::FitSignificance
32//
33// For convinience (fit) the output significance is stored in a
34// container in the parlisrt
35//
36// Version 2:
37// ----------
38// + Double_t fSignificanceExc; // significance of a known excess
39//
40// Version 3:
41// ----------
42// + TArrayD fErrors; // errors of coefficients
43//
44// Version 4:
45// ----------
46// + Double_t fErrorExcess;
47// - Double_t fSignificanceExc;
48//
49//
50//////////////////////////////////////////////////////////////////////////////
51#include "MAlphaFitter.h"
52
53#include <TF1.h>
54#include <TH1.h>
55#include <TH3.h>
56
57#include <TRandom.h>
58#include <TFeldmanCousins.h>
59
60#include <TLine.h>
61#include <TLatex.h>
62#include <TVirtualPad.h>
63
64#include "MMath.h"
65#include "MString.h"
66
67#include "MLogManip.h"
68
69ClassImp(MAlphaFitter);
70
71using namespace std;
72
73// --------------------------------------------------------------------------
74//
75// Default constructor.
76//
77MAlphaFitter::MAlphaFitter(const char *name, const char *title) : fSigInt(15),
78 fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80),
79 fPolynomOrder(2), fFitBackground(kTRUE), fFunc(0),
80 fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance)
81{
82 fName = name ? name : "MAlphaFitter";
83 fTitle = title ? title : "Fit alpha";
84
85 SetSignalFunction(kGauss);
86
87 Clear();
88}
89
90// --------------------------------------------------------------------------
91//
92// Destructor
93//
94MAlphaFitter::~MAlphaFitter()
95{
96 delete fFunc;
97}
98
99// --------------------------------------------------------------------------
100//
101// Re-initializes fFunc either according to SignalFunc_t
102//
103void MAlphaFitter::SetSignalFunction(SignalFunc_t func)
104{
105 if (gROOT->GetListOfFunctions()->FindObject(""))
106 {
107 gLog << err << "MAlphaFitter::SetSignalFunction -- '' found!" << endl;
108 return;
109 }
110
111 delete fFunc;
112 fFunc = 0;
113
114 switch (func)
115 {
116 case kGauss:
117 fFunc=new TF1("", MString::Format("gaus(0) + pol%d(3)", fPolynomOrder).Data());
118 break;
119 case kThetaSq:
120 if (fPolynomOrder>0)
121 fPolynomOrder = 1;
122 fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + expo(3)");
123 break;
124 }
125
126 fSignalFunc=func;
127
128 fFunc->SetName("Dummy");
129 gROOT->GetListOfFunctions()->Remove(fFunc);
130
131 fCoefficients.Set(3+fPolynomOrder+1);
132 fCoefficients.Reset();
133
134 fErrors.Set(3+fPolynomOrder+1);
135 fErrors.Reset();
136}
137
138// --------------------------------------------------------------------------
139//
140// Reset variables which belong to results. Reset the arrays.
141//
142void MAlphaFitter::Clear(Option_t *o)
143{
144 fSignificance=0;
145 fErrorExcess=0;
146 fEventsExcess=0;
147 fEventsSignal=0;
148 fEventsBackground=0;
149
150 fChiSqSignal=0;
151 fChiSqBg=0;
152 fIntegralMax=0;
153 fScaleFactor=1;
154
155 fCoefficients.Reset();
156 fErrors.Reset();
157}
158
159// --------------------------------------------------------------------------
160//
161// Returns fFunc->Eval(d) or 0 if fFunc==NULL
162//
163Double_t MAlphaFitter::Eval(Double_t d) const
164{
165 return fFunc ? fFunc->Eval(d) : 0;
166}
167
168// --------------------------------------------------------------------------
169//
170// This function implementes the fit to the off-data as used in Fit()
171//
172Bool_t MAlphaFitter::FitOff(TH1D &h, Int_t paint)
173{
174 if (h.GetEntries()==0)
175 return kFALSE;
176
177 // First fit a polynom in the off region
178 fFunc->FixParameter(0, 0);
179 fFunc->FixParameter(1, 0);
180 fFunc->FixParameter(2, 1);
181 fFunc->ReleaseParameter(3);
182 if (fPolynomOrder!=1)
183 fFunc->FixParameter(4, 0);
184
185 for (int i=5; i<fFunc->GetNpar(); i++)
186 if (fFitBackground)
187 fFunc->ReleaseParameter(i);
188 else
189 fFunc->SetParameter(i, 0);
190
191 if (!fFitBackground)
192 return kTRUE;
193
194 if (fSignalFunc==kThetaSq)
195 {
196 const Double_t sum = h.Integral(1, 3)/3;
197 const Double_t a = sum<=1 ? 0 : TMath::Log(sum);
198 const Double_t b = -1.7;
199
200 // Do a best-guess
201 fFunc->SetParameter(3, a);
202 fFunc->SetParameter(4, b);
203 }
204
205 // options : N do not store the function, do not draw
206 // I use integral of function in bin rather than value at bin center
207 // R use the range specified in the function range
208 // Q quiet mode
209 // E Perform better Errors estimation using Minos technique
210 if (h.Fit(fFunc, "NQI", "", fBgMin, fBgMax))
211 return kFALSE;
212
213 fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();
214
215 fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters());
216 fErrors.Set(fFunc->GetNpar());
217 for (int i=3; i<fFunc->GetNpar(); i++)
218 fErrors[i] = fFunc->GetParError(i);
219
220 // ------------------------------------
221
222 if (paint)
223 {
224 if (paint==2)
225 {
226 fFunc->SetLineColor(kBlack);
227 fFunc->SetLineWidth(1);
228 }
229 else
230 {
231 fFunc->SetRange(0, 90);
232 fFunc->SetLineColor(kRed);
233 fFunc->SetLineWidth(2);
234 }
235 fFunc->Paint("same");
236 }
237
238 return kTRUE;
239}
240
241// --------------------------------------------------------------------------
242//
243// Calculate the result of the fit and set the corresponding data members
244//
245void MAlphaFitter::FitResult(const TH1D &h)
246{
247 const Double_t alphaw = h.GetXaxis()->GetBinWidth(1);
248
249 const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999);
250
251 fIntegralMax = h.GetBinLowEdge(bin+1);
252 fEventsBackground = fFunc->Integral(0, fIntegralMax)/alphaw;
253 fEventsSignal = h.Integral(1, bin);
254 fEventsExcess = fEventsSignal-fEventsBackground;
255 fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
256 fErrorExcess = MMath::ErrorExc(fEventsSignal, fEventsBackground);
257
258 // !Finitite includes IsNaN
259 if (!TMath::Finite(fSignificance))
260 fSignificance=0;
261}
262
263// --------------------------------------------------------------------------
264//
265// This is a preliminary implementation of a alpha-fit procedure for
266// all possible source positions. It will be moved into its own
267// more powerfull class soon.
268//
269// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
270// [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
271// or
272// A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
273//
274// Parameter [1] is fixed to 0 while the alpha peak should be
275// symmetric around alpha=0.
276//
277// Parameter [4] is fixed to 0 because the first derivative at
278// alpha=0 should be 0, too.
279//
280// In a first step the background is fitted between bgmin and bgmax,
281// while the parameters [0]=0 and [2]=1 are fixed.
282//
283// In a second step the signal region (alpha<sigmax) is fittet using
284// the whole function with parameters [1], [3], [4] and [5] fixed.
285//
286// The number of excess and background events are calculated as
287// s = int(hist, 0, 1.25*sigint)
288// b = int(pol2(3), 0, 1.25*sigint)
289//
290// The Significance is calculated using the Significance() member
291// function.
292//
293Bool_t MAlphaFitter::Fit(TH1D &h, Bool_t paint)
294{
295 Clear();
296
297 // Check for the region which is not filled...
298 // if (alpha0==0)
299 // return kFALSE;
300
301 // Perform fit to the off-data
302 if (!FitOff(h, paint))
303 return kFALSE;
304
305 fFunc->ReleaseParameter(0); // It is also released by SetParLimits later on
306 //func.ReleaseParameter(1); // It is also released by SetParLimits later on
307 fFunc->ReleaseParameter(2);
308 for (int i=3; i<fFunc->GetNpar(); i++)
309 fFunc->FixParameter(i, fFunc->GetParameter(i));
310
311
312 // Do not allow signals smaller than the background
313 const Double_t alpha0 = h.GetBinContent(1);
314 const Double_t s = fSignalFunc==kGauss ? fFunc->GetParameter(3) : TMath::Exp(fFunc->GetParameter(3));
315 const Double_t A = alpha0-s;
316 //const Double_t dA = TMath::Abs(A);
317 //fFunc->SetParLimits(0, -dA*4, dA*4); // SetParLimits also releases the parameter
318 fFunc->SetParLimits(2, 0, 90); // SetParLimits also releases the parameter
319
320 // Now fit a gaus in the on region on top of the polynom
321 fFunc->SetParameter(0, A);
322 fFunc->SetParameter(2, fSigMax*0.75);
323
324 // options : N do not store the function, do not draw
325 // I use integral of function in bin rather than value at bin center
326 // R use the range specified in the function range
327 // Q quiet mode
328 // E Perform better Errors estimation using Minos technique
329 h.Fit(fFunc, "NQI", "", 0, fSigMax);
330
331 fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF();
332 fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters());
333 for (int i=0; i<3; i++)
334 fErrors[i] = fFunc->GetParError(i);
335 //const Bool_t ok = NDF>0 && chi2<2.5*NDF;
336
337 // ------------------------------------
338 if (paint)
339 {
340 fFunc->SetLineColor(kGreen);
341 fFunc->SetLineWidth(2);
342 fFunc->Paint("same");
343 }
344 // ------------------------------------
345
346 //const Double_t s = fFunc->Integral(0, fSigInt)/alphaw;
347 fFunc->SetParameter(0, 0);
348 fFunc->SetParameter(2, 1);
349 //const Double_t b = fFunc->Integral(0, fSigInt)/alphaw;
350 //fSignificance = MMath::SignificanceLiMaSigned(s, b);
351
352 // Calculate the fit result and set the corresponding data members
353 FitResult(h);
354
355 return kTRUE;
356}
357
358Double_t MAlphaFitter::DoOffFit(const TH1D &hon, const TH1D &hof, Bool_t paint)
359{
360 if (fSignalFunc!=kThetaSq)
361 return 0;
362
363 // ----------------------------------------------------------------------------
364
365 const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999);
366
367 MAlphaFitter fit(*this);
368 fit.EnableBackgroundFit();
369 fit.SetBackgroundFitMin(0);
370
371 // produce a histogram containing the off-samples from on-source and
372 // off-source in the off-source region and the on-data in the source-region
373 TH1D h(hof);
374 h.SetDirectory(0);
375 h.Add(&hon);
376
377 h.Scale(0.5);
378 for (int i=1; i<=bin+3; i++)
379 {
380 h.SetBinContent(i, hof.GetBinContent(i));
381 h.SetBinError( i, hof.GetBinError(i));
382 }
383
384 // Now fit the off-data
385 if (!fit.FitOff(h, paint?2:0)) // FIXME: Show fit!
386 return -1;
387
388 // Calculate fit-result
389 fit.FitResult(h);
390
391 // Do a gaussian error propagation to calculated the error of
392 // the background estimated from the fit
393 const Double_t ea = fit.fErrors[3];
394 const Double_t eb = fit.fErrors[4];
395 const Double_t a = fit.fCoefficients[3];
396 const Double_t b = fit.fCoefficients[4];
397
398 const Double_t t = fIntegralMax;
399
400 const Double_t ex = TMath::Exp(t*b);
401 const Double_t eab = TMath::Exp(a)/b;
402
403 const Double_t eA = ex-1;
404 const Double_t eB = t*ex - eA/b;
405
406 const Double_t w = h.GetXaxis()->GetBinWidth(1);
407
408 // Error of estimated background
409 const Double_t er = TMath::Abs(eab)*TMath::Hypot(eA*ea, eB*eb)/w;
410
411 // Calculate arbitrary scale factor from propagated error from the
412 // condition: sqrt(alpha*background) = est.background/est.error
413 // const Double_t bg = hof.Integral(1, bin);
414 // const Double_t sc = bg * er*er / (fit2.GetEventsBackground()*fit2.GetEventsBackground());
415 // Assuming that bg and fit2.GetEventsBackground() are rather identical:
416 const Double_t sc = er*er / fit.fEventsBackground;
417
418
419 /*
420 cout << MMath::SignificanceLiMaSigned(hon.Integral(1, bin), fit.GetEventsBackground()/sc, sc) << " ";
421 cout << sc << " " << fit.GetEventsBackground() << " ";
422 cout << fit.fChiSqBg << endl;
423 */
424 return sc;
425}
426
427Bool_t MAlphaFitter::Fit(const TH1D &hon, const TH1D &hof, Double_t alpha, Bool_t paint)
428{
429 TH1D h(hon);
430 h.SetDirectory(0);
431 h.Add(&hof, -1); // substracts also number of entries!
432 h.SetEntries(hon.GetEntries());
433
434 MAlphaFitter fit(*this);
435 fit.SetPolynomOrder(0);
436 if (alpha<=0 || !fit.Fit(h, paint))
437 return kFALSE;
438
439 fChiSqSignal = fit.fChiSqSignal;
440 fChiSqBg = fit.fChiSqBg;
441 fCoefficients = fit.fCoefficients;
442 fErrors = fit.fErrors;
443
444 // ----------------------------------------------------------------------------
445
446 const Double_t scale = DoOffFit(hon, hof, paint);
447 if (scale<0)
448 return kFALSE;
449
450 // ----------------------------------------------------------------------------
451
452 const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999);
453
454 fIntegralMax = hon.GetBinLowEdge(bin+1);
455 fEventsBackground = hof.Integral(1, bin);
456 fEventsSignal = hon.Integral(1, bin);
457 fEventsExcess = fEventsSignal-fEventsBackground;
458 fScaleFactor = alpha;
459 fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground/alpha, alpha);
460 fErrorExcess = MMath::ErrorExc(fEventsSignal, fEventsBackground/alpha, alpha);
461
462 // !Finitite includes IsNaN
463 if (!TMath::Finite(fSignificance))
464 fSignificance=0;
465
466 return kTRUE;
467}
468
469// --------------------------------------------------------------------------
470//
471// Calculate the upper limit for fEventsSignal number of observed events
472// and fEventsBackground number of background events.
473//
474// Therefor TFeldmanCousin is used.
475//
476// The Feldman-Cousins method as described in PRD V57 #7, p3873-3889
477//
478Double_t MAlphaFitter::CalcUpperLimit() const
479{
480 // get a FeldmanCousins calculation object with the default limits
481 // of calculating a 90% CL with the minimum signal value scanned
482 // = 0.0 and the maximum signal value scanned of 50.0
483 TFeldmanCousins f;
484 f.SetMuStep(0.05);
485 f.SetMuMax(100);
486 f.SetMuMin(0);
487 f.SetCL(90);
488
489 return f.CalculateUpperLimit(fEventsSignal, fEventsBackground);
490}
491
492void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size, Bool_t draw) const
493{
494 const Double_t w = GetGausSigma();
495 const Double_t m = fIntegralMax;
496
497 const Int_t l1 = w<=0 ? 0 : (Int_t)TMath::Ceil(-TMath::Log10(w));
498 const Int_t l2 = m<=0 ? 0 : (Int_t)TMath::Ceil(-TMath::Log10(m));
499 const TString fmt = MString::Format("\\sigma_{L/M}=%%.1f \\omega=%%.%df\\circ E=%%d B=%%d x<%%.%df \\tilde\\chi_{b}=%%.1f \\tilde\\chi_{s}=%%.1f c=%%.1f f=%%.2f",
500 l1<1?1:l1+1, l2<1?1:l2+1);
501 const TString txt = MString::Format(fmt.Data(), fSignificance, w, (int)fEventsExcess,
502 (int)fEventsBackground, m, fChiSqBg, fChiSqSignal,
503 fCoefficients[3], fScaleFactor);
504
505 // This is totaly weired but the only way to get both options
506 // working with this nonsense implementation of TLatex
507 TLatex text(x, y, txt);
508 text.SetBit(TLatex::kTextNDC);
509 text.SetTextSize(size);
510 if (draw)
511 text.DrawLatex(x, y, txt);
512 else
513 text.Paint();
514
515 TLine line;
516 line.SetLineColor(14);
517 if (draw)
518 line.DrawLine(m, gPad->GetUymin(), m, gPad->GetUymax());
519 else
520 line.PaintLine(m, gPad->GetUymin(), m, gPad->GetUymax());
521}
522
523void MAlphaFitter::Copy(TObject &o) const
524{
525 MAlphaFitter &f = static_cast<MAlphaFitter&>(o);
526
527 // Setup
528 f.fSigInt = fSigInt;
529 f.fSigMax = fSigMax;
530 f.fBgMin = fBgMin;
531 f.fBgMax = fBgMax;
532 f.fScaleMin = fScaleMin;
533 f.fScaleMax = fScaleMax;
534 f.fPolynomOrder = fPolynomOrder;
535 f.fFitBackground= fFitBackground;
536 f.fSignalFunc = fSignalFunc;
537 f.fScaleMode = fScaleMode;
538 f.fScaleUser = fScaleUser;
539 f.fStrategy = fStrategy;
540 f.fCoefficients.Set(fCoefficients.GetSize());
541 f.fCoefficients.Reset();
542 f.fErrors.Set(fCoefficients.GetSize());
543 f.fErrors.Reset();
544
545 // Result
546 f.fSignificance = fSignificance;
547 f.fErrorExcess = fErrorExcess;
548 f.fEventsExcess = fEventsExcess;
549 f.fEventsSignal = fEventsSignal;
550 f.fEventsBackground = fEventsBackground;
551 f.fChiSqSignal = fChiSqSignal;
552 f.fChiSqBg = fChiSqBg;
553 f.fIntegralMax = fIntegralMax;
554 f.fScaleFactor = fScaleFactor;
555
556 // Function
557 delete f.fFunc;
558
559 f.fFunc = new TF1(*fFunc);
560 f.fFunc->SetName("Dummy");
561 gROOT->GetListOfFunctions()->Remove(f.fFunc);
562}
563
564void MAlphaFitter::Print(Option_t *o) const
565{
566 *fLog << GetDescriptor() << ": Fitting..." << endl;
567 *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl;
568 *fLog << " ...signal function: ";
569 switch (fSignalFunc)
570 {
571 case kGauss: *fLog << "gauss(x)/pol" << fPolynomOrder; break;
572 case kThetaSq: *fLog << "gauss(sqrt(x))/expo"; break;
573 }
574 *fLog << endl;
575 if (!fFitBackground)
576 *fLog << " ...no background." << endl;
577 else
578 {
579 *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl;
580 *fLog << " ...polynom order " << fPolynomOrder << endl;
581 *fLog << " ...scale mode: ";
582 switch (fScaleMode)
583 {
584 case kNone: *fLog << "none."; break;
585 case kEntries: *fLog << "entries."; break;
586 case kIntegral: *fLog << "integral."; break;
587 case kOffRegion: *fLog << "off region (integral between " << fScaleMin << " and " << fScaleMax << ")"; break;
588 case kBackground: *fLog << "background (integral between " << fBgMin << " and " << fBgMax << ")"; break;
589 case kLeastSquare: *fLog << "least square (N/A)"; break;
590 case kUserScale: *fLog << "user def (" << fScaleUser << ")"; break;
591 }
592 *fLog << endl;
593 }
594
595 if (TString(o).Contains("result"))
596 {
597 *fLog << "Result:" << endl;
598 *fLog << " - Significance (Li/Ma) " << fSignificance << endl;
599 *fLog << " - Excess Events " << fEventsExcess << endl;
600 *fLog << " - Signal Events " << fEventsSignal << endl;
601 *fLog << " - Background Events " << fEventsBackground << endl;
602 *fLog << " - E/sqrt(B>=Alpha) " << fEventsExcess/TMath::Sqrt(TMath::Max(fEventsBackground,fScaleFactor)) << endl;
603 *fLog << " - Chi^2/ndf (Signal) " << fChiSqSignal << endl;
604 *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl;
605 *fLog << " - Signal integrated up to " << fIntegralMax << "°" << endl;
606 *fLog << " - Off Scale Alpha (Off) " << fScaleFactor << endl;
607 }
608}
609
610Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, UInt_t bin, Bool_t paint)
611{
612 const TString name(MString::Format("TempAlphaEnergy%06d", gRandom->Integer(1000000)));
613
614 TH1D *h = hon.ProjectionZ(name, 0, hon.GetNbinsX()+1, bin, bin, "E");
615 h->SetDirectory(0);
616
617 const Bool_t rc = Fit(*h, paint);
618
619 delete h;
620
621 return rc;
622}
623
624Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
625{
626 const TString name(MString::Format("TempAlphaTheta%06d", gRandom->Integer(1000000)));
627
628 TH1D *h = hon.ProjectionZ(name, bin, bin, 0, hon.GetNbinsY()+1, "E");
629 h->SetDirectory(0);
630
631 const Bool_t rc = Fit(*h, paint);
632
633 delete h;
634
635 return rc;
636}
637/*
638Bool_t MAlphaFitter::FitTime(const TH3D &hon, UInt_t bin, Bool_t paint)
639{
640 const TString name(Form("TempAlphaTime%06d", gRandom->Integer(1000000)));
641
642 hon.GetZaxis()->SetRange(bin,bin);
643 TH1D *h = (TH1D*)hon.Project3D("ye");
644 hon.GetZaxis()->SetRange(-1,-1);
645
646 h->SetDirectory(0);
647
648 const Bool_t rc = Fit(*h, paint);
649 delete h;
650 return rc;
651}
652*/
653Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint)
654{
655 const TString name(MString::Format("TempAlpha%06d", gRandom->Integer(1000000)));
656
657 TH1D *h = hon.ProjectionZ(name, 0, hon.GetNbinsX()+1, 0, hon.GetNbinsY()+1, "E");
658 h->SetDirectory(0);
659
660 const Bool_t rc = Fit(*h, paint);
661
662 delete h;
663
664 return rc;
665}
666
667Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
668{
669 const TString name1(MString::Format("TempAlpha%06d_on", gRandom->Integer(1000000)));
670 const TString name0(MString::Format("TempAlpha%06d_off", gRandom->Integer(1000000)));
671
672 TH1D *h1 = hon.ProjectionZ(name1, 0, hon.GetNbinsX()+1, bin, bin, "E");
673 h1->SetDirectory(0);
674
675 TH1D *h0 = hof.ProjectionZ(name0, 0, hof.GetNbinsX()+1, bin, bin, "E");
676 h0->SetDirectory(0);
677
678 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
679
680 delete h0;
681 delete h1;
682
683 return rc;
684}
685
686Bool_t MAlphaFitter::FitTheta(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
687{
688 const TString name1(MString::Format("TempAlpha%06d_on", gRandom->Integer(1000000)));
689 const TString name0(MString::Format("TempAlpha%06d_off", gRandom->Integer(1000000)));
690
691 TH1D *h1 = hon.ProjectionZ(name1, bin, bin, 0, hon.GetNbinsY()+1, "E");
692 h1->SetDirectory(0);
693
694 TH1D *h0 = hof.ProjectionZ(name0, bin, bin, 0, hof.GetNbinsY()+1, "E");
695 h0->SetDirectory(0);
696
697 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
698
699 delete h0;
700 delete h1;
701
702 return rc;
703}
704/*
705Bool_t MAlphaFitter::FitTime(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
706{
707 const TString name1(Form("TempAlphaTime%06d_on", gRandom->Integer(1000000)));
708 const TString name0(Form("TempAlphaTime%06d_off", gRandom->Integer(1000000)));
709
710 hon.GetZaxis()->SetRange(bin,bin);
711 TH1D *h1 = (TH1D*)hon.Project3D("ye");
712 hon.GetZaxis()->SetRange(-1,-1);
713 h1->SetDirectory(0);
714
715 hof.GetZaxis()->SetRange(bin,bin);
716 TH1D *h0 = (TH1D*)hof.Project3D("ye");
717 hof.GetZaxis()->SetRange(-1,-1);
718 h0->SetDirectory(0);
719
720 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
721
722 delete h0;
723 delete h1;
724
725 return rc;
726}
727*/
728
729Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint)
730{
731 const TString name1(MString::Format("TempAlpha%06d_on", gRandom->Integer(1000000)));
732 const TString name0(MString::Format("TempAlpha%06d_off", gRandom->Integer(1000000)));
733
734 TH1D *h1 = hon.ProjectionZ(name1, 0, hon.GetNbinsX()+1, 0, hon.GetNbinsY()+1, "E");
735 h1->SetDirectory(0);
736
737 TH1D *h0 = hof.ProjectionZ(name0, 0, hof.GetNbinsX()+1, 0, hof.GetNbinsY()+1, "E");
738 h0->SetDirectory(0);
739
740 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
741
742 delete h0;
743 delete h1;
744
745 return rc;
746}
747
748Bool_t MAlphaFitter::ApplyScaling(const TH3D &hon, TH3D &hof, UInt_t bin) const
749{
750 const TString name1(MString::Format("TempAlpha%06d_on", gRandom->Integer(1000000)));
751 const TString name0(MString::Format("TempAlpha%06d_off", gRandom->Integer(1000000)));
752
753 TH1D *h1 = hon.ProjectionZ(name1, 0, hon.GetNbinsX()+1, bin, bin, "E");
754 h1->SetDirectory(0);
755
756 TH1D *h0 = hof.ProjectionZ(name0, 0, hof.GetNbinsX()+1, bin, bin, "E");
757 h0->SetDirectory(0);
758
759 const Double_t scale = Scale(*h0, *h1);
760
761 delete h0;
762 delete h1;
763
764 for (int x=0; x<=hof.GetNbinsX()+1; x++)
765 for (int z=0; z<=hof.GetNbinsZ()+1; z++)
766 {
767 hof.SetBinContent(x, bin, z, hof.GetBinContent(x, bin, z)*scale);
768 hof.SetBinError( x, bin, z, hof.GetBinError( x, bin, z)*scale);
769 }
770
771 return scale>0;
772}
773
774Bool_t MAlphaFitter::ApplyScaling(const TH3D &hon, TH3D &hof) const
775{
776 for (int y=0; y<=hof.GetNbinsY()+1; y++)
777 ApplyScaling(hon, hof, y);
778
779 return kTRUE;
780}
781
782Double_t MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
783{
784 Float_t scaleon = 1;
785 Float_t scaleof = 1;
786 switch (fScaleMode)
787 {
788 case kNone:
789 return 1;
790
791 case kEntries:
792 scaleon = on.GetEntries();
793 scaleof = of.GetEntries();
794 break;
795
796 case kIntegral:
797 scaleon = on.Integral();
798 scaleof = of.Integral();
799 break;
800
801 case kOffRegion:
802 {
803 const Int_t min = on.GetXaxis()->FindFixBin(fScaleMin);
804 const Int_t max = on.GetXaxis()->FindFixBin(fScaleMax);
805 scaleon = on.Integral(min, max);
806 scaleof = of.Integral(min, max);
807 }
808 break;
809
810 case kBackground:
811 {
812 const Int_t min = on.GetXaxis()->FindFixBin(fBgMin);
813 const Int_t max = on.GetXaxis()->FindFixBin(fBgMax);
814 scaleon = on.Integral(min, max);
815 scaleof = of.Integral(min, max);
816 }
817 break;
818
819 case kUserScale:
820 scaleon = fScaleUser;
821 break;
822
823 // This is just to make some compiler happy
824 default:
825 return 1;
826 }
827
828 if (scaleof!=0)
829 {
830 of.Scale(scaleon/scaleof);
831 return scaleon/scaleof;
832 }
833 else
834 {
835 of.Reset();
836 return 0;
837 }
838}
839
840Double_t MAlphaFitter::GetMinimizationValue() const
841{
842 switch (fStrategy)
843 {
844 case kSignificance:
845 return -GetSignificance();
846 case kSignificanceChi2:
847 return -GetSignificance()/GetChiSqSignal();
848 case kSignificanceLogExcess:
849 if (GetEventsExcess()<1)
850 return 0;
851 return -GetSignificance()*TMath::Log10(GetEventsExcess());
852 case kSignificanceSqrtExcess:
853 if (GetEventsExcess()<1)
854 return 0;
855 return -GetSignificance()*TMath::Sqrt(GetEventsExcess());
856 case kSignificanceExcess:
857 return -GetSignificance()*GetEventsExcess();
858 case kExcess:
859 return -GetEventsExcess();
860 case kGaussSigma:
861 return GetGausSigma();
862 case kWeakSource:
863 if (GetEventsExcess()<1)
864 return 0;
865 return -GetEventsExcess()/TMath::Sqrt(TMath::Max(GetEventsBackground(), GetEventsBackground()));
866 case kWeakSourceLogExcess:
867 if (GetEventsExcess()<1)
868 return 0;
869 return -GetEventsExcess()/TMath::Sqrt(TMath::Max(GetEventsBackground(), GetEventsBackground()))*TMath::Log10(GetEventsExcess());
870 }
871 return 0;
872}
873
874Int_t MAlphaFitter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
875{
876 Bool_t rc = kFALSE;
877
878 //void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; }
879 //void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; }
880
881 if (IsEnvDefined(env, prefix, "SignalIntegralMax", print))
882 {
883 SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalIntegralMax", fSigInt));
884 rc = kTRUE;
885 }
886 if (IsEnvDefined(env, prefix, "SignalFitMax", print))
887 {
888 SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalFitMax", fSigMax));
889 rc = kTRUE;
890 }
891 if (IsEnvDefined(env, prefix, "BackgroundFitMax", print))
892 {
893 SetBackgroundFitMax(GetEnvValue(env, prefix, "BackgroundFitMax", fBgMax));
894 rc = kTRUE;
895 }
896 if (IsEnvDefined(env, prefix, "BackgroundFitMin", print))
897 {
898 SetBackgroundFitMin(GetEnvValue(env, prefix, "BackgroundFitMin", fBgMin));
899 rc = kTRUE;
900 }
901 if (IsEnvDefined(env, prefix, "ScaleMin", print))
902 {
903 SetScaleMin(GetEnvValue(env, prefix, "ScaleMin", fScaleMin));
904 rc = kTRUE;
905 }
906 if (IsEnvDefined(env, prefix, "ScaleMax", print))
907 {
908 SetScaleMax(GetEnvValue(env, prefix, "ScaleMax", fScaleMax));
909 rc = kTRUE;
910 }
911 if (IsEnvDefined(env, prefix, "PolynomOrder", print))
912 {
913 SetPolynomOrder(GetEnvValue(env, prefix, "PolynomOrder", fPolynomOrder));
914 rc = kTRUE;
915 }
916
917 if (IsEnvDefined(env, prefix, "MinimizationStrategy", print))
918 {
919 TString txt = GetEnvValue(env, prefix, "MinimizationStrategy", "");
920 txt = txt.Strip(TString::kBoth);
921 txt.ToLower();
922 if (txt==(TString)"significance")
923 fStrategy = kSignificance;
924 if (txt==(TString)"significancechi2")
925 fStrategy = kSignificanceChi2;
926 if (txt==(TString)"significanceexcess")
927 fStrategy = kSignificanceExcess;
928 if (txt==(TString)"excess")
929 fStrategy = kExcess;
930 if (txt==(TString)"gausssigma" || txt==(TString)"gaussigma")
931 fStrategy = kGaussSigma;
932 if (txt==(TString)"weaksource")
933 fStrategy = kWeakSource;
934 rc = kTRUE;
935 }
936 if (IsEnvDefined(env, prefix, "Scale", print))
937 {
938 fScaleUser = GetEnvValue(env, prefix, "Scale", fScaleUser);
939 rc = kTRUE;
940 }
941 if (IsEnvDefined(env, prefix, "ScaleMode", print))
942 {
943 TString txt = GetEnvValue(env, prefix, "ScaleMode", "");
944 txt = txt.Strip(TString::kBoth);
945 txt.ToLower();
946 if (txt==(TString)"none")
947 fScaleMode = kNone;
948 if (txt==(TString)"entries")
949 fScaleMode = kEntries;
950 if (txt==(TString)"integral")
951 fScaleMode = kIntegral;
952 if (txt==(TString)"offregion")
953 fScaleMode = kOffRegion;
954 if (txt==(TString)"background")
955 fScaleMode = kBackground;
956 if (txt==(TString)"leastsquare")
957 fScaleMode = kLeastSquare;
958 if (txt==(TString)"userscale")
959 fScaleMode = kUserScale;
960 if (txt==(TString)"fixed")
961 FixScale();
962 rc = kTRUE;
963 }
964 if (IsEnvDefined(env, prefix, "SignalFunction", print))
965 {
966 TString txt = GetEnvValue(env, prefix, "SignalFunction", "");
967 txt = txt.Strip(TString::kBoth);
968 txt.ToLower();
969 if (txt==(TString)"gauss" || txt==(TString)"gaus")
970 SetSignalFunction(kGauss);
971 if (txt==(TString)"thetasq")
972 SetSignalFunction(kThetaSq);
973 rc = kTRUE;
974 }
975
976 return rc;
977}
Note: See TracBrowser for help on using the repository browser.