Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc (revision 7999) @@ -30,119 +30,39 @@ // 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=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,6 +89,65 @@ } -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 && x1x) + x = x1; + if (rc>1 && x2>=min && x2x) + x = x2; + if (rc>2 && x3>=min && x3x) + 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,7 +155,8 @@ // 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,6 +173,4 @@ for (Int_t i=0; imax) + pos = max; + + return EvalInteg(pos-fRiseTime, pos+fFallTime); } @@ -217,10 +212,10 @@ 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,6 +228,8 @@ // 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,15 +237,16 @@ 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,9 +255,13 @@ // 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,9 +272,10 @@ { 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,16 +284,17 @@ // 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 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,7 +305,7 @@ // 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,9 +314,10 @@ { 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,65 +338,89 @@ { 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) - { - 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 +} Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (revision 7999) @@ -5,4 +5,7 @@ #include #endif + +#include +class TComplex; class MExtralgoSpline @@ -35,4 +38,6 @@ 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,12 +45,395 @@ } - 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; i0 && GetMax(i-1, xmax1, ymax1); + Bool_t rc2 = ifNum-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 = iymax) + { + 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=min && x2y2) + { + 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; iymax) + { + 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: