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

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