source: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h@ 2830

Last change on this file since 2830 was 2828, checked in by gaug, 22 years ago
*** empty log message ***
File size: 9.3 KB
Line 
1#ifndef MARS_MHCalibrationBlindPixel
2#define MARS_MHCalibrationBlindPixel
3
4#ifndef MARS_MH
5#include "MH.h"
6#endif
7
8#ifndef MARS_MHCalibrationConfig
9#include "MHCalibrationConfig.h"
10#endif
11
12#ifndef ROOT_TH1
13#include "TH1.h"
14#endif
15
16#ifndef ROOT_TH1F
17#include "TH1F.h"
18#endif
19
20#ifndef ROOT_TF1
21#include "TF1.h"
22#endif
23
24#ifndef ROOT_TPaveText
25#include "TPaveText.h"
26#endif
27
28const Double_t NoWay = 10000000000.0;
29
30class TMath;
31class MParList;
32class MHCalibrationBlindPixel : public MH
33{
34private:
35
36 TH1F* fHBlindPixelCharge; //-> Histogram with the single Phe spectrum
37 TH1I* fHBlindPixelTime; //-> Variance of summed FADC slices
38 TH1I* fHBlindPixelChargevsN; //-> Summed Charge vs. Event Nr.
39
40 TF1 *fSinglePheFit;
41 TF1 *fTimeGausFit;
42 TF1 *fSinglePhePedFit;
43
44 Axis_t fBlindPixelChargefirst;
45 Axis_t fBlindPixelChargelast;
46 Int_t fBlindPixelChargenbins;
47
48 void ResetBin(Int_t i);
49 void DrawLegend();
50
51 TPaveText *fFitLegend;
52 Bool_t fFitOK;
53
54 Double_t fLambda;
55 Double_t fMu0;
56 Double_t fMu1;
57 Double_t fSigma0;
58 Double_t fSigma1;
59
60 Double_t fLambdaErr;
61 Double_t fMu0Err;
62 Double_t fMu1Err;
63 Double_t fSigma0Err;
64 Double_t fSigma1Err;
65
66 Double_t fChisquare;
67 Double_t fProb;
68 Int_t fNdf;
69
70 Double_t fMeanTime;
71 Double_t fMeanTimeErr;
72 Double_t fSigmaTime;
73 Double_t fSigmaTimeErr;
74
75 Double_t fLambdaCheck;
76 Double_t fLambdaCheckErr;
77
78public:
79
80 MHCalibrationBlindPixel(const char *name=NULL, const char *title=NULL);
81 ~MHCalibrationBlindPixel();
82
83 typedef Double_t (*BlindPixelFitFunc)(Double_t *, Double_t *);
84
85 Bool_t FillBlindPixelCharge(Float_t q) { return fHBlindPixelCharge->Fill(q) > -1; }
86 Bool_t FillBlindPixelTime(Int_t t) { return fHBlindPixelTime->Fill(t) > -1; }
87 Bool_t FillBlindPixelChargevsN(Stat_t rq, Int_t t) { return fHBlindPixelChargevsN->Fill(t,rq) > -1; }
88
89 const Double_t GetLambda() const { return fLambda; }
90 const Double_t GetLambdaCheck() const { return fLambdaCheck; }
91 const Double_t GetMu0() const { return fMu0; }
92 const Double_t GetMu1() const { return fMu1; }
93 const Double_t GetSigma0() const { return fSigma0; }
94 const Double_t GetSigma1() const { return fSigma1; }
95
96 const Double_t GetLambdaErr() const { return fLambdaErr; }
97 const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; }
98 const Double_t GetMu0Err() const { return fMu0Err; }
99 const Double_t GetMu1Err() const { return fMu1Err; }
100 const Double_t GetSigma0Err() const { return fSigma0Err; }
101 const Double_t GetSigma1Err() const { return fSigma1Err; }
102
103 const Double_t GetChiSquare() const { return fChisquare; }
104 const Double_t GetProb() const { return fProb; }
105 const Int_t GetNdf() const { return fNdf; }
106
107 const Double_t GetMeanTime() const { return fMeanTime; }
108 const Double_t GetMeanTimeErr() const { return fMeanTimeErr; }
109 const Double_t GetSigmaTime() const { return fSigmaTime; }
110 const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; }
111
112 Bool_t SimulateSinglePhe(Double_t lambda,
113 Double_t mu0,
114 Double_t mu1,
115 Double_t sigma0,
116 Double_t sigma1);
117
118 Bool_t FitSinglePhe(Axis_t rmin=0, Axis_t rmax=0, Option_t *opt="RL0+Q");
119 Bool_t FitTime(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+Q");
120
121 void ChangeFitFunc(BlindPixelFitFunc fitfunc, Int_t par=6);
122
123 void CutAllEdges();
124
125 TObject *DrawClone(Option_t *option="") const;
126 void Draw(Option_t *option="");
127
128 Bool_t IsFitOK() { return fFitOK; }
129
130private:
131
132 BlindPixelFitFunc fgSinglePheFitFunc; //! In the beginning,
133 Int_t fgSinglePheFitNPar; //! we want to be flexible using different functions
134
135 inline static Double_t gfKto4(Double_t *x, Double_t *par)
136 {
137
138 Double_t lambda = par[0];
139
140 Double_t sum = 0.;
141 Double_t arg = 0.;
142
143 Double_t mu0 = par[1];
144 Double_t mu1 = par[2];
145
146 if (mu1 < mu0)
147 return NoWay;
148
149 Double_t sigma0 = par[3];
150 Double_t sigma1 = par[4];
151
152 if (sigma1 < sigma0)
153 return NoWay;
154
155
156 Double_t mu2 = (2.*mu1)-mu0;
157 Double_t mu3 = (3.*mu1)-(2.*mu0);
158 Double_t mu4 = (4.*mu1)-(3.*mu0);
159
160 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
161 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
162 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
163
164 Double_t lambda2 = lambda*lambda;
165 Double_t lambda3 = lambda2*lambda;
166 Double_t lambda4 = lambda3*lambda;
167
168 // k=0:
169 arg = (x[0] - mu0)/sigma0;
170 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
171
172 // k=1:
173 arg = (x[0] - mu1)/sigma1;
174 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
175
176 // k=2:
177 arg = (x[0] - mu2)/sigma2;
178 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
179
180 // k=3:
181 arg = (x[0] - mu3)/sigma3;
182 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
183
184 // k=4:
185 arg = (x[0] - mu4)/sigma4;
186 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
187
188 return TMath::Exp(-1.*lambda)*par[5]*sum;
189
190 } //->
191
192 inline static Double_t gfKto5(Double_t *x, Double_t *par)
193 {
194
195 Double_t lambda = par[0];
196
197 Double_t sum = 0.;
198 Double_t arg = 0.;
199
200 Double_t mu0 = par[1];
201 Double_t mu1 = par[2];
202
203 if (mu1 < mu0)
204 return NoWay;
205
206 Double_t sigma0 = par[3];
207 Double_t sigma1 = par[4];
208
209 if (sigma1 < sigma0)
210 return NoWay;
211
212
213 Double_t mu2 = (2.*mu1)-mu0;
214 Double_t mu3 = (3.*mu1)-(2.*mu0);
215 Double_t mu4 = (4.*mu1)-(3.*mu0);
216 Double_t mu5 = (5.*mu1)-(4.*mu0);
217
218 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
219 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
220 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
221 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
222
223 Double_t lambda2 = lambda*lambda;
224 Double_t lambda3 = lambda2*lambda;
225 Double_t lambda4 = lambda3*lambda;
226 Double_t lambda5 = lambda4*lambda;
227
228 // k=0:
229 arg = (x[0] - mu0)/sigma0;
230 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
231
232 // k=1:
233 arg = (x[0] - mu1)/sigma1;
234 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
235
236 // k=2:
237 arg = (x[0] - mu2)/sigma2;
238 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
239
240 // k=3:
241 arg = (x[0] - mu3)/sigma3;
242 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
243
244 // k=4:
245 arg = (x[0] - mu4)/sigma4;
246 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
247
248 // k=5:
249 arg = (x[0] - mu5)/sigma5;
250 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
251
252 return TMath::Exp(-1.*lambda)*par[5]*sum;
253
254 }
255
256
257 inline static Double_t gfKto6(Double_t *x, Double_t *par)
258 {
259
260 Double_t lambda = par[0];
261
262 Double_t sum = 0.;
263 Double_t arg = 0.;
264
265 Double_t mu0 = par[1];
266 Double_t mu1 = par[2];
267
268 if (mu1 < mu0)
269 return NoWay;
270
271 Double_t sigma0 = par[3];
272 Double_t sigma1 = par[4];
273
274 if (sigma1 < sigma0)
275 return NoWay;
276
277
278 Double_t mu2 = (2.*mu1)-mu0;
279 Double_t mu3 = (3.*mu1)-(2.*mu0);
280 Double_t mu4 = (4.*mu1)-(3.*mu0);
281 Double_t mu5 = (5.*mu1)-(4.*mu0);
282 Double_t mu6 = (6.*mu1)-(5.*mu0);
283
284 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
285 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
286 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
287 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
288 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
289
290 Double_t lambda2 = lambda*lambda;
291 Double_t lambda3 = lambda2*lambda;
292 Double_t lambda4 = lambda3*lambda;
293 Double_t lambda5 = lambda4*lambda;
294 Double_t lambda6 = lambda5*lambda;
295
296 // k=0:
297 arg = (x[0] - mu0)/sigma0;
298 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
299
300 // k=1:
301 arg = (x[0] - mu1)/sigma1;
302 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
303
304 // k=2:
305 arg = (x[0] - mu2)/sigma2;
306 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
307
308 // k=3:
309 arg = (x[0] - mu3)/sigma3;
310 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
311
312 // k=4:
313 arg = (x[0] - mu4)/sigma4;
314 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
315
316 // k=5:
317 arg = (x[0] - mu5)/sigma5;
318 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
319
320 // k=6:
321 arg = (x[0] - mu6)/sigma6;
322 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
323
324 return TMath::Exp(-1.*lambda)*par[5]*sum;
325
326 }
327
328 ClassDef(MHCalibrationBlindPixel, 1)
329};
330
331#endif /* MARS_MHCalibrationBlindPixel */
Note: See TracBrowser for help on using the repository browser.