source: trunk/Mars/mextralgo/MExtralgoSpline.h@ 13736

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