Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 5218) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 5221) @@ -15,5 +15,5 @@ ! * ! -! Author(s): Markus Gaug 05/2004 +! Author(s): Markus Gaug 09/2004 ! ! Copyright: MAGIC Software Development, 2002-2004 @@ -26,18 +26,120 @@ // MExtractTimeAndChargeSpline // -// Fast Spline extractor using a cubic spline algorithm of Numerical Recipes. -// It returns the integral below the interpolating spline. +// Fast Spline extractor using a cubic spline algorithm, adapted from +// 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 "y2a" get immediately divided 6. and are called here +// "fHiGainSecondDeriv" and "fLoGainSecondDeriv" although they are now not exactly +// the second derivative coefficients any more. // -// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) -// to modify the ranges. Ranges have to be an even number. In case of odd -// ranges, the last slice will be reduced by one. -// -// Defaults are: +// The calculation of the cubic-spline interpolated value "y" on a point +// "x" along the FADC-slices axis becomes: // -// fHiGainFirst = fgHiGainFirst = 3 -// fHiGainLast = fgHiGainLast = 14 -// fLoGainFirst = fgLoGainFirst = 3 -// fLoGainLast = fgLoGainLast = 14 -// +// 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 - fRiseTime) and +// (Int_t)(fAbMaxPos + fFallTime) +// (default: fRiseTime: 1.5, fFallTime: 4.5) +// 3) Sum only half the values of the edge slices +// 4) Sum 1.5*fHiGainSecondDeriv of the not-edge slices using the "natural cubic +// spline with second derivatives set to 0. at the edges. +// (Remember that fHiGainSecondDeriv had been divided by 6.) +// +// The values: fNumHiGainSamples and fNumLoGainSamples are set to: +// 1) If Charge Type: kAmplitude was chosen: 1. +// 2) If Charge Type: kIntegral was chosen: TMath::Floor(fRiseTime + fFallTime) +// +// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) +// to modify the ranges. +// +// Defaults: +// fHiGainFirst = 2 +// fHiGainLast = 14 +// fLoGainFirst = 3 +// 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: - SetTimeType(MExtractTimeAndChargeSpline::kMaximum) for extraction +// the position of the maximum (default) +// --> needs 22 function evaluations +// - SetTimeType(MExtractTimeAndChargeSpline::kHalfMaximum) for extraction +// the position of the half maximum at the rising edge. +// --> needs max. 34 function evaluations +// - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the +// computation of the amplitude at the maximum (default) +// --> no further function evaluation needed +// - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the +// computation of the integral beneith the spline between fRiseTime +// from the position of the maximum to fFallTime after the position of +// the maximum. +// --> needs one more simple summation loop over 7 slices. +// ////////////////////////////////////////////////////////////////////////////// #include "MExtractTimeAndChargeSpline.h" @@ -60,5 +162,5 @@ const Byte_t MExtractTimeAndChargeSpline::fgLoGainFirst = 3; const Byte_t MExtractTimeAndChargeSpline::fgLoGainLast = 14; -const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.001; +const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.01; const Float_t MExtractTimeAndChargeSpline::fgRiseTime = 2.0; const Float_t MExtractTimeAndChargeSpline::fgFallTime = 4.0; @@ -89,14 +191,9 @@ SetFallTime(); + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); + SetTimeType(); SetChargeType(); - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); - - fNumHiGainSamples = 1.; - fNumLoGainSamples = 1.; - fSqrtHiGainSamples = 1.; - fSqrtLoGainSamples = 1.; - } @@ -119,4 +216,72 @@ } +//------------------------------------------------------------------- +// +// Set the ranges +// In order to set the fNum...Samples variables correctly for the case, +// the integral is computed, have to overwrite this function and make an +// explicit call to SetChargeType(). +// +void MExtractTimeAndChargeSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast) +{ + + MExtractor::SetRange(hifirst, hilast, lofirst, lolast); + + if (IsExtractionType(kIntegral)) + SetChargeType(kIntegral); + +} + + +//------------------------------------------------------------------- +// +// Set the Time Extraction type. Possible are: +// - kMaximum: Search the maximum of the spline and return its position +// - kHalfMaximum: Search the half maximum left from the maximum and return +// its position +// +void MExtractTimeAndChargeSpline::SetTimeType( ExtractionType_t typ ) +{ + + CLRBIT(fFlags,kMaximum); + CLRBIT(fFlags,kHalfMaximum); + SETBIT(fFlags,typ); + +} + +//------------------------------------------------------------------- +// +// Set the Charge Extraction type. Possible are: +// - kAmplitude: Search the value of the spline at the maximum +// - kIntegral: Integral the spline from fHiGainFirst to fHiGainLast, +// by counting the edge bins only half and setting the +// second derivative to zero, there. +// +void MExtractTimeAndChargeSpline::SetChargeType( ExtractionType_t typ ) +{ + + CLRBIT(fFlags,kAmplitude); + CLRBIT(fFlags,kIntegral ); + + SETBIT(fFlags,typ); + + if (IsExtractionType(kAmplitude)) + { + fNumHiGainSamples = 1.; + fNumLoGainSamples = 1.; + fSqrtHiGainSamples = 1.; + fSqrtLoGainSamples = 1.; + } + + if (IsExtractionType(kIntegral)) + { + fNumHiGainSamples = TMath::Floor(fRiseTime + fFallTime); + fNumLoGainSamples = fNumHiGainSamples; + fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); + fSqrtLoGainSamples = fSqrtHiGainSamples; + } +} + + // -------------------------------------------------------------------------- // @@ -275,7 +440,9 @@ if (sat || maxpos < 2) { - time = IsTimeType(kMaximum) + time = IsExtractionType(kMaximum) ? (Float_t)(fHiGainFirst + maxpos) : (Float_t)(fHiGainFirst + maxpos - 1); + sum = IsExtractionType(kAmplitude) + ? fAbMax : 0.; return; } @@ -295,4 +462,5 @@ fHiGainSecondDeriv[range-1] = 0.; + for (Int_t k=range-2;k>=0;k--) fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.; @@ -338,6 +506,5 @@ // - // Test the possibility that the absolute maximum has not been found before the - // maxpos and search from maxpos to maxpos+1 in steps of 0.2 + // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2 // if (fAbMaxPos > upper-0.1) @@ -379,6 +546,6 @@ // Try a better precision. // - const Float_t up = fAbMaxPos+step-0.055; - const Float_t lo = fAbMaxPos-step+0.055; + const Float_t up = fAbMaxPos+step-0.035; + const Float_t lo = fAbMaxPos-step+0.035; const Float_t maxpossave = fAbMaxPos; @@ -387,5 +554,5 @@ b = x - lower; - step = 0.02; // step size of 42 ps + step = 0.025; // step size of 83 ps while (x range) + lastslice = range; + Int_t i = startslice; sum = 0.5*fHiGainSignal[i]; + // + // We sum 1.5 times the second deriv. coefficients because these had been + // divided by 6. already. Usually, 0.25*fHiGainSecondDeriv should be added. + // for (i=startslice+1; i