#include <CSCFindPeakTime.h>

Public Member Functions
float	averageTime (int tmax, const float *adc)
	Weighted average of time bins. More...

	CSCFindPeakTime (const edm::ParameterSet &ps)

void	fivePoleFitCharge (int tmax, const float *adc, const float &t_zero, const float &t_peak, std::vector< float > &adcsFit)

float	fivePoleFitTime (int tmax, const float *adc, float t_peak)

float	parabolaFitTime (int tmax, const float *adc)
	Parabolic fit to three time bins centered on maximum. More...

float	peakTime (int tmax, const float *adc, float t_peak)
	Basic result of this class. More...

	~CSCFindPeakTime ()

Private Attributes
bool	useAverageTime

bool	useFivePoleFit

bool	useParabolaFit

Detailed Description

This is CSCFindPeakTime

Used to provide improved estimate of SCA peak time.

Definition at line 13 of file CSCFindPeakTime.h.

Constructor & Destructor Documentation

CSCFindPeakTime::CSCFindPeakTime ( const edm::ParameterSet & ps )

explicit

Definition at line 8 of file CSCFindPeakTime.cc.

References edm::ParameterSet::getParameter(), LogTrace, useAverageTime, useFivePoleFit, and useParabolaFit.

                                                            : 
   useAverageTime(false), useParabolaFit(false), useFivePoleFit(false) {
 
   useAverageTime = ps.getParameter<bool>("UseAverageTime");
   useParabolaFit = ps.getParameter<bool>("UseParabolaFit");
   useFivePoleFit = ps.getParameter<bool>("UseFivePoleFit");
   LogTrace("CSCRecHit") << "[CSCFindPeakTime] useAverageTime=" << useAverageTime <<
     ", useParabolaFit=" << useParabolaFit << ", useFivePoleFit=" << useFivePoleFit;
 }

CSCFindPeakTime::~CSCFindPeakTime ( )

inline

Definition at line 19 of file CSCFindPeakTime.h.

References ecalMGPA::adc(), averageTime(), fivePoleFitCharge(), fivePoleFitTime(), parabolaFitTime(), peakTime(), and tmax.

19 {};

Member Function Documentation

float CSCFindPeakTime::averageTime	(	int	tmax,
		const float *	adc
	)

Weighted average of time bins.

Definition at line 34 of file CSCFindPeakTime.cc.

References objects.autophobj::float, and mps_fire::i.

Referenced by peakTime(), and ~CSCFindPeakTime().

                                                                {
    float sum  = 0.;
    float sumt = 0.;
    for (size_t i=0; i<4; ++i){
      sum  += adc[i];
      sumt += adc[i] * float( tmax - 1 + i );
 
    }
    return sumt/sum * 50.; //@@ in ns. May be some bin width offset things to handle here?
 }

void CSCFindPeakTime::fivePoleFitCharge	(	int	tmax,
		const float *	adc,
		const float &	t_zero,
		const float &	t_peak,
		std::vector< float > &	adcsFit
	)

Integrated charge after fivePoleFitTime

Definition at line 151 of file CSCFindPeakTime.cc.

References JetChargeProducer_cfi::exp, objects.autophobj::float, mps_fire::i, N, lumiQTWidget::t, x, and y.

Referenced by ~CSCFindPeakTime().

                                                                                                                                          {
 
   //@@ This code can certainly be replaced by fivePoleFitTime above, but I haven't time to do that now (Tim).
 
   float p0  = 4.f/t_peak;
   float tt0 = t_zero;
   int n_fit = 4;
   if ( tmax == 6 ) n_fit=3;
   
   float tb[4], y[4];
   for ( int t = 0; t < 4; ++t ){
     tb[t] = float(tmax + t - 1) * 50.f;
     y[t] = adc[t];
   }
 
   // Find the normalization factor for the function
   float x[4];    
   float sx2 = 0.f;
   float sxy = 0.f;
   for ( int j=0; j < n_fit; ++j ) {
     float t = tb[j];
     x[j] = ((t-tt0)*(t-tt0))*((t-tt0)*(t-tt0)) * std::exp( -p0 * (t-tt0) );
     sx2  = sx2 + x[j] * x[j];
     sxy  = sxy + x[j] * y[j];
   }
   float N = sxy / sx2;
     
 
   // Now compute charge for a given t  --> only need charges at: t_peak-50, t_peak and t_peak+50
   for ( int i = 0; i < 3; ++i ) {
     float t = t_peak + float(i - 1) * 50.f;
     float q_fitted = float(N) * ((t-tt0)*(t-tt0))*((t-tt0)*(t-tt0)) * std::exp( -p0 * (t-tt0) );
     adcsFit.push_back(q_fitted);
   }
   return;
 }

float CSCFindPeakTime::fivePoleFitTime	(	int	tmax,
		const float *	adc,
		float	t_peak
	)

Based on RecoLocalMuon/CSCStandAlone/interface/PulseTime.h by S. Durkin, and ported by D. Fortin. Comments updated by Tim Cox Apr 2009.

The SCA pulse shape should be representable by a function
N*(p0^2/256/exp(-4)) * (t-t0)^4 * exp( -p0*(t-t0) )

Rather than do a full fit with varying peak time too, assume the peak time is fixed to 133 nsec w.r.t. start time, t0, and fit for t0. The fit uses a binary search in t0, and at each step calculates the overall normalization factor between the function and the SCA pulse height as a least-squares fit over the 4 time bins tmax -1, tmax, tmax+1, tmax+2

