source: trunk/MagicSoft/Mars/mhcalib/MHPedestalPix.cc@ 6963

Last change on this file since 6963 was 6329, checked in by gaug, 20 years ago
*** empty log message ***
File size: 8.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 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 if ( TMath::IsNaN(GetMean())
139 || TMath::IsNaN(GetMeanErr())
140 || TMath::IsNaN(GetProb())
141 || TMath::IsNaN(GetSigma())
142 || TMath::IsNaN(GetSigmaErr())
143 || fProb < fProbLimit )
144 return kFALSE;
145
146 SetGausFitOK(kTRUE);
147 return kTRUE;
148}
149
150
151// -------------------------------------------------------------------------------
152//
153// Fits the histogram to a triple Gauss.
154//
155Bool_t MHPedestalPix::FitTripleGaus(const Double_t xmin, const Double_t xmax, Option_t *option)
156{
157
158 if (IsGausFitOK())
159 return kTRUE;
160
161 StripZeros(&fHGausHist,0);
162
163 TAxis *axe = fHGausHist.GetXaxis();
164 //
165 // Get the fitting ranges
166 //
167 Axis_t rmin = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetFirst()) : xmin;
168 Axis_t rmax = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetLast()) : xmax;
169
170 //
171 // First guesses for the fit (should be as close to reality as possible,
172 //
173 const Stat_t entries = fHGausHist.Integral(axe->FindBin(rmin),axe->FindBin(rmax),"width");
174 const Double_t sigma_guess = fHGausHist.GetRMS();
175 const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess;
176
177 fFGausFit = new TF1("GausFit","gaus(0)+gaus(3)+gaus(6)",rmin,rmax);
178
179 if (!fFGausFit)
180 {
181 *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit "
182 << "in: " << fName << endl;
183 return kFALSE;
184 }
185
186 //
187 // For the fits, we have to take special care since ROOT
188 // has stored the function pointer in a global list which
189 // lead to removing the object twice. We have to take out
190 // the following functions of the global list of functions
191 // as well:
192 //
193 gROOT->GetListOfFunctions()->Remove(fFGausFit);
194
195 fFGausFit->SetParameters(10.,-4.0,1.5,70.,1.5,6.,5.,7.,7.);
196 fFGausFit->SetParNames("Area_{0}","#mu_{0}","#sigma_{0}","Area_{1}","#mu_{1}","#sigma_{1}","Area_{2}","#mu_{2}","#sigma_{2}");
197 fFGausFit->SetParLimits(0,0.,area_guess*2.5);
198 fFGausFit->SetParLimits(1,-9.0,-2.2);
199 fFGausFit->SetParLimits(2,-1.0,15.);
200 fFGausFit->SetParLimits(3,0.,area_guess*10.);
201 fFGausFit->SetParLimits(4,-4.5,2.);
202 fFGausFit->SetParLimits(5,0.,(rmax-rmin)/3.);
203 fFGausFit->SetParLimits(6,0.,area_guess*5.);
204 fFGausFit->SetParLimits(7,6.,20.);
205 fFGausFit->SetParLimits(8,5.,40.);
206 fFGausFit->SetRange(rmin,rmax);
207
208 fHGausHist.Fit(fFGausFit,option);
209
210 SetMean (fFGausFit->GetParameter(4)-fFGausFit->GetParameter(1));
211 SetSigma (TMath::Sqrt(fFGausFit->GetParameter(5)*fFGausFit->GetParameter(5)
212 +fFGausFit->GetParameter(2)*fFGausFit->GetParameter(2)));
213 SetMeanErr (TMath::Sqrt(fFGausFit->GetParError(4)*fFGausFit->GetParError(4)
214 +fFGausFit->GetParError(1)*fFGausFit->GetParError(1)));
215 SetSigmaErr (TMath::Sqrt(fFGausFit->GetParError(5)*fFGausFit->GetParError(5)
216 +fFGausFit->GetParError(2)*fFGausFit->GetParError(2)));
217 SetProb (fFGausFit->GetProb());
218 //
219 // The fit result is accepted under condition:
220 // 1) The results are not nan's
221 // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
222 // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
223 //
224 if ( TMath::IsNaN(GetMean())
225 || TMath::IsNaN(GetMeanErr())
226 || TMath::IsNaN(GetProb())
227 || TMath::IsNaN(GetSigma())
228 || TMath::IsNaN(GetSigmaErr())
229 || fProb < fProbLimit )
230 return kFALSE;
231
232 SetGausFitOK(kTRUE);
233 return kTRUE;
234}
235
236
Note: See TracBrowser for help on using the repository browser.