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

Last change on this file since 2622 was 2603, checked in by gaug, 22 years ago
*** empty log message ***
File size: 12.7 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//
60// Default Constructor.
61//
62MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
63 : fPixId(-1),
64 fChargeGausFit(NULL),
65 fTimeGausFit(NULL),
66 fFitLegend(NULL),
67 fLowerFitRange(0.),
68 fFitOK(kFALSE),
69 fChargeChisquare(-1.),
70 fChargeProb(-1.),
71 fChargeNdf(-1),
72 fTimeChisquare(-1.),
73 fTimeProb(-1.),
74 fTimeNdf(-1)
75{
76
77 fName = name ? name : "MHCalibrationPixel";
78 fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits";
79
80 // Create a large number of bins, later we will rebin
81 fChargeFirst = -0.5;
82 fChargeLast = 10000. - 0.5;
83 fChargeNbins = 20000;
84
85 fHCharge = new TH1I("HCharge","Distribution of Summed FADC Slices Pixel ",
86 fChargeNbins,fChargeFirst,fChargeLast);
87 fHCharge->SetXTitle("Sum FADC Slices");
88 fHCharge->SetYTitle("Nr. of events");
89 fHCharge->Sumw2();
90
91 fHCharge->SetDirectory(NULL);
92
93 Axis_t tfirst = -0.5;
94 Axis_t tlast = 15.5;
95 Int_t ntbins = 16;
96
97 fHTime = new TH1I("HTime","Distribution of Mean Arrival Times Pixel ",
98 ntbins,tfirst,tlast);
99 fHTime->SetXTitle("Mean Arrival Times [FADC slice nr]");
100 fHTime->SetYTitle("Nr. of events");
101 fHTime->Sumw2();
102
103 fHTime->SetDirectory(NULL);
104
105 TString qvsntitle = "Sum of Charges vs. Event Number Pixel ";
106
107 // We define a reasonable number and later enlarge it if necessary
108 Int_t nqbins = 20000;
109 Axis_t nfirst = -0.5;
110 Axis_t nlast = (Axis_t)nqbins - 0.5;
111
112 fHChargevsN = new TH1I("HChargevsN",qvsntitle.Data(),
113 nqbins,nfirst,nlast);
114 fHChargevsN->SetXTitle("Event Nr.");
115 fHChargevsN->SetYTitle("Sum of FADC slices");
116
117 fHChargevsN->SetDirectory(NULL);
118
119}
120
121MHCalibrationPixel::~MHCalibrationPixel()
122{
123
124 delete fHCharge;
125 delete fHTime;
126 delete fHChargevsN;
127
128 if (fChargeGausFit)
129 delete fChargeGausFit;
130 if (fTimeGausFit)
131 delete fTimeGausFit;
132 if (fFitLegend)
133 delete fFitLegend;
134
135}
136
137
138void MHCalibrationPixel::ChangeHistId(Int_t id)
139{
140
141 fPixId = id;
142
143 TString nameQ = TString(fHCharge->GetName());
144 nameQ += id;
145 fHCharge->SetName(nameQ.Data());
146
147 TString nameT = TString(fHTime->GetName());
148 nameT += id;
149 fHTime->SetName(nameT.Data());
150
151 TString nameQvsN = TString(fHChargevsN->GetName());
152 nameQvsN += id;
153 fHChargevsN->SetName(nameQvsN.Data());
154
155 TString titleQ = TString(fHCharge->GetTitle());
156 titleQ += id;
157 fHCharge->SetTitle(titleQ.Data());
158
159 TString titleT = TString(fHTime->GetTitle());
160 titleT += id;
161 fHTime->SetTitle(titleT.Data());
162
163 TString titleQvsN = TString(fHChargevsN->GetTitle());
164 titleQvsN += id;
165 fHChargevsN->SetTitle(titleQvsN.Data());
166}
167
168
169void MHCalibrationPixel::Reset()
170{
171
172 for (Int_t i = fHCharge->FindBin(fChargeFirst);
173 i <= fHCharge->FindBin(fChargeLast); i++)
174 fHCharge->SetBinContent(i, 1.e-20);
175
176 for (Int_t i = 0; i < 16; i++)
177 fHTime->SetBinContent(i, 1.e-20);
178
179 fChargeLast = 9999.5;
180
181 fHCharge->GetXaxis()->SetRangeUser(0.,fChargeLast);
182
183 return;
184}
185
186
187// -------------------------------------------------------------------------
188//
189// Set the binnings and prepare the filling of the histograms
190//
191Bool_t MHCalibrationPixel::SetupFill(const MParList *plist)
192{
193
194 fHCharge->Reset();
195 fHTime->Reset();
196
197 return kTRUE;
198}
199
200
201
202// -------------------------------------------------------------------------
203//
204// Draw a legend with the fit results
205//
206void MHCalibrationPixel::DrawLegend()
207{
208
209 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
210
211 if (fFitOK)
212 fFitLegend->SetFillColor(80);
213 else
214 fFitLegend->SetFillColor(2);
215
216 fFitLegend->SetLabel("Results of the Gauss Fit:");
217 fFitLegend->SetTextSize(0.05);
218
219 const TString line1 =
220 Form("Mean: Q_{#mu} = %2.2f #pm %2.2f",fChargeMean,fChargeMeanErr);
221
222 fFitLegend->AddText(line1);
223
224
225 const TString line4 =
226 Form("Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fChargeSigma,fChargeSigmaErr);
227
228 fFitLegend->AddText(line4);
229
230
231 const TString line7 =
232 Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChargeChisquare,fChargeNdf);
233
234 fFitLegend->AddText(line7);
235
236
237 const TString line8 =
238 Form("Probability: %4.3f ",fChargeProb);
239
240 fFitLegend->AddText(line8);
241
242
243 if (fFitOK)
244 fFitLegend->AddText("Result of the Fit: OK");
245 else
246 fFitLegend->AddText("Result of the Fit: NOT OK");
247
248 fFitLegend->SetBit(kCanDelete);
249 fFitLegend->Draw();
250
251}
252
253
254// -------------------------------------------------------------------------
255//
256// Draw the histogram
257//
258void MHCalibrationPixel::Draw(Option_t *opt)
259{
260
261 gStyle->SetOptFit(0);
262 gStyle->SetOptStat(1111111);
263
264 TCanvas *c = MakeDefCanvas(this,600,900);
265
266 gROOT->SetSelectedPad(NULL);
267
268 c->Divide(2,2);
269
270 c->cd(1);
271 gPad->SetBorderMode(0);
272 gPad->SetLogy(1);
273 gPad->SetTicks();
274
275 fHCharge->Draw(opt);
276
277 if (fChargeGausFit)
278 {
279 if (fFitOK)
280 fChargeGausFit->SetLineColor(kGreen);
281 else
282 fChargeGausFit->SetLineColor(kRed);
283
284 fChargeGausFit->Draw("same");
285 c->Modified();
286 c->Update();
287 }
288
289
290 c->cd(2);
291 DrawLegend();
292 c->Update();
293
294 c->cd(3);
295 gStyle->SetOptStat(1111111);
296
297 gPad->SetLogy(1);
298 fHTime->Draw(opt);
299
300 if (fTimeGausFit)
301 {
302 if (fTimeChisquare > 1.)
303 fTimeGausFit->SetLineColor(kRed);
304 else
305 fTimeGausFit->SetLineColor(kGreen);
306
307 fTimeGausFit->Draw("same");
308 c->Modified();
309 c->Update();
310 }
311
312 c->Modified();
313 c->Update();
314
315 c->cd(4);
316 fHChargevsN->Draw(opt);
317}
318
319
320
321Bool_t MHCalibrationPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *option)
322{
323
324 if (fTimeGausFit)
325 return kFALSE;
326
327 rmin = (rmin != 0.) ? rmin : 4.;
328 rmax = (rmax != 0.) ? rmax : 9.;
329
330 const Stat_t entries = fHTime->GetEntries();
331 const Double_t mu_guess = fHTime->GetBinCenter(fHTime->GetMaximumBin());
332 const Double_t sigma_guess = (rmax - rmin)/2.;
333 const Double_t area_guess = entries/gkSq2Pi;
334
335 TString name = TString("GausTime");
336 name += fPixId;
337 fTimeGausFit = new TF1(name.Data(),"gaus",rmin,rmax);
338
339 if (!fTimeGausFit)
340 {
341 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
342 return kFALSE;
343 }
344
345 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
346 fTimeGausFit->SetParNames("Area","#mu","#sigma");
347 fTimeGausFit->SetParLimits(0,0.,entries);
348 fTimeGausFit->SetParLimits(1,rmin,rmax);
349 fTimeGausFit->SetParLimits(2,0.,(rmax-rmin));
350 fTimeGausFit->SetRange(rmin,rmax);
351
352 fHTime->Fit(fTimeGausFit,option);
353
354 rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2);
355 rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2);
356 fTimeGausFit->SetRange(rmin,rmax);
357
358 fHTime->Fit(fTimeGausFit,option);
359
360 fTimeChisquare = fTimeGausFit->GetChisquare();
361 fTimeNdf = fTimeGausFit->GetNDF();
362 fTimeProb = fTimeGausFit->GetProb();
363
364 fTimeMean = fTimeGausFit->GetParameter(1);
365 fTimeSigma = fTimeGausFit->GetParameter(2);
366
367 if (fTimeChisquare > 1.)
368 {
369 *fLog << warn << "Fit of the Arrival times failed ! " << endl;
370 return kFALSE;
371 }
372
373 return kTRUE;
374
375}
376
377Bool_t MHCalibrationPixel::FitCharge(Option_t *option)
378{
379
380 if (fChargeGausFit)
381 return kFALSE;
382
383 //
384 // Get the fitting ranges
385 //
386 Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirst;
387 Axis_t rmax = 0.;
388
389 //
390 // First guesses for the fit (should be as close to reality as possible,
391 // otherwise the fit goes gaga because of high number of dimensions ...
392 //
393 const Stat_t entries = fHCharge->GetEntries();
394 const Double_t area_guess = entries/gkSq2Pi;
395 const Double_t mu_guess = fHCharge->GetBinCenter(fHCharge->GetMaximumBin());
396 const Double_t sigma_guess = mu_guess/15.;
397
398 TString name = TString("ChargeGausFit");
399 name += fPixId;
400
401 fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLast);
402
403 if (!fChargeGausFit)
404 {
405 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
406 return kFALSE;
407 }
408
409 fChargeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
410 fChargeGausFit->SetParNames("Area","#mu","#sigma");
411 fChargeGausFit->SetParLimits(0,0.,entries);
412 fChargeGausFit->SetParLimits(1,rmin,fChargeLast);
413 fChargeGausFit->SetParLimits(2,0.,fChargeLast-rmin);
414 fChargeGausFit->SetRange(rmin,fChargeLast);
415
416 fHCharge->Fit(fChargeGausFit,option);
417
418 Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
419
420 rmin = (rtry < rmin ? rmin : rtry);
421 rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);
422 fChargeGausFit->SetRange(rmin,rmax);
423
424 fHCharge->Fit(fChargeGausFit,option);
425
426 // rmin = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
427 // rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);
428 // fChargeGausFit->SetRange(rmin,rmax);
429
430 // fHCharge->Fit(fChargeGausFit,option);
431
432 fChargeChisquare = fChargeGausFit->GetChisquare();
433 fChargeNdf = fChargeGausFit->GetNDF();
434 fChargeProb = fChargeGausFit->GetProb();
435 fChargeMean = fChargeGausFit->GetParameter(1);
436 fChargeMeanErr = fChargeGausFit->GetParError(1);
437 fChargeSigma = fChargeGausFit->GetParameter(2);
438 fChargeSigmaErr = fChargeGausFit->GetParError(2);
439
440 //
441 // The fit result is accepted under condition
442 // The Probability is greater than gkProbLimit (default 0.01 == 99%)
443 //
444 if (fChargeProb < gkProbLimit)
445 {
446 *fLog << warn << "Prob: " << fChargeProb << " is smaller than the allowed value: " << gkProbLimit << endl;
447 fFitOK = kFALSE;
448 return kFALSE;
449 }
450
451
452 fFitOK = kTRUE;
453
454 return kTRUE;
455}
456
457
458void MHCalibrationPixel::CutAllEdges()
459{
460
461 Int_t nbins = 50;
462
463 CutEdges(fHCharge,nbins);
464
465 fChargeFirst = fHCharge->GetBinLowEdge(fHCharge->GetXaxis()->GetFirst());
466 fChargeLast = fHCharge->GetBinLowEdge(fHCharge->GetXaxis()->GetLast())+fHCharge->GetBinWidth(0);
467 fChargeNbins = nbins;
468
469 CutEdges(fHChargevsN,0);
470
471}
472
473void MHCalibrationPixel::PrintChargeFitResult()
474{
475
476 *fLog << "Results of the Summed Charges Fit: " << endl;
477 *fLog << "Chisquare: " << fChargeChisquare << endl;
478 *fLog << "DoF: " << fChargeNdf << endl;
479 *fLog << "Probability: " << fChargeProb << endl;
480 *fLog << endl;
481
482}
483
484void MHCalibrationPixel::PrintTimeFitResult()
485{
486
487 *fLog << "Results of the Time Slices Fit: " << endl;
488 *fLog << "Chisquare: " << fTimeChisquare << endl;
489 *fLog << "Ndf: " << fTimeNdf << endl;
490 *fLog << "Probability: " << fTimeProb << endl;
491 *fLog << endl;
492
493}
Note: See TracBrowser for help on using the repository browser.