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

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