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

Last change on this file since 3110 was 3108, checked in by gaug, 23 years ago
*** empty log message ***
File size: 22.8 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 Int_t MHCalibrationPixel::fPSDNbins = 30;
64const Int_t MHCalibrationPixel::fPulserFrequency = 500;
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 if (IsUseLoGain())
271 fHPSD = MH::ProjectArray(fPSDLoGain, fPSDNbins,
272 Form("%s%d","PSDProjection",fPixId),
273 Form("%s%d","Power Spectrum Density Projection LoGain ",fPixId));
274 else
275 fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins,
276 Form("%s%d","PSDProjection",fPixId),
277 Form("%s%d","Power Spectrum Density Projection HiGain ",fPixId));
278
279 //
280 // First guesses for the fit (should be as close to reality as possible,
281 //
282 const Double_t xmax = fHPSD->GetXaxis()->GetXmax();
283
284 fPSDExpFit = new TF1(Form("%s%d","PSDExpFit",fPixId),"exp([0]-[1]*x)",0.,xmax);
285
286 const Double_t slope_guess = (TMath::Log(fHPSD->GetEntries())+1.)/xmax;
287 const Double_t offset_guess = slope_guess*xmax;
288
289 fPSDExpFit->SetParameters(offset_guess, slope_guess);
290 fPSDExpFit->SetParNames("Offset","Slope");
291 fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess);
292 fPSDExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess);
293 fPSDExpFit->SetRange(0.,xmax);
294
295 fHPSD->Fit(fPSDExpFit,"RQL0");
296
297 fPSDProb = fPSDExpFit->GetProb();
298
299 if (fPSDProb < gkProbLimit)
300 {
301 SETBIT(fFlags,kOscillating);
302 return kTRUE;
303 }
304
305 CLRBIT(fFlags,kOscillating);
306
307 return kFALSE;
308}
309
310void MHCalibrationPixel::CreatePSDXaxis(Int_t n)
311{
312
313 if (fPSDXaxis)
314 return;
315
316 fPSDXaxis = new TArrayF(n);
317
318 for (Int_t i=0;i<n;i++)
319 fPSDXaxis->AddAt((Float_t)(fPulserFrequency*i)/2./n,i);
320}
321
322void MHCalibrationPixel::CreateChargeXaxis(Int_t n)
323{
324
325 if (!fChargeXaxis)
326 {
327 fChargeXaxis = new TArrayF(n);
328 for (Int_t i=0;i<n;i++)
329 fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,i);
330 return;
331 }
332
333 if (fChargeXaxis->GetSize() == n)
334 return;
335
336 const Int_t diff = fChargeXaxis->GetSize()-n;
337
338 fChargeXaxis->Set(n);
339 if (diff < 0)
340 for (Int_t i=n-1;i>n-diff-1;i--)
341 fChargeXaxis->AddAt((Float_t)i,i);
342}
343
344void MHCalibrationPixel::CutArrayBorder(TArrayF *array)
345{
346
347 Int_t i;
348
349 for (i=array->GetSize()-1;i>=0;i--)
350 if (array->At(i) != 0)
351 {
352 array->Set(i+1);
353 break;
354 }
355}
356
357
358
359Bool_t MHCalibrationPixel::IsEmpty() const
360{
361 return !(fHChargeHiGain->GetEntries() || fHChargeLoGain->GetEntries());
362}
363
364Bool_t MHCalibrationPixel::IsUseLoGain() const
365{
366 return TESTBIT(fFlags,kUseLoGain);
367}
368
369Bool_t MHCalibrationPixel::IsChargeFitOK() const
370{
371 return TESTBIT(fFlags,kChargeFitOK);
372}
373
374Bool_t MHCalibrationPixel::IsOscillating()
375{
376
377 if (fPSDExpFit)
378 return TESTBIT(fFlags,kOscillating);
379
380 return CheckOscillations();
381
382}
383
384Bool_t MHCalibrationPixel::FillChargeLoGain(Float_t q)
385{
386 return (fHChargeLoGain->Fill(q) > -1);
387}
388
389Bool_t MHCalibrationPixel::FillAbsTimeLoGain(Float_t t)
390{
391 return (fHAbsTimeLoGain->Fill(t)> -1);
392}
393
394Bool_t MHCalibrationPixel::FillChargeHiGain(Float_t q)
395{
396 return (fHChargeHiGain->Fill(q) > -1);
397}
398
399Bool_t MHCalibrationPixel::FillAbsTimeHiGain(Float_t t)
400{
401 return (fHAbsTimeHiGain->Fill(t) > -1);
402}
403
404void MHCalibrationPixel::ChangeHistId(Int_t id)
405{
406
407 // Change only if the names have not yet been changed
408 if (fPixId == -1)
409 {
410
411 //
412 // Names Hi gain Histograms
413 //
414 fHChargeHiGain->SetName( Form("%s%d",fHChargeHiGain->GetName(), id));
415 fHAbsTimeHiGain->SetName( Form("%s%d",fHAbsTimeHiGain->GetName(), id));
416
417 //
418 // Title Hi gain Histograms
419 //
420 fHChargeHiGain->SetTitle( Form("%s%d",fHChargeHiGain->GetTitle(), id));
421 fHAbsTimeHiGain->SetTitle( Form("%s%d",fHAbsTimeHiGain->GetTitle(), id));
422
423 //
424 // Names Low Gain Histograms
425 //
426 fHChargeLoGain->SetName( Form("%s%d",fHChargeLoGain->GetName(),id));
427 fHAbsTimeLoGain->SetName( Form("%s%d",fHAbsTimeLoGain->GetName(),id));
428
429 //
430 // Titles Low Gain Histograms
431 //
432 fHChargeLoGain->SetTitle( Form("%s%d",fHChargeLoGain->GetTitle(),id));
433 fHAbsTimeLoGain->SetTitle( Form("%s%d",fHAbsTimeLoGain->GetTitle(),id));
434
435 fPixId = id;
436 }
437
438}
439
440
441Bool_t MHCalibrationPixel::UseLoGain()
442{
443
444 if (fHChargeHiGain->Integral() > fHChargeLoGain->Integral())
445 {
446 CLRBIT(fFlags,kUseLoGain);
447 return kFALSE;
448 }
449 else
450 {
451 SETBIT(fFlags,kUseLoGain);
452 return kTRUE;
453 }
454}
455
456Bool_t MHCalibrationPixel::FillGraphs(Float_t qhi,Float_t qlo)
457{
458
459 if (fTotalEntries >= fCurrentSize)
460 {
461 fCurrentSize *= 2;
462
463 fHiGains->Set(fCurrentSize);
464 fLoGains->Set(fCurrentSize);
465 }
466
467 fHiGains->AddAt(qhi,fTotalEntries);
468 fLoGains->AddAt(qlo,fTotalEntries);
469
470 fTotalEntries++;
471
472 return kTRUE;
473
474}
475
476
477// -------------------------------------------------------------------------
478//
479// Draw a legend with the fit results
480//
481void MHCalibrationPixel::DrawLegend()
482{
483
484 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
485
486 if (IsChargeFitOK())
487 fFitLegend->SetFillColor(80);
488 else
489 fFitLegend->SetFillColor(2);
490
491 fFitLegend->SetLabel("Results of the Gauss Fit:");
492 fFitLegend->SetTextSize(0.05);
493
494 const TString line1 =
495 Form("Mean: Q_{#mu} = %2.2f #pm %2.2f",fChargeMean,fChargeMeanErr);
496 TText *t1 = fFitLegend->AddText(line1);
497 t1->SetBit(kCanDelete);
498
499 const TString line4 =
500 Form("Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fChargeSigma,fChargeSigmaErr);
501 TText *t2 = fFitLegend->AddText(line4);
502 t2->SetBit(kCanDelete);
503
504 const TString line7 =
505 Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChargeChisquare,fChargeNdf);
506 TText *t3 = fFitLegend->AddText(line7);
507 t3->SetBit(kCanDelete);
508
509 const TString line8 =
510 Form("Probability: %4.3f ",fChargeProb);
511 TText *t4 = fFitLegend->AddText(line8);
512 t4->SetBit(kCanDelete);
513
514 if (IsChargeFitOK())
515 {
516 TText *t5 = fFitLegend->AddText("Result of the Fit: OK");
517 t5->SetBit(kCanDelete);
518 }
519 else
520 {
521 TText *t6 = fFitLegend->AddText("Result of the Fit: NOT OK");
522 t6->SetBit(kCanDelete);
523 }
524
525 fFitLegend->SetBit(kCanDelete);
526 fFitLegend->Draw();
527
528}
529
530
531TObject *MHCalibrationPixel::DrawClone(Option_t *option) const
532{
533
534 gROOT->SetSelectedPad(NULL);
535
536 MHCalibrationPixel *newobj = (MHCalibrationPixel*)Clone();
537
538 if (!newobj)
539 return 0;
540 newobj->SetBit(kCanDelete);
541
542 if (strlen(option))
543 newobj->Draw(option);
544 else
545 newobj->Draw(GetDrawOption());
546
547 return newobj;
548}
549
550
551
552// -------------------------------------------------------------------------
553//
554// Draw the histogram
555//
556void MHCalibrationPixel::Draw(Option_t *opt)
557{
558
559 if (!fHivsLoGain)
560 FitHiGainvsLoGain();
561
562 gStyle->SetOptFit(0);
563 gStyle->SetOptStat(111111);
564
565 gROOT->SetSelectedPad(NULL);
566
567 TCanvas *c = MH::MakeDefCanvas(this,600,900);
568 c->SetBit(kCanDelete);
569
570 c->Divide(2,5);
571
572 c->cd(1);
573 gPad->SetBorderMode(0);
574 gPad->SetTicks();
575
576 if (fHChargeHiGain->Integral() > 0)
577 gPad->SetLogy();
578
579 fHChargeHiGain->Draw(opt);
580
581 if (IsUseLoGain())
582 {
583
584 c->cd(2);
585 gPad->SetTicks();
586
587 if (fHChargeLoGain->Integral() > 0)
588 gPad->SetLogy();
589
590 fHChargeLoGain->Draw(opt);
591
592 if (fChargeGausFit)
593 {
594 fChargeGausFit->SetLineColor(IsChargeFitOK() ? kGreen : kRed);
595 fChargeGausFit->Draw("same");
596 }
597
598 c->cd(3);
599 gROOT->SetSelectedPad(NULL);
600 gStyle->SetOptFit();
601 if (fHivsLoGain)
602 fHivsLoGain->Draw("prof");
603
604 c->cd(4);
605 DrawLegend();
606 }
607 else
608 {
609 if (fChargeGausFit)
610 {
611 fChargeGausFit->SetLineColor(IsChargeFitOK() ? kGreen : kRed);
612 fChargeGausFit->Draw("same");
613 }
614
615 c->cd(2);
616 gPad->SetTicks();
617
618 if (fHChargeLoGain->Integral() > 0)
619 gPad->SetLogy(1);
620
621 fHChargeLoGain->Draw(opt);
622
623 c->cd(3);
624 DrawLegend();
625
626 c->cd(4);
627
628 gROOT->SetSelectedPad(NULL);
629 gStyle->SetOptFit();
630 if (fHivsLoGain)
631 fHivsLoGain->Draw("prof");
632 }
633
634 c->cd(5);
635 gPad->SetTicks();
636 fHAbsTimeHiGain->Draw(opt);
637
638 c->cd(6);
639 gPad->SetTicks();
640 fHAbsTimeLoGain->Draw(opt);
641
642 CreateChargeXaxis(fHiGains->GetSize());
643
644 c->cd(7);
645 gPad->SetTicks();
646 TGraph *gr1 = new TGraph(fChargeXaxis->GetSize(),
647 fChargeXaxis->GetArray(),
648 fHiGains->GetArray());
649 gr1->SetTitle("Evolution of HiGain charges with time");
650 gr1->GetXaxis()->SetTitle("Time [s]");
651 gr1->GetYaxis()->SetTitle("Sum FADC slices");
652 gr1->SetBit(kCanDelete);
653 gr1->Draw("AL");
654
655 CreateChargeXaxis(fLoGains->GetSize());
656
657 c->cd(8);
658 gPad->SetTicks();
659 TGraph *gr2 = new TGraph(fChargeXaxis->GetSize(),
660 fChargeXaxis->GetArray(),
661 fLoGains->GetArray());
662 gr2->SetTitle("Evolution of LoGain charges with time");
663 gr2->GetXaxis()->SetTitle("Time [s]");
664 gr2->GetYaxis()->SetTitle("Sum FADC slices");
665 gr2->SetBit(kCanDelete);
666 gr2->Draw("AL");
667
668 c->Modified();
669 c->Update();
670
671 c->cd(9);
672
673 TArrayF *array;
674 TString title = "Power Spectrum Density ";
675
676 if(IsUseLoGain())
677 {
678 if (!fPSDLoGain)
679 return;
680 array = fPSDLoGain;
681 title += "LoGain";
682 }
683 else
684 {
685 if (!fPSDHiGain)
686 return;
687 array = fPSDHiGain;
688 title += "HiGain";
689 }
690
691 if (!fPSDXaxis)
692 CreatePSDXaxis(array->GetSize());
693
694 TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
695 gr3->SetTitle(title.Data());
696 gr3->GetXaxis()->SetTitle("Freq. [Hz]");
697 gr3->GetYaxis()->SetTitle("P(f)");
698 gr3->SetBit(kCanDelete);
699 gr3->Draw("AL");
700
701 c->Modified();
702 c->Update();
703
704 c->cd(10);
705
706 gStyle->SetOptStat(111111);
707 gStyle->SetOptFit();
708 gPad->SetTicks();
709
710 if (fHPSD->Integral() > 0)
711 gPad->SetLogy();
712
713 fHPSD->Draw(opt);
714 fHPSD->GetXaxis()->SetTitle("P(f)");
715 fHPSD->GetYaxis()->SetTitle("Counts");
716
717 if (fPSDExpFit)
718 {
719 fPSDExpFit->SetLineColor(IsOscillating() ? kRed : kGreen);
720 fPSDExpFit->Draw("same");
721 }
722
723 c->Modified();
724 c->Update();
725
726 return;
727}
728
729void MHCalibrationPixel::FitHiGainvsLoGain()
730{
731
732 *fLog << inf << "Fit Hi Gain vs Lo Gain " << endl;
733
734 if (fHivsLoGain)
735 return;
736
737 if ((!fHiGains) || (!fLoGains) || (fHiGains->GetSize() == 0) || (fLoGains->GetSize() == 0))
738 return;
739
740 gStyle->SetOptFit();
741
742 fHivsLoGain = new TProfile("HiGainvsLoGain","Plot the High Gain vs. Low Gain",100,0.,1000.,0.,1000.);
743 fHivsLoGain->GetXaxis()->SetTitle("Sum of Charges High Gain");
744 fHivsLoGain->GetYaxis()->SetTitle("Sum of Charges Low Gain");
745 fHivsLoGain->SetName(Form("%s%d",fHivsLoGain->GetName(),fPixId));
746
747 for (Int_t i=0;i<fTotalEntries;i++)
748 fHivsLoGain->Fill(fHiGains->At(i),fLoGains->At(i),1);
749
750 fHivsLoGain->Fit("pol1","rq","",fHivsLoGain->GetMean()-2.5*fHivsLoGain->GetRMS(),
751 fHivsLoGain->GetMean()+2.5*fHivsLoGain->GetRMS());
752
753 fOffset = fHivsLoGain->GetFunction("pol1")->GetParameter(0);
754 fSlope = fHivsLoGain->GetFunction("pol1")->GetParameter(1);
755
756}
757
758
759
760Bool_t MHCalibrationPixel::FitCharge(Option_t *option)
761{
762
763 if (fChargeGausFit)
764 return kFALSE;
765
766 //
767 // Get the fitting ranges
768 //
769 Axis_t rmin = fChargeFirstHiGain;
770 Axis_t rmax = fChargeLastHiGain;
771 TH1F *hist = fHChargeHiGain;
772
773 if (TESTBIT(fFlags,kUseLoGain))
774 {
775 rmin = fChargeFirstLoGain;
776 rmax = fChargeLastLoGain;
777 hist = fHChargeLoGain;
778 }
779
780 //
781 // First guesses for the fit (should be as close to reality as possible,
782 //
783 const Stat_t entries = hist->Integral("width");
784 const Double_t mu_guess = hist->GetBinCenter(hist->GetMaximumBin());
785 const Double_t sigma_guess = (rmax-rmin)/2.;
786 const Double_t area_guess = entries/gkSq2Pi/sigma_guess;
787
788 fChargeGausFit = new TF1(Form("%s%d","ChargeGausFit",fPixId),"gaus",rmin,rmax);
789
790 if (!fChargeGausFit)
791 {
792 *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl;
793 return kFALSE;
794 }
795
796 fChargeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
797 fChargeGausFit->SetParNames("Area","#mu","#sigma");
798 fChargeGausFit->SetParLimits(0,0.,entries);
799 fChargeGausFit->SetParLimits(1,rmin,rmax);
800 fChargeGausFit->SetParLimits(2,0.,rmax-rmin);
801 fChargeGausFit->SetRange(rmin,rmax);
802
803 hist->Fit(fChargeGausFit,option);
804
805 //
806 // If we are not able to fit, try once again
807 //
808 if (fChargeGausFit->GetProb() < fProbLimit)
809 {
810
811 Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2);
812 rmin = (rtry < rmin ? rmin : rtry);
813 rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2);
814
815 fChargeGausFit->SetRange(rmin,rmax);
816 hist->Fit(fChargeGausFit,option);
817 }
818
819 fChargeChisquare = fChargeGausFit->GetChisquare();
820 fChargeNdf = fChargeGausFit->GetNDF();
821 fChargeProb = fChargeGausFit->GetProb();
822 fChargeMean = fChargeGausFit->GetParameter(1);
823 fChargeMeanErr = fChargeGausFit->GetParError(1);
824 fChargeSigma = fChargeGausFit->GetParameter(2);
825 fChargeSigmaErr = fChargeGausFit->GetParError(2);
826
827 //
828 // From the absolute time, we only take the mean and RMS
829 //
830 fAbsTimeMean = (Float_t)fHAbsTimeHiGain->GetMean();
831 fAbsTimeRms = (Float_t)fHAbsTimeHiGain->GetRMS();
832 fAbsTimeMeanErr = (Float_t)fAbsTimeRms / TMath::Sqrt(fHAbsTimeHiGain->GetEntries());
833
834 if (TESTBIT(fFlags,kUseLoGain))
835 {
836 fAbsTimeMean = fHAbsTimeLoGain->GetMean();
837 fAbsTimeRms = fHAbsTimeLoGain->GetRMS();
838 fAbsTimeMeanErr = fAbsTimeRms / TMath::Sqrt(fHAbsTimeLoGain->GetEntries());
839 }
840
841 //
842 // The fit result is accepted under condition:
843 // 1) The results are not nan's
844 // 2) The NDF is not smaller than fNDFLimit (5)
845 // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
846 //
847 // Otherwise take means and RMS of the histograms
848 //
849 if ( TMath::IsNaN(fChargeMean)
850 || TMath::IsNaN(fChargeMeanErr)
851 || TMath::IsNaN(fChargeProb)
852 || TMath::IsNaN(fChargeSigma)
853 || TMath::IsNaN(fChargeSigmaErr)
854 || (fChargeNdf < fNDFLimit)
855 || (fChargeProb < fProbLimit) )
856 {
857
858 fChargeMean = hist->GetMean();
859 fChargeMeanErr = hist->GetRMS()/hist->GetEntries();
860 fChargeSigma = hist->GetRMS();
861 fChargeSigmaErr = fChargeMeanErr/2.;
862
863 CLRBIT(fFlags,kChargeFitOK);
864 return kFALSE;
865 }
866
867 SETBIT(fFlags,kChargeFitOK);
868 return kTRUE;
869}
870
871
872void MHCalibrationPixel::CutAllEdges()
873{
874
875 Int_t nbins = 30;
876
877 CutEdges(fHChargeHiGain,nbins);
878
879 fChargeFirstHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetFirst());
880 fChargeLastHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetLast())
881 +fHChargeHiGain->GetBinWidth(0);
882
883 CutEdges(fHChargeLoGain,nbins);
884
885 fChargeFirstLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetFirst());
886 fChargeLastLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetLast())
887 +fHChargeLoGain->GetBinWidth(0);
888
889
890 CutEdges(fHAbsTimeHiGain,0);
891
892 fAbsTimeFirstHiGain = fHAbsTimeHiGain->GetBinLowEdge(fHAbsTimeHiGain->GetXaxis()->GetFirst());
893 fAbsTimeLastHiGain = fHAbsTimeHiGain->GetBinLowEdge(fHAbsTimeHiGain->GetXaxis()->GetLast())
894 +fHAbsTimeHiGain->GetBinWidth(0);
895
896 CutEdges(fHAbsTimeLoGain,0);
897
898 fAbsTimeFirstLoGain = fHAbsTimeLoGain->GetBinLowEdge(fHAbsTimeLoGain->GetXaxis()->GetFirst());
899 fAbsTimeLastLoGain = fHAbsTimeLoGain->GetBinLowEdge(fHAbsTimeLoGain->GetXaxis()->GetLast())
900 +fHAbsTimeLoGain->GetBinWidth(0);
901
902
903}
904
905void MHCalibrationPixel::PrintChargeFitResult()
906{
907
908 *fLog << all << "Results of the Summed Charges Fit: " << endl;
909 *fLog << all << "Chisquare: " << fChargeChisquare << endl;
910 *fLog << all << "DoF: " << fChargeNdf << endl;
911 *fLog << all << "Probability: " << fChargeProb << endl;
912 *fLog << all << endl;
913
914}
915
Note: See TracBrowser for help on using the repository browser.