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