source: trunk/Mars/mhcalib/MHCalibrationPix.cc@ 19844

Last change on this file since 19844 was 8907, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 6.7 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/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHCalibrationPix
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// MHCalibrationPix derives from MHGausEvents, thus all features of
33// MHGausEvents can be used by a class deriving from MHCalibrationPix
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 "MHCalibrationPix.h"
43
44#include <TH1.h>
45#include <TF1.h>
46#include <TMath.h>
47#include <TGraph.h>
48
49#include "MLog.h"
50#include "MLogManip.h"
51
52ClassImp(MHCalibrationPix);
53
54using namespace std;
55
56const Float_t MHCalibrationPix::fgBlackoutLimit = 5.;
57const Float_t MHCalibrationPix::fgPickupLimit = 5.;
58// --------------------------------------------------------------------------
59//
60// Default Constructor.
61// Sets:
62// - the default number for fPickupLimit (fgPickupLimit)
63// - the default number for fBlackoutLimit (fgBlackoutLimit)
64//
65// Initializes:
66// - all variables to 0.
67//
68MHCalibrationPix::MHCalibrationPix(const char *name, const char *title)
69{
70
71 fName = name ? name : "MHCalibrationPix";
72 fTitle = title ? title : "Calibration histogram events";
73
74 Clear();
75
76 SetBlackoutLimit();
77 SetPickupLimit();
78}
79
80// --------------------------------------------------------------------------
81//
82// Default Clear(), can be overloaded.
83//
84// Sets:
85// - all other pointers to NULL
86// - all variables to 0., except fPixId to -1
87// - all flags to kFALSE
88//
89// - all pointers
90//
91void MHCalibrationPix::Clear(Option_t *o)
92{
93
94 MHGausEvents::Clear();
95 fSaturated = 0;
96}
97
98void MHCalibrationPix::Reset()
99{
100
101 MHGausEvents::Reset();
102 fSaturated = 0;
103}
104
105
106// -----------------------------------------------------------------------------
107//
108// Bypasses the Gauss fit by taking mean and RMS from the histogram
109//
110// Errors are determined in the following way:
111// MeanErr = RMS / Sqrt(entries)
112// SigmaErr = RMS / (2.*Sqrt(entries) )
113//
114void MHCalibrationPix::BypassFit()
115{
116
117 const Stat_t entries = fHGausHist.GetEntries();
118
119 if (entries <= 0.)
120 {
121 *fLog << warn << GetDescriptor()
122 << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << GetName() << endl;
123 return;
124 }
125
126 fMean = fHGausHist.GetMean();
127 fMeanErr = fHGausHist.GetRMS() / TMath::Sqrt(entries);
128 fSigma = fHGausHist.GetRMS() ;
129 fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
130}
131
132// -------------------------------------------------------------------------------
133//
134// Return the number of "blackout" events, which are events with values higher
135// than fBlackoutLimit sigmas from the mean
136//
137//
138const Double_t MHCalibrationPix::GetBlackout() const
139{
140
141 if ((fMean == 0.) && (fSigma == 0.))
142 return -1.;
143
144 const Int_t first = fHGausHist.GetXaxis()->GetFirst();
145 const Int_t last = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma);
146
147 if (first >= last)
148 return 0.;
149
150 return fHGausHist.Integral(first, last);
151}
152
153
154// -------------------------------------------------------------------------------
155//
156// Return the number of "pickup" events, which are events with values higher
157// than fPickupLimit sigmas from the mean
158//
159//
160const Double_t MHCalibrationPix::GetPickup() const
161{
162 if (!IsValid())
163 return -1.;
164
165 const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma);
166 const Int_t last = fHGausHist.GetXaxis()->GetLast();
167
168 if (first >= last)
169 return 0.;
170
171 return fHGausHist.Integral(first, last);
172}
173
174// -----------------------------------------------------------------------------
175//
176// If flag IsGausFitOK() is set (histogram already successfully fitted),
177// returns kTRUE
178//
179// If both fMean and fSigma are still zero, call FitGaus()
180//
181// Repeats the Gauss fit in a smaller range, defined by:
182//
183// min = GetMean() - fBlackoutLimit * GetSigma();
184// max = GetMean() + fPickupLimit * GetSigma();
185//
186// The fit results are retrieved and stored in class-own variables.
187//
188// A flag IsGausFitOK() is set according to whether the fit probability
189// is smaller or bigger than fProbLimit, whether the NDF is bigger than
190// fNDFLimit and whether results are NaNs.
191//
192Bool_t MHCalibrationPix::RepeatFit(const Option_t *option)
193{
194
195 if (IsGausFitOK())
196 return kTRUE;
197
198 if (!IsValid())
199 return FitGaus();
200
201 //
202 // Get new fitting ranges
203 //
204 Axis_t rmin = GetMean() - fBlackoutLimit * GetSigma();
205 Axis_t rmax = GetMean() + fPickupLimit * GetSigma();
206
207 Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
208 Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ;
209
210 fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax);
211
212 fHGausHist.Fit(fFGausFit,option);
213
214 fMean = fFGausFit->GetParameter(1);
215 fSigma = fFGausFit->GetParameter(2);
216 fMeanErr = fFGausFit->GetParError(1) ;
217 fSigmaErr = fFGausFit->GetParError(2) ;
218 fProb = fFGausFit->GetProb() ;
219
220 //
221 // The fit result is accepted under condition:
222 // 1) The results are not nan's
223 // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
224 // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
225 //
226 if ( !TMath::Finite( fMean )
227 || !TMath::Finite( fMeanErr )
228 || !TMath::Finite( fProb )
229 || !TMath::Finite( fSigma )
230 || !TMath::Finite( fSigmaErr )
231 || fFGausFit->GetNDF() < GetNDFLimit()
232 || fProb < GetProbLimit() )
233 return kFALSE;
234
235 SetGausFitOK(kTRUE);
236 return kTRUE;
237
238}
239
Note: See TracBrowser for help on using the repository browser.