source: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc@ 3005

Last change on this file since 3005 was 2997, checked in by gaug, 22 years ago
*** empty log message ***
File size: 21.8 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 11/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHCalibrationBlindPixel
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 "MHCalibrationBlindPixel.h"
63#include "MHCalibrationConfig.h"
64
65#include <TStyle.h>
66#include <TCanvas.h>
67#include <TPaveText.h>
68
69#include <TF1.h>
70#include <TH1.h>
71#include <TRandom.h>
72
73#include "MLog.h"
74#include "MLogManip.h"
75
76ClassImp(MHCalibrationBlindPixel);
77
78using namespace std;
79
80const Int_t MHCalibrationBlindPixel::fgBlindPixelChargeNbins = 1000;
81const Int_t MHCalibrationBlindPixel::fgBlindPixelTimeNbins = 32;
82const Int_t MHCalibrationBlindPixel::fgBlindPixelChargevsNbins = 10000;
83const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeFirst = -0.25;
84const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeLast = 15.75;
85const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmp = 0.008;
86const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002;
87
88// --------------------------------------------------------------------------
89//
90// Default Constructor.
91//
92MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title)
93 : fHBlindPixelPSD(NULL),
94 fSinglePheFit(NULL),
95 fTimeGausFit(NULL),
96 fSinglePhePedFit(NULL),
97 fFitLegend(NULL)
98{
99
100 fName = name ? name : "MHCalibrationBlindPixel";
101 fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits";
102
103 // Create a large number of bins, later we will rebin
104 fBlindPixelChargefirst = -200.;
105 fBlindPixelChargelast = 800.;
106
107 fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices",
108 fgBlindPixelChargeNbins,fBlindPixelChargefirst,fBlindPixelChargelast);
109 fHBlindPixelCharge->SetXTitle("Sum FADC Slices");
110 fHBlindPixelCharge->SetYTitle("Nr. of events");
111 fHBlindPixelCharge->Sumw2();
112 fHBlindPixelCharge->SetDirectory(NULL);
113
114 fHBlindPixelTime = new TH1F("HBlindPixelTime","Distribution of Mean Arrival Times",
115 fgBlindPixelTimeNbins,fgBlindPixelTimeFirst,fgBlindPixelTimeLast);
116 fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
117 fHBlindPixelTime->SetYTitle("Nr. of events");
118 fHBlindPixelTime->Sumw2();
119 fHBlindPixelTime->SetDirectory(NULL);
120
121 fHBlindPixelChargevsN = new TH1I("HBlindPixelChargevsN","Sum of Charges vs. Event Number",
122 fgBlindPixelChargevsNbins,-0.5,(Axis_t)fgBlindPixelChargevsNbins-0.5);
123 fHBlindPixelChargevsN->SetXTitle("Event Nr.");
124 fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices");
125 fHBlindPixelChargevsN->SetDirectory(NULL);
126
127 Clear();
128}
129
130MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
131{
132
133 if (fFitLegend)
134 delete fFitLegend;
135
136 delete fHBlindPixelCharge;
137 delete fHBlindPixelTime;
138 delete fHBlindPixelChargevsN;
139
140 if (fHBlindPixelPSD)
141 delete fHBlindPixelPSD;
142 if (fSinglePheFit)
143 delete fSinglePheFit;
144 if (fTimeGausFit)
145 delete fTimeGausFit;
146 if(fSinglePhePedFit)
147 delete fSinglePhePedFit;
148
149}
150
151void MHCalibrationBlindPixel::Clear(Option_t *o)
152{
153
154 fBlindPixelChargefirst = -200.;
155 fBlindPixelChargelast = 800.;
156
157 fLambda = 0.;
158 fMu0 = 0.;
159 fMu1 = 0.;
160 fSigma0 = 0.;
161 fSigma1 = 0.;
162
163 fLambdaErr = 0.;
164 fMu0Err = 0.;
165 fMu1Err = 0.;
166 fSigma0Err = 0.;
167 fSigma1Err = 0.;
168
169 fChisquare = -1.;
170 fProb = -1.;
171 fNdf = -1;
172
173 fMeanTime = -1.;
174 fMeanTimeErr = -1.;
175 fSigmaTime = -1.;
176 fSigmaTimeErr = -1.;
177
178 fLambdaCheck = -1.;
179 fLambdaCheckErr = -1.;
180
181 fMeanPedestal = 0.;
182 fMeanPedestalErr = 0.;
183 fSigmaPedestal = 0.;
184 fSigmaPedestalErr = 0.;
185
186 fFitFunc = kEPoisson4;
187
188 if (fFitLegend)
189 delete fFitLegend;
190 if (fHBlindPixelPSD)
191 delete fHBlindPixelPSD;
192 if (fSinglePheFit)
193 delete fSinglePheFit;
194 if (fTimeGausFit)
195 delete fTimeGausFit;
196 if(fSinglePhePedFit)
197 delete fSinglePhePedFit;
198
199 return;
200}
201
202void MHCalibrationBlindPixel::Reset()
203{
204
205 Clear();
206
207 fHBlindPixelCharge->Reset();
208 fHBlindPixelTime->Reset();
209 fHBlindPixelChargevsN->Reset();
210
211}
212
213Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(Float_t q)
214{
215 return fHBlindPixelCharge->Fill(q) > -1;
216}
217
218Bool_t MHCalibrationBlindPixel::FillBlindPixelTime(Float_t t)
219{
220 return fHBlindPixelTime->Fill(t) > -1;
221}
222
223Bool_t MHCalibrationBlindPixel::FillBlindPixelChargevsN(Stat_t rq, Int_t t)
224{
225 return fHBlindPixelChargevsN->Fill(t,rq) > -1;
226}
227
228
229// -------------------------------------------------------------------------
230//
231// Draw a legend with the fit results
232//
233void MHCalibrationBlindPixel::DrawLegend()
234{
235
236 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
237
238 if (fFitOK)
239 fFitLegend->SetFillColor(80);
240 else
241 fFitLegend->SetFillColor(2);
242
243 fFitLegend->SetLabel("Results of the single PhE Fit (to k=6):");
244 fFitLegend->SetTextSize(0.05);
245
246 const TString line1 =
247 Form("Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr());
248 TText *t1 = fFitLegend->AddText(line1.Data());
249 t1->SetBit(kCanDelete);
250
251 const TString line6 =
252 Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr());
253 TText *t2 = fFitLegend->AddText(line6.Data());
254 t2->SetBit(kCanDelete);
255
256 const TString line2 =
257 Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err());
258 TText *t3 = fFitLegend->AddText(line2.Data());
259 t3->SetBit(kCanDelete);
260
261 const TString line3 =
262 Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err());
263 TText *t4 = fFitLegend->AddText(line3.Data());
264 t4->SetBit(kCanDelete);
265
266 const TString line4 =
267 Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err());
268 TText *t5 = fFitLegend->AddText(line4.Data());
269 t5->SetBit(kCanDelete);
270
271 const TString line5 =
272 Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err());
273 TText *t6 = fFitLegend->AddText(line5.Data());
274 t6->SetBit(kCanDelete);
275
276 const TString line7 =
277 Form("#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf());
278 TText *t7 = fFitLegend->AddText(line7.Data());
279 t7->SetBit(kCanDelete);
280
281 const TString line8 =
282 Form("Probability: %4.2f ",GetProb());
283 TText *t8 = fFitLegend->AddText(line8.Data());
284 t8->SetBit(kCanDelete);
285
286 if (fFitOK)
287 {
288 TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
289 t->SetBit(kCanDelete);
290 }
291 else
292 {
293 TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
294 t->SetBit(kCanDelete);
295 }
296
297 fFitLegend->SetBit(kCanDelete);
298 fFitLegend->Draw();
299
300 return;
301}
302
303TObject *MHCalibrationBlindPixel::DrawClone(Option_t *option) const
304{
305
306 gROOT->SetSelectedPad(NULL);
307
308 MHCalibrationBlindPixel *newobj = (MHCalibrationBlindPixel*)Clone();
309
310 if (!newobj)
311 return 0;
312 newobj->SetBit(kCanDelete);
313
314 if (strlen(option))
315 newobj->Draw(option);
316 else
317 newobj->Draw(GetDrawOption());
318
319 return newobj;
320}
321
322
323// -------------------------------------------------------------------------
324//
325// Draw the histogram
326//
327void MHCalibrationBlindPixel::Draw(Option_t *opt)
328{
329
330 gStyle->SetOptFit(1);
331 gStyle->SetOptStat(111111);
332
333 TCanvas *c = MakeDefCanvas(this,550,700);
334
335 c->Divide(2,3);
336
337 gROOT->SetSelectedPad(NULL);
338
339 c->cd(1);
340 gPad->SetLogx(0);
341 gPad->SetLogy(1);
342 gPad->SetTicks();
343
344 fHBlindPixelCharge->Draw(opt);
345
346 if (fFitOK)
347 fSinglePheFit->SetLineColor(kGreen);
348 else
349 fSinglePheFit->SetLineColor(kRed);
350
351 fSinglePheFit->Draw("same");
352 c->Modified();
353 c->Update();
354
355 fSinglePhePedFit->SetLineColor(kBlue);
356 fSinglePhePedFit->Draw("same");
357
358 c->cd(2);
359 DrawLegend();
360 c->Modified();
361 c->Update();
362
363 c->cd(3);
364 gPad->SetLogy(1);
365 gPad->SetBorderMode(0);
366 fHBlindPixelTime->Draw(opt);
367
368 c->cd(4);
369
370 fHBlindPixelChargevsN->Draw(opt);
371
372 c->Modified();
373 c->Update();
374
375 c->cd(5);
376
377 // fHBlindPixelPSD = (TH1F*)MFFT::PowerSpectrumDensity(fHBlindPixelChargevsN);
378 // TH1F *hist = MFFT::PowerSpectrumDensity((*this)->fHBlindPixelChargevsN);
379
380 // fHBlindPixelPSD->Draw(opt);
381 c->Modified();
382 c->Update();
383
384}
385
386
387
388Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
389{
390 gRandom->SetSeed();
391
392 if (fHBlindPixelCharge->GetIntegral() != 0)
393 {
394 *fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl;
395 *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
396 return kFALSE;
397 }
398
399 if (!InitFit(fBlindPixelChargefirst,fBlindPixelChargelast))
400 return kFALSE;
401
402 for (Int_t i=0;i<10000; i++)
403 fHBlindPixelCharge->Fill(fSinglePheFit->GetRandom());
404
405 return kTRUE;
406}
407
408Bool_t MHCalibrationBlindPixel::InitFit(Axis_t min, Axis_t max)
409{
410
411 //
412 // First guesses for the fit (should be as close to reality as possible,
413 // otherwise the fit goes gaga because of high number of dimensions ...
414 //
415 const Stat_t entries = fHBlindPixelCharge->Integral("width");
416 const Double_t lambda_guess = 0.5;
417 const Double_t maximum_bin = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
418 const Double_t norm = entries/gkSq2Pi;
419
420 //
421 // Initialize the fit function
422 //
423 switch (fFitFunc)
424 {
425 case kEPoisson4:
426 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
427 break;
428 case kEPoisson5:
429 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
430 break;
431 case kEPoisson6:
432 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
433 break;
434 case kEPolya:
435 fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8);
436 break;
437 case kEMichele:
438 break;
439
440 default:
441 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
442 return kFALSE;
443 break;
444 }
445
446 const Double_t mu_0_guess = maximum_bin;
447 const Double_t si_0_guess = 40.;
448 const Double_t mu_1_guess = mu_0_guess + 100.;
449 const Double_t si_1_guess = si_0_guess + si_0_guess;
450 // Michele
451 const Double_t lambda_1cat_guess = 0.5;
452 const Double_t lambda_1dyn_guess = 0.5;
453 const Double_t mu_1cat_guess = mu_0_guess + 50.;
454 const Double_t mu_1dyn_guess = mu_0_guess + 20.;
455 const Double_t si_1cat_guess = si_0_guess + si_0_guess;
456 const Double_t si_1dyn_guess = si_0_guess;
457 // Polya
458 const Double_t excessPoisson_guess = 0.5;
459 const Double_t delta1_guess = 8.;
460 const Double_t delta2_guess = 5.;
461 const Double_t electronicAmp_guess = fgBlindPixelElectronicAmp;
462 const Double_t electronicAmp_limit = fgBlindPixelElectronicAmpError;
463
464 //
465 // Initialize boundaries and start parameters
466 //
467 switch (fFitFunc)
468 {
469
470 case kEPoisson4:
471 if ((fMeanPedestal) && (fSigmaPedestal))
472 fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm);
473 else
474 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
475
476 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
477
478 fSinglePheFit->SetParLimits(0,0.,1.);
479
480 if ((fMeanPedestal) && (fSigmaPedestal))
481 fSinglePheFit->SetParLimits(1,
482 fMeanPedestal-1.*fMeanPedestalErr,
483 fMeanPedestal+1.*fMeanPedestalErr);
484 else
485 fSinglePheFit->SetParLimits(1,-3.,0.);
486
487 fSinglePheFit->SetParLimits(2,(max-min)/2.,max);
488
489 if ((fMeanPedestal) && (fSigmaPedestal))
490 fSinglePheFit->SetParLimits(3,
491 fSigmaPedestal-3.*fSigmaPedestalErr,
492 fSigmaPedestal+3.*fSigmaPedestalErr);
493 else
494 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
495
496 fSinglePheFit->SetParLimits(4,1.0,(max-min));
497 fSinglePheFit->SetParLimits(5,norm-0.5,norm+0.5);
498 break;
499 case kEPoisson5:
500 case kEPoisson6:
501 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
502 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
503 fSinglePheFit->SetParLimits(0,0.,1.);
504 fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
505 fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
506 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
507 fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
508 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
509 break;
510
511 case kEPolya:
512 if ((fMeanPedestal) && (fSigmaPedestal))
513 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
514 delta1_guess,delta2_guess,
515 electronicAmp_guess,
516 fSigmaPedestal,
517 norm,
518 fMeanPedestal);
519 else
520 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
521 delta1_guess,delta2_guess,
522 electronicAmp_guess,
523 si_0_guess,
524 norm, mu_0_guess);
525 fSinglePheFit->SetParNames("#lambda","b_{tot}",
526 "#delta_{1}","#delta_{2}",
527 "amp_{e}","#sigma_{0}",
528 "Area", "#mu_{0}");
529 fSinglePheFit->SetParLimits(0,0.,1.);
530 fSinglePheFit->SetParLimits(1,0.,1.);
531 fSinglePheFit->SetParLimits(2,6.,12.);
532 fSinglePheFit->SetParLimits(3,3.,8.);
533 fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
534 electronicAmp_guess+electronicAmp_limit);
535 if ((fMeanPedestal) && (fSigmaPedestal))
536 fSinglePheFit->SetParLimits(5,
537 fSigmaPedestal-3.*fSigmaPedestalErr,
538 fSigmaPedestal+3.*fSigmaPedestalErr);
539 else
540 fSinglePheFit->SetParLimits(5,min,(max-min)/1.5);
541
542 fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
543 if ((fMeanPedestal) && (fSigmaPedestal))
544 fSinglePheFit->SetParLimits(7,
545 fMeanPedestal-3.*fMeanPedestalErr,
546 fMeanPedestal+3.*fMeanPedestalErr);
547 else
548 fSinglePheFit->SetParLimits(7,-35.,15.);
549 break;
550
551 case kEMichele:
552
553
554 break;
555
556 default:
557 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
558 return kFALSE;
559 break;
560 }
561
562 return kTRUE;
563}
564
565void MHCalibrationBlindPixel::ExitFit(TF1 *f)
566{
567
568
569 //
570 // Finalize
571 //
572 switch (fFitFunc)
573 {
574
575 case kEPoisson4:
576 case kEPoisson5:
577 case kEPoisson6:
578 case kEPoisson7:
579 fLambda = fSinglePheFit->GetParameter(0);
580 fMu0 = fSinglePheFit->GetParameter(1);
581 fMu1 = fSinglePheFit->GetParameter(2);
582 fSigma0 = fSinglePheFit->GetParameter(3);
583 fSigma1 = fSinglePheFit->GetParameter(4);
584
585 fLambdaErr = fSinglePheFit->GetParError(0);
586 fMu0Err = fSinglePheFit->GetParError(1);
587 fMu1Err = fSinglePheFit->GetParError(2);
588 fSigma0Err = fSinglePheFit->GetParError(3);
589 fSigma1Err = fSinglePheFit->GetParError(4);
590 break;
591 case kEPolya:
592 fLambda = fSinglePheFit->GetParameter(0);
593 fMu0 = fSinglePheFit->GetParameter(7);
594 fMu1 = 0.;
595 fSigma0 = fSinglePheFit->GetParameter(5);
596 fSigma1 = 0.;
597
598 fLambdaErr = fSinglePheFit->GetParError(0);
599 fMu0Err = fSinglePheFit->GetParError(7);
600 fMu1Err = 0.;
601 fSigma0Err = fSinglePheFit->GetParError(5);
602 fSigma1Err = 0.;
603 default:
604 break;
605 }
606
607}
608
609
610Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt)
611{
612
613 //
614 // Get the fitting ranges
615 //
616 rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
617 rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;
618
619 if (!InitFit(rmin,rmax))
620 return kFALSE;
621
622 fHBlindPixelCharge->Fit(fSinglePheFit,opt);
623
624
625 ExitFit(fSinglePheFit);
626
627 fProb = fSinglePheFit->GetProb();
628 fChisquare = fSinglePheFit->GetChisquare();
629 fNdf = fSinglePheFit->GetNDF();
630
631 // Perform the cross-check fitting only the pedestal:
632 fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.);
633 fHBlindPixelCharge->Fit(fSinglePhePedFit,opt);
634
635 const Stat_t entries = fHBlindPixelCharge->Integral("width");
636
637 Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2);
638 fLambdaCheck = TMath::Log(entries/pedarea);
639 fLambdaCheckErr = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
640 + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);
641
642 *fLog << inf << "Results of the Blind Pixel Fit: " << endl;
643 *fLog << inf << "Chisquare: " << fChisquare << endl;
644 *fLog << inf << "DoF: " << fNdf << endl;
645 *fLog << inf << "Probability: " << fProb << endl;
646
647 //
648 // The fit result is accepted under condition that:
649 // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
650 // 2) at least 50 events are in the single Photo-electron peak
651 //
652 if (fProb < gkProbLimit)
653 {
654 *fLog << err << "ERROR: Fit Probability " << fProb
655 << " is smaller than the allowed value: " << gkProbLimit << endl;
656 fFitOK = kFALSE;
657 return kFALSE;
658 }
659
660 Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
661
662 if (contSinglePhe < 50.)
663 {
664 *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe
665 << " in the Single Photo-Electron peak " << endl;
666 fFitOK = kFALSE;
667 return kFALSE;
668 }
669 else
670 *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
671
672 fFitOK = kTRUE;
673
674 return kTRUE;
675}
676
677
678void MHCalibrationBlindPixel::CutAllEdges()
679{
680
681 Int_t nbins = 30;
682
683 CutEdges(fHBlindPixelCharge,nbins);
684
685 fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
686 fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
687
688 CutEdges(fHBlindPixelChargevsN,0);
689
690}
691
692Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt)
693{
694
695 rmin = (rmin != 0.) ? rmin : 4.;
696 rmax = (rmax != 0.) ? rmax : 9.;
697
698 const Stat_t entries = fHBlindPixelTime->Integral();
699 const Double_t mu_guess = fHBlindPixelTime->GetBinCenter(fHBlindPixelTime->GetMaximumBin());
700 const Double_t sigma_guess = (rmax - rmin)/2.;
701 const Double_t area_guess = entries/gkSq2Pi;
702
703 fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax);
704 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
705 fTimeGausFit->SetParNames("Area","#mu","#sigma");
706 fTimeGausFit->SetParLimits(0,0.,entries);
707 fTimeGausFit->SetParLimits(1,rmin,rmax);
708 fTimeGausFit->SetParLimits(2,0.,rmax-rmin);
709
710 fHBlindPixelTime->Fit(fTimeGausFit,opt);
711 rmin = fTimeGausFit->GetParameter(1) - 2.*fTimeGausFit->GetParameter(2);
712 rmax = fTimeGausFit->GetParameter(1) + 2.*fTimeGausFit->GetParameter(2);
713 fTimeGausFit->SetRange(rmin,rmax);
714
715 fHBlindPixelTime->Fit(fTimeGausFit,opt);
716
717 fMeanTime = fTimeGausFit->GetParameter(2);
718 fSigmaTime = fTimeGausFit->GetParameter(3);
719 fMeanTimeErr = fTimeGausFit->GetParError(2);
720 fSigmaTimeErr = fTimeGausFit->GetParError(3);
721
722 *fLog << inf << "Results of the Times Fit: " << endl;
723 *fLog << inf << "Chisquare: " << fTimeGausFit->GetChisquare() << endl;
724 *fLog << inf << "Ndf: " << fTimeGausFit->GetNDF() << endl;
725
726 return kTRUE;
727
728}
Note: See TracBrowser for help on using the repository browser.