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

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