Note: t0peak =4/p0 = 133 nsec, and adc[0] is arbitrarily defined a time of 0.0 nsec.

Definition at line 68 of file CSCFindPeakTime.cc.

References HiEvtPlane_cfi::chi2, constexpr, JetChargeProducer_cfi::exp, f, hcaldqm::quantity::fN, cscNeutronWriter_cfi::t0, ntuplemaker::time, and x.

Referenced by peakTime(), and ~CSCFindPeakTime().

                                                                                  {
 
   // Input is 
   // tmax   = bin# 0-7 containing max SCA pulse height  
   // adc    = 4-dim array containing SCA pulse heights in bins tmax-1 to tmax+2
   // t_peak = input estimate for SCA peak time
 
   // Returned value is improved (we hope) estimate for SCA peak time
 
   // Algorithm is to fit five-pole Semi-Gaussian function for start time of SCA pulse, t0
   // (The SCA peak is assumed to be 133 ns from t0.)
   // Note that t^4 in time domain corresponds to 1/t^5 in frequency domain (that's the 5 poles).
 
   // Initialize parameters to sensible (?) values
 
   float t0       = 0.f;
   constexpr float t0peak   = 133.f;   // this is offset of peak from start time t0
   constexpr float p0       = 4.f/t0peak;
 
   // Require that tmax is in range 2-6 of bins the eight SCA time bins 0-7
   // (Bins 0, 1 used for dynamic ped)
 
   if ( tmax < 2 || tmax > 6 ) return t_peak; //@@ Just return the input value
 
   // Set up time bins to match adc[4] input
 
   float tb[4];
   for ( int time=0; time<4; ++time ){
     tb[time] = (tmax + time -1) * 50.f;
   }
 
   // How many time bins are we fitting?
 
   int n_fit  = 4;
   if ( tmax == 6 ) n_fit = 3;
 
   float chi_min  = 1.e10f;
   float chi_last = 1.e10f;
   float tt0      = 0.f;
   float chi2     = 0.f;
   float del_t    = 100.f;
 
   float x[4];
   float sx2 = 0.f;
   float sxy = 0.f;
   float fN  = 0.f;
 
   while ( del_t > 1.f ) {
     sx2 = 0.f;
     sxy = 0.f;
         
     for ( int j=0; j < n_fit; ++j ) {
       float tdif = tb[j] - tt0;
       x[j] = (tdif * tdif) * (tdif * tdif) * std::exp( -p0 * tdif );
       sx2 += x[j] * x[j];
       sxy += x[j] * adc[j];
     }
     fN = sxy / sx2; // least squares fit over time bins i to adc[i] = fN * fivePoleFunction[i]
     
     // Compute chi^2
     chi2 = 0.0;
     for (int j=0; j < n_fit; ++j) chi2 += (adc[j] - fN * x[j]) * (adc[j] - fN * x[j]);
 
     // Test on chi^2 to decide what to do    
     if ( chi_last > chi2 ) {
       if (chi2 < chi_min ){
         t0      = tt0;
       }
       chi_last  = chi2;
       tt0       = tt0 + del_t;
     } else {
       tt0      = tt0 - 2.f * del_t;
       del_t    = 0.5*del_t;
       tt0      = tt0 + del_t;
       chi_last = 1.0e10f;
     }
   }
 
   return t0 + t0peak;
 }

float CSCFindPeakTime::parabolaFitTime	(	int	tmax,
		const float *	adc
	)

Parabolic fit to three time bins centered on maximum.

Definition at line 45 of file CSCFindPeakTime.cc.

References f, LogTrace, and tmax.

Referenced by peakTime(), and ~CSCFindPeakTime().

                                                                    {
   // 3-point parabolic fit, from Andy Kubik
  
   // We calculate offset to tmax by finding the peak of a parabola through three points
    float tpeak = tmax;
    float tcorr = 0;
 
    // By construction, input array adc is for bins tmax-1 to tmax+2
    float y1 = adc[0];
    float y2 = adc[1];
    float y3 = adc[2];
 
    // Checked and simplified... Tim Cox 08-Apr-2009
    // Denominator is not zero unless we fed in nonsense values with y2 not the peak!
    if ( (y1+y3) < 2.f*y2 ) tcorr =  0.5f * ( y1 - y3 ) / ( y1 - 2.*y2 + y3 );
    tpeak += tcorr;
 
    LogTrace("CSCFindPeakTime") << "[CSCFindPeakTime] tmax=" << tmax 
      << ", parabolic peak time is tmax+" << tcorr <<" bins, or " << tpeak*50.f << " ns";
    
    return tpeak * 50.f; // convert to ns.
 }

float CSCFindPeakTime::peakTime	(	int	tmax,
		const float *	adc,
		float	t_peak
	)

Basic result of this class.

Definition at line 18 of file CSCFindPeakTime.cc.

References averageTime(), fivePoleFitTime(), parabolaFitTime(), useAverageTime, useFivePoleFit, and useParabolaFit.

Referenced by ~CSCFindPeakTime().

                                                                         {
   if ( useAverageTime ) {
     return averageTime( tmax, adc );
   }
   else if ( useParabolaFit ) {
     return parabolaFitTime( tmax, adc );
   }
   else if ( useFivePoleFit ) {
      return fivePoleFitTime( tmax, adc, t_peak);
   }
   else {
   // return something, anyway.. may as well be average
     return averageTime( tmax, adc );
   }
 }