Ignore:
Timestamp:
10/01/06 22:19:55 (18 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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}
Note: See TracChangeset for help on using the changeset viewer.