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

Last change on this file since 2840 was 2835, checked in by gaug, 21 years ago
*** empty log message ***
File size: 10.8 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 fAna(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 Double_t mu2 = (2.*mu1)-mu0;
156 Double_t mu3 = (3.*mu1)-(2.*mu0);
157 Double_t mu4 = (4.*mu1)-(3.*mu0);
158
159 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
160 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
161 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
162
163 Double_t lambda2 = lambda*lambda;
164 Double_t lambda3 = lambda2*lambda;
165 Double_t lambda4 = lambda3*lambda;
166
167 // k=0:
168 arg = (x[0] - mu0)/sigma0;
169 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
170
171 // k=1:
172 arg = (x[0] - mu1)/sigma1;
173 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
174
175 // k=2:
176 arg = (x[0] - mu2)/sigma2;
177 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
178
179 // k=3:
180 arg = (x[0] - mu3)/sigma3;
181 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
182
183 // k=4:
184 arg = (x[0] - mu4)/sigma4;
185 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
186
187 return TMath::Exp(-1.*lambda)*par[5]*sum;
188
189 }
190
191 inline static Double_t fKto4(Double_t *x, Double_t *par)
192 {
193
194 Double_t lambda = par[0];
195
196 Double_t sum = 0.;
197 Double_t arg = 0.;
198
199 Double_t mu0 = par[1];
200 Double_t mu1 = par[2];
201
202 if (mu1 < mu0)
203 return NoWay;
204
205 Double_t sigma0 = par[3];
206 Double_t sigma1 = par[4];
207
208 if (sigma1 < sigma0)
209 return NoWay;
210
211 Double_t mu2 = (2.*mu1)-mu0;
212 Double_t mu3 = (3.*mu1)-(2.*mu0);
213 Double_t mu4 = (4.*mu1)-(3.*mu0);
214
215 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
216 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
217 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
218
219 Double_t lambda2 = lambda*lambda;
220 Double_t lambda3 = lambda2*lambda;
221 Double_t lambda4 = lambda3*lambda;
222
223 // k=0:
224 arg = (x[0] - mu0)/sigma0;
225 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
226
227 // k=1:
228 arg = (x[0] - mu1)/sigma1;
229 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
230
231 // k=2:
232 arg = (x[0] - mu2)/sigma2;
233 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
234
235 // k=3:
236 arg = (x[0] - mu3)/sigma3;
237 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
238
239 // k=4:
240 arg = (x[0] - mu4)/sigma4;
241 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
242
243 return TMath::Exp(-1.*lambda)*par[5]*sum;
244
245 }
246
247
248 inline static Double_t fKto5(Double_t *x, Double_t *par)
249 {
250
251 Double_t lambda = par[0];
252
253 Double_t sum = 0.;
254 Double_t arg = 0.;
255
256 Double_t mu0 = par[1];
257 Double_t mu1 = par[2];
258
259 if (mu1 < mu0)
260 return NoWay;
261
262 Double_t sigma0 = par[3];
263 Double_t sigma1 = par[4];
264
265 if (sigma1 < sigma0)
266 return NoWay;
267
268
269 Double_t mu2 = (2.*mu1)-mu0;
270 Double_t mu3 = (3.*mu1)-(2.*mu0);
271 Double_t mu4 = (4.*mu1)-(3.*mu0);
272 Double_t mu5 = (5.*mu1)-(4.*mu0);
273
274 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
275 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
276 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
277 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
278
279 Double_t lambda2 = lambda*lambda;
280 Double_t lambda3 = lambda2*lambda;
281 Double_t lambda4 = lambda3*lambda;
282 Double_t lambda5 = lambda4*lambda;
283
284 // k=0:
285 arg = (x[0] - mu0)/sigma0;
286 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
287
288 // k=1:
289 arg = (x[0] - mu1)/sigma1;
290 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
291
292 // k=2:
293 arg = (x[0] - mu2)/sigma2;
294 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
295
296 // k=3:
297 arg = (x[0] - mu3)/sigma3;
298 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
299
300 // k=4:
301 arg = (x[0] - mu4)/sigma4;
302 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
303
304 // k=5:
305 arg = (x[0] - mu5)/sigma5;
306 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
307
308 return TMath::Exp(-1.*lambda)*par[5]*sum;
309
310 }
311
312
313 inline static Double_t fKto6(Double_t *x, Double_t *par)
314 {
315
316 Double_t lambda = par[0];
317
318 Double_t sum = 0.;
319 Double_t arg = 0.;
320
321 Double_t mu0 = par[1];
322 Double_t mu1 = par[2];
323
324 if (mu1 < mu0)
325 return NoWay;
326
327 Double_t sigma0 = par[3];
328 Double_t sigma1 = par[4];
329
330 if (sigma1 < sigma0)
331 return NoWay;
332
333
334 Double_t mu2 = (2.*mu1)-mu0;
335 Double_t mu3 = (3.*mu1)-(2.*mu0);
336 Double_t mu4 = (4.*mu1)-(3.*mu0);
337 Double_t mu5 = (5.*mu1)-(4.*mu0);
338 Double_t mu6 = (6.*mu1)-(5.*mu0);
339
340 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
341 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
342 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
343 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
344 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
345
346 Double_t lambda2 = lambda*lambda;
347 Double_t lambda3 = lambda2*lambda;
348 Double_t lambda4 = lambda3*lambda;
349 Double_t lambda5 = lambda4*lambda;
350 Double_t lambda6 = lambda5*lambda;
351
352 // k=0:
353 arg = (x[0] - mu0)/sigma0;
354 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
355
356 // k=1:
357 arg = (x[0] - mu1)/sigma1;
358 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
359
360 // k=2:
361 arg = (x[0] - mu2)/sigma2;
362 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
363
364 // k=3:
365 arg = (x[0] - mu3)/sigma3;
366 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
367
368 // k=4:
369 arg = (x[0] - mu4)/sigma4;
370 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
371
372 // k=5:
373 arg = (x[0] - mu5)/sigma5;
374 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
375
376 // k=6:
377 arg = (x[0] - mu6)/sigma6;
378 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
379
380 return TMath::Exp(-1.*lambda)*par[5]*sum;
381
382 }
383
384 ClassDef(MHCalibrationBlindPixel, 1)
385};
386
387#endif /* MARS_MHCalibrationBlindPixel */
Note: See TracBrowser for help on using the repository browser.