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

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