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 |
|
---|
16 | class TH1D;
|
---|
17 | class TH3D;
|
---|
18 |
|
---|
19 | class MAlphaFitter : public MParContainer
|
---|
20 | {
|
---|
21 | public:
|
---|
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 | };
|
---|
39 | enum SignalFunc_t {
|
---|
40 | kGauss, kThetaSq
|
---|
41 | };
|
---|
42 |
|
---|
43 | private:
|
---|
44 | // Fitting Setup
|
---|
45 | Float_t fSigInt; // minimum of range to fit the signal
|
---|
46 | Float_t fSigMax; // maximum of range to fit the signal
|
---|
47 | Float_t fBgMin; // minimum of range to fit the background
|
---|
48 | Float_t fBgMax; // minimum of range to fit the background
|
---|
49 | Float_t fScaleMin; // minimum of range to determin the scale factor of the background
|
---|
50 | Float_t fScaleMax; // maximum of range to determin the scale factor of the background
|
---|
51 | Int_t fPolynomOrder; // order of polyom to be fitted to the background
|
---|
52 | Bool_t fFitBackground; // Backround fit: yes/no
|
---|
53 | SignalFunc_t fSignalFunc; // Type of signal function
|
---|
54 | // Result
|
---|
55 | Double_t fSignificance; // significance of signal
|
---|
56 | Double_t fEventsExcess; // calculated number of excess events (signal-bg)
|
---|
57 | Double_t fEventsSignal; // calculated number of signal events
|
---|
58 | Double_t fEventsBackground; // calculated number of bg events (fScaleFactor already applied)
|
---|
59 |
|
---|
60 | Double_t fChiSqSignal; // Reduced (chi^2/NDF) chisq of signal fit
|
---|
61 | Double_t fChiSqBg; // Reduced (chi^2/NDF) chisq of bg fit
|
---|
62 | Double_t fIntegralMax; // Calculated bin border to which it was integrated
|
---|
63 | Double_t fScaleFactor; // Scale factor determined for off-data
|
---|
64 |
|
---|
65 | TArrayD fCoefficients; // Fit result
|
---|
66 |
|
---|
67 | // Function
|
---|
68 | TF1 *fFunc; // fit function (gauss + polynom)
|
---|
69 |
|
---|
70 | // Scaling setup
|
---|
71 | ScaleMode_t fScaleMode; // scaling mode
|
---|
72 | Double_t fScaleUser; // user scale factor
|
---|
73 |
|
---|
74 | // Minimization strategy
|
---|
75 | Strategy_t fStrategy; // How to calc minimization value
|
---|
76 |
|
---|
77 | public:
|
---|
78 | // Implementing the function yourself is only about 5% faster
|
---|
79 | MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15),
|
---|
80 | fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80),
|
---|
81 | fPolynomOrder(2), fFitBackground(kTRUE), fSignalFunc(kGauss),
|
---|
82 | fCoefficients(3+fPolynomOrder+1),
|
---|
83 | fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)),
|
---|
84 | fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance)
|
---|
85 | {
|
---|
86 | fName = name ? name : "MAlphaFitter";
|
---|
87 | fTitle = title ? title : "Fit alpha";
|
---|
88 |
|
---|
89 | gROOT->GetListOfFunctions()->Remove(fFunc);
|
---|
90 | fFunc->SetName("Dummy");
|
---|
91 |
|
---|
92 | Clear();
|
---|
93 | }
|
---|
94 |
|
---|
95 | MAlphaFitter(const MAlphaFitter &f) : fFunc(0)
|
---|
96 | {
|
---|
97 | f.Copy(*this);
|
---|
98 | }
|
---|
99 | ~MAlphaFitter()
|
---|
100 | {
|
---|
101 | delete fFunc;
|
---|
102 | }
|
---|
103 |
|
---|
104 | // TObject
|
---|
105 | void Clear(Option_t *o="");
|
---|
106 | void Print(Option_t *o="") const; //*MENU*
|
---|
107 | void Copy(TObject &o) const;
|
---|
108 |
|
---|
109 | // Setter
|
---|
110 | void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; }
|
---|
111 | void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; }
|
---|
112 | void SetMinimizationStrategy(Strategy_t mode) { fStrategy = mode; }
|
---|
113 | void SetSignalIntegralMax(Float_t s) { fSigInt = s; }
|
---|
114 | void SetSignalFitMax(Float_t s) { fSigMax = s; }
|
---|
115 | void SetBackgroundFitMin(Float_t s) { fBgMin = s; }
|
---|
116 | void SetBackgroundFitMax(Float_t s) { fBgMax = s; }
|
---|
117 | void SetScaleMin(Float_t s) { fScaleMin = s; }
|
---|
118 | void SetScaleMax(Float_t s) { fScaleMax = s; }
|
---|
119 | void SetPolynomOrder(Int_t s)
|
---|
120 | {
|
---|
121 | if (s==fPolynomOrder)
|
---|
122 | return;
|
---|
123 |
|
---|
124 | fPolynomOrder = s;
|
---|
125 |
|
---|
126 | SetSignalFunction(fSignalFunc);
|
---|
127 | }
|
---|
128 | void SetSignalFunction(SignalFunc_t func)
|
---|
129 | {
|
---|
130 | delete fFunc;
|
---|
131 | switch (func)
|
---|
132 | {
|
---|
133 | case kGauss:
|
---|
134 | fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", fPolynomOrder));
|
---|
135 | break;
|
---|
136 | case kThetaSq:
|
---|
137 | fFunc=new TF1 ("", Form("[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + pol%d(3)", fPolynomOrder));
|
---|
138 | break;
|
---|
139 | }
|
---|
140 | fSignalFunc=func;
|
---|
141 | gROOT->GetListOfFunctions()->Remove(fFunc);
|
---|
142 | fFunc->SetName("Dummy");
|
---|
143 | fCoefficients.Set(3+fPolynomOrder+1);
|
---|
144 | fCoefficients.Reset();
|
---|
145 | }
|
---|
146 | void EnableBackgroundFit(Bool_t b=kTRUE) { fFitBackground=b; }
|
---|
147 |
|
---|
148 | // Getter
|
---|
149 | Double_t GetSignalIntegralMax() const { return fSigInt; }
|
---|
150 |
|
---|
151 | Double_t GetEventsExcess() const { return fEventsExcess; }
|
---|
152 | Double_t GetEventsSignal() const { return fEventsSignal; }
|
---|
153 | Double_t GetEventsBackground() const { return fEventsBackground; }
|
---|
154 |
|
---|
155 | Double_t GetSignificance() const { return fSignificance; }
|
---|
156 | Double_t GetChiSqSignal() const { return fChiSqSignal; }
|
---|
157 | Double_t GetChiSqBg() const { return fChiSqBg; }
|
---|
158 | Double_t GetScaleFactor() const { return fScaleFactor; }
|
---|
159 | Double_t GetMinimizationValue() const;
|
---|
160 |
|
---|
161 | Double_t GetGausSigma() const { return fCoefficients[2]; }
|
---|
162 | Double_t GetGausMu() const { return fCoefficients[1]; }
|
---|
163 | Double_t GetGausA() const { return fCoefficients[0]; }
|
---|
164 | Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; }
|
---|
165 | const TArrayD &GetCoefficients() const { return fCoefficients; }
|
---|
166 | Double_t Eval(Double_t d) const { return fFunc ? fFunc->Eval(d) : 0; }
|
---|
167 |
|
---|
168 | // Interface to fit
|
---|
169 | Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
|
---|
170 | Bool_t Fit(const TH1D &on, const TH1D &off, Double_t alpha, Bool_t paint=kFALSE);
|
---|
171 | Bool_t Fit(TH1D &on, TH1D *off, Double_t alpha, Bool_t paint=kFALSE)
|
---|
172 | {
|
---|
173 | return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
|
---|
174 | }
|
---|
175 | Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
|
---|
176 | {
|
---|
177 | return off ? Fit(on, *off, 1, paint) : Fit(on, paint);
|
---|
178 | }
|
---|
179 | Bool_t ScaleAndFit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
|
---|
180 | {
|
---|
181 | const Double_t alpha = off ? Scale(*off, on) : 1;
|
---|
182 | return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
|
---|
183 | }
|
---|
184 |
|
---|
185 | Bool_t FitAlpha(const TH3D &h, Bool_t paint=kFALSE);
|
---|
186 | Bool_t FitEnergy(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
|
---|
187 | Bool_t FitTheta(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
|
---|
188 | //Bool_t FitTime(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
|
---|
189 |
|
---|
190 | Bool_t FitAlpha(const TH3D &on, const TH3D &off, Bool_t paint=kFALSE);
|
---|
191 | Bool_t FitEnergy(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
|
---|
192 | Bool_t FitTheta(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
|
---|
193 | //Bool_t FitTime(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
|
---|
194 |
|
---|
195 | Bool_t FitAlpha(const TH3D &on, const TH3D *off, Bool_t paint=kFALSE)
|
---|
196 | {
|
---|
197 | return off ? FitAlpha(on, *off, paint) : FitAlpha(on, paint);
|
---|
198 | }
|
---|
199 | Bool_t FitEnergy(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
|
---|
200 | {
|
---|
201 | return off ? FitEnergy(on, *off, bin, paint) : FitEnergy(on, bin, paint);
|
---|
202 | }
|
---|
203 | Bool_t FitTheta(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
|
---|
204 | {
|
---|
205 | return off ? FitTheta(on, *off, bin, paint) : FitTheta(on, bin, paint);
|
---|
206 | }/*
|
---|
207 | Bool_t FitTime(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
|
---|
208 | {
|
---|
209 | return off ? FitTime(on, *off, bin, paint) : FitTime(on, bin, paint);
|
---|
210 | }*/
|
---|
211 |
|
---|
212 | Double_t Scale(TH1D &off, const TH1D &on) const;
|
---|
213 |
|
---|
214 | // Interface to result
|
---|
215 | void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const;
|
---|
216 |
|
---|
217 | // MTask
|
---|
218 | Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
|
---|
219 |
|
---|
220 | ClassDef(MAlphaFitter, 1)
|
---|
221 | };
|
---|
222 |
|
---|
223 | #endif
|
---|