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

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