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

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