source: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h@ 9029

Last change on this file since 9029 was 8989, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 9.4 KB
Line 
1#ifndef MARS_MAlphaFitter
2#define MARS_MAlphaFitter
3
4#ifndef MARS_MParContainer
5#include "MParContainer.h"
6#endif
7
8#ifndef ROOT_TArrayD
9#include <TArrayD.h>
10#endif
11
12#ifndef ROOT_TF1
13#include <TF1.h>
14#endif
15
16class TH1D;
17class TH3D;
18
19class MAlphaFitter : public MParContainer
20{
21public:
22 enum ScaleMode_t {
23 kNone, // No scaling
24 kEntries, // scale by the number of entries in on and off
25 kIntegral, // scale by the integral in on and off
26 kOffRegion, // scale by the integral between fScaleMin, fScaleMax in on and off
27 kBackground, // scale by the integral between fBgMin, fBgMax in on and off
28 kLeastSquare, // not yet implemented
29 kUserScale // scale by fixed factor set by SetScaleUser
30 };
31 enum Strategy_t {
32 kSignificance,
33 kSignificanceChi2,
34 kSignificanceLogExcess,
35 kSignificanceExcess,
36 kExcess,
37 kGaussSigma,
38 kWeakSource,
39 kWeakSourceLogExcess
40 };
41 enum SignalFunc_t {
42 kGauss, kThetaSq
43 };
44
45private:
46 // Fitting Setup
47 Float_t fSigInt; // minimum of range to fit the signal
48 Float_t fSigMax; // maximum of range to fit the signal
49 Float_t fBgMin; // minimum of range to fit the background
50 Float_t fBgMax; // minimum of range to fit the background
51 Float_t fScaleMin; // minimum of range to determin the scale factor of the background
52 Float_t fScaleMax; // maximum of range to determin the scale factor of the background
53 Int_t fPolynomOrder; // order of polyom to be fitted to the background
54 Bool_t fFitBackground; // Backround fit: yes/no
55 SignalFunc_t fSignalFunc; // Type of signal function
56 // Result
57 Double_t fSignificance; // significance of an unknown signal (Li/Ma 17)
58 Double_t fErrorExcess; // Simple error propagation
59 Double_t fEventsExcess; // calculated number of excess events (signal-bg)
60 Double_t fEventsSignal; // calculated number of signal events
61 Double_t fEventsBackground; // calculated number of bg events (fScaleFactor already applied)
62
63 Double_t fChiSqSignal; // Reduced (chi^2/NDF) chisq of signal fit
64 Double_t fChiSqBg; // Reduced (chi^2/NDF) chisq of bg fit
65 Double_t fIntegralMax; // Calculated bin border to which it was integrated
66 Double_t fScaleFactor; // Scale factor determined for off-data
67
68 TArrayD fCoefficients; // Fit result
69 TArrayD fErrors; // Fit errors
70
71 // Function
72 TF1 *fFunc; // fit function (gauss + polynom)
73
74 // Scaling setup
75 ScaleMode_t fScaleMode; // scaling mode
76 Double_t fScaleUser; // user scale factor
77
78 // Minimization strategy
79 Strategy_t fStrategy; // How to calc minimization value
80
81 Double_t DoOffFit(const TH1D &hon, const TH1D &hof, Bool_t paint);
82 Bool_t FitOff(TH1D &h, Int_t paint);
83 void FitResult(const TH1D &h);
84
85public:
86 // Implementing the function yourself is only about 5% faster
87 MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15),
88 fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80),
89 fPolynomOrder(2), fFitBackground(kTRUE), fSignalFunc(kGauss),
90 fCoefficients(3+fPolynomOrder+1), fErrors(3+fPolynomOrder+1),
91 fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)),
92 fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance)
93 {
94 fName = name ? name : "MAlphaFitter";
95 fTitle = title ? title : "Fit alpha";
96
97 fFunc->SetName("Dummy");
98 gROOT->GetListOfFunctions()->Remove(fFunc);
99
100 Clear();
101 }
102
103 MAlphaFitter(const MAlphaFitter &f) : fFunc(0)
104 {
105 f.Copy(*this);
106 }
107 ~MAlphaFitter()
108 {
109 delete fFunc;
110 }
111
112 // TObject
113 void Clear(Option_t *o="");
114 void Print(Option_t *o="") const; //*MENU*
115 void Copy(TObject &o) const;
116
117 // Setter
118 void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; }
119 void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; }
120 void SetMinimizationStrategy(Strategy_t mode) { fStrategy = mode; }
121 void SetSignalIntegralMax(Float_t s) { fSigInt = s; }
122 void SetSignalFitMax(Float_t s) { fSigMax = s; }
123 void SetBackgroundFitMin(Float_t s) { fBgMin = s; }
124 void SetBackgroundFitMax(Float_t s) { fBgMax = s; }
125 void SetScaleMin(Float_t s) { fScaleMin = s; }
126 void SetScaleMax(Float_t s) { fScaleMax = s; }
127 void SetPolynomOrder(Int_t s)
128 {
129 if (s==fPolynomOrder)
130 return;
131
132 fPolynomOrder = s;
133
134 SetSignalFunction(fSignalFunc);
135 }
136 void SetSignalFunction(SignalFunc_t func)
137 {
138 delete fFunc;
139 switch (func)
140 {
141 case kGauss:
142 fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", fPolynomOrder));
143 break;
144 case kThetaSq:
145// if (fPolynomOrder==0)
146// fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + pol0(3)");
147// else
148 // {
149 if (fPolynomOrder>0)
150 fPolynomOrder = 1;
151 fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + expo(3)");
152// }
153 break;
154 }
155 fSignalFunc=func;
156 fFunc->SetName("Dummy");
157 gROOT->GetListOfFunctions()->Remove(fFunc);
158
159 fCoefficients.Set(3+fPolynomOrder+1);
160 fCoefficients.Reset();
161
162 fErrors.Set(3+fPolynomOrder+1);
163 fErrors.Reset();
164 }
165 void EnableBackgroundFit(Bool_t b=kTRUE) { fFitBackground=b; }
166
167 void FixScale()
168 {
169 fScaleMode = kUserScale;
170 fScaleUser = fScaleFactor;
171 }
172
173 // Getter
174 Double_t GetSignalIntegralMax() const { return fSigInt; }
175
176 Double_t GetEventsExcess() const { return fEventsExcess; }
177 Double_t GetEventsSignal() const { return fEventsSignal; }
178 Double_t GetEventsBackground() const { return fEventsBackground; }
179
180 Double_t GetSignificance() const { return fSignificance; }
181 Double_t GetErrorExcess() const { return fErrorExcess; }
182 Double_t GetChiSqSignal() const { return fChiSqSignal; }
183 Double_t GetChiSqBg() const { return fChiSqBg; }
184 Double_t GetScaleFactor() const { return fScaleFactor; }
185 Double_t GetMinimizationValue() const;
186
187 ScaleMode_t GetScaleMode() const { return fScaleMode; }
188
189 Double_t GetGausSigma() const { return fCoefficients[2]; }
190 Double_t GetGausMu() const { return fCoefficients[1]; }
191 Double_t GetGausA() const { return fCoefficients[0]; }
192 Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; }
193 const TArrayD &GetCoefficients() const { return fCoefficients; }
194 const TArrayD &GetErrors() const { return fErrors; }
195 Double_t Eval(Double_t d) const { return fFunc ? fFunc->Eval(d) : 0; }
196
197 Double_t CalcUpperLimit() const;
198
199 // Interface to fit
200 Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
201 Bool_t Fit(const TH1D &on, const TH1D &off, Double_t alpha, Bool_t paint=kFALSE);
202 Bool_t Fit(TH1D &on, TH1D *off, Double_t alpha, Bool_t paint=kFALSE)
203 {
204 return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
205 }
206 Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
207 {
208 return off ? Fit(on, *off, 1, paint) : Fit(on, paint);
209 }
210 Bool_t ScaleAndFit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
211 {
212 const Double_t alpha = off ? Scale(*off, on) : 1;
213 return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
214 }
215
216 Bool_t FitAlpha(const TH3D &h, Bool_t paint=kFALSE);
217 Bool_t FitEnergy(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
218 Bool_t FitTheta(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
219 //Bool_t FitTime(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
220
221 Bool_t FitAlpha(const TH3D &on, const TH3D &off, Bool_t paint=kFALSE);
222 Bool_t FitEnergy(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
223 Bool_t FitTheta(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
224 //Bool_t FitTime(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
225
226 Bool_t FitAlpha(const TH3D &on, const TH3D *off, Bool_t paint=kFALSE)
227 {
228 return off ? FitAlpha(on, *off, paint) : FitAlpha(on, paint);
229 }
230 Bool_t FitEnergy(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
231 {
232 return off ? FitEnergy(on, *off, bin, paint) : FitEnergy(on, bin, paint);
233 }
234 Bool_t FitTheta(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
235 {
236 return off ? FitTheta(on, *off, bin, paint) : FitTheta(on, bin, paint);
237 }/*
238 Bool_t FitTime(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
239 {
240 return off ? FitTime(on, *off, bin, paint) : FitTime(on, bin, paint);
241 }*/
242
243 Double_t Scale(TH1D &off, const TH1D &on) const;
244
245 Bool_t ApplyScaling(const TH3D &hon, TH3D &hof, UInt_t bin) const;
246 Bool_t ApplyScaling(const TH3D &hon, TH3D &hof) const;
247
248 // Interface to result
249 void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035, Bool_t draw=kFALSE) const;
250 void DrawResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const { PaintResult(x, y, size, kTRUE); }
251
252 // MTask
253 Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
254
255 ClassDef(MAlphaFitter, 4)
256};
257
258#endif
Note: See TracBrowser for help on using the repository browser.