CSCRecHitD/src/CSCFindPeakTime.cc

0001 // This is CSCFindPeakTime
0002
0003 #include <RecoLocalMuon/CSCRecHitD/src/CSCFindPeakTime.h>
0004 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0005
0006 #include <cmath>
0007
0008 CSCFindPeakTime::CSCFindPeakTime(const edm::ParameterSet& ps)
0009     : useAverageTime(false), useParabolaFit(false), useFivePoleFit(false) {
0010   useAverageTime = ps.getParameter<bool>("UseAverageTime");
0011   useParabolaFit = ps.getParameter<bool>("UseParabolaFit");
0012   useFivePoleFit = ps.getParameter<bool>("UseFivePoleFit");
0013   LogTrace("CSCRecHit") << "[CSCFindPeakTime] useAverageTime=" << useAverageTime
0014                         << ", useParabolaFit=" << useParabolaFit << ", useFivePoleFit=" << useFivePoleFit;
0015 }
0016
0017 float CSCFindPeakTime::peakTime(int tmax, const float* adc, float t_peak) {
0018   if (useAverageTime) {
0019     return averageTime(tmax, adc);
0020   } else if (useParabolaFit) {
0021     return parabolaFitTime(tmax, adc);
0022   } else if (useFivePoleFit) {
0023     return fivePoleFitTime(tmax, adc, t_peak);
0024   } else {
0025     // return something, anyway.. may as well be average
0026     return averageTime(tmax, adc);
0027   }
0028 }
0029
0030 float CSCFindPeakTime::averageTime(int tmax, const float* adc) {
0031   float sum = 0.;
0032   float sumt = 0.;
0033   for (size_t i = 0; i < 4; ++i) {
0034     sum += adc[i];
0035     sumt += adc[i] * float(tmax - 1 + i);
0036   }
0037   return sumt / sum * 50.;  //@@ in ns. May be some bin width offset things to handle here?
0038 }
0039
0040 float CSCFindPeakTime::parabolaFitTime(int tmax, const float* adc) {
0041   // 3-point parabolic fit, from Andy Kubik
0042
0043   // We calculate offset to tmax by finding the peak of a parabola through three points
0044   float tpeak = tmax;
0045   float tcorr = 0;
0046
0047   // By construction, input array adc is for bins tmax-1 to tmax+2
0048   float y1 = adc[0];
0049   float y2 = adc[1];
0050   float y3 = adc[2];
0051
0052   // Checked and simplified... Tim Cox 08-Apr-2009
0053   // Denominator is not zero unless we fed in nonsense values with y2 not the peak!
0054   if ((y1 + y3) < 2.f * y2)
0055     tcorr = 0.5f * (y1 - y3) / (y1 - 2. * y2 + y3);
0056   tpeak += tcorr;
0057
0058   LogTrace("CSCFindPeakTime") << "[CSCFindPeakTime] tmax=" << tmax << ", parabolic peak time is tmax+" << tcorr
0059                               << " bins, or " << tpeak * 50.f << " ns";
0060
0061   return tpeak * 50.f;  // convert to ns.
0062 }
0063
0064 float CSCFindPeakTime::fivePoleFitTime(int tmax, const float* adc, float t_peak) {
0065   // Input is
0066   // tmax   = bin# 0-7 containing max SCA pulse height
0067   // adc    = 4-dim array containing SCA pulse heights in bins tmax-1 to tmax+2
0068   // t_peak = input estimate for SCA peak time
0069
0070   // Returned value is improved (we hope) estimate for SCA peak time
0071
0072   // Algorithm is to fit five-pole Semi-Gaussian function for start time of SCA pulse, t0
0073   // (The SCA peak is assumed to be 133 ns from t0.)
0074   // Note that t^4 in time domain corresponds to 1/t^5 in frequency domain (that's the 5 poles).
0075
0076   // Initialize parameters to sensible (?) values
0077
0078   float t0 = 0.f;
0079   constexpr float t0peak = 133.f;  // this is offset of peak from start time t0
0080   constexpr float p0 = 4.f / t0peak;
0081
0082   // Require that tmax is in range 2-6 of bins the eight SCA time bins 0-7
0083   // (Bins 0, 1 used for dynamic ped)
0084
0085   if (tmax < 2 || tmax > 6)
0086     return t_peak;  //@@ Just return the input value
0087
0088   // Set up time bins to match adc[4] input
0089
0090   float tb[4];
0091   for (int time = 0; time < 4; ++time) {
0092     tb[time] = (tmax + time - 1) * 50.f;
0093   }
0094
0095   // How many time bins are we fitting?
0096
0097   int n_fit = 4;
0098   if (tmax == 6)
0099     n_fit = 3;
0100
0101   float chi_min = 1.e10f;
0102   float chi_last = 1.e10f;
0103   float tt0 = 0.f;
0104   float chi2 = 0.f;
0105   float del_t = 100.f;
0106
0107   float x[4];
0108   float sx2 = 0.f;
0109   float sxy = 0.f;
0110   float fN = 0.f;
0111
0112   while (del_t > 1.f) {
0113     sx2 = 0.f;
0114     sxy = 0.f;
0115
0116     for (int j = 0; j < n_fit; ++j) {
0117       float tdif = tb[j] - tt0;
0118       x[j] = (tdif * tdif) * (tdif * tdif) * std::exp(-p0 * tdif);
0119       sx2 += x[j] * x[j];
0120       sxy += x[j] * adc[j];
0121     }
0122     fN = sxy / sx2;  // least squares fit over time bins i to adc[i] = fN * fivePoleFunction[i]
0123
0124     // Compute chi^2
0125     chi2 = 0.0;
0126     for (int j = 0; j < n_fit; ++j)
0127       chi2 += (adc[j] - fN * x[j]) * (adc[j] - fN * x[j]);
0128
0129     // Test on chi^2 to decide what to do
0130     if (chi_last > chi2) {
0131       if (chi2 < chi_min) {
0132         t0 = tt0;
0133       }
0134       chi_last = chi2;
0135       tt0 = tt0 + del_t;
0136     } else {
0137       tt0 = tt0 - 2.f * del_t;
0138       del_t = 0.5 * del_t;
0139       tt0 = tt0 + del_t;
0140       chi_last = 1.0e10f;
0141     }
0142   }
0143
0144   return t0 + t0peak;
0145 }
0146
0147 void CSCFindPeakTime::fivePoleFitCharge(
0148     int tmax, const float* adc, const float& t_zero, const float& t_peak, std::vector<float>& adcsFit) {
0149   //@@ This code can certainly be replaced by fivePoleFitTime above, but I haven't time to do that now (Tim).
0150
0151   float p0 = 4.f / t_peak;
0152   float tt0 = t_zero;
0153   int n_fit = 4;
0154   if (tmax == 6)
0155     n_fit = 3;
0156
0157   float tb[4], y[4];
0158   for (int t = 0; t < 4; ++t) {
0159     tb[t] = float(tmax + t - 1) * 50.f;
0160     y[t] = adc[t];
0161   }
0162
0163   // Find the normalization factor for the function
0164   float x[4];
0165   float sx2 = 0.f;
0166   float sxy = 0.f;
0167   for (int j = 0; j < n_fit; ++j) {
0168     float t = tb[j];
0169     x[j] = ((t - tt0) * (t - tt0)) * ((t - tt0) * (t - tt0)) * std::exp(-p0 * (t - tt0));
0170     sx2 = sx2 + x[j] * x[j];
0171     sxy = sxy + x[j] * y[j];
0172   }
0173   float N = sxy / sx2;
0174
0175   // Now compute charge for a given t  --> only need charges at: t_peak-50, t_peak and t_peak+50
0176   for (int i = 0; i < 3; ++i) {
0177     float t = t_peak + float(i - 1) * 50.f;
0178     float q_fitted = float(N) * ((t - tt0) * (t - tt0)) * ((t - tt0) * (t - tt0)) * std::exp(-p0 * (t - tt0));
0179     adcsFit.push_back(q_fitted);
0180   }
0181   return;
0182 }