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

Last change on this file since 2880 was 2880, checked in by gaug, 21 years ago
*** empty log message ***
File size: 17.4 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 // TVirtualPad *padsav = gPad;
236
237 // TObject *newobj = (TObject *)IsA()->New();
238 // Copy(*newobj);
239 // TObject *newobj = Clone();
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 // ((MHCalibrationBlindPixel*)newobj)->Draw(option);
249 // newobj->Draw(option);
250 else
251 newobj->Draw(GetDrawOption());
252 // ((MHCalibrationBlindPixel*)newobj)->Draw(GetDrawOption());
253 // newobj->Draw(GetDrawOption());
254 // if (padsav) padsav->cd();
255
256 return newobj;
257}
258
259
260// -------------------------------------------------------------------------
261//
262// Draw the histogram
263//
264void MHCalibrationBlindPixel::Draw(Option_t *opt)
265{
266
267 gStyle->SetOptFit(0);
268 gStyle->SetOptStat(111111);
269
270 TCanvas *c = MakeDefCanvas(this,550,700);
271
272 c->Divide(2,3);
273
274 gROOT->SetSelectedPad(NULL);
275
276 c->cd(1);
277 gPad->SetLogx(0);
278 gPad->SetLogy(1);
279 gPad->SetTicks();
280
281 fHBlindPixelCharge->Draw(opt);
282
283 if (fFitOK)
284 fSinglePheFit->SetLineColor(kGreen);
285 else
286 fSinglePheFit->SetLineColor(kRed);
287
288 fSinglePheFit->Draw("same");
289 c->Modified();
290 c->Update();
291
292 fSinglePhePedFit->SetLineColor(kBlue);
293 fSinglePhePedFit->Draw("same");
294
295 c->cd(2);
296 DrawLegend();
297 c->Modified();
298 c->Update();
299
300 c->cd(3);
301 gPad->SetLogy(1);
302 gPad->SetBorderMode(0);
303 fHBlindPixelTime->Draw(opt);
304
305 c->cd(4);
306
307 fHBlindPixelChargevsN->Draw(opt);
308
309 c->Modified();
310 c->Update();
311
312 c->cd(5);
313
314 // fHBlindPixelPSD = (TH1F*)MFFT::PowerSpectrumDensity(fHBlindPixelChargevsN);
315 // TH1F *hist = MFFT::PowerSpectrumDensity((*this)->fHBlindPixelChargevsN);
316
317 // fHBlindPixelPSD->Draw(opt);
318 c->Modified();
319 c->Update();
320
321}
322
323
324
325Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1)
326{
327 gRandom->SetSeed();
328
329 if (fHBlindPixelCharge->GetIntegral() != 0)
330 {
331 *fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl;
332 *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
333 return kFALSE;
334 }
335
336 if (!InitFit(fBlindPixelChargefirst,fBlindPixelChargelast))
337 return kFALSE;
338
339 for (Int_t i=0;i<10000; i++)
340 fHBlindPixelCharge->Fill(fSinglePheFit->GetRandom());
341
342 return kTRUE;
343}
344
345Bool_t MHCalibrationBlindPixel::InitFit(Axis_t min, Axis_t max)
346{
347
348 //
349 // First guesses for the fit (should be as close to reality as possible,
350 // otherwise the fit goes gaga because of high number of dimensions ...
351 //
352 const Stat_t entries = fHBlindPixelCharge->Integral("width");
353 const Double_t lambda_guess = 0.5;
354 const Double_t mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
355 const Double_t si_0_guess = 20.;
356 const Double_t mu_1_guess = mu_0_guess + 50.;
357 const Double_t si_1_guess = si_0_guess + si_0_guess;
358 const Double_t norm = entries/gkSq2Pi;
359
360 //
361 // Initialize
362 //
363 switch (fFitFunc)
364 {
365
366 case kEPoisson4:
367 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
368 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
369 fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
370 fSinglePheFit->SetParLimits(0,0.,1.);
371 fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
372 fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
373 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
374 fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
375 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
376 break;
377
378 case kEPoisson5:
379 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
380 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
381 fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
382 fSinglePheFit->SetParLimits(0,0.,1.);
383 fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
384 fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
385 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
386 fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
387 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
388 break;
389
390 case kEPoisson6:
391 fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
392 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
393 fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
394 fSinglePheFit->SetParLimits(0,0.,1.);
395 fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
396 fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
397 fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
398 fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
399 fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
400 break;
401
402 case kEPoisson7:
403 break;
404
405 case kEPolya:
406 break;
407
408 case kEMichele:
409 break;
410
411 default:
412 *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
413 return kFALSE;
414 break;
415 }
416
417 return kTRUE;
418}
419
420void MHCalibrationBlindPixel::ExitFit(TF1 *f)
421{
422
423
424 //
425 // Finalize
426 //
427 switch (fFitFunc)
428 {
429
430 case kEPoisson4:
431 case kEPoisson5:
432 case kEPoisson6:
433 case kEPoisson7:
434 fLambda = fSinglePheFit->GetParameter(0);
435 fMu0 = fSinglePheFit->GetParameter(1);
436 fMu1 = fSinglePheFit->GetParameter(2);
437 fSigma0 = fSinglePheFit->GetParameter(3);
438 fSigma1 = fSinglePheFit->GetParameter(4);
439
440 fLambdaErr = fSinglePheFit->GetParError(0);
441 fMu0Err = fSinglePheFit->GetParError(1);
442 fMu1Err = fSinglePheFit->GetParError(2);
443 fSigma0Err = fSinglePheFit->GetParError(3);
444 fSigma1Err = fSinglePheFit->GetParError(4);
445 break;
446 default:
447 break;
448 }
449
450}
451
452
453Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt)
454{
455
456 //
457 // Get the fitting ranges
458 //
459 rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst;
460 rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast;
461
462 if (!InitFit(rmin,rmax))
463 return kFALSE;
464
465 fHBlindPixelCharge->Fit(fSinglePheFit,opt);
466
467
468 ExitFit(fSinglePheFit);
469
470 fProb = fSinglePheFit->GetProb();
471 fChisquare = fSinglePheFit->GetChisquare();
472 fNdf = fSinglePheFit->GetNDF();
473
474 // Perform the cross-check fitting only the pedestal:
475 fSinglePhePedFit = new TF1("GausPed","gaus",rmin,fMu0);
476 fHBlindPixelCharge->Fit(fSinglePhePedFit,opt);
477
478 const Stat_t entries = fHBlindPixelCharge->Integral("width");
479
480 Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2);
481 fLambdaCheck = TMath::Log(entries/pedarea);
482 fLambdaCheckErr = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
483 + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);
484
485 *fLog << inf << "Results of the Blind Pixel Fit: " << endl;
486 *fLog << inf << "Chisquare: " << fChisquare << endl;
487 *fLog << inf << "DoF: " << fNdf << endl;
488 *fLog << inf << "Probability: " << fProb << endl;
489
490 //
491 // The fit result is accepted under condition that:
492 // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
493 // 2) at least 100 events are in the single Photo-electron peak
494 //
495 if (fProb < gkProbLimit)
496 {
497 *fLog << err << "ERROR: Fit Probability " << fProb
498 << " is smaller than the allowed value: " << gkProbLimit << endl;
499 fFitOK = kFALSE;
500 return kFALSE;
501 }
502
503 if (fProb < 0.01)
504 *fLog << warn << "WARNING: Fit Probability " << fProb << " is smaller than 1% " << endl;
505
506 Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
507
508 if (contSinglePhe < 100.)
509 {
510 *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe
511 << " in the Single Photo-Electron peak " << endl;
512 fFitOK = kFALSE;
513 return kFALSE;
514 }
515 else
516 *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
517
518 fFitOK = kTRUE;
519
520 return kTRUE;
521}
522
523
524void MHCalibrationBlindPixel::CutAllEdges()
525{
526
527 Int_t nbins = 30;
528
529 CutEdges(fHBlindPixelCharge,nbins);
530
531 fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
532 fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
533 fBlindPixelChargenbins = nbins;
534
535 CutEdges(fHBlindPixelChargevsN,0);
536
537}
538
539Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt)
540{
541
542 rmin = (rmin != 0.) ? rmin : 4.;
543 rmax = (rmax != 0.) ? rmax : 9.;
544
545 const Stat_t entries = fHBlindPixelTime->Integral();
546 const Double_t mu_guess = fHBlindPixelTime->GetBinCenter(fHBlindPixelTime->GetMaximumBin());
547 const Double_t sigma_guess = (rmax - rmin)/2.;
548 const Double_t area_guess = entries/gkSq2Pi;
549
550 fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax);
551 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
552 fTimeGausFit->SetParNames("Area","#mu","#sigma");
553 fTimeGausFit->SetParLimits(0,0.,entries);
554 fTimeGausFit->SetParLimits(1,rmin,rmax);
555 fTimeGausFit->SetParLimits(2,0.,rmax-rmin);
556
557 fHBlindPixelTime->Fit(fTimeGausFit,opt);
558 rmin = fTimeGausFit->GetParameter(1) - 2.*fTimeGausFit->GetParameter(2);
559 rmax = fTimeGausFit->GetParameter(1) + 2.*fTimeGausFit->GetParameter(2);
560 fTimeGausFit->SetRange(rmin,rmax);
561
562 fHBlindPixelTime->Fit(fTimeGausFit,opt);
563
564 fMeanTime = fTimeGausFit->GetParameter(2);
565 fSigmaTime = fTimeGausFit->GetParameter(3);
566 fMeanTimeErr = fTimeGausFit->GetParError(2);
567 fSigmaTimeErr = fTimeGausFit->GetParError(3);
568
569 *fLog << inf << "Results of the Times Fit: " << endl;
570 *fLog << inf << "Chisquare: " << fTimeGausFit->GetChisquare() << endl;
571 *fLog << inf << "Ndf: " << fTimeGausFit->GetNDF() << endl;
572
573 return kTRUE;
574
575}
Note: See TracBrowser for help on using the repository browser.