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

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