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

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