source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h@ 5050

Last change on this file since 5050 was 5050, checked in by gaug, 20 years ago
*** empty log message ***
File size: 16.7 KB
Line 
1#ifndef MARS_MHCalibrationChargeBlindPix
2#define MARS_MHCalibrationChargeBlindPix
3
4#ifndef MARS_MHCalibrationPix
5#include "MHCalibrationPix.h"
6#endif
7
8#ifndef ROOT_TF1
9#include <TF1.h>
10#endif
11
12#ifndef ROOT_TVector
13#include <TVector.h>
14#endif
15
16class TH1F;
17class TPaveText;
18class TText;
19class MExtractedSignalBlindPixel;
20class MRawEvtPixelIter;
21
22class MHCalibrationChargeBlindPix : public MHCalibrationPix
23{
24private:
25
26 static const Float_t fgNumSinglePheLimit; //! Default for fNumSinglePheLimit (now set to: 50)
27 static const Float_t gkSignalInitializer; //! Signal initializer (-9999.)
28
29 static const Double_t gkElectronicAmp; // Electronic Amplification after the PMT (in FADC counts/N_e)
30 static const Double_t gkElectronicAmpErr; // Error of the electronic amplification
31
32 Float_t fSinglePheCut; // Value of summed FADC slices upon which event considered as single-phe
33 Float_t fNumSinglePheLimit; // Minimum number of single-phe events
34
35 TVector fASinglePheFADCSlices; // Averaged FADC slice entries supposed single-phe events
36 TVector fAPedestalFADCSlices; // Averaged FADC slice entries supposed pedestal events
37
38 TF1 *fSinglePheFit; // Single Phe Fit (Gaussians convoluted with Poisson)
39
40 UInt_t fNumSinglePhes; // Number of entries in fASinglePheFADCSlices
41 UInt_t fNumPedestals; // Number of entries in fAPedestalFADCSlices
42
43 Double_t fLambda; // Poisson mean from Single-phe fit
44 Double_t fLambdaCheck; // Poisson mean from Pedestal fit alone
45 Double_t fMu0; // Mean of the pedestal
46 Double_t fMu1; // Mean of single-phe peak
47 Double_t fSigma0; // Sigma of the pedestal
48 Double_t fSigma1; // Sigma of single-phe peak
49 Double_t fLambdaErr; // Error of Poisson mean from Single-phe fit
50 Double_t fLambdaCheckErr; // Error of Poisson mean from Pedestal fit alone
51 Double_t fMu0Err; // Error of Mean of the pedestal
52 Double_t fMu1Err; // Error of Mean of single-phe peak
53 Double_t fSigma0Err; // Error of Sigma of the pedestal
54 Double_t fSigma1Err; // Error of Sigma of single-phe peak
55 Double_t fChisquare; // Chisquare of single-phe fit
56 Int_t fNDF; // Ndof of single-phe fit
57 Double_t fProb; // Probability of singleo-phe fit
58
59 Byte_t fFlags; // Bit-field for the flags
60 enum { kSinglePheFitOK, kPedestalFitOK }; // Possible bits to be set
61
62public:
63
64 enum FitFunc_t { kEPoisson4, kEPoisson5,
65 kEPoisson6, kEPoisson7,
66 kEPolya, kEMichele }; // Possible fit functions types
67
68private:
69
70 FitFunc_t fFitFunc;
71
72 TPaveText *fFitLegend; //! Some legend to display the fit results
73 TH1F *fHSinglePheFADCSlices; //! A histogram created and deleted only in Draw()
74 TH1F *fHPedestalFADCSlices; //! A histogram created and deleted only in Draw()
75
76 // Fit
77 Bool_t InitFit();
78 void ExitFit();
79
80 void DrawLegend(Option_t *opt="");
81
82public:
83
84 MHCalibrationChargeBlindPix(const char *name=NULL, const char *title=NULL);
85 ~MHCalibrationChargeBlindPix();
86
87 void Clear(Option_t *o="");
88 void Reset() {}
89
90 // TObject *Clone(const char *) const;
91
92 // Getters
93 const Double_t GetLambda () const { return fLambda; }
94 const Double_t GetLambdaCheck () const { return fLambdaCheck; }
95 const Double_t GetMu0 () const { return fMu0; }
96 const Double_t GetMu1 () const { return fMu1; }
97 const Double_t GetSigma0 () const { return fSigma0; }
98 const Double_t GetSigma1 () const { return fSigma1; }
99 const Double_t GetLambdaErr () const { return fLambdaErr; }
100 const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; }
101 const Double_t GetMu0Err () const { return fMu0Err; }
102 const Double_t GetMu1Err () const { return fMu1Err; }
103 const Double_t GetSigma0Err () const { return fSigma0Err; }
104 const Double_t GetSigma1Err () const { return fSigma1Err; }
105 const Float_t GetSinglePheCut () const { return fSinglePheCut; }
106
107 TVector &GetASinglePheFADCSlices() { return fASinglePheFADCSlices; }
108 const TVector &GetASinglePheFADCSlices() const { return fASinglePheFADCSlices; }
109
110 TVector &GetAPedestalFADCSlices() { return fAPedestalFADCSlices; }
111 const TVector &GetAPedestalFADCSlices() const { return fAPedestalFADCSlices; }
112
113 const Bool_t IsSinglePheFitOK() const;
114 const Bool_t IsPedestalFitOK() const;
115
116 // Setters
117 void SetFitFunc ( const FitFunc_t func ) { fFitFunc = func; }
118 void SetSinglePheCut ( const Float_t cut = 0. ) { fSinglePheCut = cut; }
119 void SetNumSinglePheLimit ( const Float_t lim =fgNumSinglePheLimit ) { fNumSinglePheLimit = lim; }
120
121 void SetSinglePheFitOK ( const Bool_t b=kTRUE);
122 void SetPedestalFitOK ( const Bool_t b=kTRUE);
123
124 // Fill histos
125 void FillSinglePheFADCSlices(const MRawEvtPixelIter &iter);
126 void FillPedestalFADCSlices( const MRawEvtPixelIter &iter);
127 void FinalizeSinglePheSpectrum();
128
129 // Draws
130 void Draw(Option_t *opt="");
131
132 // Fits
133 Bool_t FitSinglePhe (Option_t *opt="RL0+Q");
134 void FitPedestal (Option_t *opt="RL0+Q");
135
136 // Simulation
137 Bool_t SimulateSinglePhe(const Double_t lambda,
138 const Double_t mu0, const Double_t mu1,
139 const Double_t sigma0, const Double_t sigma1);
140
141private:
142
143 inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par)
144 {
145
146 Double_t lambda1cat = par[0];
147 Double_t lambda1dyn = par[1];
148 Double_t mu0 = par[2];
149 Double_t mu1cat = par[3];
150 Double_t mu1dyn = par[4];
151 Double_t sigma0 = par[5];
152 Double_t sigma1cat = par[6];
153 Double_t sigma1dyn = par[7];
154
155 Double_t sumcat = 0.;
156 Double_t sumdyn = 0.;
157 Double_t arg = 0.;
158
159 if (lambda1cat < lambda1dyn)
160 return FLT_MAX;
161
162 if (mu1cat < mu0)
163 return FLT_MAX;
164
165 if (mu1dyn < mu0)
166 return FLT_MAX;
167
168 if (mu1cat < mu1dyn)
169 return FLT_MAX;
170
171 if (sigma0 < 0.0001)
172 return FLT_MAX;
173
174 if (sigma1cat < sigma0)
175 return FLT_MAX;
176
177 if (sigma1dyn < sigma0)
178 return FLT_MAX;
179
180 Double_t mu2cat = (2.*mu1cat)-mu0;
181 Double_t mu2dyn = (2.*mu1dyn)-mu0;
182 Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
183 Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
184
185 Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0));
186 Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0));
187 Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
188 Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
189
190 Double_t lambda2cat = lambda1cat*lambda1cat;
191 Double_t lambda2dyn = lambda1dyn*lambda1dyn;
192 Double_t lambda3cat = lambda2cat*lambda1cat;
193 Double_t lambda3dyn = lambda2dyn*lambda1dyn;
194
195 // k=0:
196 arg = (x[0] - mu0)/sigma0;
197 sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
198 sumdyn = sumcat;
199
200 // k=1cat:
201 arg = (x[0] - mu1cat)/sigma1cat;
202 sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
203 // k=1dyn:
204 arg = (x[0] - mu1dyn)/sigma1dyn;
205 sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
206
207 // k=2cat:
208 arg = (x[0] - mu2cat)/sigma2cat;
209 sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
210 // k=2dyn:
211 arg = (x[0] - mu2dyn)/sigma2dyn;
212 sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
213
214
215 // k=3cat:
216 arg = (x[0] - mu3cat)/sigma3cat;
217 sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
218 // k=3dyn:
219 arg = (x[0] - mu3dyn)/sigma3dyn;
220 sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn;
221
222 sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
223 sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
224
225 return par[8]*(sumcat+sumdyn)/2.;
226
227 }
228
229 inline static Double_t fPoissonKto4(Double_t *x, Double_t *par)
230 {
231
232 Double_t lambda = par[0];
233
234 Double_t sum = 0.;
235 Double_t arg = 0.;
236
237 Double_t mu0 = par[1];
238 Double_t mu1 = par[2];
239
240 if (mu1 < mu0)
241 return FLT_MAX;
242
243 Double_t sigma0 = par[3];
244 Double_t sigma1 = par[4];
245
246 if (sigma0 < 0.0001)
247 return FLT_MAX;
248
249 if (sigma1 < sigma0)
250 return FLT_MAX;
251
252 Double_t mu2 = (2.*mu1)-mu0;
253 Double_t mu3 = (3.*mu1)-(2.*mu0);
254 Double_t mu4 = (4.*mu1)-(3.*mu0);
255
256 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
257 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
258 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
259
260 Double_t lambda2 = lambda*lambda;
261 Double_t lambda3 = lambda2*lambda;
262 Double_t lambda4 = lambda3*lambda;
263
264 // k=0:
265 arg = (x[0] - mu0)/sigma0;
266 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
267
268 // k=1:
269 arg = (x[0] - mu1)/sigma1;
270 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
271
272 // k=2:
273 arg = (x[0] - mu2)/sigma2;
274 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
275
276 // k=3:
277 arg = (x[0] - mu3)/sigma3;
278 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
279
280 // k=4:
281 arg = (x[0] - mu4)/sigma4;
282 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
283
284 return TMath::Exp(-1.*lambda)*par[5]*sum;
285
286 }
287
288
289 inline static Double_t fPoissonKto5(Double_t *x, Double_t *par)
290 {
291
292 Double_t lambda = par[0];
293
294 Double_t sum = 0.;
295 Double_t arg = 0.;
296
297 Double_t mu0 = par[1];
298 Double_t mu1 = par[2];
299
300 if (mu1 < mu0)
301 return FLT_MAX;
302
303 Double_t sigma0 = par[3];
304 Double_t sigma1 = par[4];
305
306 if (sigma0 < 0.0001)
307 return FLT_MAX;
308
309 if (sigma1 < sigma0)
310 return FLT_MAX;
311
312
313 Double_t mu2 = (2.*mu1)-mu0;
314 Double_t mu3 = (3.*mu1)-(2.*mu0);
315 Double_t mu4 = (4.*mu1)-(3.*mu0);
316 Double_t mu5 = (5.*mu1)-(4.*mu0);
317
318 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
319 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
320 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
321 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
322
323 Double_t lambda2 = lambda*lambda;
324 Double_t lambda3 = lambda2*lambda;
325 Double_t lambda4 = lambda3*lambda;
326 Double_t lambda5 = lambda4*lambda;
327
328 // k=0:
329 arg = (x[0] - mu0)/sigma0;
330 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
331
332 // k=1:
333 arg = (x[0] - mu1)/sigma1;
334 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
335
336 // k=2:
337 arg = (x[0] - mu2)/sigma2;
338 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
339
340 // k=3:
341 arg = (x[0] - mu3)/sigma3;
342 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
343
344 // k=4:
345 arg = (x[0] - mu4)/sigma4;
346 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
347
348 // k=5:
349 arg = (x[0] - mu5)/sigma5;
350 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
351
352 return TMath::Exp(-1.*lambda)*par[5]*sum;
353
354 }
355
356
357 inline static Double_t fPoissonKto6(Double_t *x, Double_t *par)
358 {
359
360 Double_t lambda = par[0];
361
362 Double_t sum = 0.;
363 Double_t arg = 0.;
364
365 Double_t mu0 = par[1];
366 Double_t mu1 = par[2];
367
368 if (mu1 < mu0)
369 return FLT_MAX;
370
371 Double_t sigma0 = par[3];
372 Double_t sigma1 = par[4];
373
374 if (sigma0 < 0.0001)
375 return FLT_MAX;
376
377 if (sigma1 < sigma0)
378 return FLT_MAX;
379
380
381 Double_t mu2 = (2.*mu1)-mu0;
382 Double_t mu3 = (3.*mu1)-(2.*mu0);
383 Double_t mu4 = (4.*mu1)-(3.*mu0);
384 Double_t mu5 = (5.*mu1)-(4.*mu0);
385 Double_t mu6 = (6.*mu1)-(5.*mu0);
386
387 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
388 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
389 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
390 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
391 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
392
393 Double_t lambda2 = lambda*lambda;
394 Double_t lambda3 = lambda2*lambda;
395 Double_t lambda4 = lambda3*lambda;
396 Double_t lambda5 = lambda4*lambda;
397 Double_t lambda6 = lambda5*lambda;
398
399 // k=0:
400 arg = (x[0] - mu0)/sigma0;
401 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
402
403 // k=1:
404 arg = (x[0] - mu1)/sigma1;
405 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
406
407 // k=2:
408 arg = (x[0] - mu2)/sigma2;
409 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
410
411 // k=3:
412 arg = (x[0] - mu3)/sigma3;
413 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
414
415 // k=4:
416 arg = (x[0] - mu4)/sigma4;
417 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
418
419 // k=5:
420 arg = (x[0] - mu5)/sigma5;
421 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
422
423 // k=6:
424 arg = (x[0] - mu6)/sigma6;
425 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
426
427 return TMath::Exp(-1.*lambda)*par[5]*sum;
428
429 }
430
431 inline static Double_t fPolya(Double_t *x, Double_t *par)
432 {
433
434 const Double_t QEcat = 0.247; // mean quantum efficiency
435 const Double_t sqrt2 = 1.4142135623731;
436 const Double_t sqrt3 = 1.7320508075689;
437 const Double_t sqrt4 = 2.;
438
439 const Double_t lambda = par[0]; // mean number of photons
440
441 const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution
442 const Double_t delta1 = par[2]; // amplification first dynode
443 const Double_t delta2 = par[3]; // amplification subsequent dynodes
444
445 const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges
446
447 const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain
448 const Double_t A = 1. + excessPoisson - QEcat
449 + 1./delta1
450 + 1./delta1/delta2
451 + 1./delta1/delta2/delta2; // variance contributions from PMT and QE
452
453 const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion
454
455 const Double_t mu0 = par[7]; // pedestal
456 const Double_t mu1 = totAmpl; // single phe position
457 const Double_t mu2 = 2*totAmpl; // double phe position
458 const Double_t mu3 = 3*totAmpl; // triple phe position
459 const Double_t mu4 = 4*totAmpl; // quadruple phe position
460
461 const Double_t sigma0 = par[5];
462 const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
463 const Double_t sigma2 = sqrt2*sigma1;
464 const Double_t sigma3 = sqrt3*sigma1;
465 const Double_t sigma4 = sqrt4*sigma1;
466
467 const Double_t lambda2 = lambda*lambda;
468 const Double_t lambda3 = lambda2*lambda;
469 const Double_t lambda4 = lambda3*lambda;
470
471 //-- calculate the area----
472 Double_t arg = (x[0] - mu0)/sigma0;
473 Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
474
475 // k=1:
476 arg = (x[0] - mu1)/sigma1;
477 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
478
479 // k=2:
480 arg = (x[0] - mu2)/sigma2;
481 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
482
483 // k=3:
484 arg = (x[0] - mu3)/sigma3;
485 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
486
487 // k=4:
488 arg = (x[0] - mu4)/sigma4;
489 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
490
491 return TMath::Exp(-1.*lambda)*par[6]*sum;
492 }
493
494 ClassDef(MHCalibrationChargeBlindPix, 1) // Histogram class for Charge Blind Pixel Calibration
495};
496
497#endif /* MARS_MHCalibrationChargeBlindPix */
Note: See TracBrowser for help on using the repository browser.