Changeset 7999 for trunk

Timestamp:

10/01/06 22:19:55 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mbase/MMath.cc (modified) (2 diffs)
• mbase/MMath.h (modified) (2 diffs)
• mextralgo/MExtralgoSpline.cc (modified) (15 diffs)
• mextralgo/MExtralgoSpline.h (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2006/10/01 Thomas Bretz
+
+   * mextralgo/MExtralgoSpline.[h,cc]:
+     - changed from the old fashined search algorithm to a completely
+       analytical approach. Still with a lot of comments containing
+       the old code
+
+   * mbase/MMath.[h,cc]:
+     - added new function to solve polynomial equations up to the
+       thirs order.
+
+
+
  2006/09/29 Thomas Bretz
 

trunk/MagicSoft/Mars/mbase/MMath.cc

@@ -39 +39 @@
 #include <TArrayD.h>
 #endif
+
+#ifndef ROOT_TComplex
+#include <TComplex.h>
+#endif
+
 // --------------------------------------------------------------------------
 //
@@ -495 +500 @@
     return rc;
 }
+
+Int_t MMath::SolvePol2(Double_t a, Double_t b, Double_t &x1, Double_t &x2)
+{
+    const Double_t r = a*a - 4*b;
+    if (r<0)
+        return 0;
+
+    if (r==0)
+    {
+        x1 = -a/2;
+        return 1;
+    }
+
+    const Double_t s = TMath::Sqrt(r);
+
+    x1 = (-a+s)/2;
+    x2 = (-a-s)/2;
+
+    return 2;
+}
+
+// --------------------------------------------------------------------------
+//
+// This is a helper function making the execution of SolverPol3 a bit faster
+//
+static inline Double_t ReMul(const TComplex &c1, const TComplex &th)
+{
+    const TComplex c2 = TComplex::Cos(th/3.);
+    return c1.Re() * c2.Re() - c1.Im() * c2.Im();
+}
+
+// --------------------------------------------------------------------------
+//
+// Solves: x^3 + ax^2 + bx + c = 0;
+// Return number of the real solutions returns as z1, z2, z3
+//
+// Algorithm adapted from http://home.att.net/~srschmitt/cubizen.heml
+// Which is based on the solution given in
+//    http://mathworld.wolfram.com/CubicEquation.html
+//
+// -------------------------------------------------------------------------
+//
+// Exact solutions of cubic polynomial equations
+// by Stephen R. Schmitt Algorithm
+//
+// An exact solution of the cubic polynomial equation:
+//
+// x^3 + a*x^2 + b*x + c = 0
+//
+// was first published by Gerolamo Cardano (1501-1576) in his treatise,
+// Ars Magna. He did not discoverer of the solution; a professor of
+// mathematics at the University of Bologna named Scipione del Ferro (ca.
+// 1465-1526) is credited as the first to find an exact solution. In the
+// years since, several improvements to the original solution have been
+// discovered. Zeno source code
+//
+// % compute real or complex roots of cubic polynomial
+// function cubic( var z1, z2, z3 : real, a, b, c : real ) : real
+// 
+//     var Q, R, D, S, T : real
+//     var im, th : real
+// 
+//     Q := (3*b - a^2)/9
+//     R := (9*b*a - 27*c - 2*a^3)/54
+//     D := Q^3 + R^2                          % polynomial discriminant
+// 
+//     if (D >= 0) then                        % complex or duplicate roots
+// 
+//         S := sgn(R + sqrt(D))*abs(R + sqrt(D))^(1/3)
+//         T := sgn(R - sqrt(D))*abs(R - sqrt(D))^(1/3)
+// 
+//         z1 := -a/3 + (S + T)               % real root
+//         z2 := -a/3 - (S + T)/2             % real part of complex root
+//         z3 := -a/3 - (S + T)/2             % real part of complex root
+//         im := abs(sqrt(3)*(S - T)/2)       % complex part of root pair
+// 
+//     else                                    % distinct real roots
+// 
+//         th := arccos(R/sqrt( -Q^3))
+//         
+//         z1 := 2*sqrt(-Q)*cos(th/3) - a/3
+//         z2 := 2*sqrt(-Q)*cos((th + 2*pi)/3) - a/3
+//         z3 := 2*sqrt(-Q)*cos((th + 4*pi)/3) - a/3
+//         im := 0
+// 
+//     end if
+// 
+//     return im                               % imaginary part
+// 
+// end function
+//
+Int_t MMath::SolvePol3(Double_t a, Double_t b, Double_t c,
+                       Double_t &x1, Double_t &x2, Double_t &x3)
+{
+    const Double_t Q = (a*a - 3*b)/9;
+    const Double_t R = (9*b*a - 27*c - 2*a*a*a)/54;
+    const Double_t D = R*R - Q*Q*Q;             // polynomial discriminant
+
+    // ----- The single-real / duplicate-roots solution -----
+
+    if (D==0)
+    {
+        const Double_t r = MMath::Sqrt3(R);
+
+        x1 = 2*r - a/3.;               // real root
+        x2 =   r - a/3.;               // real root
+
+        return 2;
+    }
+
+    if (D>0)                                    // complex or duplicate roots
+    {
+        const Double_t sqrtd = TMath::Sqrt(D);
+
+        const Double_t S = MMath::Sqrt3(R + sqrtd);
+        const Double_t T = MMath::Sqrt3(R - sqrtd);
+
+        x1 = (S+T) - a/3.;               // real root
+
+        return 1;
+
+        //z2 = (S + T)/2 - a/3.;            // real part of complex root
+        //z3 = (S + T)/2 - a/3.;            // real part of complex root
+        //im = fabs(sqrt(3)*(S - T)/2)      // complex part of root pair
+    }
+
+    // ----- The general solution with three roots ---
+
+    if (Q==0)
+        return 0;
+
+    if (Q>0) // This is here for speed reasons
+    {
+        const Double_t sqrtq = TMath::Sqrt(Q);
+        const Double_t rq    = R/TMath::Abs(Q);
+
+        const Double_t th1 = TMath::ACos(rq/sqrtq);
+        const Double_t th2 = th1 + TMath::TwoPi();
+        const Double_t th3 = th2 + TMath::TwoPi();
+
+        x1 = 2.*sqrtq * TMath::Cos(th1/3.) - a/3.;
+        x2 = 2.*sqrtq * TMath::Cos(th2/3.) - a/3.;
+        x3 = 2.*sqrtq * TMath::Cos(th3/3.) - a/3.;
+
+        return 3;
+    }
+
+    const TComplex sqrtq = TComplex::Sqrt(Q);
+    const Double_t rq    = R/TMath::Abs(Q);
+
+    const TComplex th1 = TComplex::ACos(rq/sqrtq);
+    const TComplex th2 = th1 + TMath::TwoPi();
+    const TComplex th3 = th2 + TMath::TwoPi();
+
+    // For ReMul, see bove
+    x1 = ReMul(2.*sqrtq, th1) - a/3.;
+    x2 = ReMul(2.*sqrtq, th2) - a/3.;
+    x3 = ReMul(2.*sqrtq, th3) - a/3.;
+
+    return 3;
+}

