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

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