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

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