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

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