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

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