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

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