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

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