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

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