source: branches/Corsika7500Compatibility/mhcalib/MHPedestalPix.cc@ 18728

Last change on this file since 18728 was 8020, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 8.3 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 02/2005 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHPedestalPix
28//
29// A base class for events which are believed to follow a Gaussian distribution
30// with time, e.g. calibration events, observables containing white noise, ...
31//
32// MHPedestalPix derives from MHGausEvents, thus all features of
33// MHGausEvents can be used by a class deriving from MHPedestalPix
34//
35// As an additional feature to MHGausEvents, this class offers to skip the fitting
36// to set mean, sigma and its errors directly from the histograms with the function
37// BypassFit()
38//
39// See also: MHGausEvents
40//
41//////////////////////////////////////////////////////////////////////////////
42#include "MHPedestalPix.h"
43
44#include <TH1.h>
45#include <TF1.h>
46#include <TGraph.h>
47
48#include "MLog.h"
49#include "MLogManip.h"
50
51ClassImp(MHPedestalPix);
52
53using namespace std;
54
55// --------------------------------------------------------------------------
56//
57// Default Constructor.
58//
59MHPedestalPix::MHPedestalPix(const char *name, const char *title)
60{
61
62 fName = name ? name : "MHPedestalPix";
63 fTitle = title ? title : "Pedestal histogram events";
64
65}
66
67// -------------------------------------------------------------------------------
68//
69// Fits the histogram to a double Gauss.
70//
71//
72Bool_t MHPedestalPix::FitDoubleGaus(const Double_t xmin, const Double_t xmax, Option_t *option)
73{
74
75 if (IsGausFitOK())
76 return kTRUE;
77
78 StripZeros(&fHGausHist,0);
79
80 TAxis *axe = fHGausHist.GetXaxis();
81 //
82 // Get the fitting ranges
83 //
84 Axis_t rmin = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetFirst()) : xmin;
85 Axis_t rmax = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetLast()) : xmax;
86
87 //
88 // First guesses for the fit (should be as close to reality as possible,
89 //
90 const Stat_t entries = fHGausHist.Integral(axe->FindBin(rmin),axe->FindBin(rmax),"width");
91 const Double_t sigma_guess = fHGausHist.GetRMS();
92 const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess;
93
94 fFGausFit = new TF1("GausFit","gaus(0)+gaus(3)",rmin,rmax);
95
96 if (!fFGausFit)
97 {
98 *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit "
99 << "in: " << fName << endl;
100 return kFALSE;
101 }
102
103 //
104 // For the fits, we have to take special care since ROOT
105 // has stored the function pointer in a global list which
106 // lead to removing the object twice. We have to take out
107 // the following functions of the global list of functions
108 // as well:
109 //
110 gROOT->GetListOfFunctions()->Remove(fFGausFit);
111
112 fFGausFit->SetParameters(area_guess/2.,0.,sigma_guess/2.,area_guess/2.,25.,sigma_guess/2.);
113 fFGausFit->SetParNames("Area_{0}","#mu_{0}","#sigma_{0}","Area_{1}","#mu_{1}","#sigma_{1}");
114 fFGausFit->SetParLimits(0,0.,area_guess*5.);
115 fFGausFit->SetParLimits(1,rmin,0.);
116 fFGausFit->SetParLimits(2,0.,rmax-rmin);
117 fFGausFit->SetParLimits(3,0.,area_guess*10.);
118 fFGausFit->SetParLimits(4,0.,rmax/2.);
119 fFGausFit->SetParLimits(5,0.,rmax-rmin);
120 fFGausFit->SetRange(rmin,rmax);
121
122 fHGausHist.Fit(fFGausFit,option);
123
124 SetMean (fFGausFit->GetParameter(4)-fFGausFit->GetParameter(1));
125 SetSigma (TMath::Sqrt(fFGausFit->GetParameter(5)*fFGausFit->GetParameter(5)
126 +fFGausFit->GetParameter(2)*fFGausFit->GetParameter(2)));
127 SetMeanErr (TMath::Sqrt(fFGausFit->GetParError(4)*fFGausFit->GetParError(4)
128 +fFGausFit->GetParError(1)*fFGausFit->GetParError(1)));
129 SetSigmaErr (TMath::Sqrt(fFGausFit->GetParError(5)*fFGausFit->GetParError(5)
130 +fFGausFit->GetParError(2)*fFGausFit->GetParError(2)));
131 SetProb (fFGausFit->GetProb());
132 //
133 // The fit result is accepted under condition:
134 // 1) The results are not nan's
135 // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
136 // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
137 //
138 // !Finitite means either infinite or not-a-number
139 if ( !TMath::Finite(GetMean())
140 || !TMath::Finite(GetMeanErr())
141 || !TMath::Finite(GetProb())
142 || !TMath::Finite(GetSigma())
143 || !TMath::Finite(GetSigmaErr())
144 || fProb < GetProbLimit())
145 return kFALSE;
146
147 SetGausFitOK(kTRUE);
148 return kTRUE;
149}
150
151
152// -------------------------------------------------------------------------------
153//
154// Fits the histogram to a triple Gauss.
155//
156Bool_t MHPedestalPix::FitTripleGaus(const Double_t xmin, const Double_t xmax, Option_t *option)
157{
158
159 if (IsGausFitOK())
160 return kTRUE;
161
162 StripZeros(&fHGausHist,0);
163
164 TAxis *axe = fHGausHist.GetXaxis();
165 //
166 // Get the fitting ranges
167 //
168 Axis_t rmin = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetFirst()) : xmin;
169 Axis_t rmax = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetLast()) : xmax;
170
171 //
172 // First guesses for the fit (should be as close to reality as possible,
173 //
174 const Stat_t entries = fHGausHist.Integral(axe->FindBin(rmin),axe->FindBin(rmax),"width");
175 const Double_t sigma_guess = fHGausHist.GetRMS();
176 const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess;
177
178 fFGausFit = new TF1("GausFit","gaus(0)+gaus(3)+gaus(6)",rmin,rmax);
179
180 if (!fFGausFit)
181 {
182 *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit "
183 << "in: " << fName << endl;
184 return kFALSE;
185 }
186
187 //
188 // For the fits, we have to take special care since ROOT
189 // has stored the function pointer in a global list which
190 // lead to removing the object twice. We have to take out
191 // the following functions of the global list of functions
192 // as well:
193 //
194 gROOT->GetListOfFunctions()->Remove(fFGausFit);
195
196 fFGausFit->SetParameters(10.,-4.0,1.5,70.,1.5,6.,5.,7.,7.);
197 fFGausFit->SetParNames("Area_{0}","#mu_{0}","#sigma_{0}","Area_{1}","#mu_{1}","#sigma_{1}","Area_{2}","#mu_{2}","#sigma_{2}");
198 fFGausFit->SetParLimits(0,0.,area_guess*2.5);
199 fFGausFit->SetParLimits(1,-9.0,-2.2);
200 fFGausFit->SetParLimits(2,-1.0,15.);
201 fFGausFit->SetParLimits(3,0.,area_guess*10.);
202 fFGausFit->SetParLimits(4,-4.5,2.);
203 fFGausFit->SetParLimits(5,0.,(rmax-rmin)/3.);
204 fFGausFit->SetParLimits(6,0.,area_guess*5.);
205 fFGausFit->SetParLimits(7,6.,20.);
206 fFGausFit->SetParLimits(8,5.,40.);
207 fFGausFit->SetRange(rmin,rmax);
208
209 fHGausHist.Fit(fFGausFit,option);
210
211 SetMean (fFGausFit->GetParameter(4)-fFGausFit->GetParameter(1));
212 SetSigma (TMath::Sqrt(fFGausFit->GetParameter(5)*fFGausFit->GetParameter(5)
213 +fFGausFit->GetParameter(2)*fFGausFit->GetParameter(2)));
214 SetMeanErr (TMath::Sqrt(fFGausFit->GetParError(4)*fFGausFit->GetParError(4)
215 +fFGausFit->GetParError(1)*fFGausFit->GetParError(1)));
216 SetSigmaErr (TMath::Sqrt(fFGausFit->GetParError(5)*fFGausFit->GetParError(5)
217 +fFGausFit->GetParError(2)*fFGausFit->GetParError(2)));
218 SetProb (fFGausFit->GetProb());
219 //
220 // The fit result is accepted under condition:
221 // 1) The results are not nan's
222 // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
223 // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
224 //
225 // !Finitite means either infinite or not-a-number
226 if ( !TMath::Finite(GetMean())
227 || !TMath::Finite(GetMeanErr())
228 || !TMath::Finite(GetProb())
229 || !TMath::Finite(GetSigma())
230 || !TMath::Finite(GetSigmaErr())
231 || fProb < GetProbLimit() )
232 return kFALSE;
233
234 SetGausFitOK(kTRUE);
235 return kTRUE;
236}
237
238
Note: See TracBrowser for help on using the repository browser.