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

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