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

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