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

Last change on this file since 3211 was 3191, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 28.1 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 = kEPoisson5;
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 *fLog << inf << "min: " << min << endl;
690 *fLog << "max: " << max << endl;
691
692 //
693 // First guesses for the fit (should be as close to reality as possible,
694 // otherwise the fit goes gaga because of high number of dimensions ...
695 //
696 const Stat_t entries = fHBlindPixelCharge->Integral("width");
697 const Double_t lambda_guess = 0.1;
698 const Double_t maximum_bin = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
699 const Double_t norm = entries/gkSq2Pi;
700
701 //
702 // Initialize the fit function
703 //
704 switch (fFitFunc)
705 {
706 case kEPoisson4:
707 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
708 break;
709 case kEPoisson5:
710 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
711 break;
712 case kEPoisson6:
713 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
714 break;
715 case kEPolya:
716 fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8);
717 break;
718 case kEMichele:
719 break;
720
721 default:
722 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
723 return kFALSE;
724 break;
725 }
726
727 const Double_t mu_0_guess = maximum_bin;
728 const Double_t si_0_guess = 40.;
729 const Double_t mu_1_guess = mu_0_guess + 100.;
730 const Double_t si_1_guess = si_0_guess + si_0_guess;
731 // Michele
732 const Double_t lambda_1cat_guess = 0.5;
733 const Double_t lambda_1dyn_guess = 0.5;
734 const Double_t mu_1cat_guess = mu_0_guess + 50.;
735 const Double_t mu_1dyn_guess = mu_0_guess + 20.;
736 const Double_t si_1cat_guess = si_0_guess + si_0_guess;
737 const Double_t si_1dyn_guess = si_0_guess;
738 // Polya
739 const Double_t excessPoisson_guess = 0.5;
740 const Double_t delta1_guess = 8.;
741 const Double_t delta2_guess = 5.;
742 const Double_t electronicAmp_guess = fgBlindPixelElectronicAmp;
743 const Double_t electronicAmp_limit = fgBlindPixelElectronicAmpError;
744
745 *fLog << inf << "pedestal: " << fMeanPedestal << endl;
746 *fLog << "sigma: " << fSigmaPedestal << endl;
747
748 //
749 // Initialize boundaries and start parameters
750 //
751 switch (fFitFunc)
752 {
753
754 case kEPoisson4:
755 if ((fMeanPedestal) && (fSigmaPedestal))
756 fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm);
757 else
758 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
759
760 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
761
762 fSinglePheFit->SetParLimits(0,0.,0.5);
763
764 if ((fMeanPedestal) && (fSigmaPedestal))
765 fSinglePheFit->SetParLimits(1,
766 fMeanPedestal-5.*fMeanPedestalErr,
767 fMeanPedestal+5.*fMeanPedestalErr);
768 else
769 fSinglePheFit->SetParLimits(1,-3.,0.);
770
771 fSinglePheFit->SetParLimits(2,min,max);
772
773 if ((fMeanPedestal) && (fSigmaPedestal))
774 fSinglePheFit->SetParLimits(3,
775 fSigmaPedestal-5.*fSigmaPedestalErr,
776 fSigmaPedestal+5.*fSigmaPedestalErr);
777 else
778 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
779
780 fSinglePheFit->SetParLimits(4,0.,(max-min));
781 fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.5*norm));
782 break;
783 case kEPoisson5:
784 case kEPoisson6:
785 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
786 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
787 fSinglePheFit->SetParLimits(0,0.,1.);
788 fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
789 fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
790 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
791 fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
792 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
793 break;
794
795 case kEPolya:
796 if ((fMeanPedestal) && (fSigmaPedestal))
797 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
798 delta1_guess,delta2_guess,
799 electronicAmp_guess,
800 fSigmaPedestal,
801 norm,
802 fMeanPedestal);
803 else
804 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
805 delta1_guess,delta2_guess,
806 electronicAmp_guess,
807 si_0_guess,
808 norm, mu_0_guess);
809 fSinglePheFit->SetParNames("#lambda","b_{tot}",
810 "#delta_{1}","#delta_{2}",
811 "amp_{e}","#sigma_{0}",
812 "Area", "#mu_{0}");
813 fSinglePheFit->SetParLimits(0,0.,1.);
814 fSinglePheFit->SetParLimits(1,0.,1.);
815 fSinglePheFit->SetParLimits(2,6.,12.);
816 fSinglePheFit->SetParLimits(3,3.,8.);
817 fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit,
818 electronicAmp_guess+electronicAmp_limit);
819 if ((fMeanPedestal) && (fSigmaPedestal))
820 fSinglePheFit->SetParLimits(5,
821 fSigmaPedestal-3.*fSigmaPedestalErr,
822 fSigmaPedestal+3.*fSigmaPedestalErr);
823 else
824 fSinglePheFit->SetParLimits(5,min,(max-min)/1.5);
825
826 fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
827 if ((fMeanPedestal) && (fSigmaPedestal))
828 fSinglePheFit->SetParLimits(7,
829 fMeanPedestal-3.*fMeanPedestalErr,
830 fMeanPedestal+3.*fMeanPedestalErr);
831 else
832 fSinglePheFit->SetParLimits(7,-35.,15.);
833 break;
834
835 case kEMichele:
836
837
838 break;
839
840 default:
841 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
842 return kFALSE;
843 break;
844 }
845
846 return kTRUE;
847}
848
849void MHCalibrationBlindPixel::ExitFit(TF1 *f)
850{
851
852
853 //
854 // Finalize
855 //
856 switch (fFitFunc)
857 {
858
859 case kEPoisson4:
860 case kEPoisson5:
861 case kEPoisson6:
862 case kEPoisson7:
863 fLambda = fSinglePheFit->GetParameter(0);
864 fMu0 = fSinglePheFit->GetParameter(1);
865 fMu1 = fSinglePheFit->GetParameter(2);
866 fSigma0 = fSinglePheFit->GetParameter(3);
867 fSigma1 = fSinglePheFit->GetParameter(4);
868
869 fLambdaErr = fSinglePheFit->GetParError(0);
870 fMu0Err = fSinglePheFit->GetParError(1);
871 fMu1Err = fSinglePheFit->GetParError(2);
872 fSigma0Err = fSinglePheFit->GetParError(3);
873 fSigma1Err = fSinglePheFit->GetParError(4);
874 break;
875 case kEPolya:
876 fLambda = fSinglePheFit->GetParameter(0);
877 fMu0 = fSinglePheFit->GetParameter(7);
878 fMu1 = 0.;
879 fSigma0 = fSinglePheFit->GetParameter(5);
880 fSigma1 = 0.;
881
882 fLambdaErr = fSinglePheFit->GetParError(0);
883 fMu0Err = fSinglePheFit->GetParError(7);
884 fMu1Err = 0.;
885 fSigma0Err = fSinglePheFit->GetParError(5);
886 fSigma1Err = 0.;
887 default:
888 break;
889 }
890
891}
892
893
894Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt)
895{
896
897 //
898 // Get the fitting ranges
899 //
900 rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
901 rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;
902 if (!InitFit(rmin,rmax))
903 return kFALSE;
904
905 fHBlindPixelCharge->Fit(fSinglePheFit,opt);
906
907 ExitFit(fSinglePheFit);
908
909 fProb = fSinglePheFit->GetProb();
910 fChisquare = fSinglePheFit->GetChisquare();
911 fNdf = fSinglePheFit->GetNDF();
912
913 // Perform the cross-check fitting only the pedestal:
914 fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.);
915 fHBlindPixelCharge->Fit(fSinglePhePedFit,opt);
916
917 const Stat_t entries = fHBlindPixelCharge->Integral("width");
918
919 Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2);
920 fLambdaCheck = TMath::Log(entries/pedarea);
921 fLambdaCheckErr = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
922 + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);
923
924 *fLog << inf << "Results of the Blind Pixel Fit: " << endl;
925 *fLog << "Chisquare: " << fChisquare << endl;
926 *fLog << "DoF: " << fNdf << endl;
927 *fLog << "Probability: " << fProb << endl;
928
929 //
930 // The fit result is accepted under condition that:
931 // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
932 // 2) at least 50 events are in the single Photo-electron peak
933 //
934 if (fProb < gkProbLimit)
935 {
936 *fLog << warn << "WARNING - Fit Probability " << fProb
937 << " is smaller than the allowed value: " << gkProbLimit << endl;
938 CLRBIT(fFlags,kFitOK);
939 return kFALSE;
940 }
941
942 Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
943
944 if (contSinglePhe < 50.)
945 {
946 *fLog << warn << "WARNING - Statistics is too low: Only " << contSinglePhe
947 << " in the Single Photo-Electron peak " << endl;
948 CLRBIT(fFlags,kFitOK);
949 return kFALSE;
950 }
951 else
952 *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
953
954 SETBIT(fFlags,kFitOK);
955
956 return kTRUE;
957}
958
959
960void MHCalibrationBlindPixel::CutAllEdges()
961{
962
963 Int_t nbins = 25;
964
965 StripZeros(fHBlindPixelCharge,nbins);
966
967 fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
968 fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
969
970}
971
972Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt)
973{
974
975 rmin = (rmin != 0.) ? rmin : 4.;
976 rmax = (rmax != 0.) ? rmax : 9.;
977
978 const Stat_t entries = fHBlindPixelTime->Integral();
979 const Double_t mu_guess = fHBlindPixelTime->GetBinCenter(fHBlindPixelTime->GetMaximumBin());
980 const Double_t sigma_guess = (rmax - rmin)/2.;
981 const Double_t area_guess = entries/gkSq2Pi;
982
983 fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax);
984 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
985 fTimeGausFit->SetParNames("Area","#mu","#sigma");
986 fTimeGausFit->SetParLimits(0,0.,entries);
987 fTimeGausFit->SetParLimits(1,rmin,rmax);
988 fTimeGausFit->SetParLimits(2,0.,rmax-rmin);
989
990 fHBlindPixelTime->Fit(fTimeGausFit,opt);
991 rmin = fTimeGausFit->GetParameter(1) - 2.*fTimeGausFit->GetParameter(2);
992 rmax = fTimeGausFit->GetParameter(1) + 2.*fTimeGausFit->GetParameter(2);
993 fTimeGausFit->SetRange(rmin,rmax);
994
995 fHBlindPixelTime->Fit(fTimeGausFit,opt);
996
997 fMeanTime = fTimeGausFit->GetParameter(2);
998 fSigmaTime = fTimeGausFit->GetParameter(3);
999 fMeanTimeErr = fTimeGausFit->GetParError(2);
1000 fSigmaTimeErr = fTimeGausFit->GetParError(3);
1001
1002 *fLog << inf << "Results of the Times Fit: " << endl;
1003 *fLog << "Chisquare: " << fTimeGausFit->GetChisquare() << endl;
1004 *fLog << "Ndf: " << fTimeGausFit->GetNDF() << endl;
1005
1006 return kTRUE;
1007
1008}
1009
Note: See TracBrowser for help on using the repository browser.