trunk/MagicSoft/Mars/mbase/MMath.h

@@ -47 +47 @@
     Double_t InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x);
 
+    inline Int_t SolvePol1(Double_t c, Double_t d, Double_t &x1)
+    {
+        if (c==0)
+            return 0;
+
+        x1 = -d/c;
+        return 1;
+    }
+    Int_t SolvePol2(Double_t c, Double_t d, Double_t &x1, Double_t &x2);
+    inline Int_t SolvePol2(Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2)
+    {
+        return b==0 ? SolvePol1(c, d, x1) : SolvePol2(c/b, d/b, x1, x2);
+    }
+    Int_t SolvePol3(Double_t b, Double_t c, Double_t d, Double_t &x1, Double_t &x2, Double_t &x3);
+    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)
+    {
+        return a==0 ? SolvePol2(b, c, d, x1, x2) : SolvePol3(b/a, c/a, d/a, x1, x2, x3);
+    }
+
     TArrayD LeastSqFitExpW1(Int_t n, Double_t *x, Double_t *y);
     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); }
 
+    inline Double_t Sqrt3(Double_t x) { return TMath::Sign(TMath::Power(TMath::Abs(x), 1./3), x); }
+
     inline Double_t Sgn(Double_t d) { return d<0 ? -1 : 1; }
 }

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

@@ -30 +30 @@
 //   Numerical Recipes in C++, 2nd edition, pp. 116-119.
 //   
-//   The coefficients "ya" are here denoted as "fHiGainSignal" and "fLoGainSignal"
-//   which means the FADC value subtracted by the clock-noise corrected pedestal.
+//   The coefficients "ya" are here denoted as "fVal" corresponding to
+//   the FADC value subtracted by the clock-noise corrected pedestal.
 //
 //   The coefficients "y2a" get immediately divided 6. and are called here 
-//   "fHiGainSecondDeriv" and "fLoGainSecondDeriv" although they are now not exactly 
-//   the second derivative coefficients any more. 
+//   fDer2 although they are now not exactly the second derivative
+//   coefficients any more.
 // 
 //   The calculation of the cubic-spline interpolated value "y" on a point 
