source: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc@ 9292

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