source: trunk/Mars/mhcalib/MHCalibrationChargeBlindPix.cc@ 16318

Last change on this file since 16318 was 8911, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 39.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Markus Gaug 02/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2008
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHCalibrationChargeBlindPix
28//
29// Histogram class for the charge calibration of the Blind Pixel.
30// Stores and fits the charges and stores the averaged assumed pedestal and
31// single-phe FADC slice entries. Charges are taken from MExtractedSignalBlindPix.
32// Performs the Single Photo-electron fit to extract the Poisson mean and its errors.
33//
34// Different fits can be chosen with the function SetFitFunc().
35//
36// The fit result is accepted under the condition that:
37// 1) the Probability is greater than fProbLimit (default 0.001 == 99.7%)
38// 2) at least fNumSinglePheLimit events are found in the single Photo-electron peak
39//
40// The single FADC slice entries are averaged and stored in fASinglePheFADCSlices, if
41// their sum exceeds fSinglePheCut, otherwise in fAPedestalFADCSlices.
42//
43// Used numbers are the following:
44//
45// Electronic conversion factor:
46// Assume, we have N_e electrons at the anode,
47// thus a charge of N_e*e (e = electron charge) Coulomb.
48//
49// This charge is AC coupled and runs into a R_pre = 50 Ohm resistency.
50// The corresponding current is amplified by a gain factor G_pre = 400
51// (the precision of this value still has to be checked !!!) and again AC coupled to
52// the output.
53// The corresponding signal goes through the whole transmission and
54// amplification chain and is digitized in the FADCs.
55// The conversion Signal Area to FADC counts (Conv_trans) has been measured
56// by David and Oscar to be approx. 3.9 pVs^-1
57//
58// Thus: Conversion FADC counts to Number of Electrons at Anode:
59// FADC counts = (1/Conv_tran) * G_pre * R_pre * e * N_e = 8 * 10^-4 N_e.
60//
61// Also: FADC counts = 8*10^-4 * GAIN * N_phe
62//
63// In the blind pixel, there is an additional pre-amplifier with an amplification of
64// about 10. Therefore, we have for the blind pixel:
65//
66//
67// FADC counts (Blind Pixel) = 8*10^-3 * GAIN * N_phe
68//
69//////////////////////////////////////////////////////////////////////////////
70#include "MHCalibrationChargeBlindPix.h"
71
72#include <TStyle.h>
73#include <TCanvas.h>
74#include <TPaveText.h>
75#include <TPaveStats.h>
76#include <TGraph.h>
77
78#include <TVector.h>
79#include <TF1.h>
80#include <TH1.h>
81#include <TH2D.h>
82#include <TRandom.h>
83
84#include "MLog.h"
85#include "MLogManip.h"
86
87#include "MParList.h"
88
89#include "MRawEvtPixelIter.h"
90
91ClassImp(MHCalibrationChargeBlindPix);
92
93using namespace std;
94
95const Double_t MHCalibrationChargeBlindPix::gkElectronicAmp = 0.008;
96const Double_t MHCalibrationChargeBlindPix::gkElectronicAmpErr = 0.002;
97const Float_t MHCalibrationChargeBlindPix::gkSignalInitializer = -9999.;
98const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.;
99// --------------------------------------------------------------------------
100//
101// Default Constructor.
102//
103// Sets:
104// - the default number for fNumSinglePheLimit (fgNumSinglePheLimit)
105// - the default number of bins after stripping (30)
106//
107// Initializes:
108// - all pointers to NULL
109// - fASinglePheFADCSlices(0);
110// - fAPedestalFADCSlices(0);
111//
112// Calls:
113// - Clear()
114//
115MHCalibrationChargeBlindPix::MHCalibrationChargeBlindPix(const char *name, const char *title)
116 : fSinglePheFit(NULL),
117 fFitLegend(NULL),
118 fHSinglePheFADCSlices(NULL), fHPedestalFADCSlices(NULL)
119{
120
121 fName = name ? name : "MHCalibrationChargeBlindPix";
122 fTitle = title ? title : "Statistics of the FADC sums of Blind Pixel calibration events";
123
124 //fASinglePheFADCSlices.ResizeTo(1);
125 //fAPedestalFADCSlices.ResizeTo(1);
126
127 SetNumSinglePheLimit();
128 SetProbLimit(0.001);
129 SetBinsAfterStripping(0);
130
131 Clear();
132}
133
134// --------------------------------------------------------------------------
135//
136// Default Destructor.
137//
138// Deletes (if Pointer is not NULL):
139//
140// - fFitLegend
141// - fHSinglePheFADCSlices
142// - fHPedestalFADCSlices
143// - fSinglePheFit
144//
145MHCalibrationChargeBlindPix::~MHCalibrationChargeBlindPix()
146{
147
148 //
149 // The next two lines are important for the case that
150 // the class has been stored to a file and is read again.
151 // In this case, the next two lines prevent a segm. violation
152 // in the destructor
153 //
154 gROOT->GetListOfFunctions()->Remove(fSinglePheFit);
155
156 if (fSinglePheFit)
157 delete fSinglePheFit;
158
159 if (fFitLegend)
160 delete fFitLegend;
161
162 if (fHSinglePheFADCSlices)
163 delete fHSinglePheFADCSlices;
164
165 if (fHPedestalFADCSlices)
166 delete fHPedestalFADCSlices;
167
168}
169
170// --------------------------------------------------------------------------
171//
172// Sets:
173// - all variables to 0., except the fit result variables to gkSignalInitializer
174// - all flags to kFALSE
175// - all pointers to NULL
176// - the default fit function (kEPoisson5)
177//
178// Deletes:
179// - all pointers unequal NULL
180//
181// Calls:
182// - MHCalibrationChargePix::Clear()
183//
184void MHCalibrationChargeBlindPix::Clear(Option_t *o)
185{
186
187 fLambda = gkSignalInitializer;
188 fMu0 = gkSignalInitializer;
189 fMu1 = gkSignalInitializer;
190 fSigma0 = gkSignalInitializer;
191 fSigma1 = gkSignalInitializer;
192 fLambdaErr = gkSignalInitializer;
193 fMu0Err = gkSignalInitializer;
194 fMu1Err = gkSignalInitializer;
195 fSigma0Err = gkSignalInitializer;
196 fSigma1Err = gkSignalInitializer;
197
198 fLambdaCheck = gkSignalInitializer;
199 fLambdaCheckErr = gkSignalInitializer;
200
201 fFitFunc = kEPoisson5;
202
203 fNumSinglePhes = 0;
204 fNumPedestals = 0;
205
206 fChisquare = 0.;
207 fNDF = 0 ;
208 fProb = 0.;
209
210 SetSinglePheFitOK ( kFALSE );
211 SetPedestalFitOK ( kFALSE );
212
213 if (fFitLegend)
214 {
215 delete fFitLegend;
216 fFitLegend = NULL;
217 }
218
219 if (fSinglePheFit)
220 {
221 delete fSinglePheFit;
222 fSinglePheFit = NULL;
223 }
224
225 if (fHSinglePheFADCSlices)
226 {
227 delete fHSinglePheFADCSlices;
228 fHSinglePheFADCSlices = NULL;
229 }
230
231 if (fHPedestalFADCSlices)
232 {
233 delete fHPedestalFADCSlices;
234 fHPedestalFADCSlices = NULL;
235 }
236
237 MHCalibrationPix::Clear();
238 return;
239}
240
241/*
242// --------------------------------------------------------------------------
243//
244// Our own clone function is necessary since root 3.01/06 or Mars 0.4
245// I don't know the reason.
246//
247// Creates new MHCalibrationCam
248//
249TObject *MHCalibrationChargeBlindPix::Clone(const char *) const
250{
251
252 MHCalibrationChargeBlindPix *pix = new MHCalibrationChargeBlindPix();
253 this->Copy(*pix);
254
255 this->fHGausHist.Copy(pix->fHGausHist);
256 this->fSinglePheFit->Copy(*(pix->fSinglePheFit));
257 this->fHSinglePheFADCSlices->Copy(*(pix->fHSinglePheFADCSlices));
258 this->fHPedestalFADCSlices->Copy(*(pix->fHPedestalFADCSlices));
259
260
261 return pix;
262}
263*/
264
265// --------------------------------------------------------------------------
266//
267// Set bit kSinglePheFitOK from outside
268//
269void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b )
270{
271 b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
272}
273
274// --------------------------------------------------------------------------
275//
276// Set bit kPedestalFitOK from outside
277//
278void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b)
279{
280 b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
281}
282
283// --------------------------------------------------------------------------
284//
285// Ask for status of bit kSinglePheFitOK
286//
287const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const
288{
289 return TESTBIT(fFlags,kSinglePheFitOK);
290}
291
292// --------------------------------------------------------------------------
293//
294// Ask for status of bit kPedestalFitOK
295//
296const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const
297{
298 return TESTBIT(fFlags,kPedestalFitOK);
299}
300
301/*
302// --------------------------------------------------------------------------
303//
304// Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices
305//
306void MHCalibrationChargeBlindPix::FinalizeSinglePheSpectrum()
307{
308
309 if (fNumSinglePhes > 1)
310 for (Int_t i=0;i<fASinglePheFADCSlices.GetNrows();i++)
311 fASinglePheFADCSlices[i] = fASinglePheFADCSlices[i]/fNumSinglePhes;
312 if (fNumPedestals > 1)
313 for (Int_t i=0;i<fAPedestalFADCSlices.GetNrows();i++)
314 fAPedestalFADCSlices[i] = fAPedestalFADCSlices[i]/fNumPedestals;
315}
316
317// --------------------------------------------------------------------------
318//
319// Checks again for the size and fills fASinglePheFADCSlices with the FADC slice entries
320//
321void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter)
322{
323
324 const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
325
326 if (fASinglePheFADCSlices.GetNrows() < n)
327 fASinglePheFADCSlices.ResizeTo(n);
328
329 Int_t i=0;
330
331 Byte_t *start = iter.GetHiGainSamples();
332 Byte_t *end = start + iter.GetNumHiGainSamples();
333
334 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
335 fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
336
337 start = iter.GetLoGainSamples();
338 end = start + iter.GetNumLoGainSamples();
339
340 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
341 fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
342
343 fNumSinglePhes++;
344}
345
346// --------------------------------------------------------------------------
347//
348// Checks again for the size and fills fAPedestalFADCSlices with the FADC slice entries
349//
350void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter)
351{
352
353 const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
354
355 if (fAPedestalFADCSlices.GetNrows() < n)
356 fAPedestalFADCSlices.ResizeTo(n);
357
358 Int_t i = 0;
359 Byte_t *start = iter.GetHiGainSamples();
360 Byte_t *end = start + iter.GetNumHiGainSamples();
361
362 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
363 fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
364
365 start = iter.GetLoGainSamples();
366 end = start + iter.GetNumLoGainSamples();
367
368 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
369 fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
370
371 fNumPedestals++;
372}
373*/
374
375// --------------------------------------------------------------------------
376//
377// Task to simulate single phe spectrum with the given parameters
378//
379Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
380{
381
382 gRandom->SetSeed();
383
384 if (fHGausHist.GetIntegral() != 0)
385 {
386 *fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl;
387 *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
388 return kFALSE;
389 }
390
391 if (!InitFit())
392 return kFALSE;
393
394 for (Int_t i=0;i<10000; i++)
395 fHGausHist.Fill(fSinglePheFit->GetRandom());
396
397 return kTRUE;
398}
399
400// --------------------------------------------------------------------------
401//
402// - Get the ranges from the stripped histogram
403// - choose reasonable start values for the fit
404// - initialize the fit function depending on fFitFunc
405// - initialize parameter names and limits depending on fFitFunc
406//
407Bool_t MHCalibrationChargeBlindPix::InitFit()
408{
409
410 //
411 // Get the fitting ranges
412 //
413 Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
414 Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast());
415
416 if (rmin < 0.)
417 rmin = 0.;
418
419 //
420 // First guesses for the fit (should be as close to reality as possible,
421 // otherwise the fit goes gaga because of high number of dimensions ...
422 //
423 const Stat_t entries = fHGausHist.Integral("width");
424 const Double_t lambda_guess = 0.5;
425 //const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
426 const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi());
427
428 //
429 // Initialize the fit function
430 //
431 switch (fFitFunc)
432 {
433 case kEPoisson4:
434 fSinglePheFit = new TF1("SinglePheFit",&PoissonKto4,rmin,rmax,6);
435 rmin += 6.5;
436 break;
437 case kEPoisson5:
438 fSinglePheFit = new TF1("SinglePheFit",&PoissonKto5,rmin,rmax,6);
439 rmin = 0.;
440 break;
441 case kEPoisson6:
442 fSinglePheFit = new TF1("SinglePheFit",&PoissonKto6,rmin,rmax,6);
443 break;
444 case kEPolya:
445 fSinglePheFit = new TF1("SinglePheFit",&Polya,rmin,rmax,8);
446 break;
447 case kEMichele:
448 fSinglePheFit = new TF1("SinglePheFit",&FitFuncMichele,rmin,rmax,9);
449 break;
450 default:
451 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
452 return kFALSE;
453 break;
454 }
455
456 if (!fSinglePheFit)
457 {
458 *fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl;
459 return kFALSE;
460 }
461
462 //
463 // For the fits, we have to take special care since ROOT
464 // has stored the function pointer in a global list which
465 // lead to removing the object twice. We have to take out
466 // the following functions of the global list of functions
467 // as well:
468 //
469 gROOT->GetListOfFunctions()->Remove(fSinglePheFit);
470
471 const Double_t mu_0_guess = 13.5;
472 const Double_t si_0_guess = 2.5;
473 const Double_t mu_1_guess = 30.;
474 const Double_t si_1_guess = si_0_guess + si_0_guess;
475 // Michele
476 const Double_t lambda_1cat_guess = 1.00;
477 const Double_t lambda_1dyn_guess = lambda_1cat_guess/10.;
478 const Double_t mu_1cat_guess = 50.;
479 const Double_t mu_1dyn_guess = 17.;
480 const Double_t si_1cat_guess = si_0_guess + si_0_guess;
481 const Double_t si_1dyn_guess = si_0_guess + si_0_guess/2.;
482 // Polya
483 const Double_t excessPoisson_guess = 0.5;
484 const Double_t delta1_guess = 8.;
485 const Double_t delta2_guess = 5.;
486 const Double_t electronicAmp_guess = gkElectronicAmp;
487 const Double_t electronicAmp_limit = gkElectronicAmpErr;
488
489 //
490 // Initialize boundaries and start parameters
491 //
492 switch (fFitFunc)
493 {
494
495 case kEPoisson4:
496 fSinglePheFit->SetParNames( "#lambda", "#mu_{0}", "#mu_{1}", "#sigma_{0}", "#sigma_{1}","Area");
497 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
498 fSinglePheFit->SetParLimits(0,0.,2.);
499 fSinglePheFit->SetParLimits(1,10.,17.);
500 fSinglePheFit->SetParLimits(2,17.,50.);
501 fSinglePheFit->SetParLimits(3,1.,5.);
502 fSinglePheFit->SetParLimits(4,5.,30.);
503 fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.7*norm));
504 break;
505 case kEPoisson5:
506 case kEPoisson6:
507 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
508 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,800.,si_0_guess,500.,norm);
509 fSinglePheFit->SetParLimits(0,0.,2.);
510 fSinglePheFit->SetParLimits(1,0.,100.);
511 fSinglePheFit->SetParLimits(2,300.,1500.);
512 fSinglePheFit->SetParLimits(3,30.,250.);
513 fSinglePheFit->SetParLimits(4,100.,1000.);
514 fSinglePheFit->SetParLimits(5,norm/1.5,norm*1.5);
515 break;
516
517 case kEPolya:
518 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
519 delta1_guess,delta2_guess,
520 electronicAmp_guess,
521 10.,
522 norm,
523 0.);
524 fSinglePheFit->SetParNames("#lambda","b_{tot}",
525 "#delta_{1}","#delta_{2}",
526 "amp_{e}","#sigma_{0}",
527 "Area", "#mu_{0}");
528 fSinglePheFit->SetParLimits(0,0.,1.);
529 fSinglePheFit->SetParLimits(1,0.,1.);
530 fSinglePheFit->SetParLimits(2,6.,12.);
531 fSinglePheFit->SetParLimits(3,3.,8.);
532 fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
533 electronicAmp_guess+electronicAmp_limit);
534 fSinglePheFit->SetParLimits(5,0.,40.);
535 fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
536 fSinglePheFit->SetParLimits(7,-10.,10.);
537 break;
538 case kEMichele:
539 fSinglePheFit->SetParNames("#lambda_{cat}","#lambda_{dyn}",
540 "#mu_{0}","#mu_{1cat}","#mu_{1dyn}",
541 "#sigma_{0}","#sigma_{1cat}","#sigma_{1dyn}",
542 "Area");
543 fSinglePheFit->SetParameters(lambda_1cat_guess, lambda_1dyn_guess,
544 mu_0_guess, mu_1cat_guess,mu_1dyn_guess,
545 si_0_guess, si_1cat_guess,si_1dyn_guess,
546 norm);
547 fSinglePheFit->SetParLimits(0,0.01,2.0);
548 fSinglePheFit->SetParLimits(1,0.,0.5);
549 fSinglePheFit->SetParLimits(2,10.,16.);
550 fSinglePheFit->SetParLimits(3,25.,50.);
551 fSinglePheFit->SetParLimits(4,16.,18.5);
552 fSinglePheFit->SetParLimits(5,1.,5.);
553 fSinglePheFit->SetParLimits(6,10.,50.);
554 fSinglePheFit->SetParLimits(7,5.,10.);
555 fSinglePheFit->SetParLimits(8,norm/2.,norm*2.5);
556 break;
557
558 default:
559 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
560 return kFALSE;
561 break;
562 }
563
564 fSinglePheFit->SetRange(rmin,rmax);
565
566 return kTRUE;
567}
568
569// --------------------------------------------------------------------------
570//
571// - Retrieve the parameters depending on fFitFunc
572// - Retrieve probability, Chisquare and NDF
573//
574void MHCalibrationChargeBlindPix::ExitFit()
575{
576
577
578 //
579 // Finalize
580 //
581 switch (fFitFunc)
582 {
583
584 case kEPoisson4:
585 case kEPoisson5:
586 case kEPoisson6:
587 case kEPoisson7:
588 fLambda = fSinglePheFit->GetParameter(0);
589 fMu0 = fSinglePheFit->GetParameter(1);
590 fMu1 = fSinglePheFit->GetParameter(2);
591 fSigma0 = fSinglePheFit->GetParameter(3);
592 fSigma1 = fSinglePheFit->GetParameter(4);
593
594 fLambdaErr = fSinglePheFit->GetParError(0);
595 fMu0Err = fSinglePheFit->GetParError(1);
596 fMu1Err = fSinglePheFit->GetParError(2);
597 fSigma0Err = fSinglePheFit->GetParError(3);
598 fSigma1Err = fSinglePheFit->GetParError(4);
599 break;
600 case kEPolya:
601 fLambda = fSinglePheFit->GetParameter(0);
602 fMu0 = fSinglePheFit->GetParameter(7);
603 fMu1 = 0.;
604 fSigma0 = fSinglePheFit->GetParameter(5);
605 fSigma1 = 0.;
606
607 fLambdaErr = fSinglePheFit->GetParError(0);
608 fMu0Err = fSinglePheFit->GetParError(7);
609 fMu1Err = 0.;
610 fSigma0Err = fSinglePheFit->GetParError(5);
611 fSigma1Err = 0.;
612 case kEMichele:
613 fLambda = fSinglePheFit->GetParameter(0);
614 fMu0 = fSinglePheFit->GetParameter(2);
615 fMu1 = fSinglePheFit->GetParameter(3);
616 fSigma0 = fSinglePheFit->GetParameter(5);
617 fSigma1 = fSinglePheFit->GetParameter(6);
618
619 fLambdaErr = fSinglePheFit->GetParError(0);
620 fMu0Err = fSinglePheFit->GetParError(2);
621 fMu1Err = fSinglePheFit->GetParError(3);
622 fSigma0Err = fSinglePheFit->GetParError(5);
623 fSigma1Err = fSinglePheFit->GetParError(6);
624 break;
625 default:
626 break;
627 }
628
629 fProb = fSinglePheFit->GetProb();
630 fChisquare = fSinglePheFit->GetChisquare();
631 fNDF = fSinglePheFit->GetNDF();
632
633 *fLog << all << "Results of the Blind Pixel Fit: " << endl;
634 *fLog << all << "Chisquare: " << fChisquare << endl;
635 *fLog << all << "DoF: " << fNDF << endl;
636 *fLog << all << "Probability: " << fProb << endl;
637
638}
639
640// --------------------------------------------------------------------------
641//
642// - Executes InitFit()
643// - Fits the fHGausHist with fSinglePheFit
644// - Executes ExitFit()
645//
646// The fit result is accepted under condition:
647// 1) The results are not nan's
648// 2) The NDF is not smaller than fNDFLimit (5)
649// 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
650// 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
651//
652Bool_t MHCalibrationChargeBlindPix::FitSinglePhe(Option_t *opt)
653{
654
655 if (!InitFit())
656 return kFALSE;
657
658 fHGausHist.Fit(fSinglePheFit,opt);
659
660 ExitFit();
661
662 //
663 // The fit result is accepted under condition:
664 // 1) The results are not nan's
665 // 2) The NDF is not smaller than fNDFLimit (5)
666 // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
667 // 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
668 //
669 // !Finitite means either infinite or not-a-number
670 if ( !TMath::Finite(fLambda)
671 || !TMath::Finite(fLambdaErr)
672 || !TMath::Finite(fProb)
673 || !TMath::Finite(fMu0)
674 || !TMath::Finite(fMu0Err)
675 || !TMath::Finite(fMu1)
676 || !TMath::Finite(fMu1Err)
677 || !TMath::Finite(fSigma0)
678 || !TMath::Finite(fSigma0Err)
679 || !TMath::Finite(fSigma1)
680 || !TMath::Finite(fSigma1Err)
681 || fNDF < GetNDFLimit()
682 || fProb < GetProbLimit() )
683 return kFALSE;
684
685 const Stat_t entries = fHGausHist.Integral("width");
686 const Float_t numSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
687
688 if (numSinglePhe < fNumSinglePheLimit)
689 {
690 *fLog << warn << "WARNING - Statistics is too low: Only " << numSinglePhe
691 << " in the Single Photo-Electron peak " << endl;
692 return kFALSE;
693 }
694 else
695 *fLog << all << numSinglePhe << " in Single Photo-Electron peak " << endl;
696
697 SetSinglePheFitOK();
698 return kTRUE;
699}
700
701// --------------------------------------------------------------------------
702//
703// - Retrieves limits for the fit
704// - Fits the fHGausHist with Gauss
705// - Retrieves the results to fLambdaCheck and fLambdaCheckErr
706// - Sets a flag IsPedestalFitOK()
707//
708void MHCalibrationChargeBlindPix::FitPedestal (Option_t *opt)
709{
710
711 // Perform the cross-check fitting only the pedestal:
712 const Axis_t rmin = 0.;
713 // const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
714 const Axis_t rmax = fSinglePheCut;
715
716 FitGaus(opt, rmin, rmax);
717
718 const Stat_t entries = fHGausHist.Integral("width");
719 const Double_t pedarea = fFGausFit->Integral(0.,fSinglePheCut);
720
721 fLambdaCheck = TMath::Log(entries/pedarea);
722 // estimate the error by the error of the obtained area from the Gauss-function:
723 fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0);
724
725 SetPedestalFitOK(IsGausFitOK());
726 return;
727}
728
729
730// -------------------------------------------------------------------------
731//
732// Draw a legend with the fit results
733//
734void MHCalibrationChargeBlindPix::DrawLegend(Option_t *opt)
735{
736
737 TString option(opt);
738
739 if (!fFitLegend)
740 {
741 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
742 fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (",
743 (fFitFunc == kEPoisson4) ? "Poisson(k=4))" :
744 (fFitFunc == kEPoisson5) ? "Poisson(k=5))" :
745 (fFitFunc == kEPoisson6) ? "Poisson(k=6))" :
746 (fFitFunc == kEPolya ) ? "Polya(k=4))" :
747 (fFitFunc == kEMichele ) ? "Michele)"
748 : " none )" ));
749 fFitLegend->SetTextSize(0.05);
750 }
751 else
752 fFitLegend->Clear();
753
754 const TString line1 =
755 Form("Mean: #lambda = %2.2f #pm %2.2f",fLambda,fLambdaErr);
756 TText *t1 = fFitLegend->AddText(line1.Data());
757 t1->SetBit(kCanDelete);
758
759 const TString line6 =
760 Form("Mean #lambda_{check} = %2.2f #pm %2.2f",fLambdaCheck,fLambdaCheckErr);
761 TText *t2 = fFitLegend->AddText(line6.Data());
762 t2->SetBit(kCanDelete);
763
764 if (option.Contains("datacheck"))
765 {
766 if (fLambda + 3.*fLambdaErr < fLambdaCheck - 3.*fLambdaCheckErr
767 ||
768 fLambda - 3.*fLambdaErr > fLambdaCheck + 3.*fLambdaCheckErr )
769 {
770 TText *t = fFitLegend->AddText("#lambda and #lambda_{check} more than 3#sigma apart!");
771 t->SetBit(kCanDelete);
772 }
773 }
774 else
775 {
776
777 const TString line2 =
778 Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",fMu0,fMu0Err);
779 TText *t3 = fFitLegend->AddText(line2.Data());
780 t3->SetBit(kCanDelete);
781
782 const TString line3 =
783 Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",fSigma0,fSigma0Err);
784 TText *t4 = fFitLegend->AddText(line3.Data());
785 t4->SetBit(kCanDelete);
786
787 const TString line4 =
788 Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",fMu1,fMu1Err);
789 TText *t5 = fFitLegend->AddText(line4.Data());
790 t5->SetBit(kCanDelete);
791
792 const TString line5 =
793 Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",fSigma1,fSigma1Err);
794 TText *t6 = fFitLegend->AddText(line5.Data());
795 t6->SetBit(kCanDelete);
796 }
797
798 const TString line7 =
799 Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChisquare,fNDF);
800 TText *t7 = fFitLegend->AddText(line7.Data());
801 t7->SetBit(kCanDelete);
802
803 const TString line8 =
804 Form("Probability: %6.4f ",fProb);
805 TText *t8 = fFitLegend->AddText(line8.Data());
806 t8->SetBit(kCanDelete);
807
808 if (IsSinglePheFitOK())
809 {
810 TText *t = fFitLegend->AddText("Result of the Fit: OK");
811 t->SetBit(kCanDelete);
812 }
813 else
814 {
815 TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
816 t->SetBit(kCanDelete);
817 }
818
819 fFitLegend->SetFillColor(IsSinglePheFitOK() ? 80 : 2);
820 fFitLegend->Draw();
821
822 return;
823}
824
825
826// -------------------------------------------------------------------------
827//
828// Draw the histogram
829//
830// The following options can be chosen:
831//
832// "": displays the fHGausHist, the fits, the legend and fASinglePheFADCSlices and fAPedestalFADCSlices
833// "all": executes additionally MHGausEvents::Draw(), with option "fourierevents"
834// "datacheck" display the fHGausHist, the fits and the legend
835//
836void MHCalibrationChargeBlindPix::Draw(Option_t *opt)
837{
838
839 TString option(opt);
840 option.ToLower();
841
842 Int_t win = 1;
843
844 TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600);
845 TVirtualPad *pad = NULL;
846
847 if (option.Contains("all"))
848 {
849 option.ReplaceAll("all","");
850 oldpad->Divide(2,1);
851 win = 2;
852 oldpad->cd(1);
853 TVirtualPad *newpad = gPad;
854 pad = newpad;
855 pad->Divide(2,2);
856 pad->cd(1);
857 }
858 else if (option.Contains("datacheck"))
859 {
860 pad = oldpad;
861 pad->Divide(1,2);
862 pad->cd(1);
863 fHGausHist.SetStats(0);
864 }
865 else
866 {
867 pad = oldpad;
868 pad->Divide(2,2);
869 pad->cd(1);
870 }
871
872 if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow())
873 gPad->SetLogy();
874
875 gPad->SetTicks();
876
877 fHGausHist.Draw();
878 if (fFGausFit )
879 {
880 fFGausFit->SetLineColor(kBlue);
881 fFGausFit->Draw("same");
882 if (!option.Contains("datacheck"))
883 {
884 TLine *line = new TLine(fSinglePheCut, 0., fSinglePheCut, 10.);
885 line->SetBit(kCanDelete);
886 line->SetLineColor(kBlue);
887 line->SetLineWidth(3);
888 line->DrawLine(fSinglePheCut, 0., fSinglePheCut, 2.);
889 }
890 }
891
892 if (fSinglePheFit)
893 {
894 fSinglePheFit->SetFillStyle(0);
895 fSinglePheFit->SetLineWidth(3);
896 fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed);
897 fSinglePheFit->Draw("same");
898 }
899
900 pad->cd(2);
901 DrawLegend(option.Data());
902
903 if (option.Contains("datacheck"))
904 return;
905/*
906 pad->cd(3);
907 if (fASinglePheFADCSlices.GetNrows()!=1)
908 {
909 if (fHSinglePheFADCSlices)
910 delete fHSinglePheFADCSlices;
911
912 fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices);
913 fHSinglePheFADCSlices->SetName("SinglePheFADCSlices");
914 fHSinglePheFADCSlices->SetTitle(Form("%s%4.1f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut));
915 fHSinglePheFADCSlices->SetXTitle("FADC slice number");
916 fHSinglePheFADCSlices->SetYTitle("FADC counts");
917 const Int_t nbins = fHSinglePheFADCSlices->GetNbinsX();
918 TH2D *nulls = new TH2D("Nulls",fHSinglePheFADCSlices->GetTitle(),nbins,0.,
919 fHSinglePheFADCSlices->GetXaxis()->GetBinCenter(nbins),
920 100,0.,50.);
921 nulls->SetDirectory(NULL);
922 nulls->SetBit(kCanDelete);
923 nulls->GetXaxis()->SetTitle(fHSinglePheFADCSlices->GetXaxis()->GetTitle());
924 nulls->GetYaxis()->SetTitle(fHSinglePheFADCSlices->GetYaxis()->GetTitle());
925 nulls->GetXaxis()->CenterTitle();
926 nulls->GetYaxis()->CenterTitle();
927 nulls->SetStats(0);
928 nulls->Draw();
929 fHSinglePheFADCSlices->Draw("same");
930 }
931
932 pad->cd(4);
933 if (fAPedestalFADCSlices.GetNrows()!=1)
934 {
935
936 if (fHPedestalFADCSlices)
937 delete fHPedestalFADCSlices;
938
939 fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices);
940 fHPedestalFADCSlices->SetName("PedestalFADCSlices");
941 fHPedestalFADCSlices->SetTitle(Form("%s%4.1f","Pedestal FADC Slices, Sum < ",fSinglePheCut));
942 fHPedestalFADCSlices->SetXTitle("FADC slice number");
943 fHPedestalFADCSlices->SetYTitle("FADC counts");
944 const Int_t nbins = fHPedestalFADCSlices->GetNbinsX();
945 TH2D *nullp = new TH2D("Nullp",fHPedestalFADCSlices->GetTitle(),nbins,0.,
946 fHPedestalFADCSlices->GetXaxis()->GetBinCenter(nbins),
947 100,0.,50.);
948 nullp->SetDirectory(NULL);
949 nullp->SetBit(kCanDelete);
950 nullp->GetXaxis()->SetTitle(fHPedestalFADCSlices->GetXaxis()->GetTitle());
951 nullp->GetYaxis()->SetTitle(fHPedestalFADCSlices->GetYaxis()->GetTitle());
952 nullp->GetXaxis()->CenterTitle();
953 nullp->GetYaxis()->CenterTitle();
954 nullp->SetStats(0);
955 nullp->Draw();
956 fHPedestalFADCSlices->Draw("same");
957 }
958 */
959 if (win < 2)
960 return;
961
962 oldpad->cd(2);
963 MHCalibrationPix::Draw("fourierevents");
964}
965
966Double_t MHCalibrationChargeBlindPix::FitFuncMichele(Double_t *x, Double_t *par)
967{
968 Double_t lambda1cat = par[0];
969 Double_t lambda1dyn = par[1];
970 Double_t mu0 = par[2];
971 Double_t mu1cat = par[3];
972 Double_t mu1dyn = par[4];
973 Double_t sigma0 = par[5];
974 Double_t sigma1cat = par[6];
975 Double_t sigma1dyn = par[7];
976
977 Double_t sumcat = 0.;
978 Double_t sumdyn = 0.;
979 Double_t arg = 0.;
980
981 if (lambda1cat < lambda1dyn)
982 return FLT_MAX;
983
984 if (mu1cat < mu0)
985 return FLT_MAX;
986
987 if (mu1dyn < mu0)
988 return FLT_MAX;
989
990 if (mu1cat < mu1dyn)
991 return FLT_MAX;
992
993 if (sigma0 < 0.0001)
994 return FLT_MAX;
995
996 if (sigma1cat < sigma0)
997 return FLT_MAX;
998
999 if (sigma1dyn < sigma0)
1000 return FLT_MAX;
1001
1002 Double_t mu2cat = (2.*mu1cat)-mu0;
1003 Double_t mu2dyn = (2.*mu1dyn)-mu0;
1004 Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
1005 Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
1006
1007 Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0));
1008 Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0));
1009 Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
1010 Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
1011
1012 Double_t lambda2cat = lambda1cat*lambda1cat;
1013 Double_t lambda2dyn = lambda1dyn*lambda1dyn;
1014 Double_t lambda3cat = lambda2cat*lambda1cat;
1015 Double_t lambda3dyn = lambda2dyn*lambda1dyn;
1016
1017 // k=0:
1018 arg = (x[0] - mu0)/sigma0;
1019 sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
1020 sumdyn = sumcat;
1021
1022 // k=1cat:
1023 arg = (x[0] - mu1cat)/sigma1cat;
1024 sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
1025 // k=1dyn:
1026 arg = (x[0] - mu1dyn)/sigma1dyn;
1027 sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
1028
1029 // k=2cat:
1030 arg = (x[0] - mu2cat)/sigma2cat;
1031 sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
1032 // k=2dyn:
1033 arg = (x[0] - mu2dyn)/sigma2dyn;
1034 sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
1035
1036
1037 // k=3cat:
1038 arg = (x[0] - mu3cat)/sigma3cat;
1039 sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
1040 // k=3dyn:
1041 arg = (x[0] - mu3dyn)/sigma3dyn;
1042 sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn;
1043
1044 sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
1045 sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
1046
1047 return par[8]*(sumcat+sumdyn)/2.;
1048}
1049
1050Double_t MHCalibrationChargeBlindPix::PoissonKto4(Double_t *x, Double_t *par)
1051{
1052 Double_t lambda = par[0];
1053
1054 Double_t sum = 0.;
1055 Double_t arg = 0.;
1056
1057 Double_t mu0 = par[1];
1058 Double_t mu1 = par[2];
1059
1060 if (mu1 < mu0)
1061 return FLT_MAX;
1062
1063 Double_t sigma0 = par[3];
1064 Double_t sigma1 = par[4];
1065
1066 if (sigma0 < 0.0001)
1067 return FLT_MAX;
1068
1069 if (sigma1 < sigma0)
1070 return FLT_MAX;
1071
1072 Double_t mu2 = (2.*mu1)-mu0;
1073 Double_t mu3 = (3.*mu1)-(2.*mu0);
1074 Double_t mu4 = (4.*mu1)-(3.*mu0);
1075
1076 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
1077 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
1078 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
1079
1080 Double_t lambda2 = lambda*lambda;
1081 Double_t lambda3 = lambda2*lambda;
1082 Double_t lambda4 = lambda3*lambda;
1083
1084 // k=0:
1085 arg = (x[0] - mu0)/sigma0;
1086 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
1087
1088 // k=1:
1089 arg = (x[0] - mu1)/sigma1;
1090 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
1091
1092 // k=2:
1093 arg = (x[0] - mu2)/sigma2;
1094 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
1095
1096 // k=3:
1097 arg = (x[0] - mu3)/sigma3;
1098 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
1099
1100 // k=4:
1101 arg = (x[0] - mu4)/sigma4;
1102 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
1103
1104 return TMath::Exp(-1.*lambda)*par[5]*sum;
1105}
1106
1107Double_t MHCalibrationChargeBlindPix::PoissonKto5(Double_t *x, Double_t *par)
1108{
1109 Double_t lambda = par[0];
1110
1111 Double_t sum = 0.;
1112 Double_t arg = 0.;
1113
1114 Double_t mu0 = par[1];
1115 Double_t mu1 = par[2];
1116
1117 if (mu1 < mu0)
1118 return FLT_MAX;
1119
1120 Double_t sigma0 = par[3];
1121 Double_t sigma1 = par[4];
1122
1123 if (sigma0 < 0.0001)
1124 return FLT_MAX;
1125
1126 if (sigma1 < sigma0)
1127 return FLT_MAX;
1128
1129
1130 Double_t mu2 = (2.*mu1)-mu0;
1131 Double_t mu3 = (3.*mu1)-(2.*mu0);
1132 Double_t mu4 = (4.*mu1)-(3.*mu0);
1133 Double_t mu5 = (5.*mu1)-(4.*mu0);
1134
1135 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
1136 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
1137 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
1138 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
1139
1140 Double_t lambda2 = lambda*lambda;
1141 Double_t lambda3 = lambda2*lambda;
1142 Double_t lambda4 = lambda3*lambda;
1143 Double_t lambda5 = lambda4*lambda;
1144
1145 // k=0:
1146 arg = (x[0] - mu0)/sigma0;
1147 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
1148
1149 // k=1:
1150 arg = (x[0] - mu1)/sigma1;
1151 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
1152
1153 // k=2:
1154 arg = (x[0] - mu2)/sigma2;
1155 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
1156
1157 // k=3:
1158 arg = (x[0] - mu3)/sigma3;
1159 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
1160
1161 // k=4:
1162 arg = (x[0] - mu4)/sigma4;
1163 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
1164
1165 // k=5:
1166 arg = (x[0] - mu5)/sigma5;
1167 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
1168
1169 return TMath::Exp(-1.*lambda)*par[5]*sum;
1170}
1171
1172Double_t MHCalibrationChargeBlindPix::PoissonKto6(Double_t *x, Double_t *par)
1173{
1174 Double_t lambda = par[0];
1175
1176 Double_t sum = 0.;
1177 Double_t arg = 0.;
1178
1179 Double_t mu0 = par[1];
1180 Double_t mu1 = par[2];
1181
1182 if (mu1 < mu0)
1183 return FLT_MAX;
1184
1185 Double_t sigma0 = par[3];
1186 Double_t sigma1 = par[4];
1187
1188 if (sigma0 < 0.0001)
1189 return FLT_MAX;
1190
1191 if (sigma1 < sigma0)
1192 return FLT_MAX;
1193
1194
1195 Double_t mu2 = (2.*mu1)-mu0;
1196 Double_t mu3 = (3.*mu1)-(2.*mu0);
1197 Double_t mu4 = (4.*mu1)-(3.*mu0);
1198 Double_t mu5 = (5.*mu1)-(4.*mu0);
1199 Double_t mu6 = (6.*mu1)-(5.*mu0);
1200
1201 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
1202 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
1203 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
1204 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
1205 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
1206
1207 Double_t lambda2 = lambda*lambda;
1208 Double_t lambda3 = lambda2*lambda;
1209 Double_t lambda4 = lambda3*lambda;
1210 Double_t lambda5 = lambda4*lambda;
1211 Double_t lambda6 = lambda5*lambda;
1212
1213 // k=0:
1214 arg = (x[0] - mu0)/sigma0;
1215 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
1216
1217 // k=1:
1218 arg = (x[0] - mu1)/sigma1;
1219 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
1220
1221 // k=2:
1222 arg = (x[0] - mu2)/sigma2;
1223 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
1224
1225 // k=3:
1226 arg = (x[0] - mu3)/sigma3;
1227 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
1228
1229 // k=4:
1230 arg = (x[0] - mu4)/sigma4;
1231 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
1232
1233 // k=5:
1234 arg = (x[0] - mu5)/sigma5;
1235 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
1236
1237 // k=6:
1238 arg = (x[0] - mu6)/sigma6;
1239 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
1240
1241 return TMath::Exp(-1.*lambda)*par[5]*sum;
1242}
1243
1244Double_t MHCalibrationChargeBlindPix::Polya(Double_t *x, Double_t *par)
1245{
1246 const Double_t QEcat = 0.247; // mean quantum efficiency
1247 const Double_t sqrt2 = 1.4142135623731;
1248 const Double_t sqrt3 = 1.7320508075689;
1249 const Double_t sqrt4 = 2.;
1250
1251 const Double_t lambda = par[0]; // mean number of photons
1252
1253 const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution
1254 const Double_t delta1 = par[2]; // amplification first dynode
1255 const Double_t delta2 = par[3]; // amplification subsequent dynodes
1256
1257 const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges
1258
1259 const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain
1260 const Double_t A = 1. + excessPoisson - QEcat
1261 + 1./delta1
1262 + 1./delta1/delta2
1263 + 1./delta1/delta2/delta2; // variance contributions from PMT and QE
1264
1265 const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion
1266
1267 const Double_t mu0 = par[7]; // pedestal
1268 const Double_t mu1 = totAmpl; // single phe position
1269 const Double_t mu2 = 2*totAmpl; // double phe position
1270 const Double_t mu3 = 3*totAmpl; // triple phe position
1271 const Double_t mu4 = 4*totAmpl; // quadruple phe position
1272
1273 const Double_t sigma0 = par[5];
1274 const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
1275 const Double_t sigma2 = sqrt2*sigma1;
1276 const Double_t sigma3 = sqrt3*sigma1;
1277 const Double_t sigma4 = sqrt4*sigma1;
1278
1279 const Double_t lambda2 = lambda*lambda;
1280 const Double_t lambda3 = lambda2*lambda;
1281 const Double_t lambda4 = lambda3*lambda;
1282
1283 //-- calculate the area----
1284 Double_t arg = (x[0] - mu0)/sigma0;
1285 Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
1286
1287 // k=1:
1288 arg = (x[0] - mu1)/sigma1;
1289 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
1290
1291 // k=2:
1292 arg = (x[0] - mu2)/sigma2;
1293 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
1294
1295 // k=3:
1296 arg = (x[0] - mu3)/sigma3;
1297 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
1298
1299 // k=4:
1300 arg = (x[0] - mu4)/sigma4;
1301 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
1302
1303 return TMath::Exp(-1.*lambda)*par[6]*sum;
1304}
Note: See TracBrowser for help on using the repository browser.