Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 6979) +++ /trunk/MagicSoft/Mars/Changelog (revision 6980) @@ -78,4 +78,9 @@ * mmuon/Makefile, mmuon/MuonLinkDef.h: - added MHSingleMuon + + * msignal/MExtractTimeAndchargeSpline.cc: + - fixed a bug found by Abelardo which could cause different + results with different compiler optimiztaion options due + to rounding errors. Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 6979) +++ /trunk/MagicSoft/Mars/NEWS (revision 6980) @@ -39,4 +39,8 @@ - added muon support to star + + - fixed a bug in MExtractTimeAndchargeSpline which could cause + different results with different compiler optimiztaion options + due to rounding errors. Index: /trunk/MagicSoft/Mars/callisto_Dec04Jan05.rc =================================================================== --- /trunk/MagicSoft/Mars/callisto_Dec04Jan05.rc (revision 6979) +++ /trunk/MagicSoft/Mars/callisto_Dec04Jan05.rc (revision 6980) @@ -270,5 +270,5 @@ #MJCalibration.MHCalibrationChargeCam.NumHiGainSaturationLimit: 0.02 #MJCalibration.MHCalibrationChargeCam.NumLoGainSaturationLimit: 0.005 -#MJCalibration.MHCalibrationChargeCam.ProbLimit: 0.00000001 +MJCalibration.MHCalibrationChargeCam.ProbLimit: 1.0E-18 #MJCalibration.MHCalibrationChargeCam.OverflowLimit: 0.005 #MJCalibration.MHCalibrationChargeCam.PulserFrequency: 500 Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 6979) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 6980) @@ -21,4 +21,5 @@ ! \* ======================================================================== */ + ////////////////////////////////////////////////////////////////////////////// // @@ -160,4 +161,5 @@ const Float_t MExtractTimeAndChargeSpline::fgOffsetLoGain = 1.7; // 5 ns const Float_t MExtractTimeAndChargeSpline::fgLoGainStartShift = -1.8; + // -------------------------------------------------------------------------- // @@ -1060,118 +1062,81 @@ void MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t &sum, Float_t start, Float_t last) { - - const Float_t step = 0.2; - - if (start < 0) - { - last -= start; - start = 0.; - } - - Int_t klo = int(start); - Int_t khi = klo+1; - - Float_t lo = TMath::Floor(start); - Float_t up = lo + 1.; - - const Int_t m = int((start-klo)/step); - start = step*m + klo; // Correct start for the digitization due to resolution - - Float_t x = start; - Float_t a = up-start; - Float_t b = start-lo; - - while (1) - { - - while (x last) - { - sum *= step; - return; - } - - a -= step; - b += step; - - sum += a*fHiGainSignal[klo] - + b*fHiGainSignal[khi] - + (a*a*a-a)*fHiGainSecondDeriv[klo] + const Float_t step = 0.2; + + if (start < 0) + { + last -= start; + start = 0.; + } + + const Int_t n = TMath::Nint((last-start)/step); + + sum = 0.; + for (Int_t i=1; i<=n+1; i++) + { + const Float_t x = start+i*step; + const Int_t klo = (Int_t)TMath::Floor(x); + const Int_t khi = klo + 1; + // Note: if x is close to one integer number (= a FADC sample) + // we get the same result by using that sample as klo, and the + // next one as khi, or using the sample as khi and the previous + // one as klo (the spline is of course continuous). So we do not + // expect problems from rounding issues in the argument of + // Floor() above (we have noticed differences in roundings + // depending on the compilation options). + + const Float_t a = khi - x; // Distance from x to next FADC sample + const Float_t b = x - klo; // Distance from x to previous FADC sample + + sum += a*fHiGainSignal[klo] + + b*fHiGainSignal[khi] + + (a*a*a-a)*fHiGainSecondDeriv[klo] + (b*b*b-b)*fHiGainSecondDeriv[khi]; - } - - up += 1.; - lo += 1.; - klo++; - khi++; - start += 1.; - a = 1.; - b = 0.; - } - + + // FIXME? Perhaps the integral should be done analitically + // between every two FADC slices, instead of numerically + } + sum *= step; // Transform sum in integral } + void MExtractTimeAndChargeSpline::CalcIntegralLoGain(Float_t &sum, Float_t start, Float_t last) { - - const Float_t step = 0.1; - - if (start < 0) - { - last -= start; - start = 0.; - } - - Int_t klo = int(start); - Int_t khi = klo+1; - - Float_t lo = TMath::Floor(start); - Float_t up = lo + 1.; - - const Int_t m = int((start-klo)/step); - start = step*m + klo; // Correct start for the digitization due to resolution - - Float_t x = start; - Float_t a = up-start; - Float_t b = start-lo; - - while (1) - { - - while (x last) - { - sum *= step; - return; - } - - a -= step; - b += step; - - sum += a*fLoGainSignal[klo] - + b*fLoGainSignal[khi] - + (a*a*a-a)*fLoGainSecondDeriv[klo] + const Float_t step = 0.2; + + if (start < 0) + { + last -= start; + start = 0.; + } + + const Int_t n = TMath::Nint((last-start)/step); + + sum = 0.; + for (Int_t i=1; i<=n+1; i++) + { + const Float_t x = start+i*step; + const Int_t klo = (Int_t)TMath::Floor(x); + const Int_t khi = klo + 1; + // Note: if x is close to one integer number (= a FADC sample) + // we get the same result by using that sample as klo, and the + // next one as khi, or using the sample as khi and the previous + // one as klo (the spline is of course continuous). So we do not + // expect problems from rounding issues in the argument of + // Floor() above (we have noticed differences in roundings + // depending on the compilation options). + + const Float_t a = khi - x; // Distance from x to next FADC sample + const Float_t b = x - klo; // Distance from x to previous FADC sample + + sum += a*fLoGainSignal[klo] + + b*fLoGainSignal[khi] + + (a*a*a-a)*fLoGainSecondDeriv[klo] + (b*b*b-b)*fLoGainSecondDeriv[khi]; - - } - - up += 1.; - lo += 1.; - klo++; - khi++; - start += 1.; - a = 1.; - b = 0.; - } - + + // FIXME? Perhaps the integral should be done analitically + // between every two FADC slices, instead of numerically + } + sum *= step; // Transform sum in integral } - - - // --------------------------------------------------------------------------