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

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