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

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