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

Last change on this file since 3107 was 3081, checked in by gaug, 22 years ago
*** empty log message ***
File size: 26.0 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
66#include <TStyle.h>
67#include <TCanvas.h>
68#include <TPaveText.h>
69
70#include <TGraph.h>
71#include <TF1.h>
72#include <TH1.h>
73#include <TRandom.h>
74
75#include "MFFT.h"
76#include "MLog.h"
77#include "MLogManip.h"
78
79ClassImp(MHCalibrationBlindPixel);
80
81using namespace std;
82
83const Int_t MHCalibrationBlindPixel::fgBlindPixelChargeNbins = 1000;
84const Int_t MHCalibrationBlindPixel::fgBlindPixelTimeNbins = 22;
85const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeFirst = -9.00;
86const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeLast = 12.00;
87const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmp = 0.008;
88const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002;
89
90const Int_t MHCalibrationBlindPixel::fPSDNbins = 30;
91
92// --------------------------------------------------------------------------
93//
94// Default Constructor.
95//
96MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title)
97 : fHBlindPixelPSD(NULL),
98 fSinglePheFit(NULL),
99 fTimeGausFit(NULL),
100 fSinglePhePedFit(NULL),
101 fPSDHiGain(NULL),
102 fPSDLoGain(NULL),
103 fHPSD(NULL),
104 fPSDExpFit(NULL),
105 fChargeXaxis(NULL),
106 fPSDXaxis(NULL),
107 fCurrentSize(1024),
108 fFitLegend(NULL)
109{
110
111 fName = name ? name : "MHCalibrationBlindPixel";
112 fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits";
113
114 // Create a large number of bins, later we will rebin
115 fBlindPixelChargefirst = -200.;
116 fBlindPixelChargelast = 800.;
117
118 fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices",
119 fgBlindPixelChargeNbins,fBlindPixelChargefirst,fBlindPixelChargelast);
120 fHBlindPixelCharge->SetXTitle("Sum FADC Slices");
121 fHBlindPixelCharge->SetYTitle("Nr. of events");
122 fHBlindPixelCharge->Sumw2();
123 fHBlindPixelCharge->SetDirectory(NULL);
124
125 fHBlindPixelTime = new TH1F("HBlindPixelTime","Distribution of Mean Arrival Times",
126 fgBlindPixelTimeNbins,fgBlindPixelTimeFirst,fgBlindPixelTimeLast);
127 fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
128 fHBlindPixelTime->SetYTitle("Nr. of events");
129 fHBlindPixelTime->Sumw2();
130 fHBlindPixelTime->SetDirectory(NULL);
131
132 fHiGains = new TArrayF(fCurrentSize);
133 fLoGains = new TArrayF(fCurrentSize);
134
135 Clear();
136}
137
138MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
139{
140
141 if (fFitLegend)
142 delete fFitLegend;
143
144 delete fHBlindPixelCharge;
145 delete fHBlindPixelTime;
146
147 delete fHiGains;
148 delete fLoGains;
149
150 if (fHBlindPixelPSD)
151 delete fHBlindPixelPSD;
152 if (fSinglePheFit)
153 delete fSinglePheFit;
154 if (fTimeGausFit)
155 delete fTimeGausFit;
156 if(fSinglePhePedFit)
157 delete fSinglePhePedFit;
158 if (fPSDExpFit)
159 delete fPSDExpFit;
160 if (fHPSD)
161 delete fHPSD;
162 if (fChargeXaxis)
163 delete fChargeXaxis;
164 if (fPSDXaxis)
165 delete fPSDXaxis;
166
167}
168
169void MHCalibrationBlindPixel::Clear(Option_t *o)
170{
171
172 fTotalEntries = 0;
173 fCurrentSize = 1024;
174
175 fBlindPixelChargefirst = -200.;
176 fBlindPixelChargelast = 800.;
177
178 fLambda = 0.;
179 fMu0 = 0.;
180 fMu1 = 0.;
181 fSigma0 = 0.;
182 fSigma1 = 0.;
183
184 fLambdaErr = 0.;
185 fMu0Err = 0.;
186 fMu1Err = 0.;
187 fSigma0Err = 0.;
188 fSigma1Err = 0.;
189
190 fChisquare = -1.;
191 fProb = -1.;
192 fNdf = -1;
193
194 fMeanTime = -1.;
195 fMeanTimeErr = -1.;
196 fSigmaTime = -1.;
197 fSigmaTimeErr = -1.;
198
199 fLambdaCheck = -1.;
200 fLambdaCheckErr = -1.;
201
202 fMeanPedestal = 0.;
203 fMeanPedestalErr = 0.;
204 fSigmaPedestal = 0.;
205 fSigmaPedestalErr = 0.;
206
207 fFitFunc = kEPoisson4;
208
209 if (fFitLegend)
210 delete fFitLegend;
211 if (fHBlindPixelPSD)
212 delete fHBlindPixelPSD;
213 if (fSinglePheFit)
214 delete fSinglePheFit;
215 if (fTimeGausFit)
216 delete fTimeGausFit;
217 if(fSinglePhePedFit)
218 delete fSinglePhePedFit;
219 if (fPSDExpFit)
220 delete fPSDExpFit;
221 if (fHPSD)
222 delete fHPSD;
223 if (fChargeXaxis)
224 delete fChargeXaxis;
225 if (fPSDXaxis)
226 delete fPSDXaxis;
227 if (fPSDHiGain)
228 delete fPSDHiGain;
229 if (fPSDLoGain)
230 delete fPSDLoGain;
231
232 CLRBIT(fFlags,kFitOK);
233 CLRBIT(fFlags,kOscillating);
234
235 return;
236}
237
238void MHCalibrationBlindPixel::Reset()
239{
240
241 Clear();
242
243 fHBlindPixelCharge->Reset();
244 fHBlindPixelTime->Reset();
245
246 fHiGains->Set(1024);
247 fLoGains->Set(1024);
248
249 fHiGains->Reset(0.);
250 fLoGains->Reset(0.);
251
252
253}
254
255Bool_t MHCalibrationBlindPixel::CheckOscillations()
256{
257
258 if (fPSDExpFit)
259 return IsOscillating();
260
261 MFFT fourier;
262
263 fPSDLoGain = fourier.PowerSpectrumDensity(fLoGains);
264 fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains);
265
266 fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins,
267 "PSDProjectionBlindPixel",
268 "Power Spectrum Density Projection HiGain Blind Pixel");
269 //
270 // First guesses for the fit (should be as close to reality as possible,
271 //
272 const Double_t xmax = fHPSD->GetXaxis()->GetXmax();
273
274 fPSDExpFit = new TF1("PSDExpFit","exp([0]-[1]*x)",0.,xmax);
275
276 const Double_t slope_guess = (TMath::Log(fHPSD->GetEntries())+1.)/xmax;
277 const Double_t offset_guess = slope_guess*xmax;
278
279 fPSDExpFit->SetParameters(offset_guess, slope_guess);
280 fPSDExpFit->SetParNames("Offset","Slope");
281 fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess);
282 fPSDExpFit->SetParLimits(1,slope_guess/2.,2.*slope_guess);
283 fPSDExpFit->SetRange(0.,xmax);
284
285 fHPSD->Fit(fPSDExpFit,"RQL0");
286
287 fPSDProb = fPSDExpFit->GetProb();
288
289 if (fPSDProb < gkProbLimit)
290 {
291 SETBIT(fFlags,kOscillating);
292 return kTRUE;
293 }
294
295 CLRBIT(fFlags,kOscillating);
296
297 return kFALSE;
298}
299
300void MHCalibrationBlindPixel::CreatePSDXaxis(Int_t n)
301{
302
303 if (fPSDXaxis)
304 return;
305
306 fPSDXaxis = new TArrayF(n);
307
308 for (Int_t i=0;i<n;i++)
309 fPSDXaxis->AddAt((Float_t)i,i);
310}
311
312void MHCalibrationBlindPixel::CreateChargeXaxis(Int_t n)
313{
314
315 if (!fChargeXaxis)
316 {
317 fChargeXaxis = new TArrayF(n);
318 for (Int_t i=0;i<n;i++)
319 fChargeXaxis->AddAt((Float_t)i,i);
320 return;
321 }
322
323 if (fChargeXaxis->GetSize() == n)
324 return;
325
326 const Int_t diff = fChargeXaxis->GetSize()-n;
327 fChargeXaxis->Set(n);
328 if (diff < 0)
329 for (Int_t i=n;i<n+diff;i++)
330 fChargeXaxis->AddAt((Float_t)i,i);
331}
332
333void MHCalibrationBlindPixel::CutArrayBorder(TArrayF *array)
334{
335
336 Int_t i;
337
338 for (i=array->GetSize()-1;i>=0;i--)
339 if (array->At(i) != 0)
340 {
341 array->Set(i+1);
342 break;
343 }
344}
345
346const Bool_t MHCalibrationBlindPixel::IsFitOK() const
347{
348 return TESTBIT(fFlags,kFitOK);
349}
350
351const Bool_t MHCalibrationBlindPixel::IsOscillating()
352{
353
354 if (fPSDExpFit)
355 return TESTBIT(fFlags,kOscillating);
356
357 return CheckOscillations();
358
359}
360
361Bool_t MHCalibrationBlindPixel::FillGraphs(Float_t qhi,Float_t qlo)
362{
363
364 if (fTotalEntries >= fCurrentSize)
365 {
366 fCurrentSize *= 2;
367
368 fHiGains->Set(fCurrentSize);
369 fLoGains->Set(fCurrentSize);
370 }
371
372 fHiGains->AddAt(qhi,fTotalEntries);
373 fLoGains->AddAt(qlo,fTotalEntries);
374
375 fTotalEntries++;
376
377 return kTRUE;
378
379}
380
381
382
383Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(Float_t q)
384{
385 return fHBlindPixelCharge->Fill(q) > -1;
386}
387
388Bool_t MHCalibrationBlindPixel::FillBlindPixelTime(Float_t t)
389{
390 return fHBlindPixelTime->Fill(t) > -1;
391}
392
393
394// -------------------------------------------------------------------------
395//
396// Draw a legend with the fit results
397//
398void MHCalibrationBlindPixel::DrawLegend()
399{
400
401 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
402
403 if (IsFitOK())
404 fFitLegend->SetFillColor(80);
405 else
406 fFitLegend->SetFillColor(2);
407
408 fFitLegend->SetLabel("Results of the single PhE Fit (to k=6):");
409 fFitLegend->SetTextSize(0.05);
410
411 const TString line1 =
412 Form("Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr());
413 TText *t1 = fFitLegend->AddText(line1.Data());
414 t1->SetBit(kCanDelete);
415
416 const TString line6 =
417 Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr());
418 TText *t2 = fFitLegend->AddText(line6.Data());
419 t2->SetBit(kCanDelete);
420
421 const TString line2 =
422 Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err());
423 TText *t3 = fFitLegend->AddText(line2.Data());
424 t3->SetBit(kCanDelete);
425
426 const TString line3 =
427 Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err());
428 TText *t4 = fFitLegend->AddText(line3.Data());
429 t4->SetBit(kCanDelete);
430
431 const TString line4 =
432 Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err());
433 TText *t5 = fFitLegend->AddText(line4.Data());
434 t5->SetBit(kCanDelete);
435
436 const TString line5 =
437 Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err());
438 TText *t6 = fFitLegend->AddText(line5.Data());
439 t6->SetBit(kCanDelete);
440
441 const TString line7 =
442 Form("#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf());
443 TText *t7 = fFitLegend->AddText(line7.Data());
444 t7->SetBit(kCanDelete);
445
446 const TString line8 =
447 Form("Probability: %4.2f ",GetProb());
448 TText *t8 = fFitLegend->AddText(line8.Data());
449 t8->SetBit(kCanDelete);
450
451 if (IsFitOK())
452 {
453 TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
454 t->SetBit(kCanDelete);
455 }
456 else
457 {
458 TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
459 t->SetBit(kCanDelete);
460 }
461
462 fFitLegend->SetBit(kCanDelete);
463 fFitLegend->Draw();
464
465 return;
466}
467
468TObject *MHCalibrationBlindPixel::DrawClone(Option_t *option) const
469{
470
471 gROOT->SetSelectedPad(NULL);
472
473 MHCalibrationBlindPixel *newobj = (MHCalibrationBlindPixel*)Clone();
474
475 if (!newobj)
476 return 0;
477 newobj->SetBit(kCanDelete);
478
479 if (strlen(option))
480 newobj->Draw(option);
481 else
482 newobj->Draw(GetDrawOption());
483
484 return newobj;
485}
486
487
488// -------------------------------------------------------------------------
489//
490// Draw the histogram
491//
492void MHCalibrationBlindPixel::Draw(Option_t *opt)
493{
494
495 gStyle->SetOptFit(1);
496 gStyle->SetOptStat(111111);
497
498 TCanvas *c = MakeDefCanvas(this,550,700);
499
500 c->Divide(2,4);
501
502 gROOT->SetSelectedPad(NULL);
503
504 c->cd(1);
505 gPad->SetLogx(0);
506 gPad->SetLogy(1);
507 gPad->SetTicks();
508
509 fHBlindPixelCharge->Draw(opt);
510
511 if (IsFitOK())
512 fSinglePheFit->SetLineColor(kGreen);
513 else
514 fSinglePheFit->SetLineColor(kRed);
515
516 fSinglePheFit->Draw("same");
517 c->Modified();
518 c->Update();
519
520 fSinglePhePedFit->SetLineColor(kBlue);
521 fSinglePhePedFit->Draw("same");
522
523 c->cd(2);
524 DrawLegend();
525 c->Modified();
526 c->Update();
527
528 c->cd(3);
529 gPad->SetLogy(1);
530 gPad->SetBorderMode(0);
531 fHBlindPixelTime->Draw(opt);
532
533 CutArrayBorder(fHiGains);
534 CreateChargeXaxis(fHiGains->GetSize());
535
536 c->cd(5);
537 gPad->SetTicks();
538 TGraph *gr1 = new TGraph(fChargeXaxis->GetSize(),
539 fChargeXaxis->GetArray(),
540 fHiGains->GetArray());
541 gr1->SetTitle("Evolution of HiGain charges with event number");
542 gr1->SetBit(kCanDelete);
543 gr1->Draw("AL");
544
545 CutArrayBorder(fLoGains);
546 CreateChargeXaxis(fHiGains->GetSize());
547
548 c->cd(6);
549 gPad->SetTicks();
550 TGraph *gr2 = new TGraph(fChargeXaxis->GetSize(),
551 fChargeXaxis->GetArray(),
552 fLoGains->GetArray());
553 gr2->SetTitle("Evolution of HiGain charges with event number");
554 gr2->SetBit(kCanDelete);
555 gr2->Draw("AL");
556
557 c->Modified();
558 c->Update();
559
560 c->cd(7);
561
562 TArrayF *array;
563
564 if (!fPSDHiGain)
565 return;
566 array = fPSDHiGain;
567
568 if (!fPSDXaxis)
569 CreatePSDXaxis(array->GetSize());
570
571 TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
572 gr3->SetTitle("Power Spectrum Density");
573 gr3->SetBit(kCanDelete);
574 gr3->Draw("AL");
575
576 c->Modified();
577 c->Update();
578
579 c->cd(8);
580
581 gStyle->SetOptStat(111111);
582 gPad->SetTicks();
583
584 if (fHPSD->Integral() > 0)
585 gPad->SetLogy();
586
587 fHPSD->Draw(opt);
588
589 if (fPSDExpFit)
590 {
591 fPSDExpFit->SetLineColor(IsOscillating() ? kRed : kGreen);
592 fPSDExpFit->Draw("same");
593 }
594
595 c->Modified();
596 c->Update();
597
598}
599
600
601
602Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
603{
604 gRandom->SetSeed();
605
606 if (fHBlindPixelCharge->GetIntegral() != 0)
607 {
608 *fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl;
609 *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
610 return kFALSE;
611 }
612
613 if (!InitFit(fBlindPixelChargefirst,fBlindPixelChargelast))
614 return kFALSE;
615
616 for (Int_t i=0;i<10000; i++)
617 fHBlindPixelCharge->Fill(fSinglePheFit->GetRandom());
618
619 return kTRUE;
620}
621
622Bool_t MHCalibrationBlindPixel::InitFit(Axis_t min, Axis_t max)
623{
624
625 //
626 // First guesses for the fit (should be as close to reality as possible,
627 // otherwise the fit goes gaga because of high number of dimensions ...
628 //
629 const Stat_t entries = fHBlindPixelCharge->Integral("width");
630 const Double_t lambda_guess = 0.5;
631 const Double_t maximum_bin = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
632 const Double_t norm = entries/gkSq2Pi;
633
634 //
635 // Initialize the fit function
636 //
637 switch (fFitFunc)
638 {
639 case kEPoisson4:
640 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
641 break;
642 case kEPoisson5:
643 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
644 break;
645 case kEPoisson6:
646 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
647 break;
648 case kEPolya:
649 fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8);
650 break;
651 case kEMichele:
652 break;
653
654 default:
655 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
656 return kFALSE;
657 break;
658 }
659
660 const Double_t mu_0_guess = maximum_bin;
661 const Double_t si_0_guess = 40.;
662 const Double_t mu_1_guess = mu_0_guess + 100.;
663 const Double_t si_1_guess = si_0_guess + si_0_guess;
664 // Michele
665 const Double_t lambda_1cat_guess = 0.5;
666 const Double_t lambda_1dyn_guess = 0.5;
667 const Double_t mu_1cat_guess = mu_0_guess + 50.;
668 const Double_t mu_1dyn_guess = mu_0_guess + 20.;
669 const Double_t si_1cat_guess = si_0_guess + si_0_guess;
670 const Double_t si_1dyn_guess = si_0_guess;
671 // Polya
672 const Double_t excessPoisson_guess = 0.5;
673 const Double_t delta1_guess = 8.;
674 const Double_t delta2_guess = 5.;
675 const Double_t electronicAmp_guess = fgBlindPixelElectronicAmp;
676 const Double_t electronicAmp_limit = fgBlindPixelElectronicAmpError;
677
678 //
679 // Initialize boundaries and start parameters
680 //
681 switch (fFitFunc)
682 {
683
684 case kEPoisson4:
685 if ((fMeanPedestal) && (fSigmaPedestal))
686 fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm);
687 else
688 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
689
690 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
691
692 fSinglePheFit->SetParLimits(0,0.,1.);
693
694 if ((fMeanPedestal) && (fSigmaPedestal))
695 fSinglePheFit->SetParLimits(1,
696 fMeanPedestal-1.*fMeanPedestalErr,
697 fMeanPedestal+1.*fMeanPedestalErr);
698 else
699 fSinglePheFit->SetParLimits(1,-3.,0.);
700
701 fSinglePheFit->SetParLimits(2,(max-min)/2.,max);
702
703 if ((fMeanPedestal) && (fSigmaPedestal))
704 fSinglePheFit->SetParLimits(3,
705 fSigmaPedestal-3.*fSigmaPedestalErr,
706 fSigmaPedestal+3.*fSigmaPedestalErr);
707 else
708 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
709
710 fSinglePheFit->SetParLimits(4,1.0,(max-min));
711 fSinglePheFit->SetParLimits(5,norm-0.5,norm+0.5);
712 break;
713 case kEPoisson5:
714 case kEPoisson6:
715 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
716 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
717 fSinglePheFit->SetParLimits(0,0.,1.);
718 fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
719 fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
720 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
721 fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
722 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
723 break;
724
725 case kEPolya:
726 if ((fMeanPedestal) && (fSigmaPedestal))
727 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
728 delta1_guess,delta2_guess,
729 electronicAmp_guess,
730 fSigmaPedestal,
731 norm,
732 fMeanPedestal);
733 else
734 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
735 delta1_guess,delta2_guess,
736 electronicAmp_guess,
737 si_0_guess,
738 norm, mu_0_guess);
739 fSinglePheFit->SetParNames("#lambda","b_{tot}",
740 "#delta_{1}","#delta_{2}",
741 "amp_{e}","#sigma_{0}",
742 "Area", "#mu_{0}");
743 fSinglePheFit->SetParLimits(0,0.,1.);
744 fSinglePheFit->SetParLimits(1,0.,1.);
745 fSinglePheFit->SetParLimits(2,6.,12.);
746 fSinglePheFit->SetParLimits(3,3.,8.);
747 fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
748 electronicAmp_guess+electronicAmp_limit);
749 if ((fMeanPedestal) && (fSigmaPedestal))
750 fSinglePheFit->SetParLimits(5,
751 fSigmaPedestal-3.*fSigmaPedestalErr,
752 fSigmaPedestal+3.*fSigmaPedestalErr);
753 else
754 fSinglePheFit->SetParLimits(5,min,(max-min)/1.5);
755
756 fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
757 if ((fMeanPedestal) && (fSigmaPedestal))
758 fSinglePheFit->SetParLimits(7,
759 fMeanPedestal-3.*fMeanPedestalErr,
760 fMeanPedestal+3.*fMeanPedestalErr);
761 else
762 fSinglePheFit->SetParLimits(7,-35.,15.);
763 break;
764
765 case kEMichele:
766
767
768 break;
769
770 default:
771 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
772 return kFALSE;
773 break;
774 }
775
776 return kTRUE;
777}
778
779void MHCalibrationBlindPixel::ExitFit(TF1 *f)
780{
781
782
783 //
784 // Finalize
785 //
786 switch (fFitFunc)
787 {
788
789 case kEPoisson4:
790 case kEPoisson5:
791 case kEPoisson6:
792 case kEPoisson7:
793 fLambda = fSinglePheFit->GetParameter(0);
794 fMu0 = fSinglePheFit->GetParameter(1);
795 fMu1 = fSinglePheFit->GetParameter(2);
796 fSigma0 = fSinglePheFit->GetParameter(3);
797 fSigma1 = fSinglePheFit->GetParameter(4);
798
799 fLambdaErr = fSinglePheFit->GetParError(0);
800 fMu0Err = fSinglePheFit->GetParError(1);
801 fMu1Err = fSinglePheFit->GetParError(2);
802 fSigma0Err = fSinglePheFit->GetParError(3);
803 fSigma1Err = fSinglePheFit->GetParError(4);
804 break;
805 case kEPolya:
806 fLambda = fSinglePheFit->GetParameter(0);
807 fMu0 = fSinglePheFit->GetParameter(7);
808 fMu1 = 0.;
809 fSigma0 = fSinglePheFit->GetParameter(5);
810 fSigma1 = 0.;
811
812 fLambdaErr = fSinglePheFit->GetParError(0);
813 fMu0Err = fSinglePheFit->GetParError(7);
814 fMu1Err = 0.;
815 fSigma0Err = fSinglePheFit->GetParError(5);
816 fSigma1Err = 0.;
817 default:
818 break;
819 }
820
821}
822
823
824Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt)
825{
826
827 //
828 // Get the fitting ranges
829 //
830 rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
831 rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;
832
833 if (!InitFit(rmin,rmax))
834 return kFALSE;
835
836 fHBlindPixelCharge->Fit(fSinglePheFit,opt);
837
838 ExitFit(fSinglePheFit);
839
840 fProb = fSinglePheFit->GetProb();
841 fChisquare = fSinglePheFit->GetChisquare();
842 fNdf = fSinglePheFit->GetNDF();
843
844 // Perform the cross-check fitting only the pedestal:
845 fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.);
846 fHBlindPixelCharge->Fit(fSinglePhePedFit,opt);
847
848 const Stat_t entries = fHBlindPixelCharge->Integral("width");
849
850 Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2);
851 fLambdaCheck = TMath::Log(entries/pedarea);
852 fLambdaCheckErr = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
853 + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);
854
855 *fLog << inf << "Results of the Blind Pixel Fit: " << endl;
856 *fLog << inf << "Chisquare: " << fChisquare << endl;
857 *fLog << inf << "DoF: " << fNdf << endl;
858 *fLog << inf << "Probability: " << fProb << endl;
859
860 //
861 // The fit result is accepted under condition that:
862 // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
863 // 2) at least 50 events are in the single Photo-electron peak
864 //
865 if (fProb < gkProbLimit)
866 {
867 *fLog << err << "ERROR: Fit Probability " << fProb
868 << " is smaller than the allowed value: " << gkProbLimit << endl;
869 CLRBIT(fFlags,kFitOK);
870 return kFALSE;
871 }
872
873 Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
874
875 if (contSinglePhe < 50.)
876 {
877 *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe
878 << " in the Single Photo-Electron peak " << endl;
879 CLRBIT(fFlags,kFitOK);
880 return kFALSE;
881 }
882 else
883 *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
884
885 SETBIT(fFlags,kFitOK);
886
887 return kTRUE;
888}
889
890
891void MHCalibrationBlindPixel::CutAllEdges()
892{
893
894 Int_t nbins = 20;
895
896 CutEdges(fHBlindPixelCharge,nbins);
897
898 fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
899 fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
900
901}
902
903Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt)
904{
905
906 rmin = (rmin != 0.) ? rmin : 4.;
907 rmax = (rmax != 0.) ? rmax : 9.;
908
909 const Stat_t entries = fHBlindPixelTime->Integral();
910 const Double_t mu_guess = fHBlindPixelTime->GetBinCenter(fHBlindPixelTime->GetMaximumBin());
911 const Double_t sigma_guess = (rmax - rmin)/2.;
912 const Double_t area_guess = entries/gkSq2Pi;
913
914 fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax);
915 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
916 fTimeGausFit->SetParNames("Area","#mu","#sigma");
917 fTimeGausFit->SetParLimits(0,0.,entries);
918 fTimeGausFit->SetParLimits(1,rmin,rmax);
919 fTimeGausFit->SetParLimits(2,0.,rmax-rmin);
920
921 fHBlindPixelTime->Fit(fTimeGausFit,opt);
922 rmin = fTimeGausFit->GetParameter(1) - 2.*fTimeGausFit->GetParameter(2);
923 rmax = fTimeGausFit->GetParameter(1) + 2.*fTimeGausFit->GetParameter(2);
924 fTimeGausFit->SetRange(rmin,rmax);
925
926 fHBlindPixelTime->Fit(fTimeGausFit,opt);
927
928 fMeanTime = fTimeGausFit->GetParameter(2);
929 fSigmaTime = fTimeGausFit->GetParameter(3);
930 fMeanTimeErr = fTimeGausFit->GetParError(2);
931 fSigmaTimeErr = fTimeGausFit->GetParError(3);
932
933 *fLog << inf << "Results of the Times Fit: " << endl;
934 *fLog << inf << "Chisquare: " << fTimeGausFit->GetChisquare() << endl;
935 *fLog << inf << "Ndf: " << fTimeGausFit->GetNDF() << endl;
936
937 return kTRUE;
938
939}
Note: See TracBrowser for help on using the repository browser.