source: trunk/Mars/mextralgo/MExtralgoSpline.cc@ 12392

Last change on this file since 12392 was 9226, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 9.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 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-2009
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// Solve the polynomial
123//
124// y = a*x^3 + b*x^2 + c*x + d'
125// 0 = a*x^3 + b*x^2 + c*x + d' - y
126//
127// to find y in the i-th bin. Return the result as x1, x2, x3 and the return
128// code from MMath::SolvPol3.
129//
130Int_t MExtralgoSpline::SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) 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 return MMath::SolvePol3(a, b, c, d, x1, x2, x3);
148}
149
150// --------------------------------------------------------------------------
151//
152// Returns the highest x value in [min;max[ at which the spline in
153// the bin i is equal to y
154//
155// min and max must be in the range [0;1]
156//
157// The default for min is 0, the default for max is 1
158// The default for y is 0
159//
160Double_t MExtralgoSpline::FindYdn(Int_t i, Double_t y, Double_t min, Double_t max) const
161{
162 Double_t x1, x2, x3;
163 const Int_t rc = SolvePol3(i, y, x1, x2, x3);
164
165 Double_t x = -1;
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<0 ? -2 : x+i;
175}
176
177// --------------------------------------------------------------------------
178//
179// Returns the lowest x value in [min;max[ at which the spline in
180// the bin i is equal to y
181//
182// min and max must be in the range [0;1]
183//
184// The default for min is 0, the default for max is 1
185// The default for y is 0
186//
187Double_t MExtralgoSpline::FindYup(Int_t i, Double_t y, Double_t min, Double_t max) const
188{
189 Double_t x1, x2, x3;
190 const Int_t rc = SolvePol3(i, y, x1, x2, x3);
191
192 Double_t x = 2;
193
194 if (rc>0 && x1>min && x1<=max && x1<x)
195 x = x1;
196 if (rc>1 && x2>min && x2<=max && x2<x)
197 x = x2;
198 if (rc>2 && x3>min && x3<=max && x3<x)
199 x = x3;
200
201 return x>1 ? -2 : x+i;
202}
203
204// --------------------------------------------------------------------------
205//
206// Search analytically downward for the value y of the spline, starting
207// at x, until x==0. If y is not found or out of range -2 is returned.
208//
209Double_t MExtralgoSpline::SearchYdn(Float_t x, Float_t y) const
210{
211 if (x>=fNum-1)
212 x = fNum-1.0001;
213
214 Int_t i = TMath::FloorNint(x);
215 if (i<0)
216 return -2;
217
218 Double_t rc = FindYdn(i, y, 0, x-i);
219 while (--i>=0 && rc<0)
220 rc = FindYdn(i, y);
221
222 return rc;
223}
224
225// --------------------------------------------------------------------------
226//
227// Search analytically upwards for the value y of the spline, starting
228// at x, until x==fNum-1. If y is not found or out of range -2 is returned.
229//
230Double_t MExtralgoSpline::SearchYup(Float_t x, Float_t y) const
231{
232 if (x<0)
233 x = 0.0001;
234
235 Int_t i = TMath::FloorNint(x);
236 if (i>fNum-2)
237 return -2;
238
239 Double_t rc = FindYup(i, y, x-i, 1.);
240 while (++i<fNum-1 && rc<0)
241 rc = FindYup(i, y);
242
243 return rc;
244}
245
246// --------------------------------------------------------------------------
247//
248// Do a range check an then calculate the integral from start-fRiseTime
249// to start+fFallTime. An extrapolation of 0.5 slices is allowed.
250//
251Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
252{
253 // We allow extrapolation of 1/2 slice.
254 const Float_t min = fRiseTime; //-0.5+fRiseTime;
255 const Float_t max = fNum-1-fFallTime; //fNum-0.5+fFallTime;
256
257 if (pos<min)
258 pos = min;
259 if (pos>max)
260 pos = max;
261
262 return EvalInteg(pos-fRiseTime, pos+fFallTime);
263}
264
265MArrayF MExtralgoSpline::GetIntegral(bool norm) const
266{
267 MArrayF val(fNum);
268
269 //val[0] = 0;
270
271 Double_t integ = 0;
272 for (int i=0; i<fNum-1; i++)
273 {
274 integ += EvalInteg(i);
275
276 val[i+1] = integ;
277 }
278
279 if (norm)
280 for (int i=0; i<fNum-1; i++)
281 val[i+1] /= val[fNum-1];
282
283 return val;
284}
285
286Float_t MExtralgoSpline::ExtractNoise()
287{
288 if (fNum<5)
289 return 0;
290
291 if (fExtractionType == kAmplitude)
292 {
293 const Int_t pos = gRandom->Integer(fNum-1);
294 const Float_t nsx = gRandom->Uniform();
295 return Eval(pos, nsx);
296 }
297 else
298 {
299 const Float_t pos = gRandom->Uniform(fNum-1-fRiseTime-fFallTime)+fRiseTime;
300 return CalcIntegral(pos);
301 }
302}
303
304void MExtralgoSpline::Extract(Int_t maxbin, Bool_t width)
305{
306 fSignal = 0;
307 fTime = 0;
308 fWidth = 0;
309 fSignalDev = -1;
310 fTimeDev = -1;
311 fWidthDev = -1;
312
313 if (fNum<2)
314 return;
315
316 Float_t maxpos;
317 // FIXME: Check the default if no maximum found!!!
318 GetMaxAroundI(maxbin, maxpos, fHeight);
319
320 // --- End NEW ---
321
322 if (fExtractionType == kAmplitude)
323 {
324 fTime = maxpos;
325 fTimeDev = 0;
326 fSignal = fHeight;
327 fSignalDev = 0; // means: is valid
328 return;
329 }
330
331 fSignal = CalcIntegral(maxpos);
332 fSignalDev = 0; // means: is valid
333
334 if (fExtractionType==kIntegralRel && fHeightTm<0)
335 {
336 fTime = maxpos;
337 fTimeDev = 0;
338 return;
339 }
340
341 const Float_t h = fExtractionType==kIntegralAbs ? fHeightTm : fHeight*fHeightTm;
342
343 // Search downwards for fHeight/2
344 // By doing also a search upwards we could extract the pulse width
345 fTime = SearchYdn(maxpos, h);
346 fTimeDev = 0;
347 if (width)
348 {
349 fWidth = SearchYup(maxpos, h)-fTime;
350 fWidthDev = 0;
351 }
352}
Note: See TracBrowser for help on using the repository browser.