source: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc@ 3440

Last change on this file since 3440 was 3323, checked in by gaug, 22 years ago
*** empty log message ***
File size: 24.7 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-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHCalibrationChargeBlindPix
28//
29// Performs all the Single Photo-Electron Fit to extract
30// the mean number of photons and to derive the light flux
31//
32// The fit result is accepted under condition that:
33// 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
34// 2) at least 100 events are in the single Photo-electron peak
35//
36// Used numbers are the following:
37//
38// Electronic conversion factor:
39// Assume, we have N_e electrons at the anode,
40// thus a charge of N_e*e (e = electron charge) Coulomb.
41//
42// This charge is AC coupled and runs into a R_pre = 50 Ohm resistency.
43// The corresponding current is amplified by a gain factor G_pre = 400
44// (the precision of this value still has to be checked !!!) and again AC coupled to
45// the output.
46// The corresponding signal goes through the whole transmission and
47// amplification chain and is digitized in the FADCs.
48// The conversion Signal Area to FADC counts (Conv_trans) has been measured
49// by David and Oscar to be approx. 3.9 pVs^-1
50//
51// Thus: Conversion FADC counts to Number of Electrons at Anode:
52// FADC counts = (1/Conv_tran) * G_pre * R_pre * e * N_e = 8 * 10^-4 N_e.
53//
54// Also: FADC counts = 8*10^-4 * GAIN * N_phe
55//
56// In the blind pixel, there is an additional pre-amplifier with an amplification of
57// about 10. Therefore, we have for the blind pixel:
58//
59// FADC counts (Blind Pixel) = 8*10^-3 * GAIN * N_phe
60//
61//////////////////////////////////////////////////////////////////////////////
62#include "MHCalibrationChargeBlindPix.h"
63
64#include <TStyle.h>
65#include <TCanvas.h>
66#include <TPaveText.h>
67
68#include <TVector.h>
69#include <TF1.h>
70#include <TH1.h>
71#include <TRandom.h>
72
73#include "MLog.h"
74#include "MLogManip.h"
75
76#include "MParList.h"
77
78#include "MRawEvtData.h"
79#include "MRawEvtPixelIter.h"
80
81#include "MExtractedSignalBlindPixel.h"
82#include "MCalibrationChargeBlindPix.h"
83
84ClassImp(MHCalibrationChargeBlindPix);
85
86using namespace std;
87
88const Double_t MHCalibrationChargeBlindPix::gkElectronicAmp = 0.008;
89const Double_t MHCalibrationChargeBlindPix::gkElectronicAmpErr = 0.002;
90
91const Int_t MHCalibrationChargeBlindPix::fgChargeNbins = 1400;
92const Axis_t MHCalibrationChargeBlindPix::fgChargeFirst = -200.5;
93const Axis_t MHCalibrationChargeBlindPix::fgChargeLast = 1199.5;
94
95const Float_t MHCalibrationChargeBlindPix::fgSinglePheCut = 200.;
96const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.;
97// --------------------------------------------------------------------------
98//
99// Default Constructor.
100//
101MHCalibrationChargeBlindPix::MHCalibrationChargeBlindPix(const char *name, const char *title)
102 : fBlindPix(NULL), fSignal(NULL), fRawEvt(NULL),
103 fASinglePheFADCSlices(30), fAPedestalFADCSlices(30),
104 fSinglePheFit(NULL),
105 fFitLegend(NULL),
106 fHSinglePheFADCSlices(NULL), fHPedestalFADCSlices(NULL)
107{
108
109 fName = name ? name : "MHCalibrationChargeBlindPix";
110 fTitle = title ? title : "Fill the accumulated charges and times of all Blind Pixel events and perform fits";
111
112 SetChargeNbins();
113 SetChargeFirst();
114 SetChargeLast();
115
116 SetSinglePheCut();
117 SetNumSinglePheLimit();
118
119 fHGausHist.SetName("HCalibrationChargeBlindPix");
120 fHGausHist.SetTitle("Distribution of Summed FADC slices Blind Pixel");
121 fHGausHist.SetXTitle("Sum FADC Slices");
122 fHGausHist.SetYTitle("Nr. of events");
123
124 Clear();
125}
126
127MHCalibrationChargeBlindPix::~MHCalibrationChargeBlindPix()
128{
129
130 if (fSinglePheFit)
131 delete fSinglePheFit;
132
133 if (fFitLegend)
134 delete fFitLegend;
135
136 if (fHSinglePheFADCSlices)
137 delete fHSinglePheFADCSlices;
138
139 if (fHPedestalFADCSlices)
140 delete fHPedestalFADCSlices;
141
142}
143
144void MHCalibrationChargeBlindPix::Init()
145{
146
147 fHGausHist.SetBins( fChargeNbins, fChargeFirst, fChargeLast);
148}
149
150void MHCalibrationChargeBlindPix::Clear(Option_t *o)
151{
152
153 fLambda = -999.;
154 fMu0 = -999.;
155 fMu1 = -999.;
156 fSigma0 = -999.;
157 fSigma1 = -999.;
158
159 fLambdaErr = -999.;
160 fMu0Err = -999.;
161 fMu1Err = -999.;
162 fSigma0Err = -999.;
163 fSigma1Err = -999.;
164
165 fLambdaCheck = -999.;
166 fLambdaCheckErr = -999.;
167
168 fMeanPedestal = 0.;
169 fMeanPedestalErr = 0.;
170 fSigmaPedestal = 0.;
171 fSigmaPedestalErr = 0.;
172
173 fFitFunc = kEPoisson5;
174
175 fExtractSlices = 0;
176 fNumSinglePhes = 0;
177 fNumPedestals = 0;
178
179 fNumSinglePhes = 0;
180 fNumPedestals = 0;
181
182 fChisquare = 0.;
183 fNDF = 0 ;
184 fProb = 0.;
185
186 SetSinglePheFitOK ( kFALSE );
187 SetPedestalFitOK ( kFALSE );
188
189 if (fFitLegend)
190 {
191 delete fFitLegend;
192 fFitLegend = NULL;
193 }
194
195 if (fSinglePheFit)
196 {
197 delete fSinglePheFit;
198 fSinglePheFit = NULL;
199 }
200
201 if (fHSinglePheFADCSlices)
202 {
203 delete fHSinglePheFADCSlices;
204 fHSinglePheFADCSlices = NULL;
205 }
206
207 if (fHPedestalFADCSlices)
208 {
209 delete fHPedestalFADCSlices;
210 fHPedestalFADCSlices = NULL;
211 }
212
213
214 MHCalibrationChargePix::Clear();
215 return;
216}
217
218void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b)
219{
220 b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
221}
222
223void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b)
224{
225 b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
226}
227
228const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const
229{
230 return TESTBIT(fFlags,kSinglePheFitOK);
231}
232
233const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const
234{
235 return TESTBIT(fFlags,kPedestalFitOK);
236}
237
238Bool_t MHCalibrationChargeBlindPix::SetupFill(const MParList *pList)
239{
240 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
241 if (!fRawEvt)
242 {
243 *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
244 return kFALSE;
245 }
246
247 fSignal = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel");
248 if (!fSignal)
249 {
250 *fLog << err << GetDescriptor() << ": ERROR: Could not find MExtractedSignalBlindPixel ... aborting " << endl;
251 return kFALSE;
252 }
253
254 Init();
255
256 return kTRUE;
257}
258
259Bool_t MHCalibrationChargeBlindPix::ReInit(MParList *pList)
260{
261
262 fBlindPix = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
263 if (!fBlindPix)
264 {
265 *fLog << err << GetDescriptor() << ": ERROR: Could not find MCalibrationChargeBlindPix ... aborting " << endl;
266 return kFALSE;
267 }
268
269 return kTRUE;
270}
271
272Bool_t MHCalibrationChargeBlindPix::Fill(const MParContainer *par, const Stat_t w)
273{
274
275 Float_t slices = (Float_t)fSignal->GetNumFADCSamples();
276
277 if (slices == 0.)
278 {
279 *fLog << err << "Number of used signal slices in MExtractedSignalBlindPix is zero ... abort."
280 << endl;
281 return kFALSE;
282 }
283
284 if (fExtractSlices != 0. && slices != fExtractSlices )
285 {
286 *fLog << err << "Number of used signal slices changed in MExtractedSignalCam ... abort."
287 << endl;
288 return kFALSE;
289 }
290 fExtractSlices = slices;
291
292 //
293 // Signal extraction and histogram filling
294 //
295 const Float_t signal = (Float_t)fSignal->GetExtractedSignal();
296 FillHistAndArray(signal);
297
298 //
299 // IN order to study the single-phe posistion, we extract the slices
300 //
301 MRawEvtPixelIter pixel(fRawEvt);
302 pixel.Jump(fSignal->GetBlindPixelIdx());
303
304 if (signal > fSinglePheCut)
305 FillSinglePheFADCSlices(pixel);
306 else
307 FillPedestalFADCSlices(pixel);
308
309 return kTRUE;
310}
311
312Bool_t MHCalibrationChargeBlindPix::Finalize()
313{
314
315 if (IsEmpty())
316 {
317 *fLog << err << GetDescriptor() << ": My histogram has not been filled !! " << endl;
318 return kFALSE;
319 }
320
321 CreateFourierSpectrum();
322 fBlindPix->SetOscillating ( !IsFourierSpectrumOK() );
323
324 fMeanPedestal = fSignal->GetPed();
325 fMeanPedestalErr = fSignal->GetPedErr();
326 fSigmaPedestal = fSignal->GetPedRms();
327 fSigmaPedestalErr = fSignal->GetPedRmsErr();
328
329 if (fNumSinglePhes > 1)
330 for (Int_t i=0;i<fASinglePheFADCSlices.GetNrows();i++)
331 fASinglePheFADCSlices[i] = fASinglePheFADCSlices[i]/fNumSinglePhes;
332 if (fNumPedestals > 1)
333 for (Int_t i=0;i<fAPedestalFADCSlices.GetNrows();i++)
334 fAPedestalFADCSlices[i] = fAPedestalFADCSlices[i]/fNumPedestals;
335
336 FitPedestal();
337
338 if (FitSinglePhe())
339 fBlindPix->SetSinglePheFitOK();
340
341 fBlindPix->SetLambda ( fLambda );
342 fBlindPix->SetMu0 ( fMu0 );
343 fBlindPix->SetMu0Err ( fMu0Err );
344 fBlindPix->SetMu1 ( fMu1 );
345 fBlindPix->SetMu1Err ( fMu1Err );
346 fBlindPix->SetSigma0 ( fSigma0 );
347 fBlindPix->SetSigma0Err ( fSigma0Err );
348 fBlindPix->SetSigma1 ( fSigma1 );
349 fBlindPix->SetSigma1Err ( fSigma1Err );
350 fBlindPix->SetProb ( fProb );
351
352 fBlindPix->SetLambdaCheck ( fLambdaCheck );
353 fBlindPix->SetLambdaCheckErr ( fLambdaCheckErr );
354
355 return kTRUE;
356}
357
358void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter)
359{
360
361 const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
362
363 if (fASinglePheFADCSlices.GetNrows() < n)
364 fASinglePheFADCSlices.ResizeTo(n);
365
366 Int_t i=0;
367
368 Byte_t *start = iter.GetHiGainSamples();
369 Byte_t *end = start + iter.GetNumHiGainSamples();
370
371 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
372 fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
373
374 start = iter.GetLoGainSamples();
375 end = start + iter.GetNumLoGainSamples();
376
377 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
378 fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
379
380 fNumSinglePhes++;
381}
382
383void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter)
384{
385
386 const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
387
388 if (fAPedestalFADCSlices.GetNrows() < n)
389 fAPedestalFADCSlices.ResizeTo(n);
390
391 Int_t i = 0;
392 Byte_t *start = iter.GetHiGainSamples();
393 Byte_t *end = start + iter.GetNumHiGainSamples();
394
395 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
396 fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
397
398 start = iter.GetLoGainSamples();
399 end = start + iter.GetNumLoGainSamples();
400
401 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
402 fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
403
404 fNumPedestals++;
405}
406
407
408
409Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
410{
411
412 gRandom->SetSeed();
413
414 if (fHGausHist.GetIntegral() != 0)
415 {
416 *fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl;
417 *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
418 return kFALSE;
419 }
420
421 if (!InitFit())
422 return kFALSE;
423
424 for (Int_t i=0;i<10000; i++)
425 fHGausHist.Fill(fSinglePheFit->GetRandom());
426
427 return kTRUE;
428}
429
430Bool_t MHCalibrationChargeBlindPix::InitFit()
431{
432
433 //
434 // Get the fitting ranges
435 //
436 Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
437 Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast());
438
439 if (rmin < 0.)
440 rmin = 0.;
441
442 //
443 // First guesses for the fit (should be as close to reality as possible,
444 // otherwise the fit goes gaga because of high number of dimensions ...
445 //
446 const Stat_t entries = fHGausHist.Integral("width");
447 const Double_t lambda_guess = 0.1;
448 const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
449 const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi());
450
451 //
452 // Initialize the fit function
453 //
454 switch (fFitFunc)
455 {
456 case kEPoisson4:
457 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6);
458 break;
459 case kEPoisson5:
460 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6);
461 break;
462 case kEPoisson6:
463 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6);
464 break;
465 case kEPolya:
466 fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8);
467 break;
468 case kEMichele:
469 break;
470
471 default:
472 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
473 return kFALSE;
474 break;
475 }
476
477 if (!fSinglePheFit)
478 {
479 *fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl;
480 return kFALSE;
481 }
482
483 const Double_t mu_0_guess = maximum_bin;
484 const Double_t si_0_guess = 40.;
485 const Double_t mu_1_guess = mu_0_guess + 100.;
486 const Double_t si_1_guess = si_0_guess + si_0_guess;
487 // Michele
488// const Double_t lambda_1cat_guess = 0.5;
489 // const Double_t lambda_1dyn_guess = 0.5;
490 // const Double_t mu_1cat_guess = mu_0_guess + 50.;
491 // const Double_t mu_1dyn_guess = mu_0_guess + 20.;
492 // const Double_t si_1cat_guess = si_0_guess + si_0_guess;
493 // const Double_t si_1dyn_guess = si_0_guess;
494 // Polya
495 const Double_t excessPoisson_guess = 0.5;
496 const Double_t delta1_guess = 8.;
497 const Double_t delta2_guess = 5.;
498 const Double_t electronicAmp_guess = gkElectronicAmp;
499 const Double_t electronicAmp_limit = gkElectronicAmpErr;
500
501 //
502 // Initialize boundaries and start parameters
503 //
504 switch (fFitFunc)
505 {
506
507 case kEPoisson4:
508 fSinglePheFit->SetParNames( "#lambda", "#mu_{0}", "#mu_{1}", "#sigma_{0}", "#sigma_{1}","Area");
509 fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm);
510
511 fSinglePheFit->SetParLimits(0,0.,0.5);
512 fSinglePheFit->SetParLimits(1,
513 fMeanPedestal-5.*fMeanPedestalErr,
514 fMeanPedestal+5.*fMeanPedestalErr);
515 fSinglePheFit->SetParLimits(2,rmin,rmax);
516 fSinglePheFit->SetParLimits(3,
517 fSigmaPedestal-5.*fSigmaPedestalErr,
518 fSigmaPedestal+5.*fSigmaPedestalErr);
519 fSinglePheFit->SetParLimits(4,0.,(rmax-rmin));
520 fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.5*norm));
521 break;
522 case kEPoisson5:
523 case kEPoisson6:
524 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
525 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
526 fSinglePheFit->SetParLimits(0,0.,1.);
527 fSinglePheFit->SetParLimits(1,rmin,(rmax-rmin)/1.5);
528 fSinglePheFit->SetParLimits(2,(rmax-rmin)/2.,(rmax-0.05*(rmax-rmin)));
529 fSinglePheFit->SetParLimits(3,1.0,(rmax-rmin)/2.0);
530 fSinglePheFit->SetParLimits(4,1.0,(rmax-rmin)/2.5);
531 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
532 break;
533
534 case kEPolya:
535 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
536 delta1_guess,delta2_guess,
537 electronicAmp_guess,
538 fSigmaPedestal,
539 norm,
540 fMeanPedestal);
541 fSinglePheFit->SetParNames("#lambda","b_{tot}",
542 "#delta_{1}","#delta_{2}",
543 "amp_{e}","#sigma_{0}",
544 "Area", "#mu_{0}");
545 fSinglePheFit->SetParLimits(0,0.,1.);
546 fSinglePheFit->SetParLimits(1,0.,1.);
547 fSinglePheFit->SetParLimits(2,6.,12.);
548 fSinglePheFit->SetParLimits(3,3.,8.);
549 fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
550 electronicAmp_guess+electronicAmp_limit);
551 fSinglePheFit->SetParLimits(5,
552 fSigmaPedestal-3.*fSigmaPedestalErr,
553 fSigmaPedestal+3.*fSigmaPedestalErr);
554 fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
555 fSinglePheFit->SetParLimits(7,
556 fMeanPedestal-3.*fMeanPedestalErr,
557 fMeanPedestal+3.*fMeanPedestalErr);
558 break;
559 case kEMichele:
560 break;
561
562 default:
563 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
564 return kFALSE;
565 break;
566 }
567
568 fSinglePheFit->SetRange(rmin,rmax);
569
570 return kTRUE;
571}
572
573void MHCalibrationChargeBlindPix::ExitFit()
574{
575
576
577 //
578 // Finalize
579 //
580 switch (fFitFunc)
581 {
582
583 case kEPoisson4:
584 case kEPoisson5:
585 case kEPoisson6:
586 case kEPoisson7:
587 fLambda = fSinglePheFit->GetParameter(0);
588 fMu0 = fSinglePheFit->GetParameter(1);
589 fMu1 = fSinglePheFit->GetParameter(2);
590 fSigma0 = fSinglePheFit->GetParameter(3);
591 fSigma1 = fSinglePheFit->GetParameter(4);
592
593 fLambdaErr = fSinglePheFit->GetParError(0);
594 fMu0Err = fSinglePheFit->GetParError(1);
595 fMu1Err = fSinglePheFit->GetParError(2);
596 fSigma0Err = fSinglePheFit->GetParError(3);
597 fSigma1Err = fSinglePheFit->GetParError(4);
598 break;
599 case kEPolya:
600 fLambda = fSinglePheFit->GetParameter(0);
601 fMu0 = fSinglePheFit->GetParameter(7);
602 fMu1 = 0.;
603 fSigma0 = fSinglePheFit->GetParameter(5);
604 fSigma1 = 0.;
605
606 fLambdaErr = fSinglePheFit->GetParError(0);
607 fMu0Err = fSinglePheFit->GetParError(7);
608 fMu1Err = 0.;
609 fSigma0Err = fSinglePheFit->GetParError(5);
610 fSigma1Err = 0.;
611 default:
612 break;
613 }
614
615 fProb = fSinglePheFit->GetProb();
616 fChisquare = fSinglePheFit->GetChisquare();
617 fNDF = fSinglePheFit->GetNDF();
618
619 *fLog << all << "Results of the Blind Pixel Fit: " << endl;
620 *fLog << all << "Chisquare: " << fChisquare << endl;
621 *fLog << all << "DoF: " << fNDF << endl;
622 *fLog << all << "Probability: " << fProb << endl;
623
624}
625
626
627Bool_t MHCalibrationChargeBlindPix::FitSinglePhe(Option_t *opt)
628{
629
630 if (!InitFit())
631 return kFALSE;
632
633 fHGausHist.Fit(fSinglePheFit,opt);
634
635 ExitFit();
636
637 //
638 // The fit result is accepted under condition:
639 // 1) The results are not nan's
640 // 2) The NDF is not smaller than fNDFLimit (5)
641 // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
642 // 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
643 if ( TMath::IsNaN(fLambda)
644 || TMath::IsNaN(fLambdaErr)
645 || TMath::IsNaN(fProb)
646 || TMath::IsNaN(fMu0)
647 || TMath::IsNaN(fMu0Err)
648 || TMath::IsNaN(fMu1)
649 || TMath::IsNaN(fMu1Err)
650 || TMath::IsNaN(fSigma0)
651 || TMath::IsNaN(fSigma0Err)
652 || TMath::IsNaN(fSigma1)
653 || TMath::IsNaN(fSigma1Err)
654 || fNDF < fNDFLimit
655 || fProb < fProbLimit )
656 return kFALSE;
657
658 const Stat_t entries = fHGausHist.Integral("width");
659 const Float_t numSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
660
661 if (numSinglePhe < fNumSinglePheLimit)
662 {
663 *fLog << warn << "WARNING - Statistics is too low: Only " << numSinglePhe
664 << " in the Single Photo-Electron peak " << endl;
665 return kFALSE;
666 }
667 else
668 *fLog << all << numSinglePhe << " in Single Photo-Electron peak " << endl;
669
670 SetSinglePheFitOK();
671 return kTRUE;
672}
673
674void MHCalibrationChargeBlindPix::FitPedestal (Option_t *opt)
675{
676
677 // Perform the cross-check fitting only the pedestal:
678 const Axis_t rmin = 0.;
679 const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
680
681 FitGaus(opt, rmin, rmax);
682
683 const Stat_t entries = fHGausHist.Integral("width");
684 const Double_t fitarea = fFGausFit->GetParameter(0);
685 const Double_t pedarea = fitarea * TMath::Sqrt(TMath::TwoPi()) * fFGausFit->GetParameter(2);
686
687 fLambdaCheck = TMath::Log(entries/pedarea);
688 fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0)
689 + fFGausFit->GetParError(2)/fFGausFit->GetParameter(2);
690
691
692 SetPedestalFitOK();
693 return;
694}
695
696
697// -------------------------------------------------------------------------
698//
699// Draw a legend with the fit results
700//
701void MHCalibrationChargeBlindPix::DrawLegend()
702{
703
704 if (!fFitLegend)
705 {
706 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
707 fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (",
708 (fFitFunc = kEPoisson4) ? "Poisson(k=4))" :
709 (fFitFunc = kEPoisson5) ? "Poisson(k=5))" :
710 (fFitFunc = kEPoisson6) ? "Poisson(k=4))" :
711 (fFitFunc = kEPolya ) ? "Polya(k=4))" :
712 (fFitFunc = kEMichele ) ? "Michele)" : " none )" ));
713 fFitLegend->SetTextSize(0.05);
714 }
715 else
716 fFitLegend->Clear();
717
718 const TString line1 =
719 Form("Mean: #lambda = %2.2f #pm %2.2f",fLambda,fLambdaErr);
720 TText *t1 = fFitLegend->AddText(line1.Data());
721 t1->SetBit(kCanDelete);
722
723 const TString line6 =
724 Form("Mean #lambda (check) = %2.2f #pm %2.2f",fLambdaCheck,fLambdaCheckErr);
725 TText *t2 = fFitLegend->AddText(line6.Data());
726 t2->SetBit(kCanDelete);
727
728 const TString line2 =
729 Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",fMu0,fMu0Err);
730 TText *t3 = fFitLegend->AddText(line2.Data());
731 t3->SetBit(kCanDelete);
732
733 const TString line3 =
734 Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",fSigma0,fSigma0Err);
735 TText *t4 = fFitLegend->AddText(line3.Data());
736 t4->SetBit(kCanDelete);
737
738 const TString line4 =
739 Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",fMu1,fMu1Err);
740 TText *t5 = fFitLegend->AddText(line4.Data());
741 t5->SetBit(kCanDelete);
742
743 const TString line5 =
744 Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",fSigma1,fSigma1Err);
745 TText *t6 = fFitLegend->AddText(line5.Data());
746 t6->SetBit(kCanDelete);
747
748 const TString line7 =
749 Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChisquare,fNDF);
750 TText *t7 = fFitLegend->AddText(line7.Data());
751 t7->SetBit(kCanDelete);
752
753 const TString line8 =
754 Form("Probability: %4.2f ",fProb);
755 TText *t8 = fFitLegend->AddText(line8.Data());
756 t8->SetBit(kCanDelete);
757
758 if (IsSinglePheFitOK())
759 {
760 TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
761 t->SetBit(kCanDelete);
762 }
763 else
764 {
765 TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
766 t->SetBit(kCanDelete);
767 }
768
769 fFitLegend->SetFillColor(IsSinglePheFitOK() ? 80 : 2);
770 fFitLegend->Draw();
771
772 return;
773}
774
775
776// -------------------------------------------------------------------------
777//
778// Draw the histogram
779//
780void MHCalibrationChargeBlindPix::Draw(Option_t *opt)
781{
782
783 TString option(opt);
784 option.ToLower();
785
786 Int_t win = 1;
787
788 TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600);
789 TVirtualPad *pad = NULL;
790
791 oldpad->SetBorderMode(0);
792
793 if (option.Contains("all"))
794 {
795 option.ReplaceAll("all","");
796 oldpad->Divide(2,1);
797 win = 2;
798 oldpad->cd(1);
799 TVirtualPad *newpad = gPad;
800 pad = newpad;
801 pad->Divide(2,2);
802 pad->cd(1);
803 }
804 else
805 {
806 pad = oldpad;
807 pad->Divide(2,2);
808 pad->cd(1);
809 }
810
811 if (!IsEmpty())
812 gPad->SetLogy();
813
814 gPad->SetTicks();
815
816 fHGausHist.Draw(opt);
817 if (fFGausFit)
818 {
819 fFGausFit->SetLineColor(kBlue);
820 fFGausFit->Draw("same");
821 }
822 if (fSinglePheFit)
823 {
824 fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed);
825 fSinglePheFit->Draw("same");
826 }
827
828 pad->cd(2);
829 DrawLegend();
830
831 pad->cd(3);
832 if (fHSinglePheFADCSlices)
833 delete fHSinglePheFADCSlices;
834
835 fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices);
836 fHSinglePheFADCSlices->SetName("SinglePheFADCSlices");
837 fHSinglePheFADCSlices->SetTitle(Form("%s%f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut));
838 fHSinglePheFADCSlices->SetXTitle("FADC slice number");
839 fHSinglePheFADCSlices->SetYTitle("FADC counts");
840 fHSinglePheFADCSlices->Draw(opt);
841
842 pad->cd(4);
843 if (fHPedestalFADCSlices)
844 delete fHPedestalFADCSlices;
845
846 fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices);
847 fHPedestalFADCSlices->SetName("PedestalFADCSlices");
848 fHPedestalFADCSlices->SetTitle(Form("%s%f","Pedestal FADC Slices, Sum < ",fSinglePheCut));
849 fHPedestalFADCSlices->SetXTitle("FADC slice number");
850 fHPedestalFADCSlices->SetYTitle("FADC counts");
851 fHPedestalFADCSlices->Draw(opt);
852
853 if (win < 2)
854 return;
855
856 oldpad->cd(2);
857 MHGausEvents::Draw("fourierevents");
858}
859
860
861
862
863
864
865
866
867
868
869
Note: See TracBrowser for help on using the repository browser.