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

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