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

Last change on this file since 7008 was 7005, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 18.6 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 // E Perform better Errors estimation using Minos technique
143 h.Fit(fFunc, "NQI", "", bgmin, bgmax);
144
145 fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();
146
147 // ------------------------------------
148 if (paint)
149 {
150 fFunc->SetRange(0, 90);
151 fFunc->SetLineColor(kRed);
152 fFunc->SetLineWidth(2);
153 fFunc->Paint("same");
154 }
155 // ------------------------------------
156
157 fFunc->ReleaseParameter(0);
158 //func.ReleaseParameter(1);
159 fFunc->ReleaseParameter(2);
160 fFunc->FixParameter(3, fFunc->GetParameter(3));
161 for (int i=5; i<fFunc->GetNpar(); i++)
162 fFunc->FixParameter(i, fFunc->GetParameter(i));
163
164 // Do not allow signals smaller than the background
165 const Double_t A = alpha0-fFunc->GetParameter(3);
166 const Double_t dA = TMath::Abs(A);
167 fFunc->SetParLimits(0, -dA*4, dA*4);
168 fFunc->SetParLimits(2, 0, 90);
169
170 // Now fit a gaus in the on region on top of the polynom
171 fFunc->SetParameter(0, A);
172 fFunc->SetParameter(2, sigmax*0.75);
173
174 // options : N do not store the function, do not draw
175 // I use integral of function in bin rather than value at bin center
176 // R use the range specified in the function range
177 // Q quiet mode
178 // E Perform better Errors estimation using Minos technique
179 h.Fit(fFunc, "NQI", "", 0, sigmax);
180
181 fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF();
182 fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters());
183
184 //const Bool_t ok = NDF>0 && chi2<2.5*NDF;
185
186 // ------------------------------------
187 if (paint)
188 {
189 fFunc->SetLineColor(kGreen);
190 fFunc->SetLineWidth(2);
191 fFunc->Paint("same");
192 }
193 // ------------------------------------
194
195 //const Double_t s = fFunc->Integral(0, fSigInt)/alphaw;
196 fFunc->SetParameter(0, 0);
197 fFunc->SetParameter(2, 1);
198 //const Double_t b = fFunc->Integral(0, fSigInt)/alphaw;
199 //fSignificance = MMath::SignificanceLiMaSigned(s, b);
200
201 const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999);
202
203 fIntegralMax = h.GetBinLowEdge(bin+1);
204 fEventsBackground = fFunc->Integral(0, fIntegralMax)/alphaw;
205 fEventsSignal = h.Integral(0, bin);
206 fEventsExcess = fEventsSignal-fEventsBackground;
207 fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
208
209 if (fEventsExcess<0)
210 fEventsExcess=0;
211
212 return kTRUE;
213}
214
215Bool_t MAlphaFitter::Fit(const TH1D &hon, const TH1D &hof, Double_t alpha, Bool_t paint)
216{
217 TH1D h(hon);
218 h.Add(&hof, -1); // substracts also number of entries!
219 h.SetEntries(hon.GetEntries());
220
221 MAlphaFitter fit(*this);
222 fit.SetPolynomOrder(1);
223
224 if (alpha<=0 || !fit.Fit(h, paint))
225 return kFALSE;
226
227 fChiSqSignal = fit.GetChiSqSignal();
228 fChiSqBg = fit.GetChiSqBg();
229 fCoefficients = fit.GetCoefficients();
230
231 const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999);
232
233 fIntegralMax = hon.GetBinLowEdge(bin+1);
234 fEventsBackground = hof.Integral(0, bin);
235 fEventsSignal = hon.Integral(0, bin);
236 fEventsExcess = fEventsSignal-fEventsBackground;
237 fScaleFactor = alpha;
238 fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground/alpha, alpha);
239
240 if (fEventsExcess<0)
241 fEventsExcess=0;
242/*
243 TF1 func("", "gaus(0)+pol0(3)", 0., 90.);
244
245 const Double_t A = fEventsSignal/bin;
246 const Double_t dA = TMath::Abs(A);
247 func.SetParLimits(0, -dA*4, dA*4);
248 func.SetParLimits(2, 0, 90);
249 func.SetParLimits(3, -dA, dA);
250
251 func.SetParameter(0, A);
252 func.FixParameter(1, 0);
253 func.SetParameter(2, fSigMax*0.75);
254 func.SetParameter(3, 0);
255
256 // options : N do not store the function, do not draw
257 // I use integral of function in bin rather than value at bin center
258 // R use the range specified in the function range
259 // Q quiet mode
260 // E Perform better Errors estimation using Minos technique
261 TH1D h(hon);
262 h.Add(&hof, -1);
263 h.Fit(&func, "NQI", "", 0, 90);
264
265 fChiSqSignal = func.GetChisquare()/func.GetNDF();
266
267 const Int_t bin1 = h.GetXaxis()->FindFixBin(func.GetParameter(2)*2);
268
269 fChiSqBg = 0;
270 for (int i=bin1; i<=h.GetNbinsX(); i++)
271 {
272 const Float_t val = h.GetBinContent(i);
273 fChiSqBg = val*val;
274 }
275 if (fChiSqBg>0)
276 fChiSqBg /= h.GetNbinsX()+1-bin1;
277
278 fCoefficients.Set(func.GetNpar(), func.GetParameters());
279
280 // ------------------------------------
281 if (paint)
282 {
283 func.SetLineColor(kBlue);
284 func.SetLineWidth(2);
285 func.Paint("same");
286 }
287 // ------------------------------------
288 */
289 return kTRUE;
290}
291
292void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const
293{
294 const Double_t w = GetGausSigma();
295 const Int_t l = w<=0 ? 0 : (Int_t)TMath::Ceil(-TMath::Log10(w));
296 const TString fmt = Form("\\sigma_{Li/Ma}=%%.1f \\omega=%%.%df\\circ E=%%d (\\alpha<%%.1f\\circ) (\\chi_{b}^{2}/ndf=%%.1f \\chi_{s}^{2}/ndf=%%.1f c_{0}=%%.1f)",
297 l<1?1:l);
298
299 TLatex text(x, y, Form(fmt.Data(), fSignificance, w, (int)fEventsExcess,
300 fIntegralMax, fChiSqBg, fChiSqSignal,
301 fCoefficients[3]));
302
303 text.SetBit(TLatex::kTextNDC);
304 text.SetTextSize(size);
305 text.Paint();
306}
307
308void MAlphaFitter::Copy(TObject &o) const
309{
310 MAlphaFitter &f = static_cast<MAlphaFitter&>(o);
311
312 // Setup
313 f.fSigInt = fSigInt;
314 f.fSigMax = fSigMax;
315 f.fBgMin = fBgMin;
316 f.fBgMax = fBgMax;
317 f.fScaleMin = fScaleMin;
318 f.fScaleMax = fScaleMax;
319 f.fPolynomOrder = fPolynomOrder;
320 f.fScaleMode = fScaleMode;
321 f.fScaleUser = fScaleUser;
322 f.fStrategy = fStrategy;
323 f.fCoefficients.Set(fCoefficients.GetSize());
324 f.fCoefficients.Reset();
325
326 // Result
327 f.fSignificance = fSignificance;
328 f.fEventsExcess = fEventsExcess;
329 f.fEventsSignal = fEventsSignal;
330 f.fEventsBackground = fEventsBackground;
331 f.fChiSqSignal = fChiSqSignal;
332 f.fChiSqBg = fChiSqBg;
333 f.fIntegralMax = fIntegralMax;
334 f.fScaleFactor = fScaleFactor;
335
336 // Function
337 TF1 *fcn = f.fFunc;
338 f.fFunc = new TF1(*fFunc);
339 gROOT->GetListOfFunctions()->Remove(f.fFunc);
340 f.fFunc->SetName("Dummy");
341 delete fcn;
342}
343
344void MAlphaFitter::Print(Option_t *o) const
345{
346 *fLog << GetDescriptor() << ": Fitting..." << endl;
347 *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl;
348 *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl;
349 *fLog << " ...polynom order " << fPolynomOrder << endl;
350 *fLog << " ...scale mode: ";
351 switch (fScaleMode)
352 {
353 case kNone: *fLog << "none."; break;
354 case kEntries: *fLog << "entries."; break;
355 case kIntegral: *fLog << "integral."; break;
356 case kOffRegion: *fLog << "off region."; break;
357 case kLeastSquare: *fLog << "least square."; break;
358 case kUserScale: *fLog << "user def (" << fScaleUser << ")"; break;
359 }
360 *fLog << endl;
361
362 if (TString(o).Contains("result"))
363 {
364 *fLog << "Result:" << endl;
365 *fLog << " - Significance (Li/Ma) " << fSignificance << endl;
366 *fLog << " - Excess Events " << fEventsExcess << endl;
367 *fLog << " - Signal Events " << fEventsSignal << endl;
368 *fLog << " - Background Events " << fEventsBackground << endl;
369 *fLog << " - Chi^2/ndf (Signal) " << fChiSqSignal << endl;
370 *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl;
371 *fLog << " - Signal integrated up to " << fIntegralMax << "°" << endl;
372 *fLog << " - Scale Factor (Off) " << fScaleFactor << endl;
373 }
374}
375
376Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, UInt_t bin, Bool_t paint)
377{
378 const TString name(Form("TempAlphaEnergy%06d", gRandom->Integer(1000000)));
379 TH1D *h = hon.ProjectionZ(name, -1, 9999, bin, bin, "E");
380 h->SetDirectory(0);
381
382 const Bool_t rc = Fit(*h, paint);
383 delete h;
384 return rc;
385}
386
387Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
388{
389 const TString name(Form("TempAlphaTheta%06d", gRandom->Integer(1000000)));
390 TH1D *h = hon.ProjectionZ(name, bin, bin, -1, 9999, "E");
391 h->SetDirectory(0);
392
393 const Bool_t rc = Fit(*h, paint);
394 delete h;
395 return rc;
396}
397
398Bool_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
409Bool_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
427Bool_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
445Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint)
446{
447 const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000)));
448 const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
449
450 TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, -1, 9999, "E");
451 TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, -1, 9999, "E");
452 h1->SetDirectory(0);
453 h0->SetDirectory(0);
454
455 const Bool_t rc = ScaleAndFit(*h1, h0, paint);
456
457 delete h0;
458 delete h1;
459
460 return rc;
461}
462
463Double_t MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
464{
465 Float_t scaleon = 1;
466 Float_t scaleof = 1;
467 switch (fScaleMode)
468 {
469 case kNone:
470 return 1;
471
472 case kEntries:
473 scaleon = on.GetEntries();
474 scaleof = of.GetEntries();
475 break;
476
477 case kIntegral:
478 scaleon = on.Integral();
479 scaleof = of.Integral();
480 break;
481
482 case kOffRegion:
483 {
484 const Int_t min = on.GetXaxis()->FindFixBin(fScaleMin);
485 const Int_t max = on.GetXaxis()->FindFixBin(fScaleMax);
486 scaleon = on.Integral(min, max);
487 scaleof = of.Integral(min, max);
488 }
489 break;
490
491 case kUserScale:
492 scaleon = fScaleUser;
493 break;
494
495 // This is just to make some compiler happy
496 default:
497 return 1;
498 }
499
500 if (scaleof!=0)
501 {
502 of.Scale(scaleon/scaleof);
503 return scaleon/scaleof;
504 }
505 else
506 {
507 of.Reset();
508 return 0;
509 }
510}
511
512Double_t MAlphaFitter::GetMinimizationValue() const
513{
514 switch (fStrategy)
515 {
516 case kSignificance:
517 return -GetSignificance();
518 case kSignificanceChi2:
519 return -GetSignificance()/GetChiSqSignal();
520 case kSignificanceLogExcess:
521 if (GetEventsExcess()<1)
522 return 0;
523 return -GetSignificance()*TMath::Log10(GetEventsExcess());
524 case kSignificanceExcess:
525 return -GetSignificance()*GetEventsExcess();
526 case kExcess:
527 return -GetEventsExcess();
528 }
529 return 0;
530}
531
532Int_t MAlphaFitter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
533{
534 Bool_t rc = kFALSE;
535
536 //void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; }
537 //void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; }
538
539 if (IsEnvDefined(env, prefix, "SignalIntegralMax", print))
540 {
541 SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalIntegralMax", fSigInt));
542 rc = kTRUE;
543 }
544 if (IsEnvDefined(env, prefix, "SignalFitMax", print))
545 {
546 SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalFitMax", fSigMax));
547 rc = kTRUE;
548 }
549 if (IsEnvDefined(env, prefix, "BackgroundFitMax", print))
550 {
551 SetBackgroundFitMax(GetEnvValue(env, prefix, "BackgroundFitMax", fBgMax));
552 rc = kTRUE;
553 }
554 if (IsEnvDefined(env, prefix, "BackgroundFitMin", print))
555 {
556 SetBackgroundFitMin(GetEnvValue(env, prefix, "BackgroundFitMin", fBgMin));
557 rc = kTRUE;
558 }
559 if (IsEnvDefined(env, prefix, "ScaleMin", print))
560 {
561 SetScaleMin(GetEnvValue(env, prefix, "ScaleMin", fScaleMin));
562 rc = kTRUE;
563 }
564 if (IsEnvDefined(env, prefix, "ScaleMax", print))
565 {
566 SetScaleMax(GetEnvValue(env, prefix, "ScaleMax", fScaleMax));
567 rc = kTRUE;
568 }
569 if (IsEnvDefined(env, prefix, "PolynomOrder", print))
570 {
571 SetPolynomOrder(GetEnvValue(env, prefix, "PolynomOrder", fPolynomOrder));
572 rc = kTRUE;
573 }
574
575 if (IsEnvDefined(env, prefix, "MinimizationStrategy", print))
576 {
577 TString txt = GetEnvValue(env, prefix, "MinimizationStrategy", "");
578 txt = txt.Strip(TString::kBoth);
579 txt.ToLower();
580 if (txt==(TString)"significance")
581 fStrategy = kSignificance;
582 if (txt==(TString)"significancechi2")
583 fStrategy = kSignificanceChi2;
584 if (txt==(TString)"significanceexcess")
585 fStrategy = kSignificanceExcess;
586 if (txt==(TString)"excess")
587 fStrategy = kExcess;
588 rc = kTRUE;
589 }
590
591 if (IsEnvDefined(env, prefix, "ScaleMode", print))
592 {
593 TString txt = GetEnvValue(env, prefix, "ScaleMode", "");
594 txt = txt.Strip(TString::kBoth);
595 txt.ToLower();
596 if (txt==(TString)"none")
597 fScaleMode = kNone;
598 if (txt==(TString)"entries")
599 fScaleMode = kEntries;
600 if (txt==(TString)"integral")
601 fScaleMode = kIntegral;
602 if (txt==(TString)"offregion")
603 fScaleMode = kOffRegion;
604 if (txt==(TString)"leastsquare")
605 fScaleMode = kLeastSquare;
606 if (txt==(TString)"userscale")
607 fScaleMode = kUserScale;
608 rc = kTRUE;
609 }
610 if (IsEnvDefined(env, prefix, "Scale", print))
611 {
612 fScaleUser = GetEnvValue(env, prefix, "Scale", fScaleUser);
613 rc = kTRUE;
614 }
615
616 return rc;
617}
Note: See TracBrowser for help on using the repository browser.