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

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