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

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