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

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