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

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