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