-//   "x" along the FADC-slices axis becomes:
-// 
-//   y =    a*fHiGainSignal[klo] + b*fHiGainSignal[khi] 
-//       + (a*a*a-a)*fHiGainSecondDeriv[klo] + (b*b*b-b)*fHiGainSecondDeriv[khi]
-//
-//   with:
-//   a = (khi - x)
-//   b = (x - klo)
-//
-//   and "klo" being the lower bin edge FADC index and "khi" the upper bin edge FADC index.
-//   fHiGainSignal[klo] and fHiGainSignal[khi] are the FADC values at "klo" and "khi".
-//
-//   An analogues formula is used for the low-gain values.
-//
-//   The coefficients fHiGainSecondDeriv and fLoGainSecondDeriv are calculated with the 
-//   following simplified algorithm:
-//
-//   for (Int_t i=1;i<range-1;i++) {
-//       pp                   = fHiGainSecondDeriv[i-1] + 4.;
-//       fHiGainFirstDeriv[i] = fHiGainSignal[i+1] - 2.*fHiGainSignal[i] + fHiGainSignal[i-1]
-//       fHiGainFirstDeriv[i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
-//   }
-// 
-//   for (Int_t k=range-2;k>=0;k--)
-//       fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
-// 
-//
-//   This algorithm takes advantage of the fact that the x-values are all separated by exactly 1
-//   which simplifies the Numerical Recipes algorithm.
-//   (Note that the variables "fHiGainFirstDeriv" are not real first derivative coefficients.)
-//
-//   The algorithm to search the time proceeds as follows:
-//
-//   1) Calculate all fHiGainSignal from fHiGainFirst to fHiGainLast 
-//      (note that an "overlap" to the low-gain arrays is possible: i.e. fHiGainLast>14 in the case of 
-//      the MAGIC FADCs).
-//   2) Remember the position of the slice with the highest content "fAbMax" at "fAbMaxPos".
-//   3) If one or more slices are saturated or fAbMaxPos is less than 2 slices from fHiGainFirst, 
-//      return fAbMaxPos as time and fAbMax as charge (note that the pedestal is subtracted here).
-//   4) Calculate all fHiGainSecondDeriv from the fHiGainSignal array
-//   5) Search for the maximum, starting in interval fAbMaxPos-1 in steps of 0.2 till fAbMaxPos-0.2.
-//      If no maximum is found, go to interval fAbMaxPos+1. 
-//      --> 4 function evaluations
-//   6) Search for the absolute maximum from fAbMaxPos to fAbMaxPos+1 in steps of 0.2 
-//      --> 4 function  evaluations
-//   7) Try a better precision searching from new max. position fAbMaxPos-0.2 to fAbMaxPos+0.2
-//      in steps of 0.025 (83 psec. in the case of the MAGIC FADCs).
-//      --> 14 function evaluations
-//   8) If Time Extraction Type kMaximum has been chosen, the position of the found maximum is 
-//      returned, else:
-//   9) The Half Maximum is calculated. 
-//  10) fHiGainSignal is called beginning from fAbMaxPos-1 backwards until a value smaller than fHalfMax
-//      is found at "klo". 
-//  11) Then, the spline value between "klo" and "klo"+1 is halfed by means of bisection as long as 
-//      the difference between fHalfMax and spline evaluation is less than fResolution (default: 0.01).
-//      --> maximum 12 interations.
-//   
-//  The algorithm to search the charge proceeds as follows:
-//
-//  1) If Charge Type: kAmplitude was chosen, return the Maximum of the spline, found during the 
-//     time search.
-//  2) If Charge Type: kIntegral was chosen, sum the fHiGainSignal between:
-//     (Int_t)(fAbMaxPos - fRiseTimeHiGain) and 
-//     (Int_t)(fAbMaxPos + fFallTimeHiGain)
-//     (default: fRiseTime: 1.5, fFallTime: 4.5)
-//                                           sum the fLoGainSignal between:
-//     (Int_t)(fAbMaxPos - fRiseTimeHiGain*fLoGainStretch) and 
-//     (Int_t)(fAbMaxPos + fFallTimeHiGain*fLoGainStretch)
-//     (default: fLoGainStretch: 1.5)
-//       
-//  The values: fNumHiGainSamples and fNumLoGainSamples are set to:
-//  1) If Charge Type: kAmplitude was chosen: 1.
-//  2) If Charge Type: kIntegral was chosen: fRiseTimeHiGain + fFallTimeHiGain
-//                 or: fNumHiGainSamples*fLoGainStretch in the case of the low-gain
-//
-//  Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 
-//        to modify the ranges.
-// 
-//        Defaults: 
-//        fHiGainFirst =  2 
-//        fHiGainLast  =  14
-//        fLoGainFirst =  2 
-//        fLoGainLast  =  14
-//
-//  Call: SetResolution() to define the resolution of the half-maximum search.
-//        Default: 0.01
-//
-//  Call: SetRiseTime() and SetFallTime() to define the integration ranges 
-//        for the case, the extraction type kIntegral has been chosen.
-//
-//  Call: - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the 
-//          computation of the amplitude at the maximum (default) and extraction 
-//          the position of the maximum (default)
-//          --> no further function evaluation needed
-//        - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the 
-//          computation of the integral beneith the spline between fRiseTimeHiGain
-//          from the position of the maximum to fFallTimeHiGain after the position of 
-//          the maximum. The Low Gain is computed with half a slice more at the rising
-//          edge and half a slice more at the falling edge.
-//          The time of the half maximum is returned.
-//          --> needs one function evaluations but is more precise
-//        
+//   "x" along the FADC-slices axis becomes: EvalAt(x)
+//
+//   The coefficients fDer2 are calculated with the simplified
+//   algorithm in InitDerivatives.
+//
+//   This algorithm takes advantage of the fact that the x-values are all
+//   separated by exactly 1 which simplifies the Numerical Recipes algorithm.
+//   (Note that the variables fDer are not real first derivative coefficients.)
+//
 //////////////////////////////////////////////////////////////////////////////
 #include "MExtralgoSpline.h"
 
+#include "../mbase/MMath.h"
+
 using namespace std;
 
