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

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