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

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