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

Last change on this file since 6209 was 5971, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 15.9 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// 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
53ClassImp(MAlphaFitter);
54
55using namespace std;
56
57void 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//
102Bool_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 fFunc->FixParameter(4, 0);
120 fFunc->SetParLimits(2, 0, 90);
121 fFunc->SetParLimits(3, -1, 1);
122
123 const Double_t alpha0 = h.GetBinContent(1);
124 const Double_t alphaw = h.GetXaxis()->GetBinWidth(1);
125
126 // Check for the regios which is not filled...
127 if (alpha0==0)
128 return kFALSE; //*fLog << warn << "Histogram empty." << endl;
129
130 // First fit a polynom in the off region
131 fFunc->FixParameter(0, 0);
132 fFunc->FixParameter(2, 1);
133 fFunc->ReleaseParameter(3);
134
135 for (int i=5; i<fFunc->GetNpar(); i++)
136 fFunc->ReleaseParameter(i);
137
138 // options : N do not store the function, do not draw
139 // I use integral of function in bin rather than value at bin center
140 // R use the range specified in the function range
141 // Q quiet mode
142 h.Fit(fFunc, "NQI", "", bgmin, bgmax);
143
144 fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();
145
146 // ------------------------------------
147 if (paint)
148 {
149 fFunc->SetRange(0, 90);
150 fFunc->SetLineColor(kRed);
151 fFunc->SetLineWidth(2);
152 fFunc->Paint("same");
153 }
154 // ------------------------------------
155
156 fFunc->ReleaseParameter(0);
157 //func.ReleaseParameter(1);
158 fFunc->ReleaseParameter(2);
159 fFunc->FixParameter(3, fFunc->GetParameter(3));
160 for (int i=5; i<fFunc->GetNpar(); i++)
161 fFunc->FixParameter(i, fFunc->GetParameter(i));
162
163 // Do not allow signals smaller than the background
164 const Double_t A = alpha0-fFunc->GetParameter(3);
165 const Double_t dA = TMath::Abs(A);
166 fFunc->SetParLimits(0, -dA*4, dA*4);
167 fFunc->SetParLimits(2, 0, 90);
168
169 // Now fit a gaus in the on region on top of the polynom
170 fFunc->SetParameter(0, A);
171 fFunc->SetParameter(2, sigmax*0.75);
172
173 // options : N do not store the function, do not draw
174 // I use integral of function in bin rather than value at bin center
175 // R use the range specified in the function range
176 // Q quiet mode
177 h.Fit(fFunc, "NQI", "", 0, sigmax);
178
179 fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF();
180 fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters());
181
182 //const Bool_t ok = NDF>0 && chi2<2.5*NDF;
183
184 // ------------------------------------
185 if (paint)
186 {
187 fFunc->SetLineColor(kGreen);
188 fFunc->SetLineWidth(2);
189 fFunc->Paint("same");
190 }
191 // ------------------------------------
192
193 //const Double_t s = fFunc->Integral(0, fSigInt)/alphaw;
194 fFunc->SetParameter(0, 0);
195 fFunc->SetParameter(2, 1);
196 //const Double_t b = fFunc->Integral(0, fSigInt)/alphaw;
197 //fSignificance = MMath::SignificanceLiMaSigned(s, b);
198
199 const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999);
200
201 fIntegralMax = h.GetBinLowEdge(bin+1);
202 fEventsBackground = fFunc->Integral(0, fIntegralMax)/alphaw;
203 fEventsSignal = h.Integral(0, bin);
204 fEventsExcess = fEventsSignal-fEventsBackground;
205 fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
206
207 if (fEventsExcess<0)
208 fEventsExcess=0;
209
210 return kTRUE;
211}
212
213Bool_t MAlphaFitter::Fit(TH1D &hon, TH1D &hof, Double_t alpha, Bool_t paint)
214{
215 /*
216 Clear();
217 if (hon.GetEntries()==0)
218 return kFALSE;
219 */
220
221 TH1D h(hon);
222 h.Add(&hof, -1);
223
224 MAlphaFitter fit(*this);
225 fit.SetPolynomOrder(1);
226
227 if (alpha<=0 || !fit.Fit(h, paint))
228 return kFALSE;
229
230 fChiSqSignal = fit.GetChiSqSignal();
231 fChiSqBg = fit.GetChiSqBg();
232 fCoefficients = fit.GetCoefficients();
233
234 const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999);
235
236 fIntegralMax = hon.GetBinLowEdge(bin+1);
237 fEventsBackground = hof.Integral(0, bin);
238 fEventsSignal = hon.Integral(0, bin);
239 fEventsExcess = fEventsSignal-fEventsBackground;
240 fScaleFactor = alpha;
241 fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground/alpha, alpha);
242
243 if (fEventsExcess<0)
244 fEventsExcess=0;
245/*
246 TF1 func("", "gaus(0)+pol0(3)", 0., 90.);
247
248 const Double_t A = fEventsSignal/bin;
249 const Double_t dA = TMath::Abs(A);
250 func.SetParLimits(0, -dA*4, dA*4);
251 func.SetParLimits(2, 0, 90);
252 func.SetParLimits(3, -dA, dA);
253
254 func.SetParameter(0, A);
255 func.FixParameter(1, 0);
256 func.SetParameter(2, fSigMax*0.75);
257 func.SetParameter(3, 0);
258
259 // options : N do not store the function, do not draw
260 // I use integral of function in bin rather than value at bin center
261 // R use the range specified in the function range
262 // Q quiet mode
263 TH1D h(hon);
264 h.Add(&hof, -1);
265 h.Fit(&func, "NQI", "", 0, 90);
266
267 fChiSqSignal = func.GetChisquare()/func.GetNDF();
268
269 const Int_t bin1 = h.GetXaxis()->FindFixBin(func.GetParameter(2)*2);
270
271 fChiSqBg = 0;
272 for (int i=bin1; i<=h.GetNbinsX(); i++)
273 {
274 const Float_t val = h.GetBinContent(i);
275 fChiSqBg = val*val;
276 }
277 if (fChiSqBg>0)
278 fChiSqBg /= h.GetNbinsX()+1-bin1;
279
280 fCoefficients.Set(func.GetNpar(), func.GetParameters());
281
282 // ------------------------------------
283 if (paint)
284 {
285 func.SetLineColor(kBlue);
286 func.SetLineWidth(2);
287 func.Paint("same");
288 }
289 // ------------------------------------
290 */
291 return kTRUE;
292}
293
294void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const
295{
296 TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f \\omega=%.1f\\circ E=%d (\\alpha<%.1f\\circ) (\\chi_{b}^{2}/ndf=%.1f \\chi_{s}^{2}/ndf=%.1f c_{0}=%.1f)",
297 fSignificance, GetGausSigma(),
298 (int)fEventsExcess, fIntegralMax,
299 fChiSqBg, fChiSqSignal, fCoefficients[3]));
300
301 text.SetBit(TLatex::kTextNDC);
302 text.SetTextSize(size);
303 text.Paint();
304}
305
306void MAlphaFitter::Copy(TObject &o) const
307{
308 MAlphaFitter &f = static_cast<MAlphaFitter&>(o);
309
310 f.fSigInt = fSigInt;
311 f.fSigMax = fSigMax;
312 f.fBgMin = fBgMin;
313 f.fBgMax = fBgMax;
314 f.fScaleMin = fScaleMin;
315 f.fScaleMax = fScaleMax;
316 f.fPolynomOrder = fPolynomOrder;
317 f.fScaleMode = fScaleMode;
318 f.fScaleUser = fScaleUser;
319 f.fCoefficients.Set(fCoefficients.GetSize());
320 f.fCoefficients.Reset();
321
322 TF1 *fcn = f.fFunc;
323 f.fFunc = new TF1(*fFunc);
324 gROOT->GetListOfFunctions()->Remove(f.fFunc);
325 f.fFunc->SetName("Dummy");
326 delete fcn;
327}
328
329void MAlphaFitter::Print(Option_t *o) const
330{
331 *fLog << inf << GetDescriptor() << ": Fitting..." << endl;
332 *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl;
333 *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl;
334 *fLog << " ...polynom order " << fPolynomOrder << endl;
335 *fLog << " ...scale mode: ";
336 switch (fScaleMode)
337 {
338 case kNone: *fLog << "none."; break;
339 case kEntries: *fLog << "entries."; break;
340 case kIntegral: *fLog << "integral."; break;
341 case kOffRegion: *fLog << "off region."; break;
342 case kLeastSquare: *fLog << "least square."; break;
343 case kUserScale: *fLog << "user def (" << fScaleUser << ")"; break;
344 }
345 *fLog << endl;
346
347 if (TString(o).Contains("result"))
348 {
349 *fLog << "Result:" << endl;
350 *fLog << " - Significance " << 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 " << fIntegralMax << "°" << endl;
357 }
358}
359
360Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, UInt_t bin, Bool_t paint)
361{
362 const TString name(Form("TempAlphaEnergy%06d", gRandom->Integer(1000000)));
363 TH1D *h = hon.ProjectionZ(name, -1, 9999, bin, bin, "E");
364 h->SetDirectory(0);
365
366 const Bool_t rc = Fit(*h, paint);
367 delete h;
368 return rc;
369}
370
371Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
372{
373 const TString name(Form("TempAlphaTheta%06d", gRandom->Integer(1000000)));
374 TH1D *h = hon.ProjectionZ(name, bin, bin, -1, 9999, "E");
375 h->SetDirectory(0);
376
377 const Bool_t rc = Fit(*h, paint);
378 delete h;
379 return rc;
380}
381
382Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint)
383{
384 const TString name(Form("TempAlpha%06d", gRandom->Integer(1000000)));
385 TH1D *h = hon.ProjectionZ(name, -1, 9999, -1, 9999, "E");
386 h->SetDirectory(0);
387
388 const Bool_t rc = Fit(*h, paint);
389 delete h;
390 return rc;
391}
392
393Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
394{
395 const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000)));
396 const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
397
398 TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E");
399 TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E");
400 h1->SetDirectory(0);
401 h0->SetDirectory(0);
402
403 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
404
405 delete h0;
406 delete h1;
407
408 return rc;
409}
410
411Bool_t MAlphaFitter::FitTheta(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
412{
413 const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000)));
414 const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
415
416 TH1D *h1 = hon.ProjectionZ(name1, bin, bin, -1, 9999, "E");
417 TH1D *h0 = hof.ProjectionZ(name0, bin, bin, -1, 9999, "E");
418 h1->SetDirectory(0);
419 h0->SetDirectory(0);
420
421 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
422
423 delete h0;
424 delete h1;
425
426 return rc;
427}
428
429Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint)
430{
431 const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000)));
432 const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
433
434 TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, -1, 9999, "E");
435 TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, -1, 9999, "E");
436 h1->SetDirectory(0);
437 h0->SetDirectory(0);
438
439 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
440
441 delete h0;
442 delete h1;
443
444 return rc;
445}
446
447Double_t MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
448{
449 Float_t scaleon = 1;
450 Float_t scaleof = 1;
451 switch (fScaleMode)
452 {
453 case kNone:
454 return 1;
455
456 case kEntries:
457 scaleon = on.GetEntries();
458 scaleof = of.GetEntries();
459 break;
460
461 case kIntegral:
462 scaleon = on.Integral();
463 scaleof = of.Integral();
464 break;
465
466 case kOffRegion:
467 {
468 const Int_t min = on.GetXaxis()->FindFixBin(fScaleMin);
469 const Int_t max = on.GetXaxis()->FindFixBin(fScaleMax);
470 scaleon = on.Integral(min, max);
471 scaleof = of.Integral(min, max);
472 }
473 break;
474
475 case kUserScale:
476 scaleon = fScaleUser;
477 break;
478
479 // This is just to make some compiler happy
480 default:
481 return 1;
482 }
483
484 if (scaleof!=0)
485 {
486 of.Scale(scaleon/scaleof);
487 return scaleon/scaleof;
488 }
489 else
490 {
491 of.Reset();
492 return 0;
493 }
494}
495
496Int_t MAlphaFitter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
497{
498 Bool_t rc = kFALSE;
499
500 //void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; }
501 //void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; }
502
503 if (IsEnvDefined(env, prefix, "SignalIntegralMax", print))
504 {
505 SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalIntegralMax", fSigInt));
506 rc = kTRUE;
507 }
508 if (IsEnvDefined(env, prefix, "SignalFitMax", print))
509 {
510 SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalFitMax", fSigMax));
511 rc = kTRUE;
512 }
513 if (IsEnvDefined(env, prefix, "BackgroundFitMax", print))
514 {
515 SetBackgroundFitMax(GetEnvValue(env, prefix, "BackgroundFitMax", fBgMax));
516 rc = kTRUE;
517 }
518 if (IsEnvDefined(env, prefix, "BackgroundFitMin", print))
519 {
520 SetBackgroundFitMin(GetEnvValue(env, prefix, "BackgroundFitMin", fBgMin));
521 rc = kTRUE;
522 }
523 if (IsEnvDefined(env, prefix, "ScaleMin", print))
524 {
525 SetScaleMin(GetEnvValue(env, prefix, "ScaleMin", fScaleMin));
526 rc = kTRUE;
527 }
528 if (IsEnvDefined(env, prefix, "ScaleMax", print))
529 {
530 SetScaleMax(GetEnvValue(env, prefix, "ScaleMax", fScaleMax));
531 rc = kTRUE;
532 }
533 if (IsEnvDefined(env, prefix, "PolynomOrder", print))
534 {
535 SetPolynomOrder(GetEnvValue(env, prefix, "PolynomOrder", fPolynomOrder));
536 rc = kTRUE;
537 }
538
539 return rc;
540}
Note: See TracBrowser for help on using the repository browser.