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

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