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

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