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

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