+// --------------------------------------------------------------------------
+//
+// Calculate the first and second derivative for the splie.
+//
+// The coefficients are calculated such that
+//   1) fVal[i] = Eval(i, 0)
+//   2) Eval(i-1, 1)==Eval(i, 0)
+//
+// In other words: The values with the index i describe the spline
+// between fVal[i] and fVal[i+1]
+//
 void MExtralgoSpline::InitDerivatives() const
 {
@@ -169 +89 @@
 }
 
-Float_t MExtralgoSpline::CalcIntegral(Float_t start, Float_t range) const
-{
+// --------------------------------------------------------------------------
+//
+// Returns the highest x value in [min;max[ at which the spline in
+// the bin i is equal to y
+//
+// min and max are defined to be [0;1]
+//
+// The default for min is 0, the default for max is 1
+// The defaule for y is 0
+//
+Double_t MExtralgoSpline::FindY(Int_t i, Double_t y, Double_t min, Double_t max) const
+{
+    // y = a*x^3 + b*x^2 + c*x + d'
+    // 0 = a*x^3 + b*x^2 + c*x + d' - y
+
+    // Calculate coefficients
+    const Double_t a = fDer2[i+1]-fDer2[i];
+    const Double_t b = 3*fDer2[i];
+    const Double_t c = fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1];
+    const Double_t d = fVal[i] - y;
+
+    Double_t x1, x2, x3;
+    const Int_t rc = MMath::SolvePol3(a, b, c, d, x1, x2, x3);
+
+    Double_t x = -1;
+    if (rc>0 && x1>=min && x1<max && x1>x)
+        x = x1;
+    if (rc>1 && x2>=min && x2<max && x2>x)
+        x = x2;
+    if (rc>2 && x3>=min && x3<max && x3>x)
+        x = x3;
+
+    return x<0 ? -1 : x+i;
+}
+
+// --------------------------------------------------------------------------
+//
+// Search analytically downward for the value y of the spline, starting
+// at x, until x==0. If y is not found -1 is returned.
+//
+Double_t MExtralgoSpline::SearchY(Float_t x, Float_t y) const
+{
+    if (x>=fNum-1)
+        x = fNum-1.0001;
+
+    Int_t i = TMath::FloorNint(x);
+    Double_t rc = FindY(i, y, 0, x-i);
+    while (--i>=0 && rc<0)
+        rc = FindY(i, y);
+
+    return rc;
+}
+
+// --------------------------------------------------------------------------
+//
+// Do a range check an then calculate the integral from start-fRiseTime
+// to start+fFallTime. An extrapolation of 0.5 slices is allowed.
+//
+Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
+{
+/*
     // The number of steps is calculated directly from the integration
     // 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
     // other machines!
+    const Float_t start = pos-fRiseTime;
     const Float_t step  = 0.2;
     const Float_t width = fRiseTime+fFallTime;
-    const Float_t max   = range-1 - (width+step);
+    const Float_t max   = fNum-1 - (width+step);
     const Int_t   num   = TMath::Nint(width/step);
 
@@ -193 +173 @@
     for (Int_t i=0; i<num; i++)
     {
-        const Float_t x = start+i*step;
-        const Int_t klo = (Int_t)TMath::Floor(x);
         // Note: if x is close to one integer number (= a FADC sample)
         // we get the same result by using that sample as klo, and the
@@ -203 +181 @@
         // depending on the compilation options).
 
-        sum += Eval(x, klo);
+        sum += EvalAt(start + i*step);
 
         // FIXME? Perhaps the integral should be done analitically
@@ -209 +187 @@
     }
     sum *= step; // Transform sum in integral
+
     return sum;
+    */
+
+    // In the future we will calculate the intgeral analytically.
+    // It has been tested that it gives identical results within
+    // acceptable differences.
+
+    // We allow extrapolation of 1/2 slice.
+    const Float_t min = fRiseTime;        //-0.5+fRiseTime;
+    const Float_t max = fNum-1-fFallTime; //fNum-0.5+fFallTime;
+
+    if (pos<min)
+        pos = min;
+    if (pos>max)
+        pos = max;
+
+    return EvalInteg(pos-fRiseTime, pos+fFallTime);
 }
 
@@ -217 +212 @@
 
     if (fExtractionType == kAmplitude)
-        return Eval(nsx+1, 1);
+        return Eval(1, nsx);
     else
