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

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