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

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