Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 8545) +++ /trunk/MagicSoft/Mars/NEWS (revision 8546) @@ -7,4 +7,8 @@ Also the old values seemed not exactly the PulsePos used for teh check. + + - callisto: improved calculation of spline coefficients a lot. This + leads to a further improvement of the event rate calibrating MUX + data of about 15% (175evt/s instead of 150evt/s) Index: /trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc (revision 8545) +++ /trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc (revision 8546) @@ -193,46 +193,4 @@ 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 - // numerical rounding uncertanties, because we always get the same - // 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 = fNum-1 - (width+step); - const Int_t num = TMath::Nint(width/step); - - // The order is important. In some cases (limlo-/limup-check) it can - // happen that max<0. In this case we start at 0 - if (start > max) - start = max; - if (start < 0) - start = 0; - - start += step/2; - - Double_t sum = 0.; - for (Int_t i=0; i fNum-TMath::Ceil(fFallTime)-1; - if (sat || limlo || limup) - { - fTimeDev = 1.0; - if (fExtractionType == kAmplitude) - { - fSignal = fVal[maxbin]; - fTime = maxbin; - fSignalDev = 0; // means: is valid - return; - } - - fSignal = CalcIntegral(limlo ? 0 : fNum); - fTime = maxbin - 1; - fSignalDev = 0; // means: is valid - return; - } - - // - // 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 = maxbin-1; - - Float_t maxpos = maxbin;//! Current position of the maximum of the spline - Float_t max = fVal[maxbin];//! Current maximum of the spline - - // - // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2. - // If no maximum is found, go to interval maxpos+1. - // - for (Int_t i=0; i max) - { - max = y; - maxpos = x; - } - } - - // - // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2 - // - if (maxpos > maxbin - 0.1) - { - klo = maxbin; - - for (Int_t i=0; i max) - { - max = y; - maxpos = x; - } - } - } - - // - // Now, the time, abmax and khicont and klocont are set correctly within the previous precision. - // Try a better precision. - // - 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 - - for (int i=0; i max) - { - max = y; - maxpos = x; - } - } - - // - // Second, try from time down to time-0.2 in steps of fResolution. - // - - // - // Test the possibility that the absolute maximum has not been found between - // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos - // which requires new setting of klocont and khicont - // - if (abmaxpos < klo + fResolution) - klo--; - - for (int i=TMath::Nint(TMath::Ceil((abmaxpos-lo)/step))-1; i>=0; i--) - { - const Float_t x = abmaxpos - (i+1)*step; - //const Float_t y = Eval(klo, x); - const Float_t y = Eval(klo, x-klo); - if (y > max) - { - max = y; - maxpos = x; - } - } - - fTime = maxpos; - fTimeDev = fResolution; - fSignal = CalcIntegral(maxpos); - fSignalDev = 0; // means: is valid - - return; -*/ - // --- 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; @@ -461,37 +275,3 @@ fWidthDev = 0; } - - // - // 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)maxval/2 ? -step : +step; - } - */ - - // - // Now integrate the whole thing! - // - // fTime = cnt==31 ? -1 : x; - // fTimeDev = fResolution; - // 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 8545) +++ /trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (revision 8546) @@ -77,41 +77,4 @@ Double_t SearchY(Float_t maxpos, Float_t y) const; Double_t SearchYup(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; - } - */ Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const; @@ -289,64 +252,4 @@ 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; - } - */ } @@ -359,7 +262,4 @@ return kFALSE; - // 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;