Changeset 7999


Ignore:
Timestamp:
10/01/06 22:19:55 (18 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r7998 r7999  
    1818
    1919                                                 -*-*- 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
    2033 2006/09/29 Thomas Bretz
    2134
  • trunk/MagicSoft/Mars/mbase/MMath.cc

    r7899 r7999  
    3939#include <TArrayD.h>
    4040#endif
     41
     42#ifndef ROOT_TComplex
     43#include <TComplex.h>
     44#endif
     45
    4146// --------------------------------------------------------------------------
    4247//
     
    495500    return rc;
    496501}
     502
     503Int_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//
     527static 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//
     593Int_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  
    4747    Double_t InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x);
    4848
     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
    4968    TArrayD LeastSqFitExpW1(Int_t n, Double_t *x, Double_t *y);
    5069    TArrayD LeastSqFitExp(Int_t n, Double_t *x, Double_t *y);
     
    5473    inline Int_t ModF(Double_t dbl, Double_t &frac) { Double_t rc; frac = modf(dbl, &rc); return TMath::Nint(rc); }
    5574
     75    inline Double_t Sqrt3(Double_t x) { return TMath::Sign(TMath::Power(TMath::Abs(x), 1./3), x); }
     76
    5677    inline Double_t Sgn(Double_t d) { return d<0 ? -1 : 1; }
    5778}
  • trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

    r7942 r7999  
    3030//   Numerical Recipes in C++, 2nd edition, pp. 116-119.
    3131//   
    32 //   The coefficients "ya" are here denoted as "fHiGainSignal" and "fLoGainSignal"
    33 //   which means the 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.
    3434//
    3535//   The coefficients "y2a" get immediately divided 6. and are called here
    36 //   "fHiGainSecondDeriv" and "fLoGainSecondDeriv" although they are now not exactly
    37 //   the second derivative coefficients any more.
     36//   fDer2 although they are now not exactly the second derivative
     37//   coefficients any more.
    3838//
    3939//   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//
    14249//////////////////////////////////////////////////////////////////////////////
    14350#include "MExtralgoSpline.h"
    14451
     52#include "../mbase/MMath.h"
     53
    14554using namespace std;
    14655
     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//
    14767void MExtralgoSpline::InitDerivatives() const
    14868{
     
    16989}
    17090
    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//
     101Double_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//
     131Double_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//
     149Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
     150{
     151/*
    173152    // The number of steps is calculated directly from the integration
    174153    // window. This is the only way to ensure we are not dealing with
     
    176155    // value under the same conditions -- it might still be different on
    177156    // other machines!
     157    const Float_t start = pos-fRiseTime;
    178158    const Float_t step  = 0.2;
    179159    const Float_t width = fRiseTime+fFallTime;
    180     const Float_t max   = range-1 - (width+step);
     160    const Float_t max   = fNum-1 - (width+step);
    181161    const Int_t   num   = TMath::Nint(width/step);
    182162
     
    193173    for (Int_t i=0; i<num; i++)
    194174    {
    195         const Float_t x = start+i*step;
    196         const Int_t klo = (Int_t)TMath::Floor(x);
    197175        // Note: if x is close to one integer number (= a FADC sample)
    198176        // we get the same result by using that sample as klo, and the
     
    203181        // depending on the compilation options).
    204182
    205         sum += Eval(x, klo);
     183        sum += EvalAt(start + i*step);
    206184
    207185        // FIXME? Perhaps the integral should be done analitically
     
    209187    }
    210188    sum *= step; // Transform sum in integral
     189
    211190    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);
    212207}
    213208
     
    217212
    218213    if (fExtractionType == kAmplitude)
    219         return Eval(nsx+1, 1);
     214        return Eval(1, nsx);
    220215    else
    221         return CalcIntegral(2. + nsx, fNum);
    222 }
    223 
    224 void MExtralgoSpline::Extract(Byte_t sat, Int_t maxpos)
     216        return CalcIntegral(2. + nsx);
     217}
     218
     219void MExtralgoSpline::Extract(Byte_t sat, Int_t maxbin)
    225220{
    226221    fSignal    =  0;
     
    233228    // Don't start if the maxpos is too close to the limits.
    234229    //
    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;
    237234    if (sat || limlo || limup)
    238235    {
     
    240237        if (fExtractionType == kAmplitude)
    241238        {
    242             fSignal    = fVal[maxpos];
    243             fTime      = maxpos;
     239            fSignal    = fVal[maxbin];
     240            fTime      = maxbin;
    244241            fSignalDev = 0;  // means: is valid
    245242            return;
    246243        }
    247244
    248         fSignal    = CalcIntegral(limlo ? 0 : fNum, fNum);
    249         fTime      = maxpos - 1;
     245        fSignal    = CalcIntegral(limlo ? 0 : fNum);
     246        fTime      = maxbin - 1;
    250247        fSignalDev = 0;  // means: is valid
    251248        return;
    252249    }
     250    */
    253251
    254252    fTimeDev = fResolution;
     
    257255    // Now find the maximum
    258256    //
     257
     258
     259    /*
    259260    Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
    260261
    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
    264266
    265267    //
     
    270272    {
    271273        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;
    277280        }
    278281    }
     
    281284    // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2
    282285    //
    283     if (fAbMaxPos > maxpos - 0.1)
    284     {
    285         klo = maxpos;
     286    if (maxpos > maxbin - 0.1)
     287    {
     288        klo = maxbin;
    286289
    287290        for (Int_t i=0; i<TMath::Nint(TMath::Ceil((1-0.3)/step)); i++)
    288291        {
    289292            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)
    292296            {
    293                 fAbMax    = y;
    294                 fAbMaxPos = x;
     297                max    = y;
     298                maxpos = x;
    295299            }
    296300        }
     
    301305    // Try a better precision.
    302306    //
    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;
    306310
    307311    step  = 2.*fResolution; // step size of 0.1 FADC slices
     
    310314    {
    311315        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;
    317322        }
    318323    }
     
    333338    {
    334339        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 ---
    342374
    343375    if (fExtractionType == kAmplitude)
    344376    {
    345         fTime      = fAbMaxPos;
    346         fSignal    = fAbMax;
     377        fTime      = maxpos;
     378        fSignal    = maxval;
    347379        fSignalDev = 0;  // means: is valid
    348380        return;
    349381    }
    350382
    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)
    361402            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;
    390419
    391420    //
    392421    // Now integrate the whole thing!
    393422    //
    394     fTime      = x;
    395     fSignal    = CalcIntegral(fAbMaxPos - fRiseTime, fNum);
    396     fSignalDev = 0;  // means: is valid
    397 }
     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  
    55#include <TROOT.h>
    66#endif
     7
     8#include <iostream>
     9class TComplex;
    710
    811class MExtralgoSpline
     
    3538    Float_t fSignalDev;
    3639
     40    Double_t ReMul(const TComplex &c1, const TComplex &th) const;
     41
    3742    inline Float_t Eval(Float_t val, Float_t a, Float_t deriv) const
    3843    {
     
    4045    }
    4146
    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
    47435
    48436    void InitDerivatives() const;
    49     Float_t CalcIntegral(Float_t start, Float_t range) const;
     437    Float_t CalcIntegral(Float_t start) const;
    50438
    51439public:
Note: See TracChangeset for help on using the changeset viewer.