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

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