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

Last change on this file since 3063 was 3061, checked in by gaug, 21 years ago
*** empty log message ***
File size: 15.8 KB
Line 
1#ifndef MARS_MHCalibrationBlindPixel
2#define MARS_MHCalibrationBlindPixel
3
4#ifndef MARS_MH
5#include "MH.h"
6#endif
7
8class TArrayF;
9class TH1F;
10class TH1I;
11class TF1;
12class TPaveText;
13
14class TMath;
15class MParList;
16class MHCalibrationBlindPixel : public MH
17{
18private:
19
20 static const Int_t fgBlindPixelChargeNbins;
21 static const Int_t fgBlindPixelTimeNbins;
22 static const Axis_t fgBlindPixelTimeFirst;
23 static const Axis_t fgBlindPixelTimeLast;
24 static const Double_t fgBlindPixelElectronicAmp;
25 static const Double_t fgBlindPixelElectronicAmpError;
26
27 static const Axis_t fNyquistFreq;
28 static const Axis_t fMinFreq;
29 static const Int_t fPSDNbins;
30
31 TH1F* fHBlindPixelCharge; // Histogram with the single Phe spectrum
32 TH1F* fHBlindPixelTime; // Variance of summed FADC slices
33 TH1F* fHBlindPixelPSD; // Power spectrum density of fHBlindPixelChargevsN
34
35 TF1 *fSinglePheFit;
36 TF1 *fTimeGausFit;
37 TF1 *fSinglePhePedFit;
38
39 TArrayF* fPSDHiGain; //-> Power spectrum density of fHiGains
40 TArrayF* fPSDLoGain; //-> Power spectrum density of fLoGains
41
42 TH1F* fHPSD; //->
43 TF1* fPSDExpFit; //->
44
45 TArrayF *fHiGains; //->
46 TArrayF *fLoGains; //->
47 TArrayF *fChargeXaxis; //
48 TArrayF *fPSDXaxis; //
49
50 Float_t fPSDProb;
51
52 Int_t fTotalEntries; // Number of entries
53 Int_t fCurrentSize;
54
55 Axis_t fBlindPixelChargefirst;
56 Axis_t fBlindPixelChargelast;
57
58 void DrawLegend();
59 void CreateChargeXaxis(Int_t n);
60 void CreatePSDXaxis(Int_t n);
61 void CutArrayBorder(TArrayF *array);
62
63 TPaveText *fFitLegend;
64
65 Double_t fLambda;
66 Double_t fMu0;
67 Double_t fMu1;
68 Double_t fSigma0;
69 Double_t fSigma1;
70
71 Double_t fLambdaErr;
72 Double_t fMu0Err;
73 Double_t fMu1Err;
74 Double_t fSigma0Err;
75 Double_t fSigma1Err;
76
77 Double_t fChisquare;
78 Double_t fProb;
79 Int_t fNdf;
80
81 Double_t fMeanTime;
82 Double_t fMeanTimeErr;
83 Double_t fSigmaTime;
84 Double_t fSigmaTimeErr;
85
86 Double_t fLambdaCheck;
87 Double_t fLambdaCheckErr;
88
89 Double_t fMeanPedestal;
90 Double_t fMeanPedestalErr;
91 Double_t fSigmaPedestal;
92 Double_t fSigmaPedestalErr;
93
94 Byte_t fFlags;
95
96 enum { kFitOK, kOscillating };
97
98
99public:
100
101 MHCalibrationBlindPixel(const char *name=NULL, const char *title=NULL);
102 ~MHCalibrationBlindPixel();
103
104 void Clear(Option_t *o="");
105 void Reset();
106
107 Bool_t FillBlindPixelCharge(Float_t q);
108 Bool_t FillBlindPixelTime(Float_t t);
109 Bool_t FillGraphs(Float_t qhi, Float_t qlo);
110
111 // Setters
112 void SetMeanPedestal(const Float_t f) { fMeanPedestal = f; }
113 void SetMeanPedestalErr(const Float_t f) { fMeanPedestalErr = f; }
114 void SetSigmaPedestal(const Float_t f) { fSigmaPedestal = f; }
115 void SetSigmaPedestalErr(const Float_t f) { fSigmaPedestalErr = f; }
116
117 // Getters
118 const Double_t GetLambda() const { return fLambda; }
119 const Double_t GetLambdaCheck() const { return fLambdaCheck; }
120 const Double_t GetMu0() const { return fMu0; }
121 const Double_t GetMu1() const { return fMu1; }
122 const Double_t GetSigma0() const { return fSigma0; }
123 const Double_t GetSigma1() const { return fSigma1; }
124
125 const Double_t GetLambdaErr() const { return fLambdaErr; }
126 const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; }
127 const Double_t GetMu0Err() const { return fMu0Err; }
128 const Double_t GetMu1Err() const { return fMu1Err; }
129 const Double_t GetSigma0Err() const { return fSigma0Err; }
130 const Double_t GetSigma1Err() const { return fSigma1Err; }
131
132 const Double_t GetChiSquare() const { return fChisquare; }
133 const Double_t GetProb() const { return fProb; }
134 const Int_t GetNdf() const { return fNdf; }
135
136 const Double_t GetMeanTime() const { return fMeanTime; }
137 const Double_t GetMeanTimeErr() const { return fMeanTimeErr; }
138 const Double_t GetSigmaTime() const { return fSigmaTime; }
139 const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; }
140
141 const Bool_t IsFitOK() const;
142 const Bool_t IsOscillating();
143
144 const TH1F *GetHBlindPixelPSD() const { return fHBlindPixelPSD; }
145
146 // Draws
147 TObject *DrawClone(Option_t *option="") const;
148 void Draw(Option_t *option="");
149
150
151 // Fits
152 enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele };
153
154private:
155 FitFunc_t fFitFunc;
156
157public:
158 Bool_t FitSinglePhe(Axis_t rmin=0, Axis_t rmax=0, Option_t *opt="RL0+Q");
159 Bool_t FitTime(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+Q");
160 void ChangeFitFunc(FitFunc_t func) { fFitFunc = func; }
161
162 // Simulation
163 Bool_t SimulateSinglePhe(Double_t lambda,
164 Double_t mu0,Double_t mu1,
165 Double_t sigma0,Double_t sigma1);
166
167 // Others
168 void CutAllEdges();
169 Bool_t CheckOscillations();
170
171
172private:
173
174 const static Double_t fNoWay = 10000000000.0;
175
176 Bool_t InitFit(Axis_t min, Axis_t max);
177 void ExitFit(TF1 *f);
178
179 inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par)
180 {
181
182 Double_t lambda1cat = par[0];
183 Double_t lambda1dyn = par[1];
184 Double_t mu0 = par[2];
185 Double_t mu1cat = par[3];
186 Double_t mu1dyn = par[4];
187 Double_t sigma0 = par[5];
188 Double_t sigma1cat = par[6];
189 Double_t sigma1dyn = par[7];
190
191 Double_t sumcat = 0.;
192 Double_t sumdyn = 0.;
193 Double_t arg = 0.;
194
195 if (mu1cat < mu0)
196 return fNoWay;
197
198 if (sigma1cat < sigma0)
199 return fNoWay;
200
201 // if (sigma1cat < sigma1dyn)
202 // return NoWay;
203
204 //if (mu1cat < mu1dyn)
205 // return NoWay;
206
207 // if (lambda1cat < lambda1dyn)
208 // return NoWay;
209
210 Double_t mu2cat = (2.*mu1cat)-mu0;
211 Double_t mu2dyn = (2.*mu1dyn)-mu0;
212 Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
213 Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
214
215 Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0));
216 Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0));
217 Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
218 Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
219
220 Double_t lambda2cat = lambda1cat*lambda1cat;
221 Double_t lambda2dyn = lambda1dyn*lambda1dyn;
222 Double_t lambda3cat = lambda2cat*lambda1cat;
223 Double_t lambda3dyn = lambda2dyn*lambda1dyn;
224
225 // k=0:
226 arg = (x[0] - mu0)/sigma0;
227 sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
228 sumdyn =sumcat;
229
230 // k=1cat:
231 arg = (x[0] - mu1cat)/sigma1cat;
232 sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
233 // k=1dyn:
234 arg = (x[0] - mu1dyn)/sigma1dyn;
235 sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
236
237 // k=2cat:
238 arg = (x[0] - mu2cat)/sigma2cat;
239 sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
240 // k=2dyn:
241 arg = (x[0] - mu2dyn)/sigma2dyn;
242 sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
243
244
245 // k=3cat:
246 arg = (x[0] - mu3cat)/sigma3cat;
247 sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
248 // k=3dyn:
249 arg = (x[0] - mu3dyn)/sigma3dyn;
250 sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn;
251
252 sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
253 sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
254
255 return par[8]*(sumcat+sumdyn)/2.;
256
257 }
258
259 inline static Double_t fPoissonKto4(Double_t *x, Double_t *par)
260 {
261
262 Double_t lambda = par[0];
263
264 Double_t sum = 0.;
265 Double_t arg = 0.;
266
267 Double_t mu0 = par[1];
268 Double_t mu1 = par[2];
269
270 if (mu1 < mu0)
271 return fNoWay;
272
273 Double_t sigma0 = par[3];
274 Double_t sigma1 = par[4];
275
276 if (sigma1 < sigma0)
277 return fNoWay;
278
279 Double_t mu2 = (2.*mu1)-mu0;
280 Double_t mu3 = (3.*mu1)-(2.*mu0);
281 Double_t mu4 = (4.*mu1)-(3.*mu0);
282
283 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
284 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
285 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
286
287 Double_t lambda2 = lambda*lambda;
288 Double_t lambda3 = lambda2*lambda;
289 Double_t lambda4 = lambda3*lambda;
290
291 // k=0:
292 arg = (x[0] - mu0)/sigma0;
293 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
294
295 // k=1:
296 arg = (x[0] - mu1)/sigma1;
297 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
298
299 // k=2:
300 arg = (x[0] - mu2)/sigma2;
301 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
302
303 // k=3:
304 arg = (x[0] - mu3)/sigma3;
305 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
306
307 // k=4:
308 arg = (x[0] - mu4)/sigma4;
309 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
310
311 return TMath::Exp(-1.*lambda)*par[5]*sum;
312
313 }
314
315
316 inline static Double_t fPoissonKto5(Double_t *x, Double_t *par)
317 {
318
319 Double_t lambda = par[0];
320
321 Double_t sum = 0.;
322 Double_t arg = 0.;
323
324 Double_t mu0 = par[1];
325 Double_t mu1 = par[2];
326
327 if (mu1 < mu0)
328 return fNoWay;
329
330 Double_t sigma0 = par[3];
331 Double_t sigma1 = par[4];
332
333 if (sigma1 < sigma0)
334 return fNoWay;
335
336
337 Double_t mu2 = (2.*mu1)-mu0;
338 Double_t mu3 = (3.*mu1)-(2.*mu0);
339 Double_t mu4 = (4.*mu1)-(3.*mu0);
340 Double_t mu5 = (5.*mu1)-(4.*mu0);
341
342 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
343 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
344 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
345 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
346
347 Double_t lambda2 = lambda*lambda;
348 Double_t lambda3 = lambda2*lambda;
349 Double_t lambda4 = lambda3*lambda;
350 Double_t lambda5 = lambda4*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 return TMath::Exp(-1.*lambda)*par[5]*sum;
377
378 }
379
380
381 inline static Double_t fPoissonKto6(Double_t *x, Double_t *par)
382 {
383
384 Double_t lambda = par[0];
385
386 Double_t sum = 0.;
387 Double_t arg = 0.;
388
389 Double_t mu0 = par[1];
390 Double_t mu1 = par[2];
391
392 if (mu1 < mu0)
393 return fNoWay;
394
395 Double_t sigma0 = par[3];
396 Double_t sigma1 = par[4];
397
398 if (sigma1 < sigma0)
399 return fNoWay;
400
401
402 Double_t mu2 = (2.*mu1)-mu0;
403 Double_t mu3 = (3.*mu1)-(2.*mu0);
404 Double_t mu4 = (4.*mu1)-(3.*mu0);
405 Double_t mu5 = (5.*mu1)-(4.*mu0);
406 Double_t mu6 = (6.*mu1)-(5.*mu0);
407
408 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
409 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
410 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
411 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
412 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
413
414 Double_t lambda2 = lambda*lambda;
415 Double_t lambda3 = lambda2*lambda;
416 Double_t lambda4 = lambda3*lambda;
417 Double_t lambda5 = lambda4*lambda;
418 Double_t lambda6 = lambda5*lambda;
419
420 // k=0:
421 arg = (x[0] - mu0)/sigma0;
422 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
423
424 // k=1:
425 arg = (x[0] - mu1)/sigma1;
426 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
427
428 // k=2:
429 arg = (x[0] - mu2)/sigma2;
430 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
431
432 // k=3:
433 arg = (x[0] - mu3)/sigma3;
434 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
435
436 // k=4:
437 arg = (x[0] - mu4)/sigma4;
438 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
439
440 // k=5:
441 arg = (x[0] - mu5)/sigma5;
442 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
443
444 // k=6:
445 arg = (x[0] - mu6)/sigma6;
446 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
447
448 return TMath::Exp(-1.*lambda)*par[5]*sum;
449
450 }
451
452 inline static Double_t fPolya(Double_t *x, Double_t *par)
453 {
454
455 const Double_t QEcat = 0.247; // mean quantum efficiency
456 const Double_t sqrt2 = 1.4142135623731;
457 const Double_t sqrt3 = 1.7320508075689;
458 const Double_t sqrt4 = 2.;
459
460 const Double_t lambda = par[0]; // mean number of photons
461
462 const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution
463 const Double_t delta1 = par[2]; // amplification first dynode
464 const Double_t delta2 = par[3]; // amplification subsequent dynodes
465
466 const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges
467
468 const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain
469 const Double_t A = 1. + excessPoisson - QEcat
470 + 1./delta1
471 + 1./delta1/delta2
472 + 1./delta1/delta2/delta2; // variance contributions from PMT and QE
473
474 const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion
475
476 const Double_t mu0 = par[7]; // pedestal
477 const Double_t mu1 = totAmpl; // single phe position
478 const Double_t mu2 = 2*totAmpl; // double phe position
479 const Double_t mu3 = 3*totAmpl; // triple phe position
480 const Double_t mu4 = 4*totAmpl; // quadruple phe position
481
482 const Double_t sigma0 = par[5];
483 const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
484 const Double_t sigma2 = sqrt2*sigma1;
485 const Double_t sigma3 = sqrt3*sigma1;
486 const Double_t sigma4 = sqrt4*sigma1;
487
488 const Double_t lambda2 = lambda*lambda;
489 const Double_t lambda3 = lambda2*lambda;
490 const Double_t lambda4 = lambda3*lambda;
491
492 //-- calculate the area----
493 Double_t arg = (x[0] - mu0)/sigma0;
494 Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
495
496 // k=1:
497 arg = (x[0] - mu1)/sigma1;
498 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
499
500 // k=2:
501 arg = (x[0] - mu2)/sigma2;
502 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
503
504 // k=3:
505 arg = (x[0] - mu3)/sigma3;
506 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
507
508 // k=4:
509 arg = (x[0] - mu4)/sigma4;
510 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
511
512 return TMath::Exp(-1.*lambda)*par[6]*sum;
513 }
514
515
516
517 ClassDef(MHCalibrationBlindPixel, 1) // Histograms from the Calibration Blind Pixel
518};
519
520#endif /* MARS_MHCalibrationBlindPixel */
Note: See TracBrowser for help on using the repository browser.