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

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