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

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