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

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