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