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

Last change on this file since 2879 was 2852, checked in by gaug, 21 years ago
*** empty log message ***
File size: 16.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//////////////////////////////////////////////////////////////////////////////
37#include "MHCalibrationBlindPixel.h"
38#include "MHCalibrationConfig.h"
39
40#include <TStyle.h>
41#include <TMath.h>
42#include <TPad.h>
43
44#include <TMinuit.h>
45#include <TFitter.h>
46
47#include <TF1.h>
48#include <TH2.h>
49#include <TCanvas.h>
50#include <TPaveText.h>
51#include <TRandom.h>
52
53#include <TText.h>
54
55#include "MBinning.h"
56#include "MParList.h"
57
58#include "MLog.h"
59#include "MLogManip.h"
60
61ClassImp(MHCalibrationBlindPixel);
62
63using namespace std;
64// --------------------------------------------------------------------------
65//
66// Default Constructor.
67//
68MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title)
69 : fFitLegend(NULL), fFitOK(kFALSE),
70 fLambda(0.), fMu0(0.), fMu1(0.), fSigma0(0.), fSigma1(0.),
71 fLambdaErr(0.), fMu0Err(0.), fMu1Err(0.), fSigma0Err(0.), fSigma1Err(0.),
72 fChisquare(0.), fProb(0.), fNdf(0),
73 fMeanTime(0.), fMeanTimeErr(0.), fSigmaTime(0.), fSigmaTimeErr(0.),
74 fLambdaCheck(0.), fLambdaCheckErr(0.),
75 fFitFunc(kEPoisson4)
76{
77
78 fName = name ? name : "MHCalibrationBlindPixel";
79 fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits";
80
81 // Create a large number of bins, later we will rebin
82 fBlindPixelChargefirst = -1000.;
83 fBlindPixelChargelast = gkStartBlindPixelBinNr;
84 fBlindPixelChargenbins = gkStartBlindPixelBinNr+(int)fBlindPixelChargefirst;
85
86 fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices",
87 fBlindPixelChargenbins,fBlindPixelChargefirst,fBlindPixelChargelast);
88 fHBlindPixelCharge->SetXTitle("Sum FADC Slices");
89 fHBlindPixelCharge->SetYTitle("Nr. of events");
90 fHBlindPixelCharge->Sumw2();
91 fHBlindPixelCharge->SetDirectory(NULL);
92
93 Axis_t tfirst = -0.5;
94 Axis_t tlast = 15.5;
95 Int_t nbins = 16;
96
97 fHBlindPixelTime = new TH1I("HBlindPixelTime","Distribution of Mean Arrival Times",nbins,tfirst,tlast);
98 fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
99 fHBlindPixelTime->SetYTitle("Nr. of events");
100 fHBlindPixelTime->Sumw2();
101 fHBlindPixelTime->SetDirectory(NULL);
102
103 // We define a reasonable number and later enlarge it if necessary
104 nbins = 20000;
105 Axis_t nfirst = -0.5;
106 Axis_t nlast = (Axis_t)nbins - 0.5;
107
108 fHBlindPixelChargevsN = new TH1I("HBlindPixelChargevsN","Sum of Charges vs. Event Number",nbins,nfirst,nlast);
109 fHBlindPixelChargevsN->SetXTitle("Event Nr.");
110 fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices");
111 fHBlindPixelChargevsN->SetDirectory(NULL);
112
113}
114
115MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
116{
117
118 if (fFitLegend)
119 delete fFitLegend;
120
121 delete fHBlindPixelCharge;
122 delete fHBlindPixelTime;
123 delete fHBlindPixelChargevsN;
124
125}
126
127
128
129void MHCalibrationBlindPixel::ResetBin(Int_t i)
130{
131 fHBlindPixelCharge->SetBinContent (i, 1.e-20);
132 fHBlindPixelTime->SetBinContent(i, 1.e-20);
133}
134
135
136// -------------------------------------------------------------------------
137//
138// Draw a legend with the fit results
139//
140void MHCalibrationBlindPixel::DrawLegend()
141{
142
143 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
144
145 if (fFitOK)
146 fFitLegend->SetFillColor(80);
147 else
148 fFitLegend->SetFillColor(2);
149
150 fFitLegend->SetLabel("Results of the single PhE Fit (to k=6):");
151 fFitLegend->SetTextSize(0.05);
152
153 const TString line1 =
154 Form("Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr());
155 TText *t1 = fFitLegend->AddText(line1.Data());
156 t1->SetBit(kCanDelete);
157
158 const TString line6 =
159 Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr());
160 TText *t2 = fFitLegend->AddText(line6.Data());
161 t2->SetBit(kCanDelete);
162
163 const TString line2 =
164 Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err());
165 TText *t3 = fFitLegend->AddText(line2.Data());
166 t3->SetBit(kCanDelete);
167
168 const TString line3 =
169 Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err());
170 TText *t4 = fFitLegend->AddText(line3.Data());
171 t4->SetBit(kCanDelete);
172
173 const TString line4 =
174 Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err());
175 TText *t5 = fFitLegend->AddText(line4.Data());
176 t5->SetBit(kCanDelete);
177
178 const TString line5 =
179 Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err());
180 TText *t6 = fFitLegend->AddText(line5.Data());
181 t6->SetBit(kCanDelete);
182
183 const TString line7 =
184 Form("#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf());
185 TText *t7 = fFitLegend->AddText(line7.Data());
186 t7->SetBit(kCanDelete);
187
188 const TString line8 =
189 Form("Probability: %4.2f ",GetProb());
190 TText *t8 = fFitLegend->AddText(line8.Data());
191 t8->SetBit(kCanDelete);
192
193 if (fFitOK)
194 {
195 TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
196 t->SetBit(kCanDelete);
197 }
198 else
199 {
200 TText *t = fFitLegend->AddText("Result of the Fit: NOT OK");
201 t->SetBit(kCanDelete);
202 }
203
204 fFitLegend->SetBit(kCanDelete);
205 fFitLegend->Draw();
206
207 return;
208}
209
210
211TObject *MHCalibrationBlindPixel::DrawClone(Option_t *opt) const
212{
213
214 gROOT->SetSelectedPad(NULL);
215
216 TObject *newobj = Clone();
217
218 newobj->SetBit(kCanDelete);
219 if (!newobj)
220 return 0;
221
222 // ((MHCalibrationBlindPixel*)newobj)->Draw(opt);
223 newobj->Draw(opt);
224 return newobj;
225}
226
227
228
229// -------------------------------------------------------------------------
230//
231// Draw the histogram
232//
233void MHCalibrationBlindPixel::Draw(Option_t *opt)
234{
235
236 gStyle->SetOptFit(0);
237 gStyle->SetOptStat(111111);
238
239 TCanvas *c = MakeDefCanvas(this,550,700);
240
241 c->Divide(2,2);
242
243 gROOT->SetSelectedPad(NULL);
244
245 c->cd(1);
246 gPad->SetLogy(1);
247 gPad->SetTicks();
248
249 fHBlindPixelCharge->DrawCopy(opt);
250
251 if (fFitOK)
252 fSinglePheFit.SetLineColor(kGreen);
253 else
254 fSinglePheFit.SetLineColor(kRed);
255
256 fSinglePheFit.DrawCopy("same");
257 c->Modified();
258 c->Update();
259
260 fSinglePhePedFit.SetLineColor(kBlue);
261 fSinglePhePedFit.DrawCopy("same");
262
263 c->cd(2);
264 DrawLegend();
265 c->Modified();
266 c->Update();
267
268 c->cd(3);
269 gPad->SetLogy(1);
270 gPad->SetBorderMode(0);
271 fHBlindPixelTime->DrawCopy(opt);
272
273 if (fHBlindPixelTime->GetFunction("GausTime"))
274 {
275 TF1 *tfit = fHBlindPixelTime->GetFunction("GausTime");
276 if (tfit->GetProb() < 0.01)
277 tfit->SetLineColor(kRed);
278 else
279 tfit->SetLineColor(kGreen);
280
281 tfit->DrawCopy("same");
282 c->Modified();
283 c->Update();
284 }
285
286 c->cd(4);
287
288 fHBlindPixelChargevsN->DrawCopy(opt);
289
290 c->Modified();
291 c->Update();
292}
293
294
295
296Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
297{
298 gRandom->SetSeed();
299
300 if (fHBlindPixelCharge->GetIntegral() != 0)
301 {
302 *fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl;
303 *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
304 return kFALSE;
305 }
306
307 TF1 simulateSinglePhe;
308
309 InitFit(simulateSinglePhe,fBlindPixelChargefirst,fBlindPixelChargelast);
310
311 for (Int_t i=0;i<10000; i++)
312 fHBlindPixelCharge->Fill(simulateSinglePhe.GetRandom());
313
314 return kTRUE;
315}
316
317void MHCalibrationBlindPixel::InitFit(TF1& f, Axis_t min, Axis_t max)
318{
319
320 //
321 // First guesses for the fit (should be as close to reality as possible,
322 // otherwise the fit goes gaga because of high number of dimensions ...
323 //
324 const Stat_t entries = fHBlindPixelCharge->Integral("width");
325 const Double_t lambda_guess = 0.5;
326 const Double_t mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
327 const Double_t si_0_guess = 20.;
328 const Double_t mu_1_guess = mu_0_guess + 50.;
329 const Double_t si_1_guess = si_0_guess + si_0_guess;
330 const Double_t norm = entries/gkSq2Pi;
331
332 //
333 // Initialize
334 //
335 switch (fFitFunc)
336 {
337
338 case kEPoisson4:
339 f = TF1("SinglePheFit",&fPoissonKto4,min,max,6);
340 f.SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
341 f.SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
342 f.SetParLimits(0,0.,1.);
343 f.SetParLimits(1,min,(max-min)/1.5);
344 f.SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
345 f.SetParLimits(3,1.0,(max-min)/2.0);
346 f.SetParLimits(4,1.0,(max-min)/2.5);
347 f.SetParLimits(5,norm-0.1,norm+0.1);
348 break;
349
350 case kEPoisson5:
351 f = TF1("SinglePheFit",&fPoissonKto5,min,max,6);
352 f.SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
353 f.SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
354 f.SetParLimits(0,0.,1.);
355 f.SetParLimits(1,min,(max-min)/1.5);
356 f.SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
357 f.SetParLimits(3,1.0,(max-min)/2.0);
358 f.SetParLimits(4,1.0,(max-min)/2.5);
359 f.SetParLimits(5,norm-0.1,norm+0.1);
360 break;
361
362 case kEPoisson6:
363 f = TF1("SinglePheFit",&fPoissonKto6,min,max,6);
364 f.SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
365 f.SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
366 f.SetParLimits(0,0.,1.);
367 f.SetParLimits(1,min,(max-min)/1.5);
368 f.SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
369 f.SetParLimits(3,1.0,(max-min)/2.0);
370 f.SetParLimits(4,1.0,(max-min)/2.5);
371 f.SetParLimits(5,norm-0.1,norm+0.1);
372 break;
373
374 case kEPoisson7:
375 break;
376
377 case kEPolya:
378 break;
379
380 case kEMichele:
381 break;
382
383 }
384}
385
386void MHCalibrationBlindPixel::ExitFit(TF1 &f)
387{
388
389
390 //
391 // Finalize
392 //
393 switch (fFitFunc)
394 {
395
396 case kEPoisson4:
397 case kEPoisson5:
398 case kEPoisson6:
399 case kEPoisson7:
400 fLambda = fSinglePheFit.GetParameter(0);
401 fMu0 = fSinglePheFit.GetParameter(1);
402 fMu1 = fSinglePheFit.GetParameter(2);
403 fSigma0 = fSinglePheFit.GetParameter(3);
404 fSigma1 = fSinglePheFit.GetParameter(4);
405
406 fLambdaErr = fSinglePheFit.GetParError(0);
407 fMu0Err = fSinglePheFit.GetParError(1);
408 fMu1Err = fSinglePheFit.GetParError(2);
409 fSigma0Err = fSinglePheFit.GetParError(3);
410 fSigma1Err = fSinglePheFit.GetParError(4);
411 break;
412 default:
413 break;
414 }
415
416}
417
418
419Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt)
420{
421
422 //
423 // Get the fitting ranges
424 //
425 rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
426 rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;
427
428 InitFit(fSinglePheFit,rmin,rmax);
429
430 fHBlindPixelCharge->Fit(&fSinglePheFit,opt);
431
432 ExitFit(fSinglePheFit);
433
434 fProb = fSinglePheFit.GetProb();
435 fChisquare = fSinglePheFit.GetChisquare();
436 fNdf = fSinglePheFit.GetNDF();
437
438 // Perform the cross-check fitting only the pedestal:
439 fSinglePhePedFit = TF1("GausPed","gaus",rmin,fMu0);
440 fHBlindPixelCharge->Fit(&fSinglePhePedFit,opt);
441
442 const Stat_t entries = fHBlindPixelCharge->Integral("width");
443
444 Double_t pedarea = fSinglePhePedFit.GetParameter(0)*gkSq2Pi*fSinglePhePedFit.GetParameter(2);
445 fLambdaCheck = TMath::Log(entries/pedarea);
446 fLambdaCheckErr = fSinglePhePedFit.GetParError(0)/fSinglePhePedFit.GetParameter(0)
447 + fSinglePhePedFit.GetParError(2)/fSinglePhePedFit.GetParameter(2);
448
449 *fLog << inf << "Results of the Blind Pixel Fit: " << endl;
450 *fLog << inf << "Chisquare: " << fChisquare << endl;
451 *fLog << inf << "DoF: " << fNdf << endl;
452 *fLog << inf << "Probability: " << fProb << endl;
453
454 //
455 // The fit result is accepted under condition that:
456 // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
457 // 2) at least 100 events are in the single Photo-electron peak
458 //
459 if (fProb < gkProbLimit)
460 {
461 *fLog << err << "ERROR: Fit Probability " << fProb
462 << " is smaller than the allowed value: " << gkProbLimit << endl;
463 fFitOK = kFALSE;
464 return kFALSE;
465 }
466
467 if (fProb < 0.01)
468 *fLog << warn << "WARNING: Fit Probability " << fProb << " is smaller than 1% " << endl;
469
470 Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
471
472 if (contSinglePhe < 100.)
473 {
474 *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe
475 << " in the Single Photo-Electron peak " << endl;
476 fFitOK = kFALSE;
477 return kFALSE;
478 }
479 else
480 *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
481
482 fFitOK = kTRUE;
483
484 return kTRUE;
485}
486
487
488void MHCalibrationBlindPixel::CutAllEdges()
489{
490
491 Int_t nbins = 30;
492
493 CutEdges(fHBlindPixelCharge,nbins);
494
495 fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
496 fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
497 fBlindPixelChargenbins = nbins;
498
499 CutEdges(fHBlindPixelChargevsN,0);
500
501}
502
503Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt)
504{
505
506 rmin = (rmin != 0.) ? rmin : 4.;
507 rmax = (rmax != 0.) ? rmax : 9.;
508
509 const Stat_t entries = fHBlindPixelTime->Integral();
510 const Double_t mu_guess = fHBlindPixelTime->GetBinCenter(fHBlindPixelTime->GetMaximumBin());
511 const Double_t sigma_guess = (rmax - rmin)/2.;
512 const Double_t area_guess = entries/gkSq2Pi;
513
514 fTimeGausFit = TF1("GausTime","gaus",rmin,rmax);
515 fTimeGausFit.SetParameters(area_guess,mu_guess,sigma_guess);
516 fTimeGausFit.SetParNames("Area","#mu","#sigma");
517 fTimeGausFit.SetParLimits(0,0.,entries);
518 fTimeGausFit.SetParLimits(1,rmin,rmax);
519 fTimeGausFit.SetParLimits(2,0.,rmax-rmin);
520
521 fHBlindPixelTime->Fit(&fTimeGausFit,opt);
522 rmin = fTimeGausFit.GetParameter(1) - 2.*fTimeGausFit.GetParameter(2);
523 rmax = fTimeGausFit.GetParameter(1) + 2.*fTimeGausFit.GetParameter(2);
524 fTimeGausFit.SetRange(rmin,rmax);
525
526 fHBlindPixelTime->Fit(&fTimeGausFit,opt);
527
528 fMeanTime = fTimeGausFit.GetParameter(2);
529 fSigmaTime = fTimeGausFit.GetParameter(3);
530 fMeanTimeErr = fTimeGausFit.GetParError(2);
531 fSigmaTimeErr = fTimeGausFit.GetParError(3);
532
533 *fLog << inf << "Results of the Times Fit: " << endl;
534 *fLog << inf << "Chisquare: " << fTimeGausFit.GetChisquare() << endl;
535 *fLog << inf << "Ndf: " << fTimeGausFit.GetNDF() << endl;
536
537 return kTRUE;
538
539}
Note: See TracBrowser for help on using the repository browser.