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