Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 7054) +++ /trunk/MagicSoft/Mars/Changelog (revision 7055) @@ -41,4 +41,17 @@ * mjobs/MJStar.cc: - moved processing of CC-branch to the beginning of the tasklist + + * msignal/MExtractTimeAndChargeSpline.[h,cc]: + - introduced some small changes to the validity range of + some variables + - determin the higher bound above which no search is done + analog to the lower bound using the fall-time + - CalcIntegral[Hi,Lo]Gain now returns sum. No need for a reference + - fixed calling Integral[HI,Lo]Gain in cases we are at the edge of + the valid range -- at a lot of position in the code random memory + above or below the arrays have been accessed. + - improved the numercila stability of CalcIntegral[Hi,Lo]Gain + more by calculating the number of steps from the rise and fall time. + this should at least give consistent results on the same machine! Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 7054) +++ /trunk/MagicSoft/Mars/NEWS (revision 7055) @@ -77,4 +77,9 @@ adding the step-size (which results in numerical uncertanties exspecially if multiplied with large numbers) + + A lot of fixes have been introduced which effects integrating the + spline at the edges of the valid range. In this case any memory + was randomly accessed. + ! No result obtained with the Spline before can be trusted! Due to + random memory access it might by completely random! - callisto: set new defaults in MExtractTimeAndChargeDigitalFilter: Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 7054) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 7055) @@ -432,8 +432,6 @@ } } - + fAbMax = fHiGainSignal[maxpos]; - - Float_t pp; fHiGainSecondDeriv[0] = 0.; @@ -442,7 +440,7 @@ for (Int_t i=1;i range) - { - const Int_t diff = range - int(last); - last -= diff; - start -= diff; - } - - CalcIntegralHiGain(sum, start, last); - } + sum = CalcIntegralHiGain(2. + nsx, range); + fRandomIter++; return; @@ -492,11 +479,11 @@ // - // Allow no saturated slice - // and + // Allow no saturated slice and // Don't start if the maxpos is too close to the limits. // - if (sat || maxpos < TMath::Ceil(fRiseTimeHiGain) || maxpos > range-2) - { - + const Bool_t limlo = maxpos < TMath::Ceil(fRiseTimeHiGain); + const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeHiGain)-1; + if (sat || limlo || limup) + { dtime = 1.0; if (fExtractionType == kAmplitude) @@ -506,11 +493,7 @@ return; } - - if (maxpos > range - 2) - CalcIntegralHiGain(sum, (Float_t)range - fRiseTimeHiGain - fFallTimeHiGain, (Float_t)range - 0.001); - else - CalcIntegralHiGain(sum, 0.001, fRiseTimeHiGain + fFallTimeHiGain); - - time = (Float_t)(fHiGainFirst + maxpos - 1); + + sum = CalcIntegralHiGain(limlo ? 0 : range, range); + time = (Float_t)(fHiGainFirst + maxpos - 1); return; } @@ -561,5 +544,4 @@ if (fAbMaxPos > upper-0.1) { - upper = 1. + maxpos; lower = (Float_t)maxpos; @@ -678,6 +660,5 @@ while (klo > 0) { - klo--; - if (fHiGainSignal[klo] < fHalfMax) + if (fHiGainSignal[--klo] < fHalfMax) break; } @@ -719,8 +700,5 @@ + (b*b*b-b)*fHiGainSecondDeriv[khi]; - if (y > fHalfMax) - back = kTRUE; - else - back = kFALSE; + back = y > fHalfMax; if (++cnt > maxcnt) @@ -730,21 +708,9 @@ } - time = (Float_t)fHiGainFirst + x; // // Now integrate the whole thing! // - - Float_t start = fAbMaxPos - fRiseTimeHiGain; - Float_t last = fAbMaxPos + fFallTimeHiGain; - - const Int_t diff = int(last) - range; - - if (diff > 0) - { - last -= diff; - start -= diff; - } - - CalcIntegralHiGain(sum, start, last); + time = (Float_t)fHiGainFirst + x; + sum = CalcIntegralHiGain(fAbMaxPos - fRiseTimeHiGain, range); } @@ -795,6 +761,4 @@ fAbMax = fLoGainSignal[maxpos]; - Float_t pp; - fLoGainSecondDeriv[0] = 0.; fLoGainFirstDeriv[0] = 0.; @@ -802,7 +766,7 @@ for (Int_t i=1;i range) - { - const Int_t diff = range - int(last); - last -= diff; - start -= diff; - } - - CalcIntegralLoGain(sum, start, last); - } + sum = CalcIntegralLoGain(2. + nsx, range); + fRandomIter++; return; } // - // Allow no saturated slice - // and + // Allow no saturated slice and // Don't start if the maxpos is too close to the limits. // - if (sat || maxpos < TMath::Ceil(fRiseTimeLoGain) || maxpos > range-2) - { - + const Bool_t limlo = maxpos < TMath::Ceil(fRiseTimeLoGain); + const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeLoGain)-1; + if (sat || limlo || limup) + { dtime = 1.0; if (fExtractionType == kAmplitude) @@ -864,10 +817,6 @@ return; } - - if (maxpos > range-2) - CalcIntegralLoGain(sum, (Float_t)range - fRiseTimeLoGain - fFallTimeLoGain -1., (Float_t)range - 0.001); - else - CalcIntegralLoGain(sum, 0.001, fRiseTimeLoGain + fFallTimeLoGain + 1.); - + + sum = CalcIntegralLoGain(limlo ? 0 : range, range); time = (Float_t)(fLoGainFirst + maxpos - 1); return; @@ -1080,8 +1029,5 @@ + (b*b*b-b)*fLoGainSecondDeriv[khi]; - if (y > fHalfMax) - back = kTRUE; - else - back = kFALSE; + back = y > fHalfMax; if (++cnt > maxcnt) @@ -1091,38 +1037,34 @@ } - time = x + (Int_t)fLoGainFirst; - // // Now integrate the whole thing! // - Float_t start = fAbMaxPos - fRiseTimeLoGain; - Float_t last = fAbMaxPos + fFallTimeLoGain; - - const Int_t diff = int(last) - range; - - if (diff > 0) - { - last -= diff; - start -= diff; - } - CalcIntegralLoGain(sum, start, last); + time = x + (Int_t)fLoGainFirst; + sum = CalcIntegralLoGain(fAbMaxPos - fRiseTimeLoGain, range); } -void MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t &sum, Float_t start, Float_t last) +Float_t MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t start, Float_t range) const { - const Float_t step = 0.2; - + // 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 step = 0.2; + const Float_t width = fRiseTimeHiGain+fFallTimeHiGain; + const Float_t max = range-1 - (width+step); + const Int_t num = TMath::Nint(width/step); + + // The order is important. In some cases (limlo-/limup-check) it can + // happen than max<0. In this case we start at 0 + if (start > max) + start = max; if (start < 0) - { - last -= start; - start = 0.; - } - - const Int_t n = TMath::Nint((last-start)/step); + start = 0; start += step/2; - sum = 0.; - for (Int_t i=0; i max) + start = max; if (start < 0) - { - last -= start; - start = 0.; - } - - const Int_t n = TMath::Nint((last-start)/step); + start = 0; start += step/2; - sum = 0.; - for (Int_t i=0; i