source: trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc@ 2549

Last change on this file since 2549 was 2544, checked in by gaug, 22 years ago
*** empty log message ***
File size: 10.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// MHCalibrationPixel //
28// //
29// Performs all the necessary fits to extract the mean number of photons //
30// out of the derived light flux //
31// //
32//////////////////////////////////////////////////////////////////////////////
33
34#include "MHCalibrationPixel.h"
35#include "MHCalibrationConfig.h"
36#include "MCalibrationFits.h"
37
38#include <TStyle.h>
39#include <TMath.h>
40
41#include <TFitter.h>
42
43#include <TF1.h>
44#include <TH2.h>
45#include <TCanvas.h>
46#include <TPad.h>
47#include <TPaveText.h>
48
49#include "MParList.h"
50
51#include "MLog.h"
52#include "MLogManip.h"
53
54ClassImp(MHCalibrationPixel);
55
56using namespace std;
57// --------------------------------------------------------------------------
58//
59// Default Constructor.
60//
61MHCalibrationPixel::MHCalibrationPixel(Int_t pix, const char *name, const char *title)
62 : fFitOK(kFALSE), fPixId(pix), fTGausFit(NULL), fQGausFit(NULL), fFitLegend(NULL)
63{
64
65 fName = name ? name : "MHCalibrationPixel";
66 fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits";
67
68 TString qname = "HQ";
69 TString qtitle = "Distribution of Summed FADC Slices Pixel ";
70 qname += pix;
71 qtitle += pix;
72
73 // Create a large number of bins, later we will rebin
74 fQfirst = -0.5;
75 fQlast = gkStartQlast - 0.5;
76 fQnbins = gkStartPixelBinNr;
77
78 fHQ = new TH1I( qname.Data(),qtitle.Data(),
79 fQnbins,fQfirst,fQlast);
80 fHQ->SetXTitle("Sum FADC Slices");
81 fHQ->SetYTitle("Nr. of events");
82 fHQ->Sumw2();
83
84 TString tname = "HT";
85 TString ttitle = "Distribution of Mean Arrival Times Pixel ";
86 tname += pix;
87 ttitle += pix;
88
89 Axis_t tfirst = -0.5;
90 Axis_t tlast = 15.5;
91 Int_t nbins = 16;
92
93 fHT = new TH1I(tname.Data(),ttitle.Data(),
94 nbins,tfirst,tlast);
95 fHT->SetXTitle("Mean Arrival Times [FADC slice nr]");
96 fHT->SetYTitle("Nr. of events");
97 fHT->Sumw2();
98
99 TString qvsnname = "HQvsN";
100 TString qvsntitle = "Sum of Charges vs. Event Number Pixel ";
101 qvsnname += pix;
102 qvsntitle += pix;
103
104 // We define a reasonable number and later enlarge it if necessary
105 nbins = 20000;
106 Axis_t nfirst = -0.5;
107 Axis_t nlast = (Axis_t)nbins - 0.5;
108
109 fHQvsN = new TH1I(qvsnname.Data(),qvsntitle.Data(),
110 nbins,nfirst,nlast);
111 fHQvsN->SetXTitle("Event Nr.");
112 fHQvsN->SetYTitle("Sum of FADC slices");
113
114 fQChisquare = -1.;
115 fQProb = -1.;
116 fQNdf = -1;
117
118 fTChisquare = -1.;
119 fTProb = -1.;
120 fTNdf = -1;
121
122}
123
124MHCalibrationPixel::~MHCalibrationPixel()
125{
126
127 delete fHQ;
128 delete fHT;
129 delete fHQvsN;
130
131 if (fQGausFit)
132 delete fQGausFit;
133 if (fTGausFit)
134 delete fTGausFit;
135 if (fFitLegend)
136 delete fFitLegend;
137
138}
139
140void MHCalibrationPixel::Reset()
141{
142
143 for (Int_t i = fHQ->FindBin(fQfirst); i <= fHQ->FindBin(fQlast); i++)
144 fHQ->SetBinContent(i, 1.e-20);
145
146 for (Int_t i = 0; i < 16; i++)
147 fHT->SetBinContent(i, 1.e-20);
148
149 fQlast = gkStartQlast;
150
151 fHQ->GetXaxis()->SetRangeUser(0.,fQlast);
152
153 return;
154}
155
156
157// -------------------------------------------------------------------------
158//
159// Set the binnings and prepare the filling of the histograms
160//
161Bool_t MHCalibrationPixel::SetupFill(const MParList *plist)
162{
163
164 fHQ->Reset();
165 fHT->Reset();
166
167 return kTRUE;
168}
169
170
171
172// -------------------------------------------------------------------------
173//
174// Draw a legend with the fit results
175//
176void MHCalibrationPixel::DrawLegend()
177{
178
179 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
180
181 if (fFitOK)
182 fFitLegend->SetFillColor(80);
183 else
184 fFitLegend->SetFillColor(2);
185
186 fFitLegend->SetLabel("Results of the Gauss Fit:");
187 fFitLegend->SetTextSize(0.05);
188
189 char line1[32];
190 sprintf(line1,"Mean: Q_{#mu} = %2.2f #pm %2.2f",fQMean,fQMeanErr);
191 fFitLegend->AddText(line1);
192
193 char line4[32];
194 sprintf(line4,"Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fQSigma,fQSigmaErr);
195 fFitLegend->AddText(line4);
196
197 char line7[32];
198 sprintf(line7,"#chi^{2} / N_{dof}: %4.2f / %3i",fQChisquare,fQNdf);
199 fFitLegend->AddText(line7);
200
201 char line8[32];
202 sprintf(line8,"Probability: %4.2f ",fQProb);
203 fFitLegend->AddText(line8);
204
205 if (fFitOK)
206 fFitLegend->AddText("Result of the Fit: OK");
207 else
208 fFitLegend->AddText("Result of the Fit: NOT OK");
209
210 fFitLegend->SetBit(kCanDelete);
211 fFitLegend->Draw();
212
213}
214
215
216// -------------------------------------------------------------------------
217//
218// Draw the histogram
219//
220void MHCalibrationPixel::Draw(Option_t *opt)
221{
222
223 gStyle->SetOptFit(0);
224 gStyle->SetOptStat(1111111);
225
226 TCanvas *c = MakeDefCanvas(this,600,900);
227
228 gROOT->SetSelectedPad(NULL);
229
230 c->Divide(2,2);
231
232 c->cd(1);
233 gPad->SetBorderMode(0);
234 gPad->SetLogy(1);
235 gPad->SetTicks();
236
237 fHQ->DrawCopy(opt);
238
239 if (fQGausFit)
240 {
241 if (fFitOK)
242 fQGausFit->SetLineColor(kGreen);
243 else
244 fQGausFit->SetLineColor(kRed);
245
246 fQGausFit->DrawCopy("same");
247 c->Modified();
248 c->Update();
249 }
250
251
252 c->cd(2);
253 DrawLegend();
254 c->Update();
255
256 c->cd(3);
257 gStyle->SetOptStat(1111111);
258
259 gPad->SetLogy(1);
260 fHT->DrawCopy(opt);
261
262 if (fHT->GetFunction("GausTime"))
263 {
264 TF1 *tfit = fHT->GetFunction("GausTime");
265 if (tfit->GetProb() < 0.01)
266 tfit->SetLineColor(kRed);
267 else
268 tfit->SetLineColor(kGreen);
269
270 tfit->DrawCopy("same");
271 c->Modified();
272 c->Update();
273 }
274
275 c->Modified();
276 c->Update();
277
278 c->cd(4);
279 fHQvsN->DrawCopy(opt);
280}
281
282
283
284Bool_t MHCalibrationPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *option)
285{
286
287 if (fTGausFit)
288 return kFALSE;
289
290 rmin = (rmin != 0.) ? rmin : -0.5;
291 rmax = (rmax != 0.) ? rmax : 15.5;
292
293 const Stat_t entries = fHT->GetEntries();
294 const Double_t mu_guess = fHT->GetBinCenter(fHT->GetMaximumBin());
295 const Double_t sigma_guess = (rmax - rmin)/2.;
296 const Double_t area_guess = entries/gkSq2Pi;
297
298 fTGausFit = new TF1("GausTime","gaus",rmin,rmax);
299
300 if (!fTGausFit)
301 {
302 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
303 return kFALSE;
304 }
305
306 fTGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
307 fTGausFit->SetParNames("Area","#mu","#sigma");
308 fTGausFit->SetParLimits(0,0.,entries);
309 fTGausFit->SetParLimits(1,rmin,rmax);
310 fTGausFit->SetParLimits(2,0.,rmax-rmin);
311
312 fHT->Fit("GausTime",option);
313
314 fTChisquare = fTGausFit->GetChisquare();
315 fTNdf = fTGausFit->GetNDF();
316 fTProb = fTGausFit->GetProb();
317 fTMean = fTGausFit->GetParameter(1);
318 fTSigma = fTGausFit->GetParameter(2);
319
320 if (fTProb < gkProbLimit)
321 {
322 *fLog << warn << "Fit of the Arrival times failed ! " << endl;
323 return kFALSE;
324 }
325
326 return kTRUE;
327
328}
329
330Bool_t MHCalibrationPixel::FitQ(Axis_t rmin, Axis_t rmax, Option_t *option)
331{
332
333 if (fQGausFit)
334 return kFALSE;
335
336 //
337 // Get the fitting ranges
338 //
339 rmin = (rmin != 0.) ? rmin : fQfirst;
340 rmax = (rmax != 0.) ? rmax : fQlast;
341
342 //
343 // First guesses for the fit (should be as close to reality as possible,
344 // otherwise the fit goes gaga because of high number of dimensions ...
345 //
346 const Stat_t entries = fHQ->GetEntries();
347 const Double_t ar_guess = entries/gkSq2Pi;
348 const Double_t mu_guess = fHQ->GetBinCenter(fHQ->GetMaximumBin());
349 const Double_t si_guess = mu_guess/500.;
350
351 fQGausFit = new TF1("QGausFit","gaus",rmin,rmax);
352
353 if (!fQGausFit)
354 {
355 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
356 return kFALSE;
357 }
358
359 fQGausFit->SetParameters(ar_guess,mu_guess,si_guess);
360 fQGausFit->SetParNames("Area","#mu","#sigma");
361 fQGausFit->SetParLimits(0,0.,entries);
362 fQGausFit->SetParLimits(1,rmin,rmax);
363 fQGausFit->SetParLimits(2,0.,rmax-rmin);
364
365 fHQ->Fit("QGausFit",option);
366
367 fQChisquare = fQGausFit->GetChisquare();
368 fQNdf = fQGausFit->GetNDF();
369 fQProb = fQGausFit->GetProb();
370 fQMean = fQGausFit->GetParameter(1);
371 fQMeanErr = fQGausFit->GetParError(1);
372 fQSigma = fQGausFit->GetParameter(2);
373 fQSigmaErr = fQGausFit->GetParError(2);
374
375 //
376 // The fit result is accepted under condition
377 // The Probability is greater than gkProbLimit (default 0.01 == 99%)
378 //
379 if (fQProb < gkProbLimit)
380 {
381 *fLog << warn << "Prob: " << fQProb << " is smaller than the allowed value: " << gkProbLimit << endl;
382 fFitOK = kFALSE;
383 return kFALSE;
384 }
385
386
387 fFitOK = kTRUE;
388
389 return kTRUE;
390}
391
392
393void MHCalibrationPixel::CutAllEdges()
394{
395
396 //
397 // The number 100 is necessary because it is the internal binning
398 // of ROOT functions. A call to SetNpx() does NOT help
399 // If you find another solution which WORKS!!, please tell me!!
400 //
401 Int_t nbins = 100;
402
403 CutEdges(fHQ,nbins);
404
405 fQfirst = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetFirst());
406 fQlast = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetLast())+fHQ->GetBinWidth(0);
407 fQnbins = nbins;
408
409 CutEdges(fHQvsN,0);
410
411}
412
413void MHCalibrationPixel::PrintQFitResult()
414{
415
416 *fLog << "Results of the Summed Charges Fit: " << endl;
417 *fLog << "Chisquare: " << fQChisquare << endl;
418 *fLog << "DoF: " << fQNdf << endl;
419 *fLog << "Probability: " << fQProb << endl;
420 *fLog << endl;
421
422}
423
424void MHCalibrationPixel::PrintTFitResult()
425{
426
427 *fLog << "Results of the Arrival Time Slices Fit: " << endl;
428 *fLog << "Chisquare: " << fTChisquare << endl;
429 *fLog << "Ndf: " << fTNdf << endl;
430 *fLog << "Probability: " << fTProb << endl;
431 *fLog << endl;
432
433}
Note: See TracBrowser for help on using the repository browser.