source: branches/Mars_use_drstimefiles/mextralgo/MExtralgoSpline.h@ 18558

Last change on this file since 18558 was 17879, checked in by tbretz, 10 years ago
Added direct access to integral
File size: 10.9 KB
Line 
1#ifndef MARS_MExtralgoSpline
2#define MARS_MExtralgoSpline
3
4#ifndef ROOT_TMath
5#include <TMath.h>
6#endif
7
8class MArrayF;
9class TComplex;
10
11class MExtralgoSpline
12{
13private:
14 enum
15 {
16 kIntegral = BIT(0),
17 kTimeRel = BIT(1),
18 kMaximum = BIT(2),
19 kDynWidth = BIT(3),
20 kFixedWidth = BIT(4),
21 };
22public:
23 enum ExtractionType_t
24 {
25 // For backward compatibility
26 kAmplitudeAbs = 0, // Height of maximum, absolute height leading edge
27 kAmplitudeRel = kTimeRel, // Height of maximum, relative height leading edge
28 kAmplitude = kMaximum, // Position and height of maximum
29 kIntegralAbs = kIntegral, // Integral, absolute height leading edge
30 kIntegralRel = kIntegral|kTimeRel, // Integral, relative height leading edge
31 kIntegralDyn = kTimeRel|kDynWidth, // Integrate between leading edge and falling edge
32 kIntegralFixed = kTimeRel|kFixedWidth, // Integrate between leading edge minus fRiseTime and plus fFallTime
33 };
34
35private:
36 ExtractionType_t fExtractionType;
37
38private:
39 //Bool_t fIsOwner; // Owner of derivatives....
40
41 // Input
42 Float_t const *fVal;
43 const Int_t fNum;
44
45 Float_t *fDer1;
46 Float_t *fDer2;
47
48 Float_t fRiseTime;
49 Float_t fFallTime;
50
51 Float_t fHeightTm;
52
53 // Result
54 Float_t fTime;
55 Float_t fTimeDev;
56 Float_t fWidth;
57 Float_t fWidthDev;
58 Float_t fSignal;
59 Float_t fSignalDev;
60 Float_t fHeight;
61
62 Double_t ReMul(const TComplex &c1, const TComplex &th) const;
63
64 inline Float_t Eval(Float_t val, Float_t a, Float_t deriv) const
65 {
66 return a*val + (a*a*a-a)*deriv;
67 }
68
69 // Evaluate value of spline in the interval i with x=[0;1[
70 inline Float_t Eval(const Int_t i, const Float_t x) const
71 {
72 // Eval(i,x) = (fDer2[i+1]-fDer2[i])*x*x*x + 3*fDer2[i]*x*x +
73 // (fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1])*x + fVal[i];
74
75 // x := [0; 1[
76 return Eval(fVal[i], 1-x, fDer2[i]) + Eval(fVal[i+1], x, fDer2[i+1]);
77 }
78
79 // Evaluate first derivative of spline in the interval i with x=[0;1[
80 inline Double_t EvalDeriv1(const Int_t i, const Float_t x) const
81 {
82 // x := [0; 1[
83 const Double_t difval = fVal[i+1]-fVal[i];
84 const Double_t difder = fDer2[i+1]-fDer2[i];
85
86 //return 3*difder*x*x + 6*fDer2[i]*x - 2*fDer2[i] - fDer2[i+1] + difval;
87 return 3*difder*x*x + (6*x - 2)*fDer2[i] - fDer2[i+1] + difval;
88 }
89
90 // Evaluate second derivative of spline in the interval i with x=[0;1[
91 inline Double_t EvalDeriv2(const Int_t i, const Float_t x) const
92 {
93 // x := [0; 1[
94 return 6*(fDer2[i+1]*x + fDer2[i]*(1-x));
95 }
96
97 Int_t SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) const;
98 Double_t FindYdn(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
99 Double_t FindYup(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
100 //Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
101
102 Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const;
103/*
104 inline void EvalDerivEq0(const Int_t i, Float_t &rc1, Float_t &rc2) const
105 {
106 // --- ORIGINAL CODE ---
107 Double_t sumder = fDer2[i]+fDer2[i+1];
108 Double_t difder = fDer2[i]-fDer2[i+1];
109
110 Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
111 Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
112 Double_t sqt3 = sqrt(3*sqt1 + 3*sqt2);
113 Double_t denom = -3*(fDer2[i+1]-fDer2[i]);
114
115 rc1 = (3*fDer2[i] + sqt3)/denom;
116 rc2 = (3*fDer2[i] - sqt3)/denom;
117
118 // --- NEW CODE ---
119 Double_t sumder = fDer2[i]+fDer2[i+1];
120 Double_t difder = fDer2[i]-fDer2[i+1];
121
122 Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
123 Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
124 Double_t sqt3 = sqt1+sqt2<0 ? 0 : sqrt((sqt1 + sqt2)/3);
125
126 rc1 = (fDer2[i] + sqt3)/difder;
127 rc2 = (fDer2[i] - sqt3)/difder;
128 }*/
129
130 // Calculate the "Stammfunktion" of the Eval-function
131 inline Double_t EvalPrimitive(Int_t i, Float_t x) const
132 {
133 Align(i, x);
134
135 if (x==0)
136 return -fDer2[i]/4;
137
138 if (x==1)
139 return (fVal[i+1] + fVal[i])/2 - fDer2[i+1]/4 - fDer2[i]/2;
140
141 const Double_t x2 = x*x;
142 const Double_t x4 = x2*x2;
143 const Double_t x1 = 1-x;
144 const Double_t x14 = x1*x1*x1*x1;
145
146 return x2*fVal[i+1]/2 + (x4/2-x2)*fDer2[i+1]/2 + (x-x2/2)*fVal[i] + (x2/2-x-x14/4)*fDer2[i];
147
148 }
149
150 inline void Align(Int_t &i, Float_t &x) const
151 {
152 if (i<0)
153 {
154 x += i;
155 i=0;
156 }
157 if (i>=fNum-1)
158 {
159 x += i-(fNum-2);
160 i=fNum-2;
161 }
162 }
163
164 // Calculate the intgeral of the Eval-function in
165 // bin i from 0 <= a < b < 1
166 inline Double_t EvalInteg(Int_t i, Float_t a, Float_t b) const
167 {
168 return EvalPrimitive(i, b)-EvalPrimitive(i, a);
169 }
170
171 // Identical to EvalInteg(i, 0, 1) but much faster
172 // Be carefull: NO RANGECHECK!
173 inline Double_t EvalInteg(Int_t i) const
174 {
175 return (fVal[i+1] + fVal[i])/2 - (fDer2[i+1] + fDer2[i])/4;
176 }
177
178 // Identical to sum of EvalInteg(i, 0, 1) for i=a to i=b-1,
179 // but much faster
180 // It is identical to EvalInteg(fVal[a], fVal[b])
181 // Be carefull: NO RANGECHECK!
182 inline Double_t EvalInteg(Int_t a, Int_t b) const
183 {
184 /*
185 Double_t sum = 0;
186 for (int i=a; i<b; i++)
187 sum += EvalInteg(i);
188
189 return sum;
190 */
191
192 if (a==b)
193 return 0;
194
195 Double_t sum=0;
196 for (const Float_t *ptr=fDer2+a+1; ptr<fDer2+b; ptr++)
197 sum -= *ptr;
198
199 sum -= (fDer2[a]+fDer2[b])/2;
200
201 sum /= 2;
202
203 for (const Float_t *ptr=fVal+a+1; ptr<fVal+b; ptr++)
204 sum += *ptr;
205
206 sum += (fVal[a]+fVal[b])/2;
207
208 return sum;
209 }
210
211 // Calculate the intgeral of the Eval-function betwen x0 and x1
212 inline Double_t EvalInteg(Float_t x0, Float_t x1) const
213 {
214 // RANGE CHECK MISSING!
215
216 const Int_t min = TMath::CeilNint(x0);
217 const Int_t max = TMath::FloorNint(x1);
218
219 // This happens if x0 and x1 are in the same interval
220 if (min>max)
221 return EvalInteg(max, x0-max, x1-max);
222
223 // Sum complete intervals
224 Double_t sum = EvalInteg(min, max);
225
226 // Sum the incomplete intervals at the beginning and end
227 sum += EvalInteg(min-1, 1-(min-x0), 1);
228 sum += EvalInteg(max, 0, x1-max);
229
230 // return result
231 return sum;
232 }
233
234 // We search for the maximum from x=i-1 to x=i+1
235 // (Remeber: i corresponds to the value in bin i, i+1 to the
236 // next bin and i-1 to the last bin)
237 inline void GetMaxAroundI(Int_t i, Float_t &xmax, Float_t &ymax) const
238 {
239 Float_t xmax1=0, xmax2=0;
240 Float_t ymax1=0, ymax2=0;
241
242 Bool_t rc1 = i>0 && GetMax(i-1, xmax1, ymax1);
243 Bool_t rc2 = i<fNum-1 && GetMax(i, xmax2, ymax2);
244
245 // In case the medium bin is the first or last bin
246 // take the lower or upper edge of the region into account.
247 if (i==0)
248 {
249 xmax1 = 0;
250 ymax1 = fVal[0];
251 rc1 = kTRUE;
252 }
253 if (i>=fNum-1)
254 {
255 xmax2 = fNum-1;
256 ymax2 = fVal[fNum-1];
257 rc2 = kTRUE;
258 }
259
260 // Take a default in case no maximum is found
261 // FIXME: Check THIS!!!
262 xmax=i;
263 ymax=fVal[i];
264
265 if (rc1)
266 {
267 ymax = ymax1;
268 xmax = xmax1;
269 }
270 else
271 if (rc2)
272 {
273 ymax = ymax2;
274 xmax = xmax2;
275 }
276
277 if (rc2 && ymax2>ymax)
278 {
279 ymax = ymax2;
280 xmax = xmax2;
281 }
282 }
283
284 inline Bool_t GetMax(Int_t i, Float_t &xmax, Float_t &ymax, Float_t min=0, Float_t max=1) const
285 {
286 // Find analytical maximum in the bin i in the interval [min,max[
287
288 Double_t x1=-1; // This initialisation should not really be
289 Double_t x2=-1; // necessary but makes valgriund happy.
290
291 if (!EvalDerivEq0(i, x1, x2))
292 return kFALSE;
293
294 const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(i, x1)<0;
295 const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(i, x2)<0;
296
297 if (!ismax1 && !ismax2)
298 return kFALSE;
299
300 if (ismax1 && !ismax2)
301 {
302 xmax = i+x1;
303 ymax = Eval(i, x1);
304 return kTRUE;
305 }
306
307 if (!ismax1 && ismax2)
308 {
309 xmax = i+x2;
310 ymax = Eval(i, x2);
311 return kTRUE;
312 }
313
314 // Somehting must be wrong...
315 return kFALSE;
316 }
317
318 void InitDerivatives() const;
319 Float_t CalcIntegral(Float_t start) const;
320 Float_t CalcIntegral(Float_t beg, Float_t width) const;
321
322public:
323 MExtralgoSpline(const Float_t *val, Int_t n, Float_t *der1, Float_t *der2)
324 : fExtractionType(kIntegralRel), fVal(val), fNum(n), fDer1(der1), fDer2(der2), fHeightTm(0.5), fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
325 {
326 InitDerivatives();
327 }
328
329 void SetRiseFallTime(Float_t rise, Float_t fall) { fRiseTime=rise; fFallTime=fall; }
330 void SetExtractionType(ExtractionType_t typ) { fExtractionType = typ; }
331 void SetHeightTm(Float_t h) { fHeightTm = h; }
332
333 Float_t GetTime() const { return fTime; }
334 Float_t GetWidth() const { return fWidth; }
335 Float_t GetSignal() const { return fSignal; }
336 Float_t GetHeight() const { return fHeight; }
337
338 Float_t GetTimeDev() const { return fTimeDev; }
339 Float_t GetWidthDev() const { return fWidthDev; }
340 Float_t GetSignalDev() const { return fSignalDev; }
341
342 void GetSignal(Float_t &sig, Float_t &dsig) const { sig=fSignal; dsig=fSignalDev; }
343 void GetWidth(Float_t &sig, Float_t &dsig) const { sig=fWidth; dsig=fWidthDev; }
344 void GetTime(Float_t &sig, Float_t &dsig) const { sig=fTime; dsig=fTimeDev; }
345
346 Float_t ExtractNoise(/*Int_t iter*/);
347 void Extract(Int_t maxpos, Bool_t width=kFALSE);
348
349 Float_t EvalAt(const Float_t x) const;
350 Float_t Deriv1(const Float_t x) const;
351
352 Double_t SearchYdn(Double_t maxpos, Double_t y) const;
353 Double_t SearchYup(Double_t maxpos, Double_t y) const;
354
355 Double_t SearchYdn(Float_t y) const { return SearchYdn(fNum, y); }
356 Double_t SearchYup(Float_t y) const { return SearchYup(0, y); }
357
358 MArrayF GetIntegral(bool norm=false) const;
359
360 Float_t GetIntegral(Float_t beg, Float_t end) const
361 {
362 return CalcIntegral(beg, end-beg);
363 }
364};
365
366inline Float_t MExtralgoSpline::EvalAt(const Float_t x) const
367{
368 Int_t i = TMath::FloorNint(x);
369 Float_t f = x-i;
370
371 Align(i, f);
372
373 return Eval(i, f);
374}
375
376inline Float_t MExtralgoSpline::Deriv1(const Float_t x) const
377{
378 Int_t i = TMath::FloorNint(x);
379 Float_t f = x-i;
380
381 Align(i, f);
382
383 return EvalDeriv1(i, f);
384}
385
386#endif
Note: See TracBrowser for help on using the repository browser.