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

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