source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationPix.cc@ 6515

Last change on this file since 6515 was 5688, checked in by gaug, 20 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 <TGraph.h>
47
48#include "MLog.h"
49#include "MLogManip.h"
50
51ClassImp(MHCalibrationPix);
52
53using namespace std;
54
55const Float_t MHCalibrationPix::fgBlackoutLimit = 5.;
56const Float_t MHCalibrationPix::fgPickupLimit = 5.;
57// --------------------------------------------------------------------------
58//
59// Default Constructor.
60// Sets:
61// - the default number for fPickupLimit (fgPickupLimit)
62// - the default number for fBlackoutLimit (fgBlackoutLimit)
63//
64// Initializes:
65// - all variables to 0.
66//
67MHCalibrationPix::MHCalibrationPix(const char *name, const char *title)
68{
69
70 fName = name ? name : "MHCalibrationPix";
71 fTitle = title ? title : "Calibration histogram events";
72
73 Clear();
74
75 SetBlackoutLimit();
76 SetPickupLimit();
77}
78
79// --------------------------------------------------------------------------
80//
81// Default Clear(), can be overloaded.
82//
83// Sets:
84// - all other pointers to NULL
85// - all variables to 0., except fPixId to -1
86// - all flags to kFALSE
87//
88// - all pointers
89//
90void MHCalibrationPix::Clear(Option_t *o)
91{
92
93 MHGausEvents::Clear();
94 fSaturated = 0;
95}
96
97void MHCalibrationPix::Reset()
98{
99
100 MHGausEvents::Reset();
101 fSaturated = 0;
102}
103
104
105// -----------------------------------------------------------------------------
106//
107// Bypasses the Gauss fit by taking mean and RMS from the histogram
108//
109// Errors are determined in the following way:
110// MeanErr = RMS / Sqrt(entries)
111// SigmaErr = RMS / (2.*Sqrt(entries) )
112//
113void MHCalibrationPix::BypassFit()
114{
115
116 const Stat_t entries = fHGausHist.GetEntries();
117
118 if (entries <= 0.)
119 {
120 *fLog << warn << GetDescriptor()
121 << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << GetName() << endl;
122 return;
123 }
124
125 fMean = fHGausHist.GetMean();
126 fMeanErr = fHGausHist.GetRMS() / TMath::Sqrt(entries);
127 fSigma = fHGausHist.GetRMS() ;
128 fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
129}
130
131// -------------------------------------------------------------------------------
132//
133// Return the number of "blackout" events, which are events with values higher
134// than fBlackoutLimit sigmas from the mean
135//
136//
137const Double_t MHCalibrationPix::GetBlackout() const
138{
139
140 if ((fMean == 0.) && (fSigma == 0.))
141 return -1.;
142
143 const Int_t first = fHGausHist.GetXaxis()->GetFirst();
144 const Int_t last = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma);
145
146 if (first >= last)
147 return 0.;
148
149 return fHGausHist.Integral(first, last);
150}
151
152
153// -------------------------------------------------------------------------------
154//
155// Return the number of "pickup" events, which are events with values higher
156// than fPickupLimit sigmas from the mean
157//
158//
159const Double_t MHCalibrationPix::GetPickup() const
160{
161
162 if ((fMean == 0.) && (fSigma == 0.))
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 ((fMean == 0.) && (fSigma == 0.))
199 return FitGaus();
200
201 //
202 // Get new fitting ranges
203 //
204 Axis_t rmin = fMean - fBlackoutLimit * fSigma;
205 Axis_t rmax = fMean + fPickupLimit * fSigma;
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::IsNaN ( fMean )
227 || TMath::IsNaN ( fMeanErr )
228 || TMath::IsNaN ( fProb )
229 || TMath::IsNaN ( fSigma )
230 || TMath::IsNaN ( fSigmaErr )
231 || fFGausFit->GetNDF() < fNDFLimit
232 || fProb < fProbLimit )
233 return kFALSE;
234
235 SetGausFitOK(kTRUE);
236 return kTRUE;
237
238}
239
Note: See TracBrowser for help on using the repository browser.