source: trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc@ 2798

Last change on this file since 2798 was 2791, checked in by gaug, 22 years ago
*** empty log message ***
File size: 23.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#include <TGraph.h>
43#include <TAxis.h>
44
45#include <TF1.h>
46#include <TH2.h>
47#include <TProfile.h>
48#include <TCanvas.h>
49#include <TPad.h>
50#include <TPaveText.h>
51
52#include "MParList.h"
53
54#include "MLog.h"
55#include "MLogManip.h"
56
57ClassImp(MHCalibrationPixel);
58
59using namespace std;
60
61// --------------------------------------------------------------------------
62//
63// Default Constructor.
64//
65MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
66 : fPixId(-1),
67 fTotalEntries(0),
68 fHiGains(NULL),
69 fLoGains(NULL),
70 fChargeGausFit(NULL),
71 fTimeGausFit(NULL),
72 fHivsLoGain(NULL),
73 fFitLegend(NULL),
74 fLowerFitRange(-2000.),
75 fChargeFirstHiGain(-2000.5),
76 fChargeLastHiGain(9999.5),
77 fChargeNbinsHiGain(12000),
78 fChargeFirstLoGain(-2000.5),
79 fChargeLastLoGain(9999.5),
80 fChargeNbinsLoGain(1200),
81 fChargeChisquare(-1.),
82 fChargeProb(-1.),
83 fChargeNdf(-1),
84 fTimeChisquare(-1.),
85 fTimeProb(-1.),
86 fTimeNdf(-1),
87 fTimeMean(-1.),
88 fTimeSigma(-1.),
89 fTimeLowerFitRangeHiGain(0),
90 fTimeUpperFitRangeHiGain(0),
91 fTimeLowerFitRangeLoGain(0),
92 fTimeUpperFitRangeLoGain(0),
93 fUseLoGain(kFALSE),
94 fFitOK(kFALSE),
95 fOffset(0.),
96 fSlope(0.)
97{
98
99 fName = name ? name : "MHCalibrationPixel";
100 fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits";
101
102 // Create a large number of bins, later we will rebin
103 fHChargeHiGain = new TH1F("HChargeHiGain","Distribution of Summed FADC Hi Gain Slices Pixel ",
104 fChargeNbinsHiGain,fChargeFirstHiGain,fChargeLastHiGain);
105 fHChargeHiGain->SetXTitle("Sum FADC Slices (Hi Gain)");
106 fHChargeHiGain->SetYTitle("Nr. of events");
107 fHChargeHiGain->Sumw2();
108
109 fHChargeHiGain->SetDirectory(NULL);
110
111 fHChargeLoGain = new TH1F("HChargeLoGain","Distribution of Summed FADC Lo Gain Slices Pixel ",
112 fChargeNbinsLoGain,fChargeFirstLoGain,fChargeLastLoGain);
113 fHChargeLoGain->SetXTitle("Sum FADC Slices (Lo Gain)");
114 fHChargeLoGain->SetYTitle("Nr. of events");
115 fHChargeLoGain->Sumw2();
116
117 fHChargeLoGain->SetDirectory(NULL);
118
119 Axis_t tfirst = -0.5;
120 Axis_t tlast = 15.5;
121 Int_t ntbins = 16;
122
123 fHTimeHiGain = new TH1I("HTimeHiGain","Distribution of Mean Arrival Hi Gain Times Pixel ",
124 ntbins,tfirst,tlast);
125 fHTimeHiGain->SetXTitle("Mean Arrival Times [Hi Gain FADC slice nr]");
126 fHTimeHiGain->SetYTitle("Nr. of events");
127 fHTimeHiGain->SetDirectory(NULL);
128
129 fHTimeLoGain = new TH1I("HTimeLoGain","Distribution of Mean Arrival Lo Gain Times Pixel ",
130 ntbins,tfirst,tlast);
131 fHTimeLoGain->SetXTitle("Mean Arrival Times [Lo Gain FADC slice nr]");
132 fHTimeLoGain->SetYTitle("Nr. of events");
133 fHTimeLoGain->SetDirectory(NULL);
134
135 // We define a reasonable number and later enlarge it if necessary
136 Int_t nqbins = 20000;
137 Axis_t nfirst = -0.5;
138 Axis_t nlast = (Axis_t)nqbins - 0.5;
139
140 fHChargevsNHiGain = new TH1I("HChargevsNHiGain","Sum of Hi Gain Charges vs. Event Number Pixel ",
141 nqbins,nfirst,nlast);
142 fHChargevsNHiGain->SetXTitle("Event Nr.");
143 fHChargevsNHiGain->SetYTitle("Sum of Hi Gain FADC slices");
144
145 fHChargevsNHiGain->SetDirectory(NULL);
146
147 fHChargevsNLoGain = new TH1I("HChargevsNLoGain","Sum of Lo Gain Charges vs. Event Number Pixel ",
148 nqbins,nfirst,nlast);
149 fHChargevsNLoGain->SetXTitle("Event Nr.");
150 fHChargevsNLoGain->SetYTitle("Sum of Lo Gain FADC slices");
151
152 fHChargevsNLoGain->SetDirectory(NULL);
153
154 fHiGains = new TArrayF();
155 fLoGains = new TArrayF();
156
157}
158
159MHCalibrationPixel::~MHCalibrationPixel()
160{
161
162 delete fHChargeHiGain;
163 delete fHTimeHiGain;
164 delete fHChargevsNHiGain;
165
166 delete fHChargeLoGain;
167 delete fHTimeLoGain;
168 delete fHChargevsNLoGain;
169
170 delete fHiGains;
171 delete fLoGains;
172
173 if (fChargeGausFit)
174 delete fChargeGausFit;
175 if (fTimeGausFit)
176 delete fTimeGausFit;
177 if (fFitLegend)
178 delete fFitLegend;
179 if (fHivsLoGain)
180 delete fHivsLoGain;
181
182
183}
184
185
186void MHCalibrationPixel::ChangeHistId(Int_t id)
187{
188
189 fPixId = id;
190
191 //
192 // Names Hi gain Histograms
193 //
194 TString nameQHiGain = TString(fHChargeHiGain->GetName());
195 nameQHiGain += id;
196 fHChargeHiGain->SetName(nameQHiGain.Data());
197
198 TString nameTHiGain = TString(fHTimeHiGain->GetName());
199 nameTHiGain += id;
200 fHTimeHiGain->SetName(nameTHiGain.Data());
201
202 TString nameQvsNHiGain = TString(fHChargevsNHiGain->GetName());
203 nameQvsNHiGain += id;
204 fHChargevsNHiGain->SetName(nameQvsNHiGain.Data());
205
206 //
207 // Title Hi gain Histograms
208 //
209 TString titleQHiGain = TString(fHChargeHiGain->GetTitle());
210 titleQHiGain += id;
211 fHChargeHiGain->SetTitle(titleQHiGain.Data());
212
213 TString titleTHiGain = TString(fHTimeHiGain->GetTitle());
214 titleTHiGain += id;
215 fHTimeHiGain->SetTitle(titleTHiGain.Data());
216
217 TString titleQvsNHiGain = TString(fHChargevsNHiGain->GetTitle());
218 titleQvsNHiGain += id;
219 fHChargevsNHiGain->SetTitle(titleQvsNHiGain.Data());
220
221 //
222 // Names Low Gain Histograms
223 //
224 TString nameQLoGain = TString(fHChargeLoGain->GetName());
225 nameQLoGain += id;
226 fHChargeLoGain->SetName(nameQLoGain.Data());
227
228 TString nameTLoGain = TString(fHTimeLoGain->GetName());
229 nameTLoGain += id;
230 fHTimeLoGain->SetName(nameTLoGain.Data());
231
232 TString nameQvsNLoGain = TString(fHChargevsNLoGain->GetName());
233 nameQvsNLoGain += id;
234 fHChargevsNLoGain->SetName(nameQvsNLoGain.Data());
235
236 //
237 // Titles Low Gain Histograms
238 //
239 TString titleQLoGain = TString(fHChargeLoGain->GetTitle());
240 titleQLoGain += id;
241 fHChargeLoGain->SetTitle(titleQLoGain.Data());
242
243 TString titleTLoGain = TString(fHTimeLoGain->GetTitle());
244 titleTLoGain += id;
245 fHTimeLoGain->SetTitle(titleTLoGain.Data());
246
247 TString titleQvsNLoGain = TString(fHChargevsNLoGain->GetTitle());
248 titleQvsNLoGain += id;
249 fHChargevsNLoGain->SetTitle(titleQvsNLoGain.Data());
250}
251
252
253void MHCalibrationPixel::Reset()
254{
255
256 for (Int_t i = fHChargeHiGain->FindBin(fChargeFirstHiGain);
257 i <= fHChargeHiGain->FindBin(fChargeLastHiGain); i++)
258 fHChargeHiGain->SetBinContent(i, 1.e-20);
259
260 for (Int_t i = 0; i < 16; i++)
261 fHTimeHiGain->SetBinContent(i, 1.e-20);
262
263 fChargeLastHiGain = 9999.5;
264
265 fHChargeHiGain->GetXaxis()->SetRangeUser(0.,fChargeLastHiGain);
266
267 for (Int_t i = fHChargeLoGain->FindBin(fChargeFirstLoGain);
268 i <= fHChargeLoGain->FindBin(fChargeLastLoGain); i++)
269 fHChargeLoGain->SetBinContent(i, 1.e-20);
270
271 for (Int_t i = 0; i < 16; i++)
272 fHTimeLoGain->SetBinContent(i, 1.e-20);
273
274 fChargeLastLoGain = 9999.5;
275
276 fHChargeLoGain->GetXaxis()->SetRangeUser(0.,fChargeLastLoGain);
277
278 return;
279}
280
281
282// -------------------------------------------------------------------------
283//
284// Set the binnings and prepare the filling of the histograms
285//
286Bool_t MHCalibrationPixel::SetupFill(const MParList *plist)
287{
288
289 fHChargeHiGain->Reset();
290 fHTimeHiGain->Reset();
291 fHChargeLoGain->Reset();
292 fHTimeLoGain->Reset();
293
294 return kTRUE;
295}
296
297
298Bool_t MHCalibrationPixel::UseLoGain()
299{
300
301 if (fHChargeHiGain->Integral() > fHChargeLoGain->Integral())
302 {
303 fUseLoGain = kFALSE;
304 return kFALSE;
305 }
306 else
307 {
308 fUseLoGain = kTRUE;
309 return kTRUE;
310 }
311}
312
313void MHCalibrationPixel::SetPointInGraph(Float_t qhi,Float_t qlo)
314{
315
316 fHiGains->Set(++fTotalEntries);
317 fLoGains->Set(fTotalEntries);
318
319 fHiGains->AddAt(qhi,fTotalEntries-1);
320 fLoGains->AddAt(qlo,fTotalEntries-1);
321
322}
323
324
325// -------------------------------------------------------------------------
326//
327// Draw a legend with the fit results
328//
329void MHCalibrationPixel::DrawLegend()
330{
331
332 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
333
334 if (fFitOK)
335 fFitLegend->SetFillColor(80);
336 else
337 fFitLegend->SetFillColor(2);
338
339 fFitLegend->SetLabel("Results of the Gauss Fit:");
340 fFitLegend->SetTextSize(0.05);
341
342 const TString line1 =
343 Form("Mean: Q_{#mu} = %2.2f #pm %2.2f",fChargeMean,fChargeMeanErr);
344
345 fFitLegend->AddText(line1);
346
347
348 const TString line4 =
349 Form("Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fChargeSigma,fChargeSigmaErr);
350
351 fFitLegend->AddText(line4);
352
353
354 const TString line7 =
355 Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChargeChisquare,fChargeNdf);
356
357 fFitLegend->AddText(line7);
358
359
360 const TString line8 =
361 Form("Probability: %4.3f ",fChargeProb);
362
363 fFitLegend->AddText(line8);
364
365 if (fFitOK)
366 fFitLegend->AddText("Result of the Fit: OK");
367 else
368 fFitLegend->AddText("Result of the Fit: NOT OK");
369
370 fFitLegend->SetBit(kCanDelete);
371 fFitLegend->Draw();
372
373}
374
375
376// -------------------------------------------------------------------------
377//
378// Draw the histogram
379//
380void MHCalibrationPixel::Draw(Option_t *opt)
381{
382
383 if (!fHivsLoGain)
384 FitHiGainvsLoGain();
385
386 gStyle->SetOptFit(0);
387 gStyle->SetOptStat(1111111);
388
389 TCanvas *c = MakeDefCanvas(this,600,900);
390
391 gROOT->SetSelectedPad(NULL);
392
393 c->Divide(2,4);
394
395 c->cd(1);
396 gPad->SetBorderMode(0);
397 gPad->SetTicks();
398
399 if (fHChargeHiGain->Integral() > 0)
400 gPad->SetLogy(1);
401 else
402 gPad->SetLogy(0);
403
404 fHChargeHiGain->DrawCopy(opt);
405
406 c->Modified();
407 c->Update();
408
409 if (fUseLoGain)
410 {
411
412 c->cd(2);
413 gPad->SetTicks();
414
415 if (fHChargeLoGain->Integral() > 0)
416 gPad->SetLogy(1);
417 else
418 gPad->SetLogy(0);
419
420 fHChargeLoGain->DrawCopy(opt);
421
422 if (fChargeGausFit)
423 {
424 if (fFitOK)
425 fChargeGausFit->SetLineColor(kGreen);
426 else
427 fChargeGausFit->SetLineColor(kRed);
428
429 fChargeGausFit->DrawCopy("same");
430 }
431 c->Modified();
432 c->Update();
433
434 c->cd(3);
435 gROOT->SetSelectedPad(NULL);
436 gStyle->SetOptFit();
437 fHivsLoGain->Draw("prof");
438
439 gPad->Modified();
440 gPad->Update();
441
442 c->cd(4);
443 DrawLegend();
444
445 }
446 else
447 {
448 if (fChargeGausFit)
449 {
450 if (fFitOK)
451 fChargeGausFit->SetLineColor(kGreen);
452 else
453 fChargeGausFit->SetLineColor(kRed);
454
455 fChargeGausFit->DrawCopy("same");
456 }
457
458 c->cd(2);
459 gPad->SetTicks();
460
461 if (fHChargeLoGain->Integral() > 0)
462 gPad->SetLogy(1);
463 else
464 gPad->SetLogy(0);
465
466 fHChargeLoGain->DrawCopy(opt);
467 c->Modified();
468 c->Update();
469
470 c->cd(3);
471 DrawLegend();
472
473 c->cd(4);
474
475 gROOT->SetSelectedPad(NULL);
476 gStyle->SetOptFit();
477 fHivsLoGain->Draw("prof");
478 gPad->Modified();
479 gPad->Update();
480 }
481
482 c->Modified();
483 c->Update();
484
485 c->cd(5);
486 gStyle->SetOptStat(1111111);
487
488 gPad->SetTicks();
489 gPad->SetLogy(0);
490 fHTimeHiGain->DrawCopy(opt);
491 c->Modified();
492 c->Update();
493
494 if (fUseLoGain)
495 {
496
497 c->cd(6);
498 gPad->SetTicks();
499 gPad->SetLogy(0);
500 fHTimeLoGain->DrawCopy(opt);
501 c->Modified();
502 c->Update();
503
504 if (fTimeGausFit)
505 {
506 if (fTimeChisquare > 20.)
507 fTimeGausFit->SetLineColor(kRed);
508 else
509 fTimeGausFit->SetLineColor(kGreen);
510
511 fTimeGausFit->DrawCopy("same");
512 c->Modified();
513 c->Update();
514 }
515 }
516 else
517 {
518 if (fTimeGausFit)
519 {
520 if (fTimeChisquare > 20.)
521 fTimeGausFit->SetLineColor(kRed);
522 else
523 fTimeGausFit->SetLineColor(kGreen);
524
525 fTimeGausFit->DrawCopy("same");
526 c->Modified();
527 c->Update();
528 }
529
530 c->cd(6);
531 gPad->SetTicks();
532 gPad->SetLogy(0);
533 fHTimeLoGain->DrawCopy(opt);
534 c->Modified();
535 c->Update();
536
537 }
538 c->Modified();
539 c->Update();
540
541 c->cd(7);
542 gPad->SetTicks();
543 fHChargevsNHiGain->DrawCopy(opt);
544 c->Modified();
545 c->Update();
546
547 c->cd(8);
548 gPad->SetTicks();
549 fHChargevsNLoGain->DrawCopy(opt);
550 c->Modified();
551 c->Update();
552
553
554}
555
556void MHCalibrationPixel::FitHiGainvsLoGain()
557{
558
559 if (fHivsLoGain)
560 return;
561
562 gStyle->SetOptFit();
563
564 fHivsLoGain = new TProfile("HiGainvsLoGain","Plot the High Gain vs. Low Gain",100,0.,1000.,0.,1000.);
565 fHivsLoGain->GetXaxis()->SetTitle("Sum of Charges High Gain");
566 fHivsLoGain->GetYaxis()->SetTitle("Sum of Charges Low Gain");
567
568 TString titleHiLoGain = TString(fHivsLoGain->GetName());
569 titleHiLoGain += fPixId;
570 fHivsLoGain->SetName(titleHiLoGain.Data());
571
572
573 for (Int_t i=0;i<fTotalEntries;i++)
574 fHivsLoGain->Fill(fHiGains->At(i),fLoGains->At(i),1);
575
576 fHivsLoGain->Fit("pol1","rq","",fHivsLoGain->GetMean()-2.5*fHivsLoGain->GetRMS(),fHivsLoGain->GetMean()+2.*fHivsLoGain->GetRMS());
577
578 fOffset = fHivsLoGain->GetFunction("pol1")->GetParameter(0);
579 fSlope = fHivsLoGain->GetFunction("pol1")->GetParameter(1);
580
581}
582
583
584Bool_t MHCalibrationPixel::FitTimeHiGain(Axis_t rmin, Axis_t rmax, Option_t *option)
585{
586
587 if (fTimeGausFit)
588 return kFALSE;
589
590 rmin = (rmin != 0.) ? rmin : (Axis_t)fTimeLowerFitRangeHiGain;
591 rmax = (rmax != 0.) ? rmax : (Axis_t)fTimeUpperFitRangeHiGain;
592
593 const Stat_t entries = fHTimeHiGain->Integral();
594 const Double_t mu_guess = fHTimeHiGain->GetBinCenter(fHTimeHiGain->GetMaximumBin());
595 const Double_t sigma_guess = (rmax - rmin)/2.;
596 const Double_t area_guess = entries/gkSq2Pi;
597
598 TString name = TString("GausTime");
599 name += fPixId;
600 fTimeGausFit = new TF1(name.Data(),"gaus",rmin,rmax);
601
602 if (!fTimeGausFit)
603 {
604 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
605 return kFALSE;
606 }
607
608 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
609 fTimeGausFit->SetParNames("Area","#mu","#sigma");
610 fTimeGausFit->SetParLimits(0,0.,entries);
611 fTimeGausFit->SetParLimits(1,rmin,rmax);
612 fTimeGausFit->SetParLimits(2,0.,(rmax-rmin));
613 fTimeGausFit->SetRange(rmin,rmax);
614
615 fHTimeHiGain->Fit(fTimeGausFit,option);
616
617 rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2);
618 rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2);
619 fTimeGausFit->SetRange(rmin,rmax);
620
621 fHTimeHiGain->Fit(fTimeGausFit,option);
622
623 fTimeChisquare = fTimeGausFit->GetChisquare();
624 fTimeNdf = fTimeGausFit->GetNDF();
625 fTimeProb = fTimeGausFit->GetProb();
626
627 fTimeMean = fTimeGausFit->GetParameter(1);
628 fTimeSigma = fTimeGausFit->GetParameter(2);
629
630 if (fTimeChisquare > 30.) // Cannot use Probability because Ndf is sometimes < 1
631 {
632 *fLog << warn << "Fit of the Arrival times failed ! " << endl;
633 return kFALSE;
634 }
635
636 return kTRUE;
637
638}
639
640Bool_t MHCalibrationPixel::FitTimeLoGain(Axis_t rmin, Axis_t rmax, Option_t *option)
641{
642
643 if (fTimeGausFit)
644 return kFALSE;
645
646 rmin = (rmin != 0.) ? rmin : (Axis_t)fTimeLowerFitRangeLoGain;
647 rmax = (rmax != 0.) ? rmax : (Axis_t)fTimeUpperFitRangeLoGain;
648
649 const Stat_t entries = fHTimeLoGain->Integral();
650 const Double_t mu_guess = fHTimeLoGain->GetBinCenter(fHTimeLoGain->GetMaximumBin());
651 const Double_t sigma_guess = (rmax - rmin)/2.;
652 const Double_t area_guess = entries/gkSq2Pi;
653
654 TString name = TString("GausTime");
655 name += fPixId;
656 fTimeGausFit = new TF1(name.Data(),"gaus",rmin,rmax);
657
658 if (!fTimeGausFit)
659 {
660 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
661 return kFALSE;
662 }
663
664 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
665 fTimeGausFit->SetParNames("Area","#mu","#sigma");
666 fTimeGausFit->SetParLimits(0,0.,entries);
667 fTimeGausFit->SetParLimits(1,rmin,rmax);
668 fTimeGausFit->SetParLimits(2,0.,(rmax-rmin));
669 fTimeGausFit->SetRange(rmin,rmax);
670
671 fHTimeLoGain->Fit(fTimeGausFit,option);
672
673 rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2);
674 rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2);
675 fTimeGausFit->SetRange(rmin,rmax);
676
677 fHTimeLoGain->Fit(fTimeGausFit,option);
678
679 fTimeChisquare = fTimeGausFit->GetChisquare();
680 fTimeNdf = fTimeGausFit->GetNDF();
681 fTimeProb = fTimeGausFit->GetProb();
682
683 fTimeMean = fTimeGausFit->GetParameter(1);
684 fTimeSigma = fTimeGausFit->GetParameter(2);
685
686 if (fTimeChisquare > 30.) // Cannot use Probability because Ndf is sometimes < 1
687 {
688 *fLog << warn << "Fit of the Arrival times failed ! " << endl;
689 return kFALSE;
690 }
691
692 return kTRUE;
693
694}
695
696Bool_t MHCalibrationPixel::FitChargeHiGain(Option_t *option)
697{
698
699 if (fChargeGausFit)
700 return kFALSE;
701
702 //
703 // Get the fitting ranges
704 //
705 Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirstHiGain;
706 Axis_t rmax = 0.;
707
708 //
709 // First guesses for the fit (should be as close to reality as possible,
710 // otherwise the fit goes gaga because of high number of dimensions ...
711 //
712 const Stat_t entries = fHChargeHiGain->Integral();
713 const Double_t area_guess = entries/gkSq2Pi;
714 const Double_t mu_guess = fHChargeHiGain->GetBinCenter(fHChargeHiGain->GetMaximumBin());
715 const Double_t sigma_guess = mu_guess/15.;
716
717 TString name = TString("ChargeGausFit");
718 name += fPixId;
719
720 fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLastHiGain);
721
722 if (!fChargeGausFit)
723 {
724 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
725 return kFALSE;
726 }
727
728 fChargeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
729 fChargeGausFit->SetParNames("Area","#mu","#sigma");
730 fChargeGausFit->SetParLimits(0,0.,entries);
731 fChargeGausFit->SetParLimits(1,rmin,fChargeLastHiGain);
732 fChargeGausFit->SetParLimits(2,0.,fChargeLastHiGain-rmin);
733 fChargeGausFit->SetRange(rmin,fChargeLastHiGain);
734
735 fHChargeHiGain->Fit(fChargeGausFit,option);
736
737 Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.0*fChargeGausFit->GetParameter(2);
738
739 rmin = (rtry < rmin ? rmin : rtry);
740 rmax = fChargeGausFit->GetParameter(1) + 3.5*fChargeGausFit->GetParameter(2);
741 fChargeGausFit->SetRange(rmin,rmax);
742
743 fHChargeHiGain->Fit(fChargeGausFit,option);
744
745 fChargeChisquare = fChargeGausFit->GetChisquare();
746 fChargeNdf = fChargeGausFit->GetNDF();
747 fChargeProb = fChargeGausFit->GetProb();
748 fChargeMean = fChargeGausFit->GetParameter(1);
749 fChargeMeanErr = fChargeGausFit->GetParError(1);
750 fChargeSigma = fChargeGausFit->GetParameter(2);
751 fChargeSigmaErr = fChargeGausFit->GetParError(2);
752
753 //
754 // The fit result is accepted under condition
755 // The Probability is greater than gkProbLimit (default 0.001 == 99.9%)
756 //
757 if (fChargeProb < gkProbLimit)
758 {
759 *fLog << warn << "Prob: " << fChargeProb
760 << " is smaller than the allowed value: " << gkProbLimit << endl;
761 fFitOK = kFALSE;
762 return kFALSE;
763 }
764
765 fFitOK = kTRUE;
766 return kTRUE;
767}
768
769
770Bool_t MHCalibrationPixel::FitChargeLoGain(Option_t *option)
771{
772
773 if (fChargeGausFit)
774 return kFALSE;
775
776 //
777 // Get the fitting ranges
778 //
779 Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirstLoGain;
780 Axis_t rmax = 0.;
781
782 //
783 // First guesses for the fit (should be as close to reality as possible,
784 // otherwise the fit goes gaga because of high number of dimensions ...
785 //
786 const Stat_t entries = fHChargeLoGain->Integral();
787 const Double_t area_guess = entries/gkSq2Pi;
788 const Double_t mu_guess = fHChargeLoGain->GetBinCenter(fHChargeLoGain->GetMaximumBin());
789 const Double_t sigma_guess = mu_guess/15.;
790
791 TString name = TString("ChargeGausFit");
792 name += fPixId;
793
794 fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLastLoGain);
795
796 if (!fChargeGausFit)
797 {
798 *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl;
799 return kFALSE;
800 }
801
802 fChargeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
803 fChargeGausFit->SetParNames("Area","#mu","#sigma");
804 fChargeGausFit->SetParLimits(0,0.,entries);
805 fChargeGausFit->SetParLimits(1,rmin,fChargeLastLoGain);
806 fChargeGausFit->SetParLimits(2,0.,fChargeLastLoGain-rmin);
807 fChargeGausFit->SetRange(rmin,fChargeLastLoGain);
808
809 fHChargeLoGain->Fit(fChargeGausFit,option);
810
811 Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.*fChargeGausFit->GetParameter(2);
812
813 rmin = (rtry < rmin ? rmin : rtry);
814 rmax = fChargeGausFit->GetParameter(1) + 3.5*fChargeGausFit->GetParameter(2);
815 fChargeGausFit->SetRange(rmin,rmax);
816
817 fHChargeLoGain->Fit(fChargeGausFit,option);
818
819 // rmin = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
820 // rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);
821 // fChargeGausFit->SetRange(rmin,rmax);
822
823 // fHChargeLoGain->Fit(fChargeGausFit,option);
824
825 fChargeChisquare = fChargeGausFit->GetChisquare();
826 fChargeNdf = fChargeGausFit->GetNDF();
827 fChargeProb = fChargeGausFit->GetProb();
828 fChargeMean = fChargeGausFit->GetParameter(1);
829 fChargeMeanErr = fChargeGausFit->GetParError(1);
830 fChargeSigma = fChargeGausFit->GetParameter(2);
831 fChargeSigmaErr = fChargeGausFit->GetParError(2);
832
833 //
834 // The fit result is accepted under condition
835 // The Probability is greater than gkProbLimit (default 0.01 == 99%)
836 //
837 if (fChargeProb < gkProbLimit)
838 {
839 *fLog << warn << "Prob: " << fChargeProb << " is smaller than the allowed value: " << gkProbLimit << endl;
840 fFitOK = kFALSE;
841 return kFALSE;
842 }
843
844
845 fFitOK = kTRUE;
846
847 return kTRUE;
848}
849
850
851void MHCalibrationPixel::CutAllEdges()
852{
853
854 Int_t nbins = 30;
855
856 CutEdges(fHChargeHiGain,nbins);
857
858 fChargeFirstHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetFirst());
859 fChargeLastHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetLast())
860 +fHChargeHiGain->GetBinWidth(0);
861 fChargeNbinsHiGain = nbins;
862
863 CutEdges(fHChargeLoGain,nbins);
864
865 fChargeFirstLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetFirst());
866 fChargeLastLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetLast())
867 +fHChargeLoGain->GetBinWidth(0);
868 fChargeNbinsLoGain = nbins;
869
870 CutEdges(fHChargevsNHiGain,0);
871 CutEdges(fHChargevsNLoGain,0);
872
873}
874
875void MHCalibrationPixel::PrintChargeFitResult()
876{
877
878 *fLog << "Results of the Summed Charges Fit: " << endl;
879 *fLog << "Chisquare: " << fChargeChisquare << endl;
880 *fLog << "DoF: " << fChargeNdf << endl;
881 *fLog << "Probability: " << fChargeProb << endl;
882 *fLog << endl;
883
884}
885
886void MHCalibrationPixel::PrintTimeFitResult()
887{
888
889 *fLog << "Results of the Time Slices Fit: " << endl;
890 *fLog << "Chisquare: " << fTimeChisquare << endl;
891 *fLog << "Ndf: " << fTimeNdf << endl;
892 *fLog << "Probability: " << fTimeProb << endl;
893 *fLog << endl;
894
895}
Note: See TracBrowser for help on using the repository browser.