source: trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc@ 2668

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