source: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc@ 9223

Last change on this file since 9223 was 9212, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 8.6 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 analyzing 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! Author(s): Thomas Bretz <mailto:tbretz@astro.uni-wuerzbrug.de>
18! Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2002-2008
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MExtralgoSpline
28//
29// Fast Spline extractor using a cubic spline algorithm, adapted from
30// Numerical Recipes in C++, 2nd edition, pp. 116-119.
31//
32// The coefficients "ya" are here denoted as "fVal" corresponding to
33// the FADC value subtracted by the clock-noise corrected pedestal.
34//
35// The coefficients "y2a" get immediately divided 6. and are called here
36// fDer2 although they are now not exactly the second derivative
37// coefficients any more.
38//
39// The calculation of the cubic-spline interpolated value "y" on a point
40// "x" along the FADC-slices axis becomes: EvalAt(x)
41//
42// The coefficients fDer2 are calculated with the simplified
43// algorithm in InitDerivatives.
44//
45// This algorithm takes advantage of the fact that the x-values are all
46// separated by exactly 1 which simplifies the Numerical Recipes algorithm.
47// (Note that the variables fDer are not real first derivative coefficients.)
48//
49// Note, this spline is not optimized to be evaluated many many times, but
50// it is optimized to be initialized very fast with new values again and
51// again.
52//
53//////////////////////////////////////////////////////////////////////////////
54#include "MExtralgoSpline.h"
55
56#include <TRandom.h>
57
58#include "../mbase/MMath.h"
59#include "../mbase/MArrayF.h"
60
61using namespace std;
62
63// --------------------------------------------------------------------------
64//
65// Calculate the first and second derivative for the splie.
66//
67// The coefficients are calculated such that
68// 1) fVal[i] = Eval(i, 0)
69// 2) Eval(i-1, 1)==Eval(i, 0)
70//
71// In other words: The values with the index i describe the spline
72// between fVal[i] and fVal[i+1]
73//
74void MExtralgoSpline::InitDerivatives() const
75{
76 if (fNum<2)
77 return;
78
79 // Look up table for coefficients
80 static MArrayF lut;
81
82 // If the lut is not yet large enough: resize and reclaculate
83 if (fNum>(Int_t)lut.GetSize())
84 {
85 lut.Set(fNum);
86
87 lut[0] = 0.;
88 for (Int_t i=1; i<fNum-1; i++)
89 lut[i] = -1.0/(lut[i-1] + 4);
90 }
91
92 // Calculate the coefficients used to get reproduce the first and
93 // second derivative.
94 fDer1[0] = 0.;
95 for (Int_t i=1; i<fNum-1; i++)
96 {
97 const Float_t d1 = fVal[i+1] - 2*fVal[i] + fVal[i-1];
98 fDer1[i] = (fDer1[i-1]-d1)*lut[i];
99 }
100
101 fDer2[fNum-1] = 0.;
102 for (Int_t k=fNum-2; k>=0; k--)
103 fDer2[k] = lut[k]*fDer2[k+1] + fDer1[k];
104}
105
106// --------------------------------------------------------------------------
107//
108// Return the two results x1 and x2 of f'(x)=0 for the third order
109// polynomial (spline) in the interval i. Return the number of results.
110// (0 if the fist derivative does not have a null-point)
111//
112Int_t MExtralgoSpline::EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const
113{
114 const Double_t difder = fDer2[i+1]-fDer2[i];
115 const Double_t difval = fVal[i+1] -fVal[i];
116
117 return MMath::SolvePol2(3*difder, 6*fDer2[i], difval-2*fDer2[i]-fDer2[i+1], x1, x2);
118}
119
120// --------------------------------------------------------------------------
121//
122// Returns the highest x value in [min;max[ at which the spline in
123// the bin i is equal to y
124//
125// min and max are defined to be [0;1]
126//
127// The default for min is 0, the default for max is 1
128// The defaule for y is 0
129//
130Double_t MExtralgoSpline::FindY(Int_t i, Bool_t downwards, Double_t y, Double_t min, Double_t max) const
131{
132 // y = a*x^3 + b*x^2 + c*x + d'
133 // 0 = a*x^3 + b*x^2 + c*x + d' - y
134
135 // Calculate coefficients
136 const Double_t a = fDer2[i+1]-fDer2[i];
137 const Double_t b = 3*fDer2[i];
138 const Double_t c = fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1];
139 const Double_t d = fVal[i] - y;
140
141 // If the first derivative is nowhere==0 and it is increasing
142 // in one point, and the value we search is outside of the
143 // y-interval... it cannot be there
144 // if (c>0 && (d>0 || fVal[i+1]<y) && b*b<3*c*a)
145 // return -2;
146
147 Double_t x1, x2, x3;
148 const Int_t rc = MMath::SolvePol3(a, b, c, d, x1, x2, x3);
149
150 if (downwards==kTRUE)
151 {
152 Double_t x = -1;
153
154 if (rc>0 && x1>=min && x1<max && x1>x)
155 x = x1;
156 if (rc>1 && x2>=min && x2<max && x2>x)
157 x = x2;
158 if (rc>2 && x3>=min && x3<max && x3>x)
159 x = x3;
160
161 return x<0 ? -2 : x+i;
162 }
163 else
164 {
165 Double_t x = 2;
166
167 if (rc>0 && x1>min && x1<=max && x1<x)
168 x = x1;
169 if (rc>1 && x2>min && x2<=max && x2<x)
170 x = x2;
171 if (rc>2 && x3>min && x3<=max && x3<x)
172 x = x3;
173
174 return x>1 ? -2 : x+i;
175 }
176
177 return -2;
178}
179
180// --------------------------------------------------------------------------
181//
182// Search analytically downward for the value y of the spline, starting
183// at x, until x==0. If y is not found -2 is returned.
184//
185Double_t MExtralgoSpline::SearchYdn(Float_t x, Float_t y) const
186{
187 if (x>=fNum-1)
188 x = fNum-1.0001;
189
190 Int_t i = TMath::FloorNint(x);
191 if (i<0)
192 return -2;
193
194 Double_t rc = FindY(i, kTRUE, y, 0, x-i);
195 while (--i>=0 && rc<0)
196 rc = FindY(i, kTRUE, y);
197
198 return rc;
199}
200
201Double_t MExtralgoSpline::SearchYup(Float_t x, Float_t y) const
202{
203 if (x<0)
204 x = 0.0001;
205
206 Int_t i = TMath::FloorNint(x);
207 if (i>fNum-2)
208 return -2;
209
210 Double_t rc = FindY(i, kFALSE, y, x-i, 1.);
211 while (++i<fNum-1 && rc<0)
212 rc = FindY(i, kFALSE, y);
213
214 return rc;
215}
216
217// --------------------------------------------------------------------------
218//
219// Do a range check an then calculate the integral from start-fRiseTime
220// to start+fFallTime. An extrapolation of 0.5 slices is allowed.
221//
222Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
223{
224 // In the future we will calculate the intgeral analytically.
225 // It has been tested that it gives identical results within
226 // acceptable differences.
227
228 // We allow extrapolation of 1/2 slice.
229 const Float_t min = fRiseTime; //-0.5+fRiseTime;
230 const Float_t max = fNum-1-fFallTime; //fNum-0.5+fFallTime;
231
232 if (pos<min)
233 pos = min;
234 if (pos>max)
235 pos = max;
236
237 return EvalInteg(pos-fRiseTime, pos+fFallTime);
238}
239
240Float_t MExtralgoSpline::ExtractNoise()
241{
242 if (fNum<5)
243 return 0;
244
245 if (fExtractionType == kAmplitude)
246 {
247 const Int_t pos = gRandom->Integer(fNum-1);
248 const Float_t nsx = gRandom->Uniform();
249 return Eval(pos, nsx);
250 }
251 else
252 {
253 const Float_t pos = gRandom->Uniform(fNum-1-fRiseTime-fFallTime)+fRiseTime;
254 return CalcIntegral(pos);
255 }
256}
257
258void MExtralgoSpline::Extract(Int_t maxbin, Bool_t width)
259{
260 fSignal = 0;
261 fTime = 0;
262 fWidth = 0;
263 fSignalDev = -1;
264 fTimeDev = -1;
265 fWidthDev = -1;
266
267 if (fNum<2)
268 return;
269
270 Float_t maxpos;
271 // FIXME: Check the default if no maximum found!!!
272 GetMaxAroundI(maxbin, maxpos, fHeight);
273
274 // --- End NEW ---
275
276 if (fExtractionType == kAmplitude)
277 {
278 fTime = maxpos;
279 fTimeDev = 0;
280 fSignal = fHeight;
281 fSignalDev = 0; // means: is valid
282 return;
283 }
284
285 fSignal = CalcIntegral(maxpos);
286 fSignalDev = 0; // means: is valid
287
288 if (fExtractionType==kIntegralRel && fHeightTm<0)
289 {
290 fTime = maxpos;
291 fTimeDev = 0;
292 return;
293 }
294
295 const Float_t h = fExtractionType==kIntegralAbs ? fHeightTm : fHeight*fHeightTm;
296
297 // Search downwards for fHeight/2
298 // By doing also a search upwards we could extract the pulse width
299 fTime = SearchYdn(maxpos, h);
300 fTimeDev = 0;
301 if (width)
302 {
303 fWidth = SearchYup(maxpos, h)-fTime;
304 fWidthDev = 0;
305 }
306}
Note: See TracBrowser for help on using the repository browser.