source: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h@ 9223

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