-        return CalcIntegral(2. + nsx, fNum);
-}
-
-void MExtralgoSpline::Extract(Byte_t sat, Int_t maxpos)
+        return CalcIntegral(2. + nsx);
+}
+
+void MExtralgoSpline::Extract(Byte_t sat, Int_t maxbin)
 {
     fSignal    =  0;
@@ -233 +228 @@
     // Don't start if the maxpos is too close to the limits.
     //
-    const Bool_t limlo = maxpos <      TMath::Ceil(fRiseTime);
-    const Bool_t limup = maxpos > fNum-TMath::Ceil(fFallTime)-1;
+
+    /*
+    const Bool_t limlo = maxbin <      TMath::Ceil(fRiseTime);
+    const Bool_t limup = maxbin > fNum-TMath::Ceil(fFallTime)-1;
     if (sat || limlo || limup)
     {
@@ -240 +237 @@
         if (fExtractionType == kAmplitude)
         {
-            fSignal    = fVal[maxpos];
-            fTime      = maxpos;
+            fSignal    = fVal[maxbin];
+            fTime      = maxbin;
             fSignalDev = 0;  // means: is valid
             return;
         }
 
-        fSignal    = CalcIntegral(limlo ? 0 : fNum, fNum);
-        fTime      = maxpos - 1;
+        fSignal    = CalcIntegral(limlo ? 0 : fNum);
+        fTime      = maxbin - 1;
         fSignalDev = 0;  // means: is valid
         return;
     }
+    */
 
     fTimeDev = fResolution;
@@ -257 +255 @@
     // Now find the maximum
     //
+
+
+    /*
     Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
 
-    Int_t klo = maxpos-1;
-    Float_t fAbMaxPos = maxpos;//! Current position of the maximum of the spline
-    Float_t fAbMax = fVal[maxpos];//! Current maximum of the spline
+    Int_t klo = maxbin-1;
+
+    Float_t maxpos = maxbin;//! Current position of the maximum of the spline
+    Float_t max = fVal[maxbin];//! Current maximum of the spline
 
     //
@@ -270 +272 @@
     {
         const Float_t x = klo + step*(i+1);
-        const Float_t y = Eval(x, klo);
-        if (y > fAbMax)
-        {
-            fAbMax    = y;
-            fAbMaxPos = x;
+        //const Float_t y = Eval(klo, x);
+        const Float_t y = Eval(klo, x-klo);
+        if (y > max)
+        {
+            max    = y;
+            maxpos = x;
         }
     }
@@ -281 +284 @@
     // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2
     //
-    if (fAbMaxPos > maxpos - 0.1)
-    {
-        klo = maxpos;
+    if (maxpos > maxbin - 0.1)
+    {
+        klo = maxbin;
 
         for (Int_t i=0; i<TMath::Nint(TMath::Ceil((1-0.3)/step)); i++)
         {
             const Float_t x = klo + step*(i+1);
-            const Float_t y = Eval(x, klo);
-            if (y > fAbMax)
+            //const Float_t y = Eval(klo, x);
+            const Float_t y = Eval(klo, x-klo);
+            if (y > max)
             {
-                fAbMax    = y;
-                fAbMaxPos = x;
+                max    = y;
+                maxpos = x;
             }
         }
@@ -301 +305 @@
     // Try a better precision.
     //
-    const Float_t up = fAbMaxPos+step - 3.0*fResolution;
-    const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
-    const Float_t abmaxpos = fAbMaxPos;
+    const Float_t up = maxpos+step - 3.0*fResolution;
+    const Float_t lo = maxpos-step + 3.0*fResolution;
+    const Float_t abmaxpos = maxpos;
 
     step  = 2.*fResolution; // step size of 0.1 FADC slices
@@ -310 +314 @@
     {
         const Float_t x = abmaxpos + (i+1)*step;
-        const Float_t y = Eval(x, klo);
-        if (y > fAbMax)
-        {
-            fAbMax    = y;
-            fAbMaxPos = x;
+        //const Float_t y = Eval(klo, x);
+        const Float_t y = Eval(klo, x-klo);
+        if (y > max)
+        {
+            max    = y;
+            maxpos = x;
         }
     }
@@ -333 +338 @@
     {
         const Float_t x = abmaxpos - (i+1)*step;
-        const Float_t y = Eval(x, klo);
-        if (y > fAbMax)
-        {
-            fAbMax    = y;
-            fAbMaxPos = x;
-        }
-    }
+        //const Float_t y = Eval(klo, x);
+        const Float_t y = Eval(klo, x-klo);
+        if (y > max)
+        {
+            max    = y;
+            maxpos = x;
+        }
+    }
+  */
+
+    // --- Start NEW ---
+
+    // This block extracts values very similar to the old algorithm...
+    // for max>10
+    /* Most accurate (old identical) version:
+
+    Float_t xmax=maxpos, ymax=Eval(maxpos-1, 1);
+    Int_t rc = GetMaxPos(maxpos-1, xmax, ymax);
+    if (xmax==maxpos)
+    {
+        GetMaxPos(maxpos, xmax, ymax);
+
+        Float_t y = Eval(maxpos, 1);
+        if (y>ymax)
+        {
+            ymax = y;
+            xmax = maxpos+1;
+        }
+    }*/
+
+    Float_t maxpos, maxval;
+    GetMaxAroundI(maxbin, maxpos, maxval);
+
+    // --- End NEW ---
 
     if (fExtractionType == kAmplitude)
     {
-        fTime      = fAbMaxPos;
-        fSignal    = fAbMax;
+        fTime      = maxpos;
+        fSignal    = maxval;
         fSignalDev = 0;  // means: is valid
         return;
     }
 
-    Float_t fHalfMax = fAbMax/2.;//! Current half maximum of the spline
-
-    //
-    // Now, loop from the maximum bin leftward down in order to find the
-    // position of the half maximum. First, find the right FADC slice:
-    //
-    klo = maxpos;
-    while (klo > 0)
-    {
-        if (fVal[--klo] < fHalfMax)
+    //
+    // Loop from the beginning of the slice upwards to reach the maxhalf:
+    // With means of bisection:
+    //
+    /*
+    static const Float_t sqrt2 = TMath::Sqrt(2.);
+
+    step = sqrt2*3*0.061;//fRiseTime;
+    Float_t x = maxpos-0.86-3*0.061;//fRiseTime*1.25;
+
+//    step = sqrt2*0.5;//fRiseTime;
+//    Float_t x = maxpos-1.25;//fRiseTime*1.25;
+
+    Int_t  cnt  =0;
+    while (cnt++<30)
+    {
+        const Float_t y=EvalAt(x);
+
+        if (TMath::Abs(y-maxval/2)<fResolution)
             break;
-    }
-
-    //
-    // Loop from the beginning of the slice upwards to reach the fHalfMax:
-    // With means of bisection:
-    //
-    step = 0.5;
-
-    Int_t maxcnt = 20;
-    Int_t cnt    = 0;
-
-    Float_t y    = Eval(klo, klo); // FIXME: IS THIS CORRECT???????
-    Float_t x    = klo;
-    Bool_t back  = kFALSE;
-
-    while (TMath::Abs(y-fHalfMax) > fResolution)
-    {
-        x += back ? -step : +step;
-
-        const Float_t y = Eval(x, klo);
-
-        back = y > fHalfMax;
-
-        if (++cnt > maxcnt)
-            break;
-
-        step /= 2.;
-    }
+
+        step /= sqrt2; // /2
+        x += y>maxval/2 ? -step : +step;
+    }
+    */
+
+    // Search downwards for maxval/2
+    // By doing also a search upwards we could extract the pulse width
+    const Double_t x1 = SearchY(maxpos, maxval/2);
+
+    fTime      = x1;
+    fSignal    = CalcIntegral(maxpos);
+    fSignalDev = 0;  // means: is valid
+
+    //if (fSignal>100)
+    //    cout << "V=" << maxpos-x1 << endl;
 
     //
     // Now integrate the whole thing!
     //
-    fTime      = x;
-    fSignal    = CalcIntegral(fAbMaxPos - fRiseTime, fNum);
-    fSignalDev = 0;  // means: is valid
-}
+    //   fTime      = cnt==31 ? -1 : x;
+    //   fSignal    = cnt==31 ? CalcIntegral(x) : CalcIntegral(maxpos);
+    //   fSignalDev = 0;  // means: is valid
+}

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

