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

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