Changeset 7999
- Timestamp:
- 10/01/06 22:19:55 (18 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r7998 r7999 18 18 19 19 -*-*- END OF LINE -*-*- 20 2006/10/01 Thomas Bretz 21 22 * mextralgo/MExtralgoSpline.[h,cc]: 23 - changed from the old fashined search algorithm to a completely 24 analytical approach. Still with a lot of comments containing 25 the old code 26 27 * mbase/MMath.[h,cc]: 28 - added new function to solve polynomial equations up to the 29 thirs order. 30 31 32 20 33 2006/09/29 Thomas Bretz 21 34 -
trunk/MagicSoft/Mars/mbase/MMath.cc
r7899 r7999 39 39 #include <TArrayD.h> 40 40 #endif 41 42 #ifndef ROOT_TComplex 43 #include <TComplex.h> 44 #endif 45 41 46 // -------------------------------------------------------------------------- 42 47 // … … 495 500 return rc; 496 501 } 502 503 Int_t MMath::SolvePol2(Double_t a, Double_t b, Double_t &x1, Double_t &x2) 504 { 505 const Double_t r = a*a - 4*b; 506 if (r<0) 507 return 0; 508 509 if (r==0) 510 { 511 x1 = -a/2; 512 return 1; 513 } 514 515 const Double_t s = TMath::Sqrt(r); 516 517 x1 = (-a+s)/2; 518 x2 = (-a-s)/2; 519 520 return 2; 521 } 522 523 // -------------------------------------------------------------------------- 524 // 525 // This is a helper function making the execution of SolverPol3 a bit faster 526 // 527 static inline Double_t ReMul(const TComplex &c1, const TComplex &th) 528 { 529 const TComplex c2 = TComplex::Cos(th/3.); 530 return c1.Re() * c2.Re() - c1.Im() * c2.Im(); 531 } 532 533 // -------------------------------------------------------------------------- 534 // 535 // Solves: x^3 + ax^2 + bx + c = 0; 536 // Return number of the real solutions returns as z1, z2, z3 537 // 538 // Algorithm adapted from http://home.att.net/~srschmitt/cubizen.heml 539 // Which is based on the solution given in 540 // http://mathworld.wolfram.com/CubicEquation.html 541 // 542 // ------------------------------------------------------------------------- 543 // 544 // Exact solutions of cubic polynomial equations 545 // by Stephen R. Schmitt Algorithm 546 // 547 // An exact solution of the cubic polynomial equation: 548 // 549 // x^3 + a*x^2 + b*x + c = 0 550 // 551 // was first published by Gerolamo Cardano (1501-1576) in his treatise, 552 // Ars Magna. He did not discoverer of the solution; a professor of 553 // mathematics at the University of Bologna named Scipione del Ferro (ca. 554 // 1465-1526) is credited as the first to find an exact solution. In the 555 // years since, several improvements to the original solution have been 556 // discovered. Zeno source code 557 // 558 // % compute real or complex roots of cubic polynomial 559 // function cubic( var z1, z2, z3 : real, a, b, c : real ) : real 560 // 561 // var Q, R, D, S, T : real 562 // var im, th : real 563 // 564 // Q := (3*b - a^2)/9 565 // R := (9*b*a - 27*c - 2*a^3)/54 566 // D := Q^3 + R^2 % polynomial discriminant 567 // 568 // if (D >= 0) then % complex or duplicate roots 569 // 570 // S := sgn(R + sqrt(D))*abs(R + sqrt(D))^(1/3) 571 // T := sgn(R - sqrt(D))*abs(R - sqrt(D))^(1/3) 572 // 573 // z1 := -a/3 + (S + T) % real root 574 // z2 := -a/3 - (S + T)/2 % real part of complex root 575 // z3 := -a/3 - (S + T)/2 % real part of complex root 576 // im := abs(sqrt(3)*(S - T)/2) % complex part of root pair 577 // 578 // else % distinct real roots 579 // 580 // th := arccos(R/sqrt( -Q^3)) 581 // 582 // z1 := 2*sqrt(-Q)*cos(th/3) - a/3 583 // z2 := 2*sqrt(-Q)*cos((th + 2*pi)/3) - a/3 584 // z3 := 2*sqrt(-Q)*cos((th + 4*pi)/3) - a/3 585 // im := 0 586 // 587 // end if 588 // 589 // return im % imaginary part 590 // 591 // end function 592 // 593 Int_t MMath::SolvePol3(Double_t a, Double_t b, Double_t c, 594 Double_t &x1, Double_t &x2, Double_t &x3) 595 { 596 const Double_t Q = (a*a - 3*b)/9; 597 const Double_t R = (9*b*a - 27*c - 2*a*a*a)/54; 598 const Double_t D = R*R - Q*Q*Q; // polynomial discriminant 599 600 // ----- The single-real / duplicate-roots solution ----- 601 602 if (D==0) 603 { 604 const Double_t r = MMath::Sqrt3(R); 605 606 x1 = 2*r - a/3.; // real root 607 x2 = r - a/3.; // real root 608 609 return 2; 610 } 611 612 if (D>0) // complex or duplicate roots 613 { 614 const Double_t sqrtd = TMath::Sqrt(D); 615 616 const Double_t S = MMath::Sqrt3(R + sqrtd); 617 const Double_t T = MMath::Sqrt3(R - sqrtd); 618 619 x1 = (S+T) - a/3.; // real root 620 621 return 1; 622 623 //z2 = (S + T)/2 - a/3.; // real part of complex root 624 //z3 = (S + T)/2 - a/3.; // real part of complex root 625 //im = fabs(sqrt(3)*(S - T)/2) // complex part of root pair 626 } 627 628 // ----- The general solution with three roots --- 629 630 if (Q==0) 631 return 0; 632 633 if (Q>0) // This is here for speed reasons 634 { 635 const Double_t sqrtq = TMath::Sqrt(Q); 636 const Double_t rq = R/TMath::Abs(Q); 637 638 const Double_t th1 = TMath::ACos(rq/sqrtq); 639 const Double_t th2 = th1 + TMath::TwoPi(); 640 const Double_t th3 = th2 + TMath::TwoPi(); 641 642 x1 = 2.*sqrtq * TMath::Cos(th1/3.) - a/3.; 643 x2 = 2.*sqrtq * TMath::Cos(th2/3.) - a/3.; 644 x3 = 2.*sqrtq * TMath::Cos(th3/3.) - a/3.; 645 646 return 3; 647 } 648 649 const TComplex sqrtq = TComplex::Sqrt(Q); 650 const Double_t rq = R/TMath::Abs(Q); 651 652 const TComplex th1 = TComplex::ACos(rq/sqrtq); 653 const TComplex th2 = th1 + TMath::TwoPi(); 654 const TComplex th3 = th2 + TMath::TwoPi(); 655 656 // For ReMul, see bove 657 x1 = ReMul(2.*sqrtq, th1) - a/3.; 658 x2 = ReMul(2.*sqrtq, th2) - a/3.; 659 x3 = ReMul(2.*sqrtq, th3) - a/3.; 660 661 return 3; 662 } -
trunk/MagicSoft/Mars/mbase/MMath.h
r7971 r7999 47 47 Double_t InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x); 48 48 49 inline Int_t SolvePol1(Double_t c, Double_t d, Double_t &x1) 50 { 51 if (c==0) 52 return 0; 53 54 x1 = -d/c; 55 return 1; 56 } 57 Int_t SolvePol2(Double_t c, Double_t d, Double_t &x1, Double_t &x2); 58 inline Int_t SolvePol2(Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2) 59 { 60 return b==0 ? SolvePol1(c, d, x1) : SolvePol2(c/b, d/b, x1, x2); 61 } 62 Int_t SolvePol3(Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2, Double_t &x3); 63 inline Int_t SolvePol3(Double_t a, Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2, Double_t &x3) 64 { 65 return a==0 ? SolvePol2(b, c, d, x1, x2) : SolvePol3(b/a, c/a, d/a, x1, x2, x3); 66 } 67 49 68 TArrayD LeastSqFitExpW1(Int_t n, Double_t *x, Double_t *y); 50 69 TArrayD LeastSqFitExp(Int_t n, Double_t *x, Double_t *y); … … 54 73 inline Int_t ModF(Double_t dbl, Double_t &frac) { Double_t rc; frac = modf(dbl, &rc); return TMath::Nint(rc); } 55 74 75 inline Double_t Sqrt3(Double_t x) { return TMath::Sign(TMath::Power(TMath::Abs(x), 1./3), x); } 76 56 77 inline Double_t Sgn(Double_t d) { return d<0 ? -1 : 1; } 57 78 } -
trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc
r7942 r7999 30 30 // Numerical Recipes in C++, 2nd edition, pp. 116-119. 31 31 // 32 // The coefficients "ya" are here denoted as "f HiGainSignal" and "fLoGainSignal"33 // which meansthe FADC value subtracted by the clock-noise corrected pedestal.32 // The coefficients "ya" are here denoted as "fVal" corresponding to 33 // the FADC value subtracted by the clock-noise corrected pedestal. 34 34 // 35 35 // The coefficients "y2a" get immediately divided 6. and are called here 36 // "fHiGainSecondDeriv" and "fLoGainSecondDeriv" although they are now not exactly37 // the second derivative coefficients any more.36 // fDer2 although they are now not exactly the second derivative 37 // coefficients any more. 38 38 // 39 39 // The calculation of the cubic-spline interpolated value "y" on a point 40 // "x" along the FADC-slices axis becomes: 41 // 42 // y = a*fHiGainSignal[klo] + b*fHiGainSignal[khi] 43 // + (a*a*a-a)*fHiGainSecondDeriv[klo] + (b*b*b-b)*fHiGainSecondDeriv[khi] 44 // 45 // with: 46 // a = (khi - x) 47 // b = (x - klo) 48 // 49 // and "klo" being the lower bin edge FADC index and "khi" the upper bin edge FADC index. 50 // fHiGainSignal[klo] and fHiGainSignal[khi] are the FADC values at "klo" and "khi". 51 // 52 // An analogues formula is used for the low-gain values. 53 // 54 // The coefficients fHiGainSecondDeriv and fLoGainSecondDeriv are calculated with the 55 // following simplified algorithm: 56 // 57 // for (Int_t i=1;i<range-1;i++) { 58 // pp = fHiGainSecondDeriv[i-1] + 4.; 59 // fHiGainFirstDeriv[i] = fHiGainSignal[i+1] - 2.*fHiGainSignal[i] + fHiGainSignal[i-1] 60 // fHiGainFirstDeriv[i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp; 61 // } 62 // 63 // for (Int_t k=range-2;k>=0;k--) 64 // fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.; 65 // 66 // 67 // This algorithm takes advantage of the fact that the x-values are all separated by exactly 1 68 // which simplifies the Numerical Recipes algorithm. 69 // (Note that the variables "fHiGainFirstDeriv" are not real first derivative coefficients.) 70 // 71 // The algorithm to search the time proceeds as follows: 72 // 73 // 1) Calculate all fHiGainSignal from fHiGainFirst to fHiGainLast 74 // (note that an "overlap" to the low-gain arrays is possible: i.e. fHiGainLast>14 in the case of 75 // the MAGIC FADCs). 76 // 2) Remember the position of the slice with the highest content "fAbMax" at "fAbMaxPos". 77 // 3) If one or more slices are saturated or fAbMaxPos is less than 2 slices from fHiGainFirst, 78 // return fAbMaxPos as time and fAbMax as charge (note that the pedestal is subtracted here). 79 // 4) Calculate all fHiGainSecondDeriv from the fHiGainSignal array 80 // 5) Search for the maximum, starting in interval fAbMaxPos-1 in steps of 0.2 till fAbMaxPos-0.2. 81 // If no maximum is found, go to interval fAbMaxPos+1. 82 // --> 4 function evaluations 83 // 6) Search for the absolute maximum from fAbMaxPos to fAbMaxPos+1 in steps of 0.2 84 // --> 4 function evaluations 85 // 7) Try a better precision searching from new max. position fAbMaxPos-0.2 to fAbMaxPos+0.2 86 // in steps of 0.025 (83 psec. in the case of the MAGIC FADCs). 87 // --> 14 function evaluations 88 // 8) If Time Extraction Type kMaximum has been chosen, the position of the found maximum is 89 // returned, else: 90 // 9) The Half Maximum is calculated. 91 // 10) fHiGainSignal is called beginning from fAbMaxPos-1 backwards until a value smaller than fHalfMax 92 // is found at "klo". 93 // 11) Then, the spline value between "klo" and "klo"+1 is halfed by means of bisection as long as 94 // the difference between fHalfMax and spline evaluation is less than fResolution (default: 0.01). 95 // --> maximum 12 interations. 96 // 97 // The algorithm to search the charge proceeds as follows: 98 // 99 // 1) If Charge Type: kAmplitude was chosen, return the Maximum of the spline, found during the 100 // time search. 101 // 2) If Charge Type: kIntegral was chosen, sum the fHiGainSignal between: 102 // (Int_t)(fAbMaxPos - fRiseTimeHiGain) and 103 // (Int_t)(fAbMaxPos + fFallTimeHiGain) 104 // (default: fRiseTime: 1.5, fFallTime: 4.5) 105 // sum the fLoGainSignal between: 106 // (Int_t)(fAbMaxPos - fRiseTimeHiGain*fLoGainStretch) and 107 // (Int_t)(fAbMaxPos + fFallTimeHiGain*fLoGainStretch) 108 // (default: fLoGainStretch: 1.5) 109 // 110 // The values: fNumHiGainSamples and fNumLoGainSamples are set to: 111 // 1) If Charge Type: kAmplitude was chosen: 1. 112 // 2) If Charge Type: kIntegral was chosen: fRiseTimeHiGain + fFallTimeHiGain 113 // or: fNumHiGainSamples*fLoGainStretch in the case of the low-gain 114 // 115 // Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 116 // to modify the ranges. 117 // 118 // Defaults: 119 // fHiGainFirst = 2 120 // fHiGainLast = 14 121 // fLoGainFirst = 2 122 // fLoGainLast = 14 123 // 124 // Call: SetResolution() to define the resolution of the half-maximum search. 125 // Default: 0.01 126 // 127 // Call: SetRiseTime() and SetFallTime() to define the integration ranges 128 // for the case, the extraction type kIntegral has been chosen. 129 // 130 // Call: - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the 131 // computation of the amplitude at the maximum (default) and extraction 132 // the position of the maximum (default) 133 // --> no further function evaluation needed 134 // - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the 135 // computation of the integral beneith the spline between fRiseTimeHiGain 136 // from the position of the maximum to fFallTimeHiGain after the position of 137 // the maximum. The Low Gain is computed with half a slice more at the rising 138 // edge and half a slice more at the falling edge. 139 // The time of the half maximum is returned. 140 // --> needs one function evaluations but is more precise 141 // 40 // "x" along the FADC-slices axis becomes: EvalAt(x) 41 // 42 // The coefficients fDer2 are calculated with the simplified 43 // algorithm in InitDerivatives. 44 // 45 // This algorithm takes advantage of the fact that the x-values are all 46 // separated by exactly 1 which simplifies the Numerical Recipes algorithm. 47 // (Note that the variables fDer are not real first derivative coefficients.) 48 // 142 49 ////////////////////////////////////////////////////////////////////////////// 143 50 #include "MExtralgoSpline.h" 144 51 52 #include "../mbase/MMath.h" 53 145 54 using namespace std; 146 55 56 // -------------------------------------------------------------------------- 57 // 58 // Calculate the first and second derivative for the splie. 59 // 60 // The coefficients are calculated such that 61 // 1) fVal[i] = Eval(i, 0) 62 // 2) Eval(i-1, 1)==Eval(i, 0) 63 // 64 // In other words: The values with the index i describe the spline 65 // between fVal[i] and fVal[i+1] 66 // 147 67 void MExtralgoSpline::InitDerivatives() const 148 68 { … … 169 89 } 170 90 171 Float_t MExtralgoSpline::CalcIntegral(Float_t start, Float_t range) const 172 { 91 // -------------------------------------------------------------------------- 92 // 93 // Returns the highest x value in [min;max[ at which the spline in 94 // the bin i is equal to y 95 // 96 // min and max are defined to be [0;1] 97 // 98 // The default for min is 0, the default for max is 1 99 // The defaule for y is 0 100 // 101 Double_t MExtralgoSpline::FindY(Int_t i, Double_t y, Double_t min, Double_t max) const 102 { 103 // y = a*x^3 + b*x^2 + c*x + d' 104 // 0 = a*x^3 + b*x^2 + c*x + d' - y 105 106 // Calculate coefficients 107 const Double_t a = fDer2[i+1]-fDer2[i]; 108 const Double_t b = 3*fDer2[i]; 109 const Double_t c = fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1]; 110 const Double_t d = fVal[i] - y; 111 112 Double_t x1, x2, x3; 113 const Int_t rc = MMath::SolvePol3(a, b, c, d, x1, x2, x3); 114 115 Double_t x = -1; 116 if (rc>0 && x1>=min && x1<max && x1>x) 117 x = x1; 118 if (rc>1 && x2>=min && x2<max && x2>x) 119 x = x2; 120 if (rc>2 && x3>=min && x3<max && x3>x) 121 x = x3; 122 123 return x<0 ? -1 : x+i; 124 } 125 126 // -------------------------------------------------------------------------- 127 // 128 // Search analytically downward for the value y of the spline, starting 129 // at x, until x==0. If y is not found -1 is returned. 130 // 131 Double_t MExtralgoSpline::SearchY(Float_t x, Float_t y) const 132 { 133 if (x>=fNum-1) 134 x = fNum-1.0001; 135 136 Int_t i = TMath::FloorNint(x); 137 Double_t rc = FindY(i, y, 0, x-i); 138 while (--i>=0 && rc<0) 139 rc = FindY(i, y); 140 141 return rc; 142 } 143 144 // -------------------------------------------------------------------------- 145 // 146 // Do a range check an then calculate the integral from start-fRiseTime 147 // to start+fFallTime. An extrapolation of 0.5 slices is allowed. 148 // 149 Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const 150 { 151 /* 173 152 // The number of steps is calculated directly from the integration 174 153 // window. This is the only way to ensure we are not dealing with … … 176 155 // value under the same conditions -- it might still be different on 177 156 // other machines! 157 const Float_t start = pos-fRiseTime; 178 158 const Float_t step = 0.2; 179 159 const Float_t width = fRiseTime+fFallTime; 180 const Float_t max = range-1 - (width+step);160 const Float_t max = fNum-1 - (width+step); 181 161 const Int_t num = TMath::Nint(width/step); 182 162 … … 193 173 for (Int_t i=0; i<num; i++) 194 174 { 195 const Float_t x = start+i*step;196 const Int_t klo = (Int_t)TMath::Floor(x);197 175 // Note: if x is close to one integer number (= a FADC sample) 198 176 // we get the same result by using that sample as klo, and the … … 203 181 // depending on the compilation options). 204 182 205 sum += Eval (x, klo);183 sum += EvalAt(start + i*step); 206 184 207 185 // FIXME? Perhaps the integral should be done analitically … … 209 187 } 210 188 sum *= step; // Transform sum in integral 189 211 190 return sum; 191 */ 192 193 // In the future we will calculate the intgeral analytically. 194 // It has been tested that it gives identical results within 195 // acceptable differences. 196 197 // We allow extrapolation of 1/2 slice. 198 const Float_t min = fRiseTime; //-0.5+fRiseTime; 199 const Float_t max = fNum-1-fFallTime; //fNum-0.5+fFallTime; 200 201 if (pos<min) 202 pos = min; 203 if (pos>max) 204 pos = max; 205 206 return EvalInteg(pos-fRiseTime, pos+fFallTime); 212 207 } 213 208 … … 217 212 218 213 if (fExtractionType == kAmplitude) 219 return Eval( nsx+1, 1);214 return Eval(1, nsx); 220 215 else 221 return CalcIntegral(2. + nsx , fNum);222 } 223 224 void MExtralgoSpline::Extract(Byte_t sat, Int_t max pos)216 return CalcIntegral(2. + nsx); 217 } 218 219 void MExtralgoSpline::Extract(Byte_t sat, Int_t maxbin) 225 220 { 226 221 fSignal = 0; … … 233 228 // Don't start if the maxpos is too close to the limits. 234 229 // 235 const Bool_t limlo = maxpos < TMath::Ceil(fRiseTime); 236 const Bool_t limup = maxpos > fNum-TMath::Ceil(fFallTime)-1; 230 231 /* 232 const Bool_t limlo = maxbin < TMath::Ceil(fRiseTime); 233 const Bool_t limup = maxbin > fNum-TMath::Ceil(fFallTime)-1; 237 234 if (sat || limlo || limup) 238 235 { … … 240 237 if (fExtractionType == kAmplitude) 241 238 { 242 fSignal = fVal[max pos];243 fTime = max pos;239 fSignal = fVal[maxbin]; 240 fTime = maxbin; 244 241 fSignalDev = 0; // means: is valid 245 242 return; 246 243 } 247 244 248 fSignal = CalcIntegral(limlo ? 0 : fNum , fNum);249 fTime = max pos- 1;245 fSignal = CalcIntegral(limlo ? 0 : fNum); 246 fTime = maxbin - 1; 250 247 fSignalDev = 0; // means: is valid 251 248 return; 252 249 } 250 */ 253 251 254 252 fTimeDev = fResolution; … … 257 255 // Now find the maximum 258 256 // 257 258 259 /* 259 260 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one 260 261 261 Int_t klo = maxpos-1; 262 Float_t fAbMaxPos = maxpos;//! Current position of the maximum of the spline 263 Float_t fAbMax = fVal[maxpos];//! Current maximum of the spline 262 Int_t klo = maxbin-1; 263 264 Float_t maxpos = maxbin;//! Current position of the maximum of the spline 265 Float_t max = fVal[maxbin];//! Current maximum of the spline 264 266 265 267 // … … 270 272 { 271 273 const Float_t x = klo + step*(i+1); 272 const Float_t y = Eval(x, klo); 273 if (y > fAbMax) 274 { 275 fAbMax = y; 276 fAbMaxPos = x; 274 //const Float_t y = Eval(klo, x); 275 const Float_t y = Eval(klo, x-klo); 276 if (y > max) 277 { 278 max = y; 279 maxpos = x; 277 280 } 278 281 } … … 281 284 // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2 282 285 // 283 if ( fAbMaxPos > maxpos- 0.1)284 { 285 klo = max pos;286 if (maxpos > maxbin - 0.1) 287 { 288 klo = maxbin; 286 289 287 290 for (Int_t i=0; i<TMath::Nint(TMath::Ceil((1-0.3)/step)); i++) 288 291 { 289 292 const Float_t x = klo + step*(i+1); 290 const Float_t y = Eval(x, klo); 291 if (y > fAbMax) 293 //const Float_t y = Eval(klo, x); 294 const Float_t y = Eval(klo, x-klo); 295 if (y > max) 292 296 { 293 fAbMax = y;294 fAbMaxPos = x;297 max = y; 298 maxpos = x; 295 299 } 296 300 } … … 301 305 // Try a better precision. 302 306 // 303 const Float_t up = fAbMaxPos+step - 3.0*fResolution;304 const Float_t lo = fAbMaxPos-step + 3.0*fResolution;305 const Float_t abmaxpos = fAbMaxPos;307 const Float_t up = maxpos+step - 3.0*fResolution; 308 const Float_t lo = maxpos-step + 3.0*fResolution; 309 const Float_t abmaxpos = maxpos; 306 310 307 311 step = 2.*fResolution; // step size of 0.1 FADC slices … … 310 314 { 311 315 const Float_t x = abmaxpos + (i+1)*step; 312 const Float_t y = Eval(x, klo); 313 if (y > fAbMax) 314 { 315 fAbMax = y; 316 fAbMaxPos = x; 316 //const Float_t y = Eval(klo, x); 317 const Float_t y = Eval(klo, x-klo); 318 if (y > max) 319 { 320 max = y; 321 maxpos = x; 317 322 } 318 323 } … … 333 338 { 334 339 const Float_t x = abmaxpos - (i+1)*step; 335 const Float_t y = Eval(x, klo); 336 if (y > fAbMax) 337 { 338 fAbMax = y; 339 fAbMaxPos = x; 340 } 341 } 340 //const Float_t y = Eval(klo, x); 341 const Float_t y = Eval(klo, x-klo); 342 if (y > max) 343 { 344 max = y; 345 maxpos = x; 346 } 347 } 348 */ 349 350 // --- Start NEW --- 351 352 // This block extracts values very similar to the old algorithm... 353 // for max>10 354 /* Most accurate (old identical) version: 355 356 Float_t xmax=maxpos, ymax=Eval(maxpos-1, 1); 357 Int_t rc = GetMaxPos(maxpos-1, xmax, ymax); 358 if (xmax==maxpos) 359 { 360 GetMaxPos(maxpos, xmax, ymax); 361 362 Float_t y = Eval(maxpos, 1); 363 if (y>ymax) 364 { 365 ymax = y; 366 xmax = maxpos+1; 367 } 368 }*/ 369 370 Float_t maxpos, maxval; 371 GetMaxAroundI(maxbin, maxpos, maxval); 372 373 // --- End NEW --- 342 374 343 375 if (fExtractionType == kAmplitude) 344 376 { 345 fTime = fAbMaxPos;346 fSignal = fAbMax;377 fTime = maxpos; 378 fSignal = maxval; 347 379 fSignalDev = 0; // means: is valid 348 380 return; 349 381 } 350 382 351 Float_t fHalfMax = fAbMax/2.;//! Current half maximum of the spline 352 353 // 354 // Now, loop from the maximum bin leftward down in order to find the 355 // position of the half maximum. First, find the right FADC slice: 356 // 357 klo = maxpos; 358 while (klo > 0) 359 { 360 if (fVal[--klo] < fHalfMax) 383 // 384 // Loop from the beginning of the slice upwards to reach the maxhalf: 385 // With means of bisection: 386 // 387 /* 388 static const Float_t sqrt2 = TMath::Sqrt(2.); 389 390 step = sqrt2*3*0.061;//fRiseTime; 391 Float_t x = maxpos-0.86-3*0.061;//fRiseTime*1.25; 392 393 // step = sqrt2*0.5;//fRiseTime; 394 // Float_t x = maxpos-1.25;//fRiseTime*1.25; 395 396 Int_t cnt =0; 397 while (cnt++<30) 398 { 399 const Float_t y=EvalAt(x); 400 401 if (TMath::Abs(y-maxval/2)<fResolution) 361 402 break; 362 } 363 364 // 365 // Loop from the beginning of the slice upwards to reach the fHalfMax: 366 // With means of bisection: 367 // 368 step = 0.5; 369 370 Int_t maxcnt = 20; 371 Int_t cnt = 0; 372 373 Float_t y = Eval(klo, klo); // FIXME: IS THIS CORRECT??????? 374 Float_t x = klo; 375 Bool_t back = kFALSE; 376 377 while (TMath::Abs(y-fHalfMax) > fResolution) 378 { 379 x += back ? -step : +step; 380 381 const Float_t y = Eval(x, klo); 382 383 back = y > fHalfMax; 384 385 if (++cnt > maxcnt) 386 break; 387 388 step /= 2.; 389 } 403 404 step /= sqrt2; // /2 405 x += y>maxval/2 ? -step : +step; 406 } 407 */ 408 409 // Search downwards for maxval/2 410 // By doing also a search upwards we could extract the pulse width 411 const Double_t x1 = SearchY(maxpos, maxval/2); 412 413 fTime = x1; 414 fSignal = CalcIntegral(maxpos); 415 fSignalDev = 0; // means: is valid 416 417 //if (fSignal>100) 418 // cout << "V=" << maxpos-x1 << endl; 390 419 391 420 // 392 421 // Now integrate the whole thing! 393 422 // 394 fTime =x;395 fSignal = CalcIntegral(fAbMaxPos - fRiseTime, fNum);396 fSignalDev = 0; // means: is valid397 } 423 // fTime = cnt==31 ? -1 : x; 424 // fSignal = cnt==31 ? CalcIntegral(x) : CalcIntegral(maxpos); 425 // fSignalDev = 0; // means: is valid 426 } -
trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h
r7942 r7999 5 5 #include <TROOT.h> 6 6 #endif 7 8 #include <iostream> 9 class TComplex; 7 10 8 11 class MExtralgoSpline … … 35 38 Float_t fSignalDev; 36 39 40 Double_t ReMul(const TComplex &c1, const TComplex &th) const; 41 37 42 inline Float_t Eval(Float_t val, Float_t a, Float_t deriv) const 38 43 { … … 40 45 } 41 46 42 inline Float_t Eval(const Float_t x, const Int_t i) const 43 { 44 const Float_t b = x-i; 45 return Eval(fVal[i], 1-b, fDer2[i]) + Eval(fVal[i+1], b, fDer2[i+1]); 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 xmax=i; 241 ymax=fVal[i]; 242 243 if (rc1) 244 { 245 ymax = ymax1; 246 xmax = xmax1; 247 } 248 else 249 if (rc2) 250 { 251 ymax = ymax2; 252 xmax = xmax2; 253 } 254 255 if (rc2 && ymax2>ymax) 256 { 257 ymax = ymax2; 258 xmax = xmax2; 259 } 260 /* 261 // Search real maximum in [i-0.5, i+1.5] 262 Float_t xmax1, xmax2, xmax3; 263 Float_t ymax1, ymax2, ymax3; 264 265 Bool_t rc1 = i>0 && GetMax(i-1, xmax1, ymax1, 0.5, 1.0); 266 Bool_t rc2 = GetMax(i, xmax2, ymax2, 0.0, 1.0); 267 Bool_t rc3 = i<fNum-1 && GetMax(i+1, xmax3, ymax3, 0.0, 0.5); 268 269 // In case the medium bin is the first or last bin 270 // take the lower or upper edge of the region into account. 271 if (i==0) 272 { 273 xmax1 = 0; 274 ymax1 = Eval(0, 0); 275 rc1 = kTRUE; 276 } 277 if (i==fNum-1) 278 { 279 xmax3 = fNum-1e-5; 280 ymax3 = Eval(fNum-1, 1); 281 rc3 = kTRUE; 282 } 283 284 // Take a real default in case no maximum is found 285 xmax=i+0.5; 286 ymax=Eval(i, 0.5); 287 288 //if (!rc1 && !rc2 && !rc3) 289 // cout << "!!!!!!!!!!!!!!!" << endl; 290 291 if (rc1) 292 { 293 ymax = ymax1; 294 xmax = xmax1; 295 } 296 else 297 if (rc2) 298 { 299 ymax = ymax2; 300 xmax = xmax2; 301 } 302 else 303 if (rc3) 304 { 305 ymax = ymax3; 306 xmax = xmax3; 307 } 308 309 if (rc2 && ymax2>ymax) 310 { 311 ymax = ymax2; 312 xmax = xmax2; 313 } 314 if (rc3 && ymax3>ymax) 315 { 316 ymax = ymax3; 317 xmax = xmax3; 318 } 319 */ } 320 321 inline Bool_t GetMax(Int_t i, Float_t &xmax, Float_t &ymax, Float_t min=0, Float_t max=1) const 322 { 323 // Find analytical maximum in the bin i in the interval [min,max[ 324 325 Float_t x1, x2; 326 EvalDerivEq0(i, x1, x2); 327 // const Float_t x1 = EvalDerivEq0S1(i); 328 // const Float_t x2 = EvalDerivEq0S2(i); 329 330 const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(x1, i)<0; 331 const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(x2, i)<0; 332 333 if (!ismax1 && !ismax2) 334 return kFALSE; 335 336 if (ismax1 && !ismax2) 337 { 338 xmax = i+x1; 339 ymax = Eval(i, x1); 340 return kTRUE; 341 } 342 343 if (!ismax1 && ismax2) 344 { 345 xmax = i+x2; 346 ymax = Eval(i, x2); 347 return kTRUE; 348 } 349 350 // Somehting must be wrong... 351 return kFALSE; 352 /* 353 std::cout << "?????????????" << std::endl; 354 355 const Double_t y1 = Eval(i, x1); 356 const Double_t y2 = Eval(i, x2); 357 358 if (y1>y2) 359 { 360 xmax = i+x1; 361 ymax = Eval(i, x1); 362 return kTRUE; 363 } 364 else 365 { 366 xmax = i+x2; 367 ymax = Eval(i, x2); 368 return kTRUE; 369 } 370 371 return kFALSE;*/ 372 } 373 /* 374 inline Int_t GetMaxPos(Int_t i, Float_t &xmax, Float_t &ymax) const 375 { 376 Double_t x[3]; 377 378 x[0] = 0; 379 // x[1] = 1; // This means we miss a possible maximum at the 380 // upper edge of the last interval... 381 382 x[1] = EvalDerivEq0S1(i); 383 x[2] = EvalDerivEq0S2(i); 384 385 //y[0] = Eval(i, x[0]); 386 //y[1] = Eval(i, x[1]); 387 //y[1] = Eval(i, x[1]); 388 //y[2] = Eval(i, x[2]); 389 390 Int_t rc = 0; 391 Double_t max = Eval(i, x[0]); 392 393 for (Int_t j=1; j<3; j++) 394 { 395 if (x[j]<=0 || x[j]>=1) 396 continue; 397 398 const Float_t y = Eval(i, x[j]); 399 if (y>max) 400 { 401 max = y; 402 rc = j; 403 } 404 } 405 406 if (max>ymax) 407 { 408 xmax = x[rc]+i; 409 ymax = max; 410 } 411 412 return rc; 413 } 414 415 inline void GetMaxPos(Int_t min, Int_t max, Float_t &xmax, Float_t &ymax) const 416 { 417 Float_t xmax=-1; 418 Float_t ymax=-FLT_MAX; 419 420 for (int i=min; i<max; i++) 421 GetMaxPos(i, xmax, ymax); 422 423 for (int i=min+1; i<max; i++) 424 { 425 Float_t y = Eval(i, 0); 426 if (y>ymax) 427 { 428 ymax = y; 429 xmax = i; 430 } 431 } 432 433 }*/ 434 47 435 48 436 void InitDerivatives() const; 49 Float_t CalcIntegral(Float_t start , Float_t range) const;437 Float_t CalcIntegral(Float_t start) const; 50 438 51 439 public:
Note:
See TracChangeset
for help on using the changeset viewer.