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

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