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

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