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

Last change on this file since 8339 was 8316, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 12.5 KB
Line 
1#ifndef MARS_MExtralgoSpline
2#define MARS_MExtralgoSpline
3
4#ifndef ROOT_TROOT
5#include <TROOT.h>
6#endif
7
8#include <iostream>
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// Float_t fResolution;
35
36 // Result
37 Float_t fTime;
38 Float_t fTimeDev;
39 Float_t fWidth;
40 Float_t fWidthDev;
41 Float_t fSignal;
42 Float_t fSignalDev;
43 Float_t fHeight;
44
45 Double_t ReMul(const TComplex &c1, const TComplex &th) const;
46
47 inline Float_t Eval(Float_t val, Float_t a, Float_t deriv) const
48 {
49 return a*val + (a*a*a-a)*deriv;
50 }
51
52 // Evaluate value of spline in the interval i with x=[0;1[
53 inline Float_t Eval(const Int_t i, const Float_t x) const
54 {
55 // Eval(i,x) = (fDer2[i+1]-fDer2[i])*x*x*x + 3*fDer2[i]*x*x +
56 // (fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1])*x + fVal[i];
57
58 // x := [0; 1[
59 return Eval(fVal[i], 1-x, fDer2[i]) + Eval(fVal[i+1], x, fDer2[i+1]);
60 }
61
62 /*
63 inline Float_t EvalAt(const Float_t x) const
64 {
65 Int_t i = TMath::FloorNint(x);
66
67 // handle under- and overflow of the array-range by extrapolation
68 if (i<0)
69 i=0;
70 if (i>fNum-2)
71 i = fNum-2;
72
73 return Eval(i, x-i);
74 }
75 */
76
77 // Evaluate first derivative of spline in the interval i with x=[0;1[
78 inline Double_t EvalDeriv1(const Float_t x, const Int_t i) const
79 {
80 // x := [0; 1[
81 const Double_t difval = fVal[i+1]-fVal[i];
82 const Double_t difder = fDer2[i+1]-fDer2[i];
83
84 return 3*difder*x*x + 6*fDer2[i]*x - 2*fDer2[i] - fDer2[i+1] + difval;
85 }
86
87 // Evaluate second derivative of spline in the interval i with x=[0;1[
88 inline Double_t EvalDeriv2(const Float_t x, const Int_t i) const
89 {
90 // x := [0; 1[
91 return 6*(fDer2[i+1]*x + fDer2[i]*(1-x));
92 }
93
94 Double_t FindY(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
95 Double_t SearchY(Float_t maxpos, Float_t y) const;
96 Double_t SearchYup(Float_t maxpos, Float_t y) const;
97/*
98 // Evaluate first solution for a possible maximum (x|first deriv==0)
99 inline Double_t EvalDerivEq0S1(const Int_t i) const
100 {
101 // return the x value [0;1[ at which the derivative is zero (solution1)
102
103 Double_t sumder = fDer2[i]+fDer2[i+1];
104 Double_t difder = fDer2[i]-fDer2[i+1];
105
106 Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
107 Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
108
109 Double_t x = 3*fDer2[i] - sqrt(3*sqt1 + 3*sqt2);
110
111 Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
112
113 return -x/denom;
114 }
115
116 // Evaluate second solution for a possible maximum (x|first deriv==0)
117 inline Double_t EvalDerivEq0S2(const Int_t i) const
118 {
119 // return the x value [0;1[ at which the derivative is zero (solution2)
120
121 Double_t sumder = fDer2[i]+fDer2[i+1];
122 Double_t difder = fDer2[i]-fDer2[i+1];
123
124 Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
125 Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
126
127 Double_t x = 3*fDer2[i] + sqrt(3*sqt1 + 3*sqt2);
128
129 Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
130
131 return -x/denom;
132 }
133 */
134
135 inline void EvalDerivEq0(const Int_t i, Float_t &rc1, Float_t &rc2) const
136 {
137 Double_t sumder = fDer2[i]+fDer2[i+1];
138 Double_t difder = fDer2[i]-fDer2[i+1];
139
140 Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
141 Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
142 Double_t sqt3 = sqrt(3*sqt1 + 3*sqt2);
143 Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
144
145 rc1 = -(3*fDer2[i] + sqt3)/denom;
146 rc2 = -(3*fDer2[i] - sqt3)/denom;
147 }
148
149 // Calculate the "Stammfunktion" of the Eval-function
150 inline Double_t EvalPrimitive(Int_t i, Float_t x) const
151 {
152 /* TO BE CHECKED!
153 if (x==0)
154 return 0;
155
156 if (x==1)
157 return (fVal[i+1]+fVal[i])/2 - fDer2[i+1]/4;
158 */
159 Align(i, x);
160
161 const Double_t x2 = x*x;
162 const Double_t x4 = x2*x2;
163 const Double_t x1 = 1-x;
164 const Double_t x14 = x1*x1*x1*x1;
165
166 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];
167 }
168
169 inline void Align(Int_t &i, Float_t &x) const
170 {
171 if (i<0)
172 {
173 x += i;
174 i=0;
175 }
176 if (i>=fNum-1)
177 {
178 x += i-(fNum-2);
179 i=fNum-2;
180 }
181 }
182
183 // Calculate the intgeral of the Eval-function in
184 // bin i from a=[0;1[ to b=[0;1[
185 inline Double_t EvalInteg(Int_t i, Float_t a=0, Float_t b=1) const
186 {
187 // This is to make sure that we never access invalid
188 // memory, even if this should never happen.
189 // If it happens anyhow we extraolate the spline
190 Align(i, a);
191 Align(i, b);
192
193 return EvalPrimitive(i, b)-EvalPrimitive(i, a);
194 }
195
196 // Calculate the intgeral of the Eval-function betwen x0 and x1
197 inline Double_t EvalInteg(Float_t x0, Float_t x1) const
198 {
199 const Int_t min = TMath::CeilNint(x0);
200 const Int_t max = TMath::FloorNint(x1);
201
202 // This happens if x0 and x1 are in the same interval
203 if (min>max)
204 return EvalInteg(max, x0-max, x1-max);
205
206 // Sum complete intervals
207 Double_t sum = 0;
208 for (int i=min; i<max; i++)
209 sum += EvalInteg(i);
210
211 // Sum the incomplete intervals at the beginning and end
212 sum += EvalInteg(min-1, 1-(min-x0), 1);
213 sum += EvalInteg(max, 0, x1-max);
214
215 // return result
216 return sum;
217 }
218
219 // We search for the maximum from x=i-1 to x=i+1
220 // (Remeber: i corresponds to the value in bin i, i+1 to the
221 // next bin and i-1 to the last bin)
222 inline void GetMaxAroundI(Int_t i, Float_t &xmax, Float_t &ymax) const
223 {
224 Float_t xmax1=0, xmax2=0;
225 Float_t ymax1=0, ymax2=0;
226
227 Bool_t rc1 = i>0 && GetMax(i-1, xmax1, ymax1);
228 Bool_t rc2 = i<fNum-1 && GetMax(i, xmax2, ymax2);
229
230 // In case the medium bin is the first or last bin
231 // take the lower or upper edge of the region into account.
232 if (i==0)
233 {
234 xmax1 = 0;
235 ymax1 = fVal[0];
236 rc1 = kTRUE;
237 }
238 if (i>=fNum-1)
239 {
240 xmax2 = fNum-1;
241 ymax2 = fVal[fNum-1];
242 rc2 = kTRUE;
243 }
244
245 // Take a default in case no maximum is found
246 // FIXME: Check THIS!!!
247 xmax=i;
248 ymax=fVal[i];
249
250 if (rc1)
251 {
252 ymax = ymax1;
253 xmax = xmax1;
254 }
255 else
256 if (rc2)
257 {
258 ymax = ymax2;
259 xmax = xmax2;
260 }
261
262 if (rc2 && ymax2>ymax)
263 {
264 ymax = ymax2;
265 xmax = xmax2;
266 }
267 /*
268 // Search real maximum in [i-0.5, i+1.5]
269 Float_t xmax1, xmax2, xmax3;
270 Float_t ymax1, ymax2, ymax3;
271
272 Bool_t rc1 = i>0 && GetMax(i-1, xmax1, ymax1, 0.5, 1.0);
273 Bool_t rc2 = GetMax(i, xmax2, ymax2, 0.0, 1.0);
274 Bool_t rc3 = i<fNum-1 && GetMax(i+1, xmax3, ymax3, 0.0, 0.5);
275
276 // In case the medium bin is the first or last bin
277 // take the lower or upper edge of the region into account.
278 if (i==0)
279 {
280 xmax1 = 0;
281 ymax1 = Eval(0, 0);
282 rc1 = kTRUE;
283 }
284 if (i==fNum-1)
285 {
286 xmax3 = fNum-1e-5;
287 ymax3 = Eval(fNum-1, 1);
288 rc3 = kTRUE;
289 }
290
291 // Take a real default in case no maximum is found
292 xmax=i+0.5;
293 ymax=Eval(i, 0.5);
294
295 //if (!rc1 && !rc2 && !rc3)
296 // cout << "!!!!!!!!!!!!!!!" << endl;
297
298 if (rc1)
299 {
300 ymax = ymax1;
301 xmax = xmax1;
302 }
303 else
304 if (rc2)
305 {
306 ymax = ymax2;
307 xmax = xmax2;
308 }
309 else
310 if (rc3)
311 {
312 ymax = ymax3;
313 xmax = xmax3;
314 }
315
316 if (rc2 && ymax2>ymax)
317 {
318 ymax = ymax2;
319 xmax = xmax2;
320 }
321 if (rc3 && ymax3>ymax)
322 {
323 ymax = ymax3;
324 xmax = xmax3;
325 }
326*/ }
327
328 inline Bool_t GetMax(Int_t i, Float_t &xmax, Float_t &ymax, Float_t min=0, Float_t max=1) const
329 {
330 // Find analytical maximum in the bin i in the interval [min,max[
331
332 Float_t x1, x2;
333 EvalDerivEq0(i, x1, x2);
334 // const Float_t x1 = EvalDerivEq0S1(i);
335 // const Float_t x2 = EvalDerivEq0S2(i);
336
337 const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(x1, i)<0;
338 const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(x2, i)<0;
339
340 if (!ismax1 && !ismax2)
341 return kFALSE;
342
343 if (ismax1 && !ismax2)
344 {
345 xmax = i+x1;
346 ymax = Eval(i, x1);
347 return kTRUE;
348 }
349
350 if (!ismax1 && ismax2)
351 {
352 xmax = i+x2;
353 ymax = Eval(i, x2);
354 return kTRUE;
355 }
356
357 // Somehting must be wrong...
358 return kFALSE;
359 /*
360 std::cout << "?????????????" << std::endl;
361
362 const Double_t y1 = Eval(i, x1);
363 const Double_t y2 = Eval(i, x2);
364
365 if (y1>y2)
366 {
367 xmax = i+x1;
368 ymax = Eval(i, x1);
369 return kTRUE;
370 }
371 else
372 {
373 xmax = i+x2;
374 ymax = Eval(i, x2);
375 return kTRUE;
376 }
377
378 return kFALSE;*/
379 }
380/*
381 inline Int_t GetMaxPos(Int_t i, Float_t &xmax, Float_t &ymax) const
382 {
383 Double_t x[3];
384
385 x[0] = 0;
386 // x[1] = 1; // This means we miss a possible maximum at the
387 // upper edge of the last interval...
388
389 x[1] = EvalDerivEq0S1(i);
390 x[2] = EvalDerivEq0S2(i);
391
392 //y[0] = Eval(i, x[0]);
393 //y[1] = Eval(i, x[1]);
394 //y[1] = Eval(i, x[1]);
395 //y[2] = Eval(i, x[2]);
396
397 Int_t rc = 0;
398 Double_t max = Eval(i, x[0]);
399
400 for (Int_t j=1; j<3; j++)
401 {
402 if (x[j]<=0 || x[j]>=1)
403 continue;
404
405 const Float_t y = Eval(i, x[j]);
406 if (y>max)
407 {
408 max = y;
409 rc = j;
410 }
411 }
412
413 if (max>ymax)
414 {
415 xmax = x[rc]+i;
416 ymax = max;
417 }
418
419 return rc;
420 }
421
422 inline void GetMaxPos(Int_t min, Int_t max, Float_t &xmax, Float_t &ymax) const
423 {
424 Float_t xmax=-1;
425 Float_t ymax=-FLT_MAX;
426
427 for (int i=min; i<max; i++)
428 GetMaxPos(i, xmax, ymax);
429
430 for (int i=min+1; i<max; i++)
431 {
432 Float_t y = Eval(i, 0);
433 if (y>ymax)
434 {
435 ymax = y;
436 xmax = i;
437 }
438 }
439
440 }*/
441
442
443 void InitDerivatives() const;
444 Float_t CalcIntegral(Float_t start) const;
445
446public:
447 MExtralgoSpline(const Float_t *val, Int_t n, Float_t *der1, Float_t *der2)
448 : fExtractionType(kIntegralRel), fVal(val), fNum(n), fDer1(der1), fDer2(der2), fHeightTm(0.5), fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
449 {
450 InitDerivatives();
451 }
452
453 void SetRiseFallTime(Float_t rise, Float_t fall) { fRiseTime=rise; fFallTime=fall; }
454 void SetExtractionType(ExtractionType_t typ) { fExtractionType = typ; }
455 void SetHeightTm(Float_t h) { fHeightTm = h; }
456 // void SetResolution(Float_t res) { fResolution=res; }
457
458 Float_t GetTime() const { return fTime; }
459 Float_t GetWidth() const { return fWidth; }
460 Float_t GetSignal() const { return fSignal; }
461 Float_t GetHeight() const { return fHeight; }
462
463 Float_t GetTimeDev() const { return fTimeDev; }
464 Float_t GetWidthDev() const { return fWidthDev; }
465 Float_t GetSignalDev() const { return fSignalDev; }
466
467 void GetSignal(Float_t &sig, Float_t &dsig) const { sig=fSignal; dsig=fSignalDev; }
468 void GetWidth(Float_t &sig, Float_t &dsig) const { sig=fWidth; dsig=fWidthDev; }
469 void GetTime(Float_t &sig, Float_t &dsig) const { sig=fTime; dsig=fTimeDev; }
470
471 Float_t ExtractNoise(/*Int_t iter*/);
472 void Extract(Byte_t sat, Int_t maxpos, Bool_t width=kFALSE);
473};
474
475#endif
Note: See TracBrowser for help on using the repository browser.