source: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h@ 4603

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