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

Last change on this file since 4464 was 4399, 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();
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 Double_t offset = par[9];
180
181 Double_t sumcat = 0.;
182 Double_t sumdyn = 0.;
183 Double_t arg = 0.;
184
185 if (lambda1cat < lambda1dyn)
186 return FLT_MAX;
187
188 if (mu1cat < mu0)
189 return FLT_MAX;
190
191 if (mu1dyn < mu0)
192 return FLT_MAX;
193
194 if (mu1cat < mu1dyn)
195 return FLT_MAX;
196
197 if (sigma0 < 0.0001)
198 return FLT_MAX;
199
200 if (sigma1cat < sigma0)
201 return FLT_MAX;
202
203 if (sigma1dyn < sigma0)
204 return FLT_MAX;
205
206 Double_t mu2cat = (2.*mu1cat)-mu0;
207 Double_t mu2dyn = (2.*mu1dyn)-mu0;
208 Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
209 Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
210
211 Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0));
212 Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0));
213 Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
214 Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
215
216 Double_t lambda2cat = lambda1cat*lambda1cat;
217 Double_t lambda2dyn = lambda1dyn*lambda1dyn;
218 Double_t lambda3cat = lambda2cat*lambda1cat;
219 Double_t lambda3dyn = lambda2dyn*lambda1dyn;
220
221 // k=0:
222 arg = (x[0] - mu0)/sigma0;
223 sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
224 sumdyn =sumcat;
225
226 // k=1cat:
227 arg = (x[0] - mu1cat)/sigma1cat;
228 sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
229 // k=1dyn:
230 arg = (x[0] - mu1dyn)/sigma1dyn;
231 sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
232
233 // k=2cat:
234 arg = (x[0] - mu2cat)/sigma2cat;
235 sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
236 // k=2dyn:
237 arg = (x[0] - mu2dyn)/sigma2dyn;
238 sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
239
240
241 // k=3cat:
242 arg = (x[0] - mu3cat)/sigma3cat;
243 sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
244 // k=3dyn:
245 arg = (x[0] - mu3dyn)/sigma3dyn;
246 sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn;
247
248 sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
249 sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
250
251 return par[8]*(sumcat+sumdyn)/2. + offset;
252
253 }
254
255 inline static Double_t fPoissonKto4(Double_t *x, Double_t *par)
256 {
257
258 Double_t lambda = par[0];
259
260 Double_t sum = 0.;
261 Double_t arg = 0.;
262
263 Double_t mu0 = par[1];
264 Double_t mu1 = par[2];
265
266 if (mu1 < mu0)
267 return FLT_MAX;
268
269 Double_t sigma0 = par[3];
270 Double_t sigma1 = par[4];
271
272 if (sigma0 < 0.0001)
273 return FLT_MAX;
274
275 if (sigma1 < sigma0)
276 return FLT_MAX;
277
278 Double_t mu2 = (2.*mu1)-mu0;
279 Double_t mu3 = (3.*mu1)-(2.*mu0);
280 Double_t mu4 = (4.*mu1)-(3.*mu0);
281
282 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
283 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
284 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
285
286 Double_t lambda2 = lambda*lambda;
287 Double_t lambda3 = lambda2*lambda;
288 Double_t lambda4 = lambda3*lambda;
289
290 // k=0:
291 arg = (x[0] - mu0)/sigma0;
292 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
293
294 // k=1:
295 arg = (x[0] - mu1)/sigma1;
296 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
297
298 // k=2:
299 arg = (x[0] - mu2)/sigma2;
300 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
301
302 // k=3:
303 arg = (x[0] - mu3)/sigma3;
304 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
305
306 // k=4:
307 arg = (x[0] - mu4)/sigma4;
308 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
309
310 return TMath::Exp(-1.*lambda)*par[5]*sum;
311
312 }
313
314
315 inline static Double_t fPoissonKto5(Double_t *x, Double_t *par)
316 {
317
318 Double_t lambda = par[0];
319
320 Double_t sum = 0.;
321 Double_t arg = 0.;
322
323 Double_t mu0 = par[1];
324 Double_t mu1 = par[2];
325
326 if (mu1 < mu0)
327 return FLT_MAX;
328
329 Double_t sigma0 = par[3];
330 Double_t sigma1 = par[4];
331
332 if (sigma0 < 0.0001)
333 return FLT_MAX;
334
335 if (sigma1 < sigma0)
336 return FLT_MAX;
337
338
339 Double_t mu2 = (2.*mu1)-mu0;
340 Double_t mu3 = (3.*mu1)-(2.*mu0);
341 Double_t mu4 = (4.*mu1)-(3.*mu0);
342 Double_t mu5 = (5.*mu1)-(4.*mu0);
343
344 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
345 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
346 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
347 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
348
349 Double_t lambda2 = lambda*lambda;
350 Double_t lambda3 = lambda2*lambda;
351 Double_t lambda4 = lambda3*lambda;
352 Double_t lambda5 = lambda4*lambda;
353
354 // k=0:
355 arg = (x[0] - mu0)/sigma0;
356 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
357
358 // k=1:
359 arg = (x[0] - mu1)/sigma1;
360 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
361
362 // k=2:
363 arg = (x[0] - mu2)/sigma2;
364 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
365
366 // k=3:
367 arg = (x[0] - mu3)/sigma3;
368 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
369
370 // k=4:
371 arg = (x[0] - mu4)/sigma4;
372 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
373
374 // k=5:
375 arg = (x[0] - mu5)/sigma5;
376 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
377
378 return TMath::Exp(-1.*lambda)*par[5]*sum;
379
380 }
381
382
383 inline static Double_t fPoissonKto6(Double_t *x, Double_t *par)
384 {
385
386 Double_t lambda = par[0];
387
388 Double_t sum = 0.;
389 Double_t arg = 0.;
390
391 Double_t mu0 = par[1];
392 Double_t mu1 = par[2];
393
394 if (mu1 < mu0)
395 return FLT_MAX;
396
397 Double_t sigma0 = par[3];
398 Double_t sigma1 = par[4];
399
400 if (sigma0 < 0.0001)
401 return FLT_MAX;
402
403 if (sigma1 < sigma0)
404 return FLT_MAX;
405
406
407 Double_t mu2 = (2.*mu1)-mu0;
408 Double_t mu3 = (3.*mu1)-(2.*mu0);
409 Double_t mu4 = (4.*mu1)-(3.*mu0);
410 Double_t mu5 = (5.*mu1)-(4.*mu0);
411 Double_t mu6 = (6.*mu1)-(5.*mu0);
412
413 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
414 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
415 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
416 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
417 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
418
419 Double_t lambda2 = lambda*lambda;
420 Double_t lambda3 = lambda2*lambda;
421 Double_t lambda4 = lambda3*lambda;
422 Double_t lambda5 = lambda4*lambda;
423 Double_t lambda6 = lambda5*lambda;
424
425 // k=0:
426 arg = (x[0] - mu0)/sigma0;
427 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
428
429 // k=1:
430 arg = (x[0] - mu1)/sigma1;
431 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
432
433 // k=2:
434 arg = (x[0] - mu2)/sigma2;
435 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
436
437 // k=3:
438 arg = (x[0] - mu3)/sigma3;
439 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
440
441 // k=4:
442 arg = (x[0] - mu4)/sigma4;
443 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
444
445 // k=5:
446 arg = (x[0] - mu5)/sigma5;
447 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
448
449 // k=6:
450 arg = (x[0] - mu6)/sigma6;
451 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
452
453 return TMath::Exp(-1.*lambda)*par[5]*sum;
454
455 }
456
457 inline static Double_t fPolya(Double_t *x, Double_t *par)
458 {
459
460 const Double_t QEcat = 0.247; // mean quantum efficiency
461 const Double_t sqrt2 = 1.4142135623731;
462 const Double_t sqrt3 = 1.7320508075689;
463 const Double_t sqrt4 = 2.;
464
465 const Double_t lambda = par[0]; // mean number of photons
466
467 const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution
468 const Double_t delta1 = par[2]; // amplification first dynode
469 const Double_t delta2 = par[3]; // amplification subsequent dynodes
470
471 const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges
472
473 const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain
474 const Double_t A = 1. + excessPoisson - QEcat
475 + 1./delta1
476 + 1./delta1/delta2
477 + 1./delta1/delta2/delta2; // variance contributions from PMT and QE
478
479 const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion
480
481 const Double_t mu0 = par[7]; // pedestal
482 const Double_t mu1 = totAmpl; // single phe position
483 const Double_t mu2 = 2*totAmpl; // double phe position
484 const Double_t mu3 = 3*totAmpl; // triple phe position
485 const Double_t mu4 = 4*totAmpl; // quadruple phe position
486
487 const Double_t sigma0 = par[5];
488 const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
489 const Double_t sigma2 = sqrt2*sigma1;
490 const Double_t sigma3 = sqrt3*sigma1;
491 const Double_t sigma4 = sqrt4*sigma1;
492
493 const Double_t lambda2 = lambda*lambda;
494 const Double_t lambda3 = lambda2*lambda;
495 const Double_t lambda4 = lambda3*lambda;
496
497 //-- calculate the area----
498 Double_t arg = (x[0] - mu0)/sigma0;
499 Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
500
501 // k=1:
502 arg = (x[0] - mu1)/sigma1;
503 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
504
505 // k=2:
506 arg = (x[0] - mu2)/sigma2;
507 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
508
509 // k=3:
510 arg = (x[0] - mu3)/sigma3;
511 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
512
513 // k=4:
514 arg = (x[0] - mu4)/sigma4;
515 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
516
517 return TMath::Exp(-1.*lambda)*par[6]*sum;
518 }
519
520
521
522 ClassDef(MHCalibrationChargeBlindPix, 1) // Histogram class for Charge Blind Pixel Calibration
523};
524
525#endif /* MARS_MHCalibrationChargeBlindPix */
Note: See TracBrowser for help on using the repository browser.