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

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