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

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