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

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