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

Last change on this file since 4220 was 4220, checked in by gaug, 20 years ago
*** empty log message ***
File size: 32.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// 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 ChangeFitFunc().
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// FADC counts (Blind Pixel) = 8*10^-3 * GAIN * N_phe
67//
68//////////////////////////////////////////////////////////////////////////////
69#include "MHCalibrationChargeBlindPix.h"
70
71#include <TStyle.h>
72#include <TCanvas.h>
73#include <TPaveText.h>
74
75#include <TVector.h>
76#include <TF1.h>
77#include <TH1.h>
78#include <TRandom.h>
79
80#include "MLog.h"
81#include "MLogManip.h"
82
83#include "MParList.h"
84
85#include "MRawEvtData.h"
86#include "MRawEvtPixelIter.h"
87
88#include "MExtractedSignalBlindPixel.h"
89#include "MCalibrationChargeBlindPix.h"
90
91ClassImp(MHCalibrationChargeBlindPix);
92
93using namespace std;
94
95const Double_t MHCalibrationChargeBlindPix::gkElectronicAmp = 0.008;
96const Double_t MHCalibrationChargeBlindPix::gkElectronicAmpErr = 0.002;
97
98const Int_t MHCalibrationChargeBlindPix::fgChargeNbins = 1000;
99const Axis_t MHCalibrationChargeBlindPix::fgChargeFirst = -0.5;
100const Axis_t MHCalibrationChargeBlindPix::fgChargeLast = 5999.5;
101const Float_t MHCalibrationChargeBlindPix::fgSinglePheCut = 600.;
102const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.;
103// --------------------------------------------------------------------------
104//
105// Default Constructor.
106//
107// Sets:
108// - the default number for fNbins (fgChargeNbins)
109// - the default number for fFirst (fgChargeFirst)
110// - the default number for fLast (fgChargeLast)
111// - the default number for fSinglePheCut (fgSingePheCut)
112// - the default number for fNumSinglePheLimit (fgNumSinglePheLimit)
113// - the default number of bins after stripping (30)
114//
115// - the default name of the fHGausHist ("HCalibrationChargeBlindPix")
116// - the default title of the fHGausHist ("Distribution of Summed FADC slices Blind Pixel ")
117// - the default x-axis title for fHGausHist ("Sum FADC Slices")
118// - the default y-axis title for fHGausHist ("Nr. of events")
119//
120// Initializes:
121// - all pointers to NULL
122// - fASinglePheFADCSlices(0);
123// - fAPedestalFADCSlices(0);
124//
125// Calls:
126// - Clear()
127//
128MHCalibrationChargeBlindPix::MHCalibrationChargeBlindPix(const char *name, const char *title)
129 : fBlindPix(NULL), fSignal(NULL), fRawEvt(NULL),
130 fSinglePheFit(NULL),
131 fFitLegend(NULL),
132 fHSinglePheFADCSlices(NULL), fHPedestalFADCSlices(NULL)
133{
134
135 fName = name ? name : "MHCalibrationChargeBlindPix";
136 fTitle = title ? title : "Statistics of the FADC sums of Blind Pixel calibration events";
137
138 SetNbins( fgChargeNbins );
139 SetFirst( fgChargeFirst );
140 SetLast ( fgChargeLast );
141
142 fASinglePheFADCSlices.ResizeTo(1);
143 fAPedestalFADCSlices.ResizeTo(1);
144
145 SetSinglePheCut();
146 SetNumSinglePheLimit();
147 SetProbLimit(0.00001);
148 SetBinsAfterStripping(50);
149
150 fHGausHist.SetName("HCalibrationChargeBlindPix");
151 fHGausHist.SetTitle("Distribution of Summed FADC slices Blind Pixel");
152 fHGausHist.SetXTitle("Sum FADC Slices");
153 fHGausHist.SetYTitle("Nr. of events");
154
155 Clear();
156}
157
158// --------------------------------------------------------------------------
159//
160// Default Destructor.
161//
162// Deletes (if Pointer is not NULL):
163//
164// - fSinglePheFit
165// - fFitLegend
166// - fHSinglePheFADCSlices
167// - fHPedestalFADCSlices
168//
169MHCalibrationChargeBlindPix::~MHCalibrationChargeBlindPix()
170{
171
172 if (fSinglePheFit)
173 delete fSinglePheFit;
174
175 if (fFitLegend)
176 delete fFitLegend;
177
178 if (fHSinglePheFADCSlices)
179 delete fHSinglePheFADCSlices;
180
181 if (fHPedestalFADCSlices)
182 delete fHPedestalFADCSlices;
183
184}
185
186// --------------------------------------------------------------------------
187//
188// Sets:
189// - all variables to 0., except the fit result variables to -999.
190// - all flags to kFALSE
191// - all pointers to NULL
192// - the default fit function (kEPoisson5)
193//
194// Deletes:
195// - all pointers unequal NULL
196//
197// Calls:
198// - MHCalibrationChargePix::Clear()
199//
200void MHCalibrationChargeBlindPix::Clear(Option_t *o)
201{
202
203 fLambda = -999.;
204 fMu0 = -999.;
205 fMu1 = -999.;
206 fSigma0 = -999.;
207 fSigma1 = -999.;
208
209 fLambdaErr = -999.;
210 fMu0Err = -999.;
211 fMu1Err = -999.;
212 fSigma0Err = -999.;
213 fSigma1Err = -999.;
214
215 fLambdaCheck = -999.;
216 fLambdaCheckErr = -999.;
217
218 fFitFunc = kEMichele;
219 // fFitFunc = kEPoisson4;
220
221 fNumSinglePhes = 0;
222 fNumPedestals = 0;
223
224 fChisquare = 0.;
225 fNDF = 0 ;
226 fProb = 0.;
227
228 SetSinglePheFitOK ( kFALSE );
229 SetPedestalFitOK ( kFALSE );
230
231 if (fFitLegend)
232 {
233 delete fFitLegend;
234 fFitLegend = NULL;
235 }
236
237 if (fSinglePheFit)
238 {
239 delete fSinglePheFit;
240 fSinglePheFit = NULL;
241 }
242
243 if (fHSinglePheFADCSlices)
244 {
245 delete fHSinglePheFADCSlices;
246 fHSinglePheFADCSlices = NULL;
247 }
248
249 if (fHPedestalFADCSlices)
250 {
251 delete fHPedestalFADCSlices;
252 fHPedestalFADCSlices = NULL;
253 }
254
255
256 MHCalibrationChargePix::Clear();
257 return;
258}
259
260// --------------------------------------------------------------------------
261//
262// Set bit kSinglePheFitOK from outside
263//
264void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b)
265{
266 b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
267}
268
269// --------------------------------------------------------------------------
270//
271// Set bit kPedestalFitOK from outside
272//
273void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b)
274{
275 b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
276}
277
278// --------------------------------------------------------------------------
279//
280// Ask for status of bit kSinglePheFitOK
281//
282const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const
283{
284 return TESTBIT(fFlags,kSinglePheFitOK);
285}
286
287// --------------------------------------------------------------------------
288//
289// Ask for status of bit kPedestalFitOK
290//
291const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const
292{
293 return TESTBIT(fFlags,kPedestalFitOK);
294}
295
296// --------------------------------------------------------------------------
297//
298// Gets the pointers to:
299// - MRawEvtData
300// - MExtractedSignalBlindPixel
301//
302// Calls:
303// - MHGausHist::InitBins()
304//
305Bool_t MHCalibrationChargeBlindPix::SetupFill(const MParList *pList)
306{
307
308 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
309 if (!fRawEvt)
310 {
311 *fLog << err << "MRawEvtData not found... aborting." << endl;
312 return kFALSE;
313 }
314
315 fSignal = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel");
316 if (!fSignal)
317 {
318 *fLog << err << "MExtractedSignalBlindPixel not found... aborting " << endl;
319 return kFALSE;
320 }
321
322 InitBins();
323
324 return kTRUE;
325}
326
327// --------------------------------------------------------------------------
328//
329// Gets or creates the pointers to:
330// - MCalibrationChargeBlindPix
331//
332Bool_t MHCalibrationChargeBlindPix::ReInit(MParList *pList)
333{
334
335 fBlindPix = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
336 if (!fBlindPix)
337 return kFALSE;
338
339 return kTRUE;
340}
341
342// --------------------------------------------------------------------------
343//
344// Retrieves from MExtractedSignalBlindPixel:
345// - number of FADC samples
346// - extracted signal
347// - blind Pixel ID
348//
349// Resizes (if necessary):
350// - fASinglePheFADCSlices to sum of HiGain and LoGain samples
351// - fAPedestalFADCSlices to sum of HiGain and LoGain samples
352//
353// Fills the following histograms:
354// - MHGausEvents::FillHistAndArray(signal)
355//
356// Creates MRawEvtPixelIter, jumps to blind pixel ID,
357// fills the vectors fASinglePheFADCSlices and fAPedestalFADCSlices
358// with the full FADC slices, depending on the size of the signal w.r.t. fSinglePheCut
359//
360Bool_t MHCalibrationChargeBlindPix::Fill(const MParContainer *par, const Stat_t w)
361{
362
363 const Int_t samples = (Int_t)fRawEvt->GetNumHiGainSamples()+(Int_t)fRawEvt->GetNumLoGainSamples();
364
365 if (!fASinglePheFADCSlices.IsValid())
366 {
367 fASinglePheFADCSlices.ResizeTo(samples);
368 fAPedestalFADCSlices.ResizeTo(samples);
369 }
370
371 if (fASinglePheFADCSlices.GetNrows() != samples)
372 {
373 fASinglePheFADCSlices.ResizeTo(samples);
374 fAPedestalFADCSlices.ResizeTo(samples);
375 }
376
377 Float_t slices = (Float_t)fSignal->GetNumFADCSamples();
378
379 if (slices == 0.)
380 {
381 *fLog << err << "Number of used signal slices in MExtractedSignalBlindPix is zero ... abort."
382 << endl;
383 return kFALSE;
384 }
385
386 //
387 // Signal extraction and histogram filling
388 //
389 const Float_t signal = (Float_t)fSignal->GetExtractedSignal();
390 if (signal > -0.5)
391 FillHistAndArray(signal);
392 else
393 return kTRUE;
394
395 //
396 // IN order to study the single-phe posistion, we extract the slices
397 //
398 const Int_t blindpixIdx = fSignal->GetBlindPixelIdx();
399
400 MRawEvtPixelIter pixel(fRawEvt);
401 pixel.Jump(blindpixIdx);
402
403 if (signal > fSinglePheCut)
404 FillSinglePheFADCSlices(pixel);
405 else
406 FillPedestalFADCSlices(pixel);
407
408 return kTRUE;
409}
410
411// --------------------------------------------------------------------------
412//
413// Returns kFALSE, if empty
414//
415// - Creates the fourier spectrum and sets bit MHGausEvents::IsFourierSpectrumOK()
416// - Retrieves the pedestals from MExtractedSignalBlindPixel
417// - Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices
418// - Executes FitPedestal()
419// - Executes FitSinglePhe()
420// - Retrieves fit results and stores them in MCalibrationChargeBlindPix
421//
422Bool_t MHCalibrationChargeBlindPix::Finalize()
423{
424
425 if (IsEmpty())
426 {
427 *fLog << err << GetDescriptor() << ": My histogram has not been filled !! " << endl;
428 return kFALSE;
429 }
430
431 CreateFourierSpectrum();
432 fBlindPix->SetOscillating ( !IsFourierSpectrumOK() );
433
434 fBlindPix->SetValid(kTRUE);
435
436 fMeanPedestal = fSignal->GetPed();
437 fMeanPedestalErr = fSignal->GetPedErr();
438 fSigmaPedestal = fSignal->GetPedRms();
439 fSigmaPedestalErr = fSignal->GetPedRmsErr();
440
441 if (fNumSinglePhes > 1)
442 for (Int_t i=0;i<fASinglePheFADCSlices.GetNrows();i++)
443 fASinglePheFADCSlices[i] = fASinglePheFADCSlices[i]/fNumSinglePhes;
444 if (fNumPedestals > 1)
445 for (Int_t i=0;i<fAPedestalFADCSlices.GetNrows();i++)
446 fAPedestalFADCSlices[i] = fAPedestalFADCSlices[i]/fNumPedestals;
447
448 FitPedestal();
449
450 if (FitSinglePhe())
451 fBlindPix->SetSinglePheFitOK();
452 else
453 fBlindPix->SetValid(IsPedestalFitOK());
454
455 fBlindPix->SetLambda ( fLambdaCheck );
456 fBlindPix->SetLambdaVar ( fLambdaCheckErr*fLambdaCheckErr );
457 fBlindPix->SetMu0 ( fMu0 );
458 fBlindPix->SetMu0Err ( fMu0Err );
459 fBlindPix->SetMu1 ( fMu1 );
460 fBlindPix->SetMu1Err ( fMu1Err );
461 fBlindPix->SetSigma0 ( fSigma0 );
462 fBlindPix->SetSigma0Err ( fSigma0Err );
463 fBlindPix->SetSigma1 ( fSigma1 );
464 fBlindPix->SetSigma1Err ( fSigma1Err );
465 fBlindPix->SetProb ( fProb );
466
467 fBlindPix->SetLambdaCheck ( fLambda );
468 fBlindPix->SetLambdaCheckErr ( fLambdaErr );
469
470 return kTRUE;
471}
472
473
474// --------------------------------------------------------------------------
475//
476// Checks again for the size and fills fASinglePheFADCSlices with the FADC slice entries
477//
478void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter)
479{
480
481 const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
482
483 if (fASinglePheFADCSlices.GetNrows() < n)
484 fASinglePheFADCSlices.ResizeTo(n);
485
486 Int_t i=0;
487
488 Byte_t *start = iter.GetHiGainSamples();
489 Byte_t *end = start + iter.GetNumHiGainSamples();
490
491 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
492 fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
493
494 start = iter.GetLoGainSamples();
495 end = start + iter.GetNumLoGainSamples();
496
497 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
498 fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr;
499
500 fNumSinglePhes++;
501}
502
503// --------------------------------------------------------------------------
504//
505// Checks again for the size and fills fAPedestalFADCSlices with the FADC slice entries
506//
507void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter)
508{
509
510 const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples();
511
512 if (fAPedestalFADCSlices.GetNrows() < n)
513 fAPedestalFADCSlices.ResizeTo(n);
514
515 Int_t i = 0;
516 Byte_t *start = iter.GetHiGainSamples();
517 Byte_t *end = start + iter.GetNumHiGainSamples();
518
519 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
520 fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
521
522 start = iter.GetLoGainSamples();
523 end = start + iter.GetNumLoGainSamples();
524
525 for (Byte_t *ptr = start; ptr < end; ptr++, i++)
526 fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr;
527
528 fNumPedestals++;
529}
530
531
532// --------------------------------------------------------------------------
533//
534// Task to simulate single phe spectrum with the given parameters
535//
536Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
537{
538
539 gRandom->SetSeed();
540
541 if (fHGausHist.GetIntegral() != 0)
542 {
543 *fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl;
544 *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
545 return kFALSE;
546 }
547
548 if (!InitFit())
549 return kFALSE;
550
551 for (Int_t i=0;i<10000; i++)
552 fHGausHist.Fill(fSinglePheFit->GetRandom());
553
554 return kTRUE;
555}
556
557// --------------------------------------------------------------------------
558//
559// - Get the ranges from the stripped histogram
560// - choose reasonable start values for the fit
561// - initialize the fit function depending on fFitFunc
562// - initialize parameter names and limits depending on fFitFunc
563//
564Bool_t MHCalibrationChargeBlindPix::InitFit()
565{
566
567 //
568 // Get the fitting ranges
569 //
570 Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
571 Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast());
572
573 if (rmin < 0.)
574 rmin = 0.;
575
576 //
577 // First guesses for the fit (should be as close to reality as possible,
578 // otherwise the fit goes gaga because of high number of dimensions ...
579 //
580 const Stat_t entries = fHGausHist.Integral("width");
581 const Double_t lambda_guess = 0.1;
582 const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
583 const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi());
584
585 //
586 // Initialize the fit function
587 //
588 switch (fFitFunc)
589 {
590 case kEPoisson4:
591 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6);
592 break;
593 case kEPoisson5:
594 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6);
595 break;
596 case kEPoisson6:
597 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6);
598 break;
599 case kEPolya:
600 fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8);
601 break;
602 case kEMichele:
603 fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,10);
604 break;
605 default:
606 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
607 return kFALSE;
608 break;
609 }
610
611 if (!fSinglePheFit)
612 {
613 *fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl;
614 return kFALSE;
615 }
616
617 const Double_t mu_0_guess = maximum_bin;
618 const Double_t si_0_guess = 40.;
619 const Double_t mu_1_guess = mu_0_guess + 4000.;
620 const Double_t si_1_guess = si_0_guess + si_0_guess;
621 // Michele
622 const Double_t lambda_1cat_guess = 0.05;
623 const Double_t lambda_1dyn_guess = lambda_1cat_guess/10.;
624 const Double_t mu_1cat_guess = 1000.;
625 const Double_t mu_1dyn_guess = 2500.;
626 const Double_t si_1cat_guess = si_0_guess+ 500.;
627 const Double_t si_1dyn_guess = si_0_guess+ 1000.;
628 const Double_t offset_guess = 0.5;
629 // Polya
630 const Double_t excessPoisson_guess = 0.5;
631 const Double_t delta1_guess = 8.;
632 const Double_t delta2_guess = 5.;
633 const Double_t electronicAmp_guess = gkElectronicAmp;
634 const Double_t electronicAmp_limit = gkElectronicAmpErr;
635
636 //
637 // Initialize boundaries and start parameters
638 //
639 switch (fFitFunc)
640 {
641
642 case kEPoisson4:
643 fSinglePheFit->SetParNames( "#lambda", "#mu_{0}", "#mu_{1}", "#sigma_{0}", "#sigma_{1}","Area");
644 fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm);
645
646 fSinglePheFit->SetParLimits(0,0.,0.5);
647 fSinglePheFit->SetParLimits(1,
648 fMeanPedestal-5.*fMeanPedestalErr,
649 fMeanPedestal+5.*fMeanPedestalErr);
650 fSinglePheFit->SetParLimits(2,rmin,rmax);
651 fSinglePheFit->SetParLimits(3,
652 fSigmaPedestal-5.*fSigmaPedestalErr,
653 fSigmaPedestal+5.*fSigmaPedestalErr);
654 fSinglePheFit->SetParLimits(4,0.,(rmax-rmin));
655 fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.5*norm));
656 break;
657 case kEPoisson5:
658 case kEPoisson6:
659 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
660 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
661 fSinglePheFit->SetParLimits(0,0.,1.);
662 fSinglePheFit->SetParLimits(1,rmin,(rmax-rmin)/1.5);
663 fSinglePheFit->SetParLimits(2,(rmax-rmin)/2.,(rmax-0.05*(rmax-rmin)));
664 fSinglePheFit->SetParLimits(3,1.0,(rmax-rmin)/2.0);
665 fSinglePheFit->SetParLimits(4,1.0,(rmax-rmin)/2.5);
666 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
667 break;
668
669 case kEPolya:
670 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
671 delta1_guess,delta2_guess,
672 electronicAmp_guess,
673 fSigmaPedestal,
674 norm,
675 fMeanPedestal);
676 fSinglePheFit->SetParNames("#lambda","b_{tot}",
677 "#delta_{1}","#delta_{2}",
678 "amp_{e}","#sigma_{0}",
679 "Area", "#mu_{0}");
680 fSinglePheFit->SetParLimits(0,0.,1.);
681 fSinglePheFit->SetParLimits(1,0.,1.);
682 fSinglePheFit->SetParLimits(2,6.,12.);
683 fSinglePheFit->SetParLimits(3,3.,8.);
684 fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
685 electronicAmp_guess+electronicAmp_limit);
686 fSinglePheFit->SetParLimits(5,
687 fSigmaPedestal-3.*fSigmaPedestalErr,
688 fSigmaPedestal+3.*fSigmaPedestalErr);
689 fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
690 fSinglePheFit->SetParLimits(7,
691 fMeanPedestal-3.*fMeanPedestalErr,
692 fMeanPedestal+3.*fMeanPedestalErr);
693 break;
694 case kEMichele:
695 fSinglePheFit->SetParameters(lambda_1cat_guess, lambda_1dyn_guess,
696 20., mu_1cat_guess,mu_1dyn_guess,
697 si_0_guess, si_1cat_guess,si_1dyn_guess,
698 norm, offset_guess);
699 fSinglePheFit->SetParNames("#lambda_{cat}","#lambda_{dyn}",
700 "#mu_{0}","#mu_{1cat}","#mu_{1dyn}",
701 "#sigma_{0}","#sigma_{1cat}","#sigma_{1dyn}",
702 "Area","offset");
703 fSinglePheFit->SetParLimits(0,0.,0.5);
704 fSinglePheFit->SetParLimits(1,0.,0.05);
705 fSinglePheFit->SetParLimits(2,0.,fSinglePheCut);
706 fSinglePheFit->SetParLimits(3,fSinglePheCut,5500.);
707 fSinglePheFit->SetParLimits(4,1500.,3500.);
708 fSinglePheFit->SetParLimits(5,0.,fSinglePheCut);
709 fSinglePheFit->SetParLimits(6,500.,1000.);
710 fSinglePheFit->SetParLimits(7,300.,1500.);
711 fSinglePheFit->SetParLimits(8,norm/1.1,norm*1.1);
712 fSinglePheFit->SetParLimits(9,0.,1.);
713 break;
714
715 default:
716 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
717 return kFALSE;
718 break;
719 }
720
721 fSinglePheFit->SetRange(rmin,rmax);
722
723 return kTRUE;
724}
725
726// --------------------------------------------------------------------------
727//
728// - Retrieve the parameters depending on fFitFunc
729// - Retrieve probability, Chisquare and NDF
730//
731void MHCalibrationChargeBlindPix::ExitFit()
732{
733
734
735 //
736 // Finalize
737 //
738 switch (fFitFunc)
739 {
740
741 case kEPoisson4:
742 case kEPoisson5:
743 case kEPoisson6:
744 case kEPoisson7:
745 fLambda = fSinglePheFit->GetParameter(0);
746 fMu0 = fSinglePheFit->GetParameter(1);
747 fMu1 = fSinglePheFit->GetParameter(2);
748 fSigma0 = fSinglePheFit->GetParameter(3);
749 fSigma1 = fSinglePheFit->GetParameter(4);
750
751 fLambdaErr = fSinglePheFit->GetParError(0);
752 fMu0Err = fSinglePheFit->GetParError(1);
753 fMu1Err = fSinglePheFit->GetParError(2);
754 fSigma0Err = fSinglePheFit->GetParError(3);
755 fSigma1Err = fSinglePheFit->GetParError(4);
756 break;
757 case kEPolya:
758 fLambda = fSinglePheFit->GetParameter(0);
759 fMu0 = fSinglePheFit->GetParameter(7);
760 fMu1 = 0.;
761 fSigma0 = fSinglePheFit->GetParameter(5);
762 fSigma1 = 0.;
763
764 fLambdaErr = fSinglePheFit->GetParError(0);
765 fMu0Err = fSinglePheFit->GetParError(7);
766 fMu1Err = 0.;
767 fSigma0Err = fSinglePheFit->GetParError(5);
768 fSigma1Err = 0.;
769 case kEMichele:
770 fLambda = fSinglePheFit->GetParameter(0);
771 fMu0 = fSinglePheFit->GetParameter(2);
772 fMu1 = fSinglePheFit->GetParameter(3);
773 fSigma0 = fSinglePheFit->GetParameter(5);
774 fSigma1 = fSinglePheFit->GetParameter(6);
775
776 fLambdaErr = fSinglePheFit->GetParError(0);
777 fMu0Err = fSinglePheFit->GetParError(2);
778 fMu1Err = fSinglePheFit->GetParError(3);
779 fSigma0Err = fSinglePheFit->GetParError(5);
780 fSigma1Err = fSinglePheFit->GetParError(6);
781 break;
782 default:
783 break;
784 }
785
786 fProb = fSinglePheFit->GetProb();
787 fChisquare = fSinglePheFit->GetChisquare();
788 fNDF = fSinglePheFit->GetNDF();
789
790 *fLog << all << "Results of the Blind Pixel Fit: " << endl;
791 *fLog << all << "Chisquare: " << fChisquare << endl;
792 *fLog << all << "DoF: " << fNDF << endl;
793 *fLog << all << "Probability: " << fProb << endl;
794
795}
796
797// --------------------------------------------------------------------------
798//
799// - Executes InitFit()
800// - Fits the fHGausHist with fSinglePheFit
801// - Executes ExitFit()
802//
803// The fit result is accepted under condition:
804// 1) The results are not nan's
805// 2) The NDF is not smaller than fNDFLimit (5)
806// 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
807// 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
808//
809Bool_t MHCalibrationChargeBlindPix::FitSinglePhe(Option_t *opt)
810{
811
812 if (!InitFit())
813 return kFALSE;
814
815 fHGausHist.Fit(fSinglePheFit,opt);
816
817 ExitFit();
818
819 //
820 // The fit result is accepted under condition:
821 // 1) The results are not nan's
822 // 2) The NDF is not smaller than fNDFLimit (5)
823 // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
824 // 4) at least fNumSinglePheLimit events are in the single Photo-electron peak
825 //
826 if ( TMath::IsNaN(fLambda)
827 || TMath::IsNaN(fLambdaErr)
828 || TMath::IsNaN(fProb)
829 || TMath::IsNaN(fMu0)
830 || TMath::IsNaN(fMu0Err)
831 || TMath::IsNaN(fMu1)
832 || TMath::IsNaN(fMu1Err)
833 || TMath::IsNaN(fSigma0)
834 || TMath::IsNaN(fSigma0Err)
835 || TMath::IsNaN(fSigma1)
836 || TMath::IsNaN(fSigma1Err)
837 || fNDF < fNDFLimit
838 || fProb < fProbLimit )
839 return kFALSE;
840
841 const Stat_t entries = fHGausHist.Integral("width");
842 const Float_t numSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
843
844 if (numSinglePhe < fNumSinglePheLimit)
845 {
846 *fLog << warn << "WARNING - Statistics is too low: Only " << numSinglePhe
847 << " in the Single Photo-Electron peak " << endl;
848 return kFALSE;
849 }
850 else
851 *fLog << all << numSinglePhe << " in Single Photo-Electron peak " << endl;
852
853 SetSinglePheFitOK();
854 return kTRUE;
855}
856
857// --------------------------------------------------------------------------
858//
859// - Retrieves limits for the fit
860// - Fits the fHGausHist with Gauss
861// - Retrieves the results to fLambdaCheck and fLambdaCheckErr
862// - Sets a flag IsPedestalFitOK()
863//
864void MHCalibrationChargeBlindPix::FitPedestal (Option_t *opt)
865{
866
867 // Perform the cross-check fitting only the pedestal:
868 const Axis_t rmin = 0.;
869 // const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
870 const Axis_t rmax = fSinglePheCut;
871
872 FitGaus(opt, rmin, rmax);
873
874 const Stat_t entries = fHGausHist.Integral("width");
875 const Double_t pedarea = fFGausFit->Integral(0.,fSinglePheCut);
876
877 fLambdaCheck = TMath::Log(entries/pedarea);
878 // estimate the error by the error of the obtained area from the Gauss-function:
879 fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0);
880
881 SetPedestalFitOK(IsGausFitOK());
882 return;
883}
884
885
886// -------------------------------------------------------------------------
887//
888// Draw a legend with the fit results
889//
890void MHCalibrationChargeBlindPix::DrawLegend()
891{
892
893 if (!fFitLegend)
894 {
895 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
896 fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (",
897 (fFitFunc == kEPoisson4) ? "Poisson(k=4))" :
898 (fFitFunc == kEPoisson5) ? "Poisson(k=5))" :
899 (fFitFunc == kEPoisson6) ? "Poisson(k=6))" :
900 (fFitFunc == kEPolya ) ? "Polya(k=4))" :
901 (fFitFunc == kEMichele ) ? "Michele)" : " none )" ));
902 fFitLegend->SetTextSize(0.05);
903 }
904 else
905 fFitLegend->Clear();
906
907 const TString line1 =
908 Form("Mean: #lambda = %2.2f #pm %2.2f",fLambda,fLambdaErr);
909 TText *t1 = fFitLegend->AddText(line1.Data());
910 t1->SetBit(kCanDelete);
911
912 const TString line6 =
913 Form("Mean #lambda (check) = %2.2f #pm %2.2f",fLambdaCheck,fLambdaCheckErr);
914 TText *t2 = fFitLegend->AddText(line6.Data());
915 t2->SetBit(kCanDelete);
916
917 const TString line2 =
918 Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",fMu0,fMu0Err);
919 TText *t3 = fFitLegend->AddText(line2.Data());
920 t3->SetBit(kCanDelete);
921
922 const TString line3 =
923 Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",fSigma0,fSigma0Err);
924 TText *t4 = fFitLegend->AddText(line3.Data());
925 t4->SetBit(kCanDelete);
926
927 const TString line4 =
928 Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",fMu1,fMu1Err);
929 TText *t5 = fFitLegend->AddText(line4.Data());
930 t5->SetBit(kCanDelete);
931
932 const TString line5 =
933 Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",fSigma1,fSigma1Err);
934 TText *t6 = fFitLegend->AddText(line5.Data());
935 t6->SetBit(kCanDelete);
936
937 const TString line7 =
938 Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChisquare,fNDF);
939 TText *t7 = fFitLegend->AddText(line7.Data());
940 t7->SetBit(kCanDelete);
941
942 const TString line8 =
943 Form("Probability: %4.2f ",fProb);
944 TText *t8 = fFitLegend->AddText(line8.Data());
945 t8->SetBit(kCanDelete);
946
947 if (IsSinglePheFitOK())
948 {
949 TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
950 t->SetBit(kCanDelete);
951 }
952 else
953 {
954 TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
955 t->SetBit(kCanDelete);
956 }
957
958 fFitLegend->SetFillColor(IsSinglePheFitOK() ? 80 : 2);
959 fFitLegend->Draw();
960
961 return;
962}
963
964
965// -------------------------------------------------------------------------
966//
967// Draw the histogram
968//
969// The following options can be chosen:
970//
971// "": displays the fHGausHist, the legend and fASinglePheFADCSlices and fAPedestalFADCSlices
972// "all": executes additionally MHGausEvents::Draw(), with option "fourierevents"
973//
974void MHCalibrationChargeBlindPix::Draw(Option_t *opt)
975{
976
977 TString option(opt);
978 option.ToLower();
979
980 Int_t win = 1;
981
982 TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600);
983 TVirtualPad *pad = NULL;
984
985 oldpad->SetBorderMode(0);
986
987 if (option.Contains("all"))
988 {
989 option.ReplaceAll("all","");
990 oldpad->Divide(2,1);
991 win = 2;
992 oldpad->cd(1);
993 TVirtualPad *newpad = gPad;
994 pad = newpad;
995 pad->Divide(2,2);
996 pad->cd(1);
997 }
998 else
999 {
1000 pad = oldpad;
1001 pad->Divide(2,2);
1002 pad->cd(1);
1003 }
1004
1005 if (!IsEmpty())
1006 gPad->SetLogy();
1007
1008 gPad->SetTicks();
1009
1010 fHGausHist.Draw(opt);
1011 if (fFGausFit)
1012 {
1013 fFGausFit->SetLineColor(kBlue);
1014 fFGausFit->Draw("same");
1015 TLine *line = new TLine(fSinglePheCut, 0., fSinglePheCut, 10.);
1016 line->SetBit(kCanDelete);
1017 line->SetLineColor(kBlue);
1018 line->SetLineWidth(3);
1019 line->DrawLine(fSinglePheCut, 0., fSinglePheCut, 2.);
1020 }
1021 if (fSinglePheFit)
1022 {
1023 fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed);
1024 fSinglePheFit->Draw("same");
1025 }
1026
1027
1028 pad->cd(2);
1029 DrawLegend();
1030
1031 pad->cd(3);
1032
1033 if (fASinglePheFADCSlices.GetNrows()!=1)
1034 {
1035 if (fHSinglePheFADCSlices)
1036 delete fHSinglePheFADCSlices;
1037 fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices);
1038 fHSinglePheFADCSlices->SetName("SinglePheFADCSlices");
1039 fHSinglePheFADCSlices->SetTitle(Form("%s%f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut));
1040 fHSinglePheFADCSlices->SetXTitle("FADC slice number");
1041 fHSinglePheFADCSlices->SetYTitle("FADC counts");
1042 fHSinglePheFADCSlices->Draw(opt);
1043 }
1044
1045 pad->cd(4);
1046 if (fAPedestalFADCSlices.GetNrows()!=1)
1047 {
1048
1049 if (fHPedestalFADCSlices)
1050 delete fHPedestalFADCSlices;
1051
1052 fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices);
1053 fHPedestalFADCSlices->SetName("PedestalFADCSlices");
1054 fHPedestalFADCSlices->SetTitle(Form("%s%f","Pedestal FADC Slices, Sum < ",fSinglePheCut));
1055 fHPedestalFADCSlices->SetXTitle("FADC slice number");
1056 fHPedestalFADCSlices->SetYTitle("FADC counts");
1057 fHPedestalFADCSlices->Draw(opt);
1058 }
1059
1060 if (win < 2)
1061 return;
1062
1063 oldpad->cd(2);
1064 MHGausEvents::Draw("fourierevents");
1065}
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
Note: See TracBrowser for help on using the repository browser.