@@ -5 +5 @@
 #include <TROOT.h>
 #endif
+
+#include <iostream>
+class TComplex;
 
 class MExtralgoSpline
@@ -35 +38 @@
     Float_t fSignalDev;
 
+    Double_t ReMul(const TComplex &c1, const TComplex &th) const;
+
     inline Float_t Eval(Float_t val, Float_t a, Float_t deriv) const
     {
@@ -40 +45 @@
     }
 
-    inline Float_t Eval(const Float_t x, const Int_t i) const
-    {
-        const Float_t b = x-i;
-        return Eval(fVal[i], 1-b, fDer2[i]) + Eval(fVal[i+1], b, fDer2[i+1]);
-    }
+    // Evaluate value of spline in the interval i with x=[0;1[
+    inline Float_t Eval(const Int_t i, const Float_t x) const
+    {
+        // Eval(i,x) =  (fDer2[i+1]-fDer2[i])*x*x*x + 3*fDer2[i]*x*x +
+        //              (fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1])*x + fVal[i];
+
+        // x := [0; 1[
+        return Eval(fVal[i], 1-x, fDer2[i]) + Eval(fVal[i+1], x, fDer2[i+1]);
+    }
+
+    /*
+    inline Float_t EvalAt(const Float_t x) const
+    {
+        Int_t i = TMath::FloorNint(x);
+
+        // handle under- and overflow of the array-range by extrapolation
+        if (i<0)
+            i=0;
+        if (i>fNum-2)
+            i = fNum-2;
+
+        return Eval(i, x-i);
+    }
+    */
+
+    // Evaluate first derivative of spline in the interval i with x=[0;1[
+    inline Double_t EvalDeriv1(const Float_t x, const Int_t i) const
+    {
+        // x := [0; 1[
+        const Double_t difval = fVal[i+1]-fVal[i];
+        const Double_t difder = fDer2[i+1]-fDer2[i];
+
+        return 3*difder*x*x + 6*fDer2[i]*x - 2*fDer2[i] - fDer2[i+1] + difval;
+    }
+
+    // Evaluate second derivative of spline in the interval i with x=[0;1[
+    inline Double_t EvalDeriv2(const Float_t x, const Int_t i) const
+    {
+        // x := [0; 1[
+        return 6*(fDer2[i+1]*x + fDer2[i]*(1-x));
+    }
+
+    Double_t FindY(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
+    Double_t SearchY(Float_t maxpos, Float_t y) const;
+/*
+    // Evaluate first solution for a possible maximum (x|first deriv==0)
+    inline Double_t EvalDerivEq0S1(const Int_t i) const
+    {
+        // return the x value [0;1[ at which the derivative is zero (solution1)
+
+        Double_t sumder = fDer2[i]+fDer2[i+1];
+        Double_t difder = fDer2[i]-fDer2[i+1];
+
+        Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
+        Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
+
+        Double_t x = 3*fDer2[i] - sqrt(3*sqt1 + 3*sqt2);
+
+        Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
+
+        return -x/denom;
+    }
+
+    // Evaluate second solution for a possible maximum (x|first deriv==0)
+    inline Double_t EvalDerivEq0S2(const Int_t i) const
+    {
+        // return the x value [0;1[ at which the derivative is zero (solution2)
+
+        Double_t sumder = fDer2[i]+fDer2[i+1];
+        Double_t difder = fDer2[i]-fDer2[i+1];
+
+        Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
+        Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
+
+        Double_t x = 3*fDer2[i] + sqrt(3*sqt1 + 3*sqt2);
+
+        Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
+
+        return -x/denom;
+    }
+    */
+
+    inline void EvalDerivEq0(const Int_t i, Float_t &rc1, Float_t &rc2) const
+    {
+        Double_t sumder = fDer2[i]+fDer2[i+1];
+        Double_t difder = fDer2[i]-fDer2[i+1];
+
+        Double_t sqt1  = sumder*sumder - fDer2[i]*fDer2[i+1];
+        Double_t sqt2  = difder*(fVal[i+1]-fVal[i]);
+        Double_t sqt3  = sqrt(3*sqt1 + 3*sqt2);
+        Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
+
+        rc1 = -(3*fDer2[i] + sqt3)/denom;
+        rc2 = -(3*fDer2[i] - sqt3)/denom;
+    }
+
+    // Calculate the "Stammfunktion" of the Eval-function
+    inline Double_t EvalPrimitive(Int_t i, Float_t x) const
+    {
+        /* TO BE CHECKED!
+        if (x==0)
+            return 0;
+
+        if (x==1)
+            return (fVal[i+1]+fVal[i])/2 - fDer2[i+1]/4;
+            */
+        Align(i, x);
+
+        const Double_t x2  = x*x;
+        const Double_t x4  = x2*x2;
+        const Double_t x1  = 1-x;
+        const Double_t x14 = x1*x1*x1*x1;
+
+        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];
+    }
+
+    inline void Align(Int_t &i, Float_t &x) const
+    {
+        if (i<0)
+        {
+            x += i;
+            i=0;
+        }
+        if (i>=fNum-1)
+        {
+            x += i-(fNum-2);
+            i=fNum-2;
+        }
+    }
+
+    // Calculate the intgeral of the Eval-function in
+    // bin i from a=[0;1[ to b=[0;1[
+    inline Double_t EvalInteg(Int_t i, Float_t a=0, Float_t b=1) const
+    {
+        // This is to make sure that we never access invalid
+        // memory, even if this should never happen.
+        // If it happens anyhow we extraolate the spline
+        Align(i, a);
+        Align(i, b);
+
+        return EvalPrimitive(i, b)-EvalPrimitive(i, a);
+    }
+
+    // Calculate the intgeral of the Eval-function betwen x0 and x1
+    inline Double_t EvalInteg(Float_t x0, Float_t x1) const
+    {
+        const Int_t min = TMath::CeilNint(x0);
+        const Int_t max = TMath::FloorNint(x1);
+
+        // This happens if x0 and x1 are in the same interval
+        if (min>max)
+            return EvalInteg(max, x0-max, x1-max);
+
+        // Sum complete intervals
+        Double_t sum = 0;
+        for (int i=min; i<max; i++)
+            sum += EvalInteg(i);
+
+        // Sum the incomplete intervals at the beginning and end
+        sum += EvalInteg(min-1, 1-(min-x0), 1);
+        sum += EvalInteg(max,   0, x1-max);
+
+        // return result
+        return sum;
+    }
+
+    // We search for the maximum from x=i-1 to x=i+1
+    // (Remeber: i corresponds to the value in bin i, i+1 to the
+    //  next bin and i-1 to the last bin)
+    inline void GetMaxAroundI(Int_t i, Float_t &xmax, Float_t &ymax) const
+    {
+        Float_t xmax1, xmax2;
+        Float_t ymax1, ymax2;
+
+        Bool_t rc1 = i>0      && GetMax(i-1, xmax1, ymax1);
+        Bool_t rc2 = i<fNum-1 && GetMax(i,   xmax2, ymax2);
+
+        // In case the medium bin is the first or last bin
+        // take the lower or upper edge of the region into account.
+        if (i==0)
+        {
+            xmax1 = 0;
+            ymax1 = fVal[0];
+            rc1 = kTRUE;
+        }
+        if (i>fNum-2)
+        {
+            xmax2 = fNum-1;
+            ymax2 = fVal[fNum-1];
+            rc2 = kTRUE;
+        }
+
+        // Take a default in case no maximum is found
+        xmax=i;
+        ymax=fVal[i];
+
+        if (rc1)
+        {
+            ymax = ymax1;
+            xmax = xmax1;
+        }
+        else
+            if (rc2)
+            {
+                ymax = ymax2;
+                xmax = xmax2;
+            }
+
+        if (rc2 && ymax2>ymax)
+        {
+            ymax = ymax2;
+            xmax = xmax2;
+        }
+  /*
+        // Search real maximum in [i-0.5, i+1.5]
+        Float_t xmax1, xmax2, xmax3;
+        Float_t ymax1, ymax2, ymax3;
+
+        Bool_t rc1 = i>0 && GetMax(i-1, xmax1, ymax1, 0.5, 1.0);
+        Bool_t rc2 = GetMax(i,   xmax2, ymax2, 0.0, 1.0);
+        Bool_t rc3 = i<fNum-1 && GetMax(i+1, xmax3, ymax3, 0.0, 0.5);
+
+        // In case the medium bin is the first or last bin
+        // take the lower or upper edge of the region into account.
+        if (i==0)
+        {
+            xmax1 = 0;
+            ymax1 = Eval(0, 0);
+            rc1 = kTRUE;
+        }
+        if (i==fNum-1)
+        {
+            xmax3 = fNum-1e-5;
+            ymax3 = Eval(fNum-1, 1);
+            rc3 = kTRUE;
+        }
+
+        // Take a real default in case no maximum is found
+        xmax=i+0.5;
+        ymax=Eval(i, 0.5);
+
+        //if (!rc1 && !rc2 && !rc3)
+        //    cout << "!!!!!!!!!!!!!!!" << endl;
+
+        if (rc1)
+        {
+            ymax = ymax1;
+            xmax = xmax1;
+        }
+        else
+            if (rc2)
+            {
+                ymax = ymax2;
+                xmax = xmax2;
+            }
+            else
+                if (rc3)
+                {
+                    ymax = ymax3;
+                    xmax = xmax3;
+                }
+
+        if (rc2 && ymax2>ymax)
+        {
+            ymax = ymax2;
+            xmax = xmax2;
+        }
+        if (rc3 && ymax3>ymax)
+        {
+            ymax = ymax3;
+            xmax = xmax3;
+        }
+*/    }
+
+    inline Bool_t GetMax(Int_t i, Float_t &xmax, Float_t &ymax, Float_t min=0, Float_t max=1) const
+    {
+        // Find analytical maximum in the bin i in the interval [min,max[
+
+        Float_t x1, x2;
+        EvalDerivEq0(i, x1, x2);
+        // const Float_t x1 = EvalDerivEq0S1(i);
+        // const Float_t x2 = EvalDerivEq0S2(i);
+
+        const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(x1, i)<0;
+        const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(x2, i)<0;
+
+        if (!ismax1 && !ismax2)
+            return kFALSE;
+
+        if (ismax1 && !ismax2)
+        {
+            xmax = i+x1;
+            ymax = Eval(i, x1);
+            return kTRUE;
+        }
+
+        if (!ismax1 && ismax2)
+        {
+            xmax = i+x2;
+            ymax = Eval(i, x2);
+            return kTRUE;
+        }
+
+        // Somehting must be wrong...
+        return kFALSE;
+        /*
+        std::cout << "?????????????" << std::endl;
+
+        const Double_t y1 = Eval(i, x1);
+        const Double_t y2 = Eval(i, x2);
+
+        if (y1>y2)
+        {
+            xmax = i+x1;
+            ymax = Eval(i, x1);
+            return kTRUE;
+        }
+        else
+        {
+            xmax = i+x2;
+            ymax = Eval(i, x2);
+            return kTRUE;
+        }
+
+        return kFALSE;*/
+    }
+/*
+    inline Int_t GetMaxPos(Int_t i, Float_t &xmax, Float_t &ymax) const
+    {
+        Double_t x[3];
+
+        x[0] = 0;
+        // x[1] = 1; // This means we miss a possible maximum at the
+                            // upper edge of the last interval...
+
+        x[1] = EvalDerivEq0S1(i);
+        x[2] = EvalDerivEq0S2(i);
+
+        //y[0] = Eval(i, x[0]);
+        //y[1] = Eval(i, x[1]);
+        //y[1] = Eval(i, x[1]);
+        //y[2] = Eval(i, x[2]);
+
+        Int_t rc = 0;
+        Double_t max = Eval(i, x[0]);
+
+        for (Int_t j=1; j<3; j++)
+        {
+            if (x[j]<=0 || x[j]>=1)
+                continue;
+
+            const Float_t y = Eval(i, x[j]);
+            if (y>max)
+            {
+                max = y;
+                rc = j;
+            }
+        }
+
+        if (max>ymax)
+        {
+            xmax = x[rc]+i;
+            ymax = max;
+        }
+
+        return rc;
+    }
+
+    inline void GetMaxPos(Int_t min, Int_t max, Float_t &xmax, Float_t &ymax) const
+    {
+        Float_t xmax=-1;
+        Float_t ymax=-FLT_MAX;
+
+        for (int i=min; i<max; i++)
+            GetMaxPos(i, xmax, ymax);
+
+        for (int i=min+1; i<max; i++)
+        {
+            Float_t y = Eval(i, 0);
+            if (y>ymax)
+            {
+                ymax = y;
+                xmax = i;
+            }
+        }
+
+    }*/
+
 
     void InitDerivatives() const;
-    Float_t CalcIntegral(Float_t start, Float_t range) const;
+    Float_t CalcIntegral(Float_t start) const;
 
 public: