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

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