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

Last change on this file since 2969 was 2950, checked in by gaug, 21 years ago
*** empty log message ***
File size: 22.8 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 Reset();
334
335 return kTRUE;
336}
337
338
339Bool_t MHCalibrationPixel::UseLoGain()
340{
341
342 if (fHChargeHiGain->Integral() > fHChargeLoGain->Integral())
343 {
344 CLRBIT(fFlags,kUseLoGain);
345 return kFALSE;
346 }
347 else
348 {
349 SETBIT(fFlags,kUseLoGain);
350 return kTRUE;
351 }
352}
353
354Bool_t MHCalibrationPixel::FillPointInGraph(Float_t qhi,Float_t qlo)
355{
356
357 fHiGains->Set(++fTotalEntries);
358 fLoGains->Set(fTotalEntries);
359
360 fHiGains->AddAt(qhi,fTotalEntries-1);
361 fLoGains->AddAt(qlo,fTotalEntries-1);
362
363 return kTRUE;
364
365}
366
367
368// -------------------------------------------------------------------------
369//
370// Draw a legend with the fit results
371//
372void MHCalibrationPixel::DrawLegend()
373{
374
375 fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
376
377 if (IsFitOK())
378 fFitLegend->SetFillColor(80);
379 else
380 fFitLegend->SetFillColor(2);
381
382 fFitLegend->SetLabel("Results of the Gauss Fit:");
383 fFitLegend->SetTextSize(0.05);
384
385 const TString line1 =
386 Form("Mean: Q_{#mu} = %2.2f #pm %2.2f",fChargeMean,fChargeMeanErr);
387 TText *t1 = fFitLegend->AddText(line1);
388 t1->SetBit(kCanDelete);
389
390 const TString line4 =
391 Form("Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fChargeSigma,fChargeSigmaErr);
392 TText *t2 = fFitLegend->AddText(line4);
393 t2->SetBit(kCanDelete);
394
395 const TString line7 =
396 Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChargeChisquare,fChargeNdf);
397 TText *t3 = fFitLegend->AddText(line7);
398 t3->SetBit(kCanDelete);
399
400 const TString line8 =
401 Form("Probability: %4.3f ",fChargeProb);
402 TText *t4 = fFitLegend->AddText(line8);
403 t4->SetBit(kCanDelete);
404
405 if (IsFitOK())
406 {
407 TText *t5 = fFitLegend->AddText("Result of the Fit: OK");
408 t5->SetBit(kCanDelete);
409 }
410 else
411 {
412 TText *t6 = fFitLegend->AddText("Result of the Fit: NOT OK");
413 t6->SetBit(kCanDelete);
414 }
415
416 fFitLegend->SetBit(kCanDelete);
417 fFitLegend->Draw();
418
419}
420
421
422TObject *MHCalibrationPixel::DrawClone(Option_t *option) const
423{
424
425 gROOT->SetSelectedPad(NULL);
426
427 MHCalibrationPixel *newobj = (MHCalibrationPixel*)Clone();
428
429 if (!newobj)
430 return 0;
431 newobj->SetBit(kCanDelete);
432
433
434 if (strlen(option))
435 newobj->Draw(option);
436 else
437 newobj->Draw(GetDrawOption());
438
439 return newobj;
440}
441
442
443
444// -------------------------------------------------------------------------
445//
446// Draw the histogram
447//
448void MHCalibrationPixel::Draw(Option_t *opt)
449{
450
451 if (!fHivsLoGain)
452 FitHiGainvsLoGain();
453
454 gStyle->SetOptFit(0);
455 gStyle->SetOptStat(111111);
456
457 gROOT->SetSelectedPad(NULL);
458
459 TCanvas *c = MakeDefCanvas(this,600,900);
460
461 c->Divide(2,4);
462
463 c->cd(1);
464 gPad->SetBorderMode(0);
465 gPad->SetTicks();
466
467 if (fHChargeHiGain->Integral() > 0)
468 gPad->SetLogy(1);
469 else
470 gPad->SetLogy(0);
471
472 fHChargeHiGain->Draw(opt);
473
474 c->Modified();
475 c->Update();
476
477 if (IsUseLoGain())
478 {
479
480 c->cd(2);
481 gPad->SetTicks();
482
483 if (fHChargeLoGain->Integral() > 0)
484 gPad->SetLogy(1);
485 else
486 gPad->SetLogy(0);
487
488 fHChargeLoGain->Draw(opt);
489
490 if (fChargeGausFit)
491 {
492 if (IsFitOK())
493 fChargeGausFit->SetLineColor(kGreen);
494 else
495 fChargeGausFit->SetLineColor(kRed);
496
497 fChargeGausFit->Draw("same");
498 }
499 c->Modified();
500 c->Update();
501
502 c->cd(3);
503 gROOT->SetSelectedPad(NULL);
504 gStyle->SetOptFit();
505 if (fHivsLoGain)
506 fHivsLoGain->Draw("prof");
507 gPad->Modified();
508 gPad->Update();
509
510 c->cd(4);
511 DrawLegend();
512
513 }
514 else
515 {
516 if (fChargeGausFit)
517 {
518 if (IsFitOK())
519 fChargeGausFit->SetLineColor(kGreen);
520 else
521 fChargeGausFit->SetLineColor(kRed);
522
523 fChargeGausFit->Draw("same");
524 }
525
526 c->cd(2);
527 gPad->SetTicks();
528
529 if (fHChargeLoGain->Integral() > 0)
530 gPad->SetLogy(1);
531 else
532 gPad->SetLogy(0);
533
534 fHChargeLoGain->Draw(opt);
535 c->Modified();
536 c->Update();
537
538 c->cd(3);
539 DrawLegend();
540
541 c->cd(4);
542
543 gROOT->SetSelectedPad(NULL);
544 gStyle->SetOptFit();
545 if (fHivsLoGain)
546 fHivsLoGain->Draw("prof");
547 gPad->Modified();
548 gPad->Update();
549 }
550
551 c->Modified();
552 c->Update();
553
554 c->cd(5);
555 gStyle->SetOptStat(1111111);
556
557 gPad->SetTicks();
558 gPad->SetLogy(0);
559 fHTimeHiGain->Draw(opt);
560 c->Modified();
561 c->Update();
562
563 if (IsUseLoGain())
564 {
565
566 c->cd(6);
567 gPad->SetTicks();
568 gPad->SetLogy(0);
569 fHTimeLoGain->Draw(opt);
570 c->Modified();
571 c->Update();
572
573 if (fTimeGausFit)
574 {
575 if (fTimeChisquare > 20.)
576 fTimeGausFit->SetLineColor(kRed);
577 else
578 fTimeGausFit->SetLineColor(kGreen);
579
580 fTimeGausFit->Draw("same");
581 c->Modified();
582 c->Update();
583 }
584 }
585 else
586 {
587 if (fTimeGausFit)
588 {
589 if (fTimeChisquare > 20.)
590 fTimeGausFit->SetLineColor(kRed);
591 else
592 fTimeGausFit->SetLineColor(kGreen);
593
594 fTimeGausFit->Draw("same");
595 c->Modified();
596 c->Update();
597 }
598
599 c->cd(6);
600 gPad->SetTicks();
601 gPad->SetLogy(0);
602 fHTimeLoGain->Draw(opt);
603 c->Modified();
604 c->Update();
605
606 }
607 c->Modified();
608 c->Update();
609
610 c->cd(7);
611 gPad->SetTicks();
612 fHChargevsNHiGain->Draw(opt);
613 c->Modified();
614 c->Update();
615
616 c->cd(8);
617 gPad->SetTicks();
618 fHChargevsNLoGain->Draw(opt);
619 c->Modified();
620 c->Update();
621
622 return;
623}
624
625void MHCalibrationPixel::FitHiGainvsLoGain()
626{
627
628 *fLog << inf << "Fit Hi Gain vs Lo Gain " << endl;
629
630 if (fHivsLoGain)
631 return;
632
633 if ((!fHiGains) || (!fLoGains) || (fHiGains->GetSize() == 0) || (fLoGains->GetSize() == 0))
634 return;
635
636 gStyle->SetOptFit();
637
638 fHivsLoGain = new TProfile("HiGainvsLoGain","Plot the High Gain vs. Low Gain",100,0.,1000.,0.,1000.);
639 fHivsLoGain->GetXaxis()->SetTitle("Sum of Charges High Gain");
640 fHivsLoGain->GetYaxis()->SetTitle("Sum of Charges Low Gain");
641
642 TString titleHiLoGain = TString(fHivsLoGain->GetName());
643 titleHiLoGain += fPixId;
644 fHivsLoGain->SetName(titleHiLoGain.Data());
645
646 for (Int_t i=0;i<fTotalEntries;i++)
647 fHivsLoGain->Fill(fHiGains->At(i),fLoGains->At(i),1);
648
649 fHivsLoGain->Fit("pol1","rq","",fHivsLoGain->GetMean()-2.5*fHivsLoGain->GetRMS(),
650 fHivsLoGain->GetMean()+2.0*fHivsLoGain->GetRMS());
651
652 fOffset = fHivsLoGain->GetFunction("pol1")->GetParameter(0);
653 fSlope = fHivsLoGain->GetFunction("pol1")->GetParameter(1);
654
655}
656
657
658Bool_t MHCalibrationPixel::FitTimeHiGain(Axis_t rmin, Axis_t rmax, Option_t *option)
659{
660
661 if (fTimeGausFit)
662 return kFALSE;
663
664 rmin = (rmin != 0.) ? rmin : (Axis_t)fTimeLowerFitRangeHiGain;
665 rmax = (rmax != 0.) ? rmax : (Axis_t)fTimeUpperFitRangeHiGain;
666
667 const Stat_t entries = fHTimeHiGain->Integral();
668 const Double_t mu_guess = fHTimeHiGain->GetBinCenter(fHTimeHiGain->GetMaximumBin());
669 const Double_t sigma_guess = (rmax - rmin)/2.;
670 const Double_t area_guess = entries/gkSq2Pi;
671
672 TString name = TString("GausTime");
673 name += fPixId;
674 fTimeGausFit = new TF1(name.Data(),"gaus",rmin,rmax);
675
676 if (!fTimeGausFit)
677 {
678 *fLog << warn << dbginf << "WARNING: Could not create fit function for Time fit" << endl;
679 return kFALSE;
680 }
681
682 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
683 fTimeGausFit->SetParNames("Area","#mu","#sigma");
684 fTimeGausFit->SetParLimits(0,0.,entries);
685 fTimeGausFit->SetParLimits(1,rmin,rmax);
686 fTimeGausFit->SetParLimits(2,0.,(rmax-rmin));
687 fTimeGausFit->SetRange(rmin,rmax);
688
689 fHTimeHiGain->Fit(fTimeGausFit,option);
690
691 rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2);
692 rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2);
693 fTimeGausFit->SetRange(rmin,rmax);
694
695 fHTimeHiGain->Fit(fTimeGausFit,option);
696
697 fTimeChisquare = fTimeGausFit->GetChisquare();
698 fTimeNdf = fTimeGausFit->GetNDF();
699 fTimeProb = fTimeGausFit->GetProb();
700
701 fTimeMean = fTimeGausFit->GetParameter(1);
702 fTimeSigma = fTimeGausFit->GetParameter(2);
703
704 if (TMath::IsNaN(fTimeMean) || TMath::IsNaN(fTimeSigma))
705 return kFALSE;
706
707 if (TMath::IsNaN(fTimeChisquare) || fTimeChisquare > 20.) // Cannot use Probability because Ndf is sometimes < 1
708 return kFALSE;
709
710 return kTRUE;
711
712}
713
714Bool_t MHCalibrationPixel::FitTimeLoGain(Axis_t rmin, Axis_t rmax, Option_t *option)
715{
716
717 if (fTimeGausFit)
718 return kFALSE;
719
720 rmin = (rmin != 0.) ? rmin : (Axis_t)fTimeLowerFitRangeLoGain;
721 rmax = (rmax != 0.) ? rmax : (Axis_t)fTimeUpperFitRangeLoGain;
722
723 const Stat_t entries = fHTimeLoGain->Integral();
724 const Double_t mu_guess = fHTimeLoGain->GetBinCenter(fHTimeLoGain->GetMaximumBin());
725 const Double_t sigma_guess = (rmax - rmin)/2.;
726 const Double_t area_guess = entries/gkSq2Pi;
727
728 TString name = TString("GausTime");
729 name += fPixId;
730 fTimeGausFit = new TF1(name.Data(),"gaus",rmin,rmax);
731
732 if (!fTimeGausFit)
733 {
734 *fLog << warn << dbginf << "WARNING: Could not create fit function for Time fit" << endl;
735 return kFALSE;
736 }
737
738 fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
739 fTimeGausFit->SetParNames("Area","#mu","#sigma");
740 fTimeGausFit->SetParLimits(0,0.,entries);
741 fTimeGausFit->SetParLimits(1,rmin,rmax);
742 fTimeGausFit->SetParLimits(2,0.,(rmax-rmin));
743 fTimeGausFit->SetRange(rmin,rmax);
744
745 fHTimeLoGain->Fit(fTimeGausFit,option);
746
747 rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2);
748 rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2);
749 fTimeGausFit->SetRange(rmin,rmax);
750
751 fHTimeLoGain->Fit(fTimeGausFit,option);
752
753 fTimeChisquare = fTimeGausFit->GetChisquare();
754 fTimeNdf = fTimeGausFit->GetNDF();
755 fTimeProb = fTimeGausFit->GetProb();
756
757 fTimeMean = fTimeGausFit->GetParameter(1);
758 fTimeSigma = fTimeGausFit->GetParameter(2);
759
760 if (fTimeChisquare > 20.) // Cannot use Probability because Ndf is sometimes < 1
761 {
762 *fLog << warn << "WARNING: Fit of the Arrival times failed ! " << endl;
763 return kFALSE;
764 }
765
766 return kTRUE;
767
768}
769
770Bool_t MHCalibrationPixel::FitCharge(Option_t *option)
771{
772
773 if (fChargeGausFit)
774 return kFALSE;
775
776 //
777 // Get the fitting ranges
778 //
779 Axis_t rmin = fChargeFirstHiGain;
780 if (TESTBIT(fFlags,kUseLoGain))
781 rmin = fChargeFirstLoGain;
782
783 Axis_t rmax = fChargeLastHiGain;
784 if (TESTBIT(fFlags,kUseLoGain))
785 rmin = fChargeFirstLoGain;
786
787 TH1F *hist = fHChargeHiGain;
788 if (TESTBIT(fFlags,kUseLoGain))
789 hist = fHChargeLoGain;
790
791 //
792 // First guesses for the fit (should be as close to reality as possible,
793 // otherwise the fit goes gaga because of high number of dimensions ...
794 //
795 const Stat_t entries = hist->Integral();
796 const Double_t area_guess = entries/gkSq2Pi;
797 const Double_t mu_guess = hist->GetBinCenter(hist->GetMaximumBin());
798 const Double_t sigma_guess = mu_guess/15.;
799
800 TString name = TString("ChargeGausFit");
801 name += fPixId;
802
803 fChargeGausFit = new TF1(name.Data(),"gaus",rmin,rmax);
804
805 if (!fChargeGausFit)
806 {
807 *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl;
808 return kFALSE;
809 }
810
811 fChargeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
812 fChargeGausFit->SetParNames("Area","#mu","#sigma");
813 fChargeGausFit->SetParLimits(0,0.,entries);
814 fChargeGausFit->SetParLimits(1,rmin,rmax);
815 fChargeGausFit->SetParLimits(2,0.,rmax-rmin);
816 fChargeGausFit->SetRange(rmin,rmax);
817
818 hist->Fit(fChargeGausFit,option);
819
820 //
821 // If we are not able to fit, try once again
822 //
823 if (fChargeGausFit->GetProb() < gkProbLimit)
824 {
825
826 Axis_t rtry = fChargeGausFit->GetParameter(1) - 3.0*fChargeGausFit->GetParameter(2);
827 rmin = (rtry < rmin ? rmin : rtry);
828 rmax = fChargeGausFit->GetParameter(1) + 3.0*fChargeGausFit->GetParameter(2);
829 fChargeGausFit->SetRange(rmin,rmax);
830
831 fHChargeHiGain->Fit(fChargeGausFit,option);
832 }
833
834 fChargeChisquare = fChargeGausFit->GetChisquare();
835 fChargeNdf = fChargeGausFit->GetNDF();
836 fChargeProb = fChargeGausFit->GetProb();
837 fChargeMean = fChargeGausFit->GetParameter(1);
838 fChargeMeanErr = fChargeGausFit->GetParError(1);
839 fChargeSigma = fChargeGausFit->GetParameter(2);
840 fChargeSigmaErr = fChargeGausFit->GetParError(2);
841
842 //
843 // The fit result is accepted under condition:
844 // The Results are not nan's
845 // The Probability is greater than gkProbLimit (default 0.001 == 99.9%)
846 //
847 if (TMath::IsNaN(fChargeMean) || TMath::IsNaN(fChargeMeanErr))
848 {
849 CLRBIT(fFlags,kFitOK);
850 return kFALSE;
851 }
852
853 if ((fChargeProb < gkProbLimit) || (TMath::IsNaN(fChargeProb)))
854 {
855 CLRBIT(fFlags,kFitOK);
856 return kFALSE;
857 }
858
859 SETBIT(fFlags,kFitOK);
860 return kTRUE;
861}
862
863
864
865
866void MHCalibrationPixel::CutAllEdges()
867{
868
869 Int_t nbins = 30;
870
871 CutEdges(fHChargeHiGain,nbins);
872
873 fChargeFirstHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetFirst());
874 fChargeLastHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetLast())
875 +fHChargeHiGain->GetBinWidth(0);
876
877 CutEdges(fHChargeLoGain,nbins);
878
879 fChargeFirstLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetFirst());
880 fChargeLastLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetLast())
881 +fHChargeLoGain->GetBinWidth(0);
882
883 CutEdges(fHChargevsNHiGain,0);
884 CutEdges(fHChargevsNLoGain,0);
885
886}
887
888void MHCalibrationPixel::PrintChargeFitResult()
889{
890
891 *fLog << all << "Results of the Summed Charges Fit: " << endl;
892 *fLog << all << "Chisquare: " << fChargeChisquare << endl;
893 *fLog << all << "DoF: " << fChargeNdf << endl;
894 *fLog << all << "Probability: " << fChargeProb << endl;
895 *fLog << all << endl;
896
897}
898
899void MHCalibrationPixel::PrintTimeFitResult()
900{
901
902 *fLog << all << "Results of the Time Slices Fit: " << endl;
903 *fLog << all << "Chisquare: " << fTimeChisquare << endl;
904 *fLog << all << "Ndf: " << fTimeNdf << endl;
905 *fLog << all << "Probability: " << fTimeProb << endl;
906 *fLog << all << endl;
907
908}
Note: See TracBrowser for help on using the repository browser.