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

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