EcalUncalibRecHitRatioMethodAlgo< C > Class Template Reference

#include <EcalUncalibRecHitRatioMethodAlgo.h>

Classes
struct	CalculatedRecHit
struct	Ratio
struct	Tmax

Public Member Functions
void	computeAmplitude (std::vector< double > &amplitudeFitParameters)
void	computeTime (std::vector< double > &timeFitParameters, std::pair< double, double > &timeFitLimits, std::vector< double > &amplitudeFitParameters)
bool	fixMGPAslew (const C &dataFrame)
CalculatedRecHit	getCalculatedRecHit ()
void	init (const C &dataFrame, const double *pedestals, const double *pedestalRMSes, const double *gainRatios)
virtual EcalUncalibratedRecHit	makeRecHit (const C &dataFrame, const double *pedestals, const double *pedestalRMSes, const double *gainRatios, std::vector< double > &timeFitParameters, std::vector< double > &amplitudeFitParameters, std::pair< double, double > &timeFitLimits)
virtual	~EcalUncalibRecHitRatioMethodAlgo ()

Protected Attributes
std::vector< double >	amplitudeErrors_
std::vector< double >	amplitudes_
double	ampMaxError_
CalculatedRecHit	calculatedRechit_
int	num_
double	pedestal_
std::vector< Ratio >	ratios_
std::vector< Tmax >	times_
std::vector< Tmax >	timesAB_

Detailed Description

template<class C>
class EcalUncalibRecHitRatioMethodAlgo< C >

Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse using a ratio method

Id:

EcalUncalibRecHitRatioMethodAlgo.h,v 1.8 2010/06/10 10:13:05 franzoni Exp

Date:

2010/06/10 10:13:05

Revision:

1.8

Author

A. Ledovskoy (Design) - M. Balazs (Implementation)

Definition at line 20 of file EcalUncalibRecHitRatioMethodAlgo.h.

Constructor & Destructor Documentation

template<class C>

virtual EcalUncalibRecHitRatioMethodAlgo< C >::~EcalUncalibRecHitRatioMethodAlgo ( )

inlinevirtual

Definition at line 43 of file EcalUncalibRecHitRatioMethodAlgo.h.

{ };

Member Function Documentation

template<class C >

void EcalUncalibRecHitRatioMethodAlgo< C >::computeAmplitude

(

std::vector< double > &

amplitudeFitParameters

)

Definition at line 491 of file EcalUncalibRecHitRatioMethodAlgo.h.

References beta, funct::exp(), f, i, and funct::log().

Referenced by EcalUncalibRecHitWorkerRatio::run(), and EcalUncalibRecHitWorkerGlobal::run().

{
    //                                                            //
    //             CALCULATE AMPLITUDE                            //
    //                                                            //


    double alpha = amplitudeFitParameters[0];
    double beta = amplitudeFitParameters[1];

    // calculate pedestal, again

    double pedestalLimit = calculatedRechit_.timeMax - (alpha * beta) - 1.0;
    double sumA = 0;
    double sumF = 0;
    double sumAF = 0;
    double sumFF = 0;
    double sum1 = 0;
    for (unsigned int i = 0; i < amplitudes_.size(); i++) {
            double err2 = amplitudeErrors_[i]*amplitudeErrors_[i];
        double f = 0;
        double termOne = 1 + (i - calculatedRechit_.timeMax) / (alpha * beta);
        if (termOne > 1.e-5) f = exp(alpha * log(termOne)) * exp(-(i - calculatedRechit_.timeMax) / beta);

        // apply range of interesting samples

        if ( (i < pedestalLimit)
             || (f > 0.6 && i <= calculatedRechit_.timeMax)
             || (f > 0.4 && i >= calculatedRechit_.timeMax)) {
              sum1  += 1/err2;
              sumA  += amplitudes_[i]/err2;
              sumF  += f/err2;
              sumAF += f*amplitudes_[i]/err2;
              sumFF += f*f/err2;
        }
    }

    calculatedRechit_.amplitudeMax = 0;
    if(sum1 > 0){
      double denom = sumFF*sum1 - sumF*sumF;
      if(fabs(denom)>1.0e-20){
        calculatedRechit_.amplitudeMax = (sumAF*sum1 - sumA*sumF)/denom;
      }
    }
}

template<class C >

void EcalUncalibRecHitRatioMethodAlgo< C >::computeTime	(	std::vector< double > &	timeFitParameters,
		std::pair< double, double > &	timeFitLimits,
		std::vector< double > &	amplitudeFitParameters
	)

Definition at line 187 of file EcalUncalibRecHitRatioMethodAlgo.h.

References alpha, beta, funct::exp(), f, i, EcalUncalibRecHitRatioMethodAlgo< C >::Tmax::index, getHLTprescales::index, j, gen::k, funct::log(), evf::evtn::offset(), Gflash::Rmax, Gflash::Rmin, mathSSE::sqrt(), launcher::step, tmax, and relativeConstraints::value.

Referenced by EcalUncalibRecHitWorkerRatio::run(), and EcalUncalibRecHitWorkerGlobal::run().

{
  //                                                          //
  //              RATIO METHOD FOR TIME STARTS HERE           //
  //                                                          //
  double ampMaxAlphaBeta = 0;
  double tMaxAlphaBeta = 5;
  double tMaxErrorAlphaBeta = 999;
  double tMaxRatio = 5;
  double tMaxErrorRatio = 999;

  double sumAA = 0;
  double sumA  = 0;
  double sum1  = 0;
  double sum0  = 0;
  double sumAf = 0;
  double sumff = 0;
  double NullChi2 = 0;

  // null hypothesis = no pulse, pedestal only
  for(unsigned int i = 0; i < amplitudes_.size(); i++){
    double err2 = amplitudeErrors_[i]*amplitudeErrors_[i];
    sum0  += 1;
    sum1  += 1/err2;
    sumA  += amplitudes_[i]/err2;
    sumAA += amplitudes_[i]*amplitudes_[i]/err2;
  }
  if(sum0>0){
    NullChi2 = (sumAA - sumA*sumA/sum1)/sum0;
  }else{
    // not enough samples to reconstruct the pulse
    return;
  }

  // Make all possible Ratio's based on any pair of samples i and j
  // (j>i) with positive amplitudes_
  //
  //       Ratio[k] = Amp[i]/Amp[j]
  //       where Amp[i] is pedestal subtracted ADC value in a time sample [i]
  //
  double alphabeta = amplitudeFitParameters[0]*amplitudeFitParameters[1];
  double alpha = amplitudeFitParameters[0];
  double beta = amplitudeFitParameters[1];

  for(unsigned int i = 0; i < amplitudes_.size()-1; i++){
    for(unsigned int j = i+1; j < amplitudes_.size(); j++){

      if(amplitudes_[i]>1 && amplitudes_[j]>1){

    // ratio
    double Rtmp = amplitudes_[i]/amplitudes_[j];

    // error^2 due to stat fluctuations of time samples
    // (uncorrelated for both samples)

    double err1 = Rtmp*Rtmp*( (amplitudeErrors_[i]*amplitudeErrors_[i]/(amplitudes_[i]*amplitudes_[i])) + (amplitudeErrors_[j]*amplitudeErrors_[j]/(amplitudes_[j]*amplitudes_[j])) );

    // error due to fluctuations of pedestal (common to both samples)
    double stat;
    if(num_>0) stat = num_;      // num presampeles used to compute pedestal
    else       stat = 1;         // pedestal from db
    double err2 = amplitudeErrors_[j]*(amplitudes_[i]-amplitudes_[j])/(amplitudes_[j]*amplitudes_[j])/sqrt(stat);

    //error due to integer round-down. It is relevant to low
    //amplitudes_ in gainID=1 and negligible otherwise.
        double err3 = 0.289/amplitudes_[j];

    double totalError = sqrt(err1 + err2*err2 +err3*err3);


    // don't include useless ratios
    if(totalError < 1.0
       && Rtmp>0.001
       && Rtmp<exp(double(j-i)/beta)-0.001
       ){
      Ratio currentRatio = { i, (j-i), Rtmp, totalError };
      ratios_.push_back(currentRatio);
    }

      }

    }
  }

  // No useful ratios, return zero amplitude and no time measurement
  if(!ratios_.size() >0)
    return;

  // make a vector of Tmax measurements that correspond to each ratio
  // and based on Alpha-Beta parameterization of the pulse shape

  for(unsigned int i = 0; i < ratios_.size(); i++){

    double stepOverBeta = double(ratios_[i].step)/beta;
    double offset = double(ratios_[i].index) + alphabeta;

    double Rmin = ratios_[i].value - ratios_[i].error;
    if(Rmin<0.001) Rmin=0.001;

    double Rmax = ratios_[i].value + ratios_[i].error;
    double RLimit = exp(stepOverBeta)-0.001;
    if( Rmax > RLimit ) Rmax = RLimit;

    double time1 = offset - ratios_[i].step/(exp((stepOverBeta-log(Rmin))/alpha)-1.0);
    double time2 = offset - ratios_[i].step/(exp((stepOverBeta-log(Rmax))/alpha)-1.0);

    // this is the time measurement based on the ratios[i]
    double tmax = 0.5 * (time1 + time2);
    double tmaxerr = 0.5 * sqrt( (time1 - time2)*(time1 - time2) );

    // calculate chi2
    sumAf = 0;
    sumff = 0;
    for(unsigned int it = 0; it < amplitudes_.size(); it++){
      double err2 = amplitudeErrors_[it]*amplitudeErrors_[it];
      double offset = (double(it) - tmax)/alphabeta;
      double term1 = 1.0 + offset;
      if(term1>1e-6){
    double f = exp( alpha*(log(1.0+offset) - offset) );
    sumAf += amplitudes_[it]*f/err2;
    sumff += f*f/err2;
      }
    }

    double chi2 = sumAA;
    double amp = 0;
    if( sumff > 0 ){
      chi2 = sumAA - sumAf*sumAf/sumff;
      amp = sumAf/sumff;
    }
    chi2 /= sum0;

    // choose reasonable measurements. One might argue what is
    // reasonable and what is not.
    if(chi2 > 0 && tmaxerr > 0 && tmax > 0){
      Tmax currentTmax={ ratios_[i].index, ratios_[i].step, tmax, tmaxerr, amp, chi2 };
      timesAB_.push_back(currentTmax);
    }
  }

  // no reasonable time measurements!
  if( !(timesAB_.size()> 0))
    return;

  // find minimum chi2
  double chi2min = 1.0e+9;
  double timeMinimum = 5;
  double errorMinimum = 999;
  for(unsigned int i = 0; i < timesAB_.size(); i++){
    if( timesAB_[i].chi2 <= chi2min ){
      chi2min = timesAB_[i].chi2;
      timeMinimum = timesAB_[i].value;
      errorMinimum = timesAB_[i].error;
    }
  }

  // make a weighted average of tmax measurements with "small" chi2
  // (within 1 sigma of statistical uncertainty :-)
  double chi2Limit = chi2min + 1.0;
  double time_max = 0;
  double time_wgt = 0;
  for(unsigned int i = 0; i < timesAB_.size(); i++){
    if(  timesAB_[i].chi2 < chi2Limit  ){
      double inverseSigmaSquared = 1.0/(timesAB_[i].error*timesAB_[i].error);
      time_wgt += inverseSigmaSquared;
      time_max += timesAB_[i].value*inverseSigmaSquared;
    }
  }

  tMaxAlphaBeta =  time_max/time_wgt;
  tMaxErrorAlphaBeta = 1.0/sqrt(time_wgt);

  // find amplitude and chi2
  sumAf = 0;
  sumff = 0;
  for(unsigned int i = 0; i < amplitudes_.size(); i++){
    double err2 = amplitudeErrors_[i]*amplitudeErrors_[i];
    double offset = (double(i) - tMaxAlphaBeta)/alphabeta;
    double term1 = 1.0 + offset;
    if(term1>1e-6){
      double f = exp( alpha*(log(1.0+offset) - offset) );
      sumAf += amplitudes_[i]*f/err2;
      sumff += f*f/err2;
    }
  }

  if( sumff > 0 ){
    ampMaxAlphaBeta  = sumAf/sumff;
    double chi2AlphaBeta = (sumAA - sumAf*sumAf/sumff)/sum0;
    if(chi2AlphaBeta > NullChi2){
      // null hypothesis is better
      return;
    }

  }else{
    // no visible pulse here
    return;
  }

  // if we got to this point, we have a reconstructied Tmax
  // using RatioAlphaBeta Method. To summarize:
  //
  //     tMaxAlphaBeta      - Tmax value
  //     tMaxErrorAlphaBeta - error on Tmax, but I would not trust it
  //     ampMaxAlphaBeta    - amplitude of the pulse
  //     ampMaxError_        - uncertainty of the time sample with max amplitude
  //



  // Do Ratio's Method with "large" pulses
  if( ampMaxAlphaBeta/ampMaxError_ > 5.0 ){

    // make a vector of Tmax measurements that correspond to each
    // ratio. Each measurement have it's value and the error

    double time_max = 0;
    double time_wgt = 0;


    for (unsigned int i = 0; i < ratios_.size(); i++) {

          if(ratios_[i].step == 1
              && ratios_[i].value >= timeFitLimits.first
              && ratios_[i].value <= timeFitLimits.second
            ){

        double time_max_i = ratios_[i].index;

        // calculate polynomial for Tmax

        double u = timeFitParameters[timeFitParameters.size() - 1];
        for (int k = timeFitParameters.size() - 2; k >= 0; k--) {
            u = u * ratios_[i].value + timeFitParameters[k];
        }

        // calculate derivative for Tmax error
        double du =
            (timeFitParameters.size() -
             1) * timeFitParameters[timeFitParameters.size() - 1];
        for (int k = timeFitParameters.size() - 2; k >= 1; k--) {
            du = du * ratios_[i].value + k * timeFitParameters[k];
        }


        // running sums for weighted average
        double errorsquared =
            ratios_[i].error * ratios_[i].error * du * du;
        if (errorsquared > 0) {

            time_max += (time_max_i - u) / errorsquared;
            time_wgt += 1.0 / errorsquared;
            Tmax currentTmax =
                { ratios_[i].index, 1, (time_max_i - u),
             sqrt(errorsquared),0,1 };
            times_.push_back(currentTmax);

        }
          }
        }


    // calculate weighted average of all Tmax measurements
    if (time_wgt > 0) {
          tMaxRatio = time_max/time_wgt;
          tMaxErrorRatio = 1.0/sqrt(time_wgt);

          // combine RatioAlphaBeta and Ratio Methods

          if( ampMaxAlphaBeta/ampMaxError_ > 10.0 ){

            // use pure Ratio Method
            calculatedRechit_.timeMax = tMaxRatio;
            calculatedRechit_.timeError = tMaxErrorRatio;

          }else{

            // combine two methods
            calculatedRechit_.timeMax = ( tMaxAlphaBeta*(10.0-ampMaxAlphaBeta/ampMaxError_) + tMaxRatio*(ampMaxAlphaBeta/ampMaxError_ - 5.0) )/5.0;
            calculatedRechit_.timeError = ( tMaxErrorAlphaBeta*(10.0-ampMaxAlphaBeta/ampMaxError_) + tMaxErrorRatio*(ampMaxAlphaBeta/ampMaxError_ - 5.0) )/5.0;

          }

        }else{

          // use RatioAlphaBeta Method
          calculatedRechit_.timeMax = tMaxAlphaBeta;
          calculatedRechit_.timeError = tMaxErrorAlphaBeta;

        }

  }else{

    // use RatioAlphaBeta Method
    calculatedRechit_.timeMax = tMaxAlphaBeta;
    calculatedRechit_.timeError = tMaxErrorAlphaBeta;

  }
}

template<class C>

bool EcalUncalibRecHitRatioMethodAlgo< C >::fixMGPAslew

(

const C &

dataFrame

)

Definition at line 156 of file EcalUncalibRecHitRatioMethodAlgo.h.

References query::result.

Referenced by EcalUncalibRecHitWorkerRatio::run(), and EcalUncalibRecHitWorkerGlobal::run().

{

  // This fuction finds sample(s) preceeding gain switching and
  // inflates errors on this sample, therefore, making this sample
  // invisible for Ratio Method. Only gain switching DOWN is
  // considered Only gainID=1,2,3 are considered. In case of the
  // saturation (gainID=0), we keep "distorted" sample because it is
  // the only chance to make time measurement; the qualilty of it will
  // be bad anyway.

  bool result = false;

  int GainIdPrev;
  int GainIdNext;
  for (int iSample = 1; iSample < C::MAXSAMPLES; iSample++) {
    GainIdPrev = dataFrame.sample(iSample-1).gainId();
    GainIdNext = dataFrame.sample(iSample).gainId();
    if( GainIdPrev>=1 && GainIdPrev<=3 && 
        GainIdNext>=1 && GainIdNext<=3 && 
        GainIdPrev<GainIdNext ){
      amplitudes_[iSample-1]=1e-9;
      amplitudeErrors_[iSample-1]=1e+9;
      result = true;      
    }
  }
  return result;

}

template<class C>

CalculatedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::getCalculatedRecHit ( )

inline

Definition at line 57 of file EcalUncalibRecHitRatioMethodAlgo.h.

Referenced by EcalUncalibRecHitWorkerRatio::run(), and EcalUncalibRecHitWorkerGlobal::run().

{ return calculatedRechit_; };

template<class C>

void EcalUncalibRecHitRatioMethodAlgo< C >::init	(	const C &	dataFrame,
		const double *	pedestals,
		const double *	pedestalRMSes,
		const double *	gainRatios
	)

Definition at line 75 of file EcalUncalibRecHitRatioMethodAlgo.h.

Referenced by EcalUncalibRecHitWorkerRatio::run(), and EcalUncalibRecHitWorkerGlobal::run().

{
    calculatedRechit_.timeMax = 5;
    calculatedRechit_.amplitudeMax = 0;
    calculatedRechit_.timeError = -999;
    amplitudes_.clear();
    amplitudes_.reserve(C::MAXSAMPLES);
    amplitudeErrors_.clear();
    amplitudeErrors_.reserve(C::MAXSAMPLES);
    ratios_.clear();
    ratios_.reserve(C::MAXSAMPLES*(C::MAXSAMPLES-1)/2);
    times_.clear();
    times_.reserve(C::MAXSAMPLES*(C::MAXSAMPLES-1)/2);
    timesAB_.clear();
    timesAB_.reserve(C::MAXSAMPLES*(C::MAXSAMPLES-1)/2);
    
    // to obtain gain 12 pedestal:
    // -> if it's in gain 12, use first sample
    // --> average it with second sample if in gain 12 and 3-sigma-noise compatible (better LF noise cancellation)
    // -> else use pedestal from database
    pedestal_ = 0;
    num_      = 0;
    if (dataFrame.sample(0).gainId() == 1) {
        pedestal_ += double (dataFrame.sample(0).adc());
        num_++;
    }
    if (num_!=0 &&
        dataFrame.sample(1).gainId() == 1 && 
        fabs(dataFrame.sample(1).adc()-dataFrame.sample(0).adc())<3*pedestalRMSes[0]) {
            pedestal_ += double (dataFrame.sample(1).adc());
            num_++;
    }
    if (num_ != 0)
        pedestal_ /= num_;
    else
        pedestal_ = pedestals[0];

        // fill vector of amplitudes, pedestal subtracted and vector
        // of amplitude uncertainties Also, find the uncertainty of a
        // sample with max amplitude. We will use it later.

        ampMaxError_ = 0;
        double ampMaxValue = -1000;

    // ped-subtracted and gain-renormalized samples. It is VERY
    // IMPORTANT to have samples one clock apart which means to
    // have vector size equal to MAXSAMPLES
    double sample;
        double sampleError;
    int GainId;
        for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
          GainId = dataFrame.sample(iSample).gainId();

          if (GainId == 1) {
            sample      = double (dataFrame.sample(iSample).adc() - pedestal_);
            sampleError = pedestalRMSes[0];
          } else if (GainId == 2 || GainId == 3){
            sample      = (double (dataFrame.sample(iSample).adc() - pedestals[GainId - 1])) *gainRatios[GainId - 1];
            sampleError = pedestalRMSes[GainId-1]*gainRatios[GainId-1];
          } else {
        sample      = 1e-9;  // GainId=0 case falls here, from saturation
        sampleError = 1e+9;  // inflate error so won't generate ratio considered for the measurement 
      }


          if(sampleError>0){
            amplitudes_.push_back(sample);
            amplitudeErrors_.push_back(sampleError);
            if(ampMaxValue < sample){
              ampMaxValue = sample;
              ampMaxError_ = sampleError;
            }
          }else{
        // inflate error for useless samples
            amplitudes_.push_back(sample);
            amplitudeErrors_.push_back(1e+9);
      }
        }
}

template<class C>

EcalUncalibratedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::makeRecHit ( const C &  dataFrame, const double *  pedestals, const double *  pedestalRMSes, const double *  gainRatios, std::vector< double > &  timeFitParameters, std::vector< double > &  amplitudeFitParameters, std::pair< double, double > &  timeFitLimits  )

virtual

Definition at line 542 of file EcalUncalibRecHitRatioMethodAlgo.h.

References init.

{

        init( dataFrame, pedestals, pedestalRMSes, gainRatios );
        computeTime( timeFitParameters, timeFitLimits, amplitudeFitParameters );
        computeAmplitude( amplitudeFitParameters );

    // 1st parameters is id
    //
    // 2nd parameters is amplitude. It is calculated by this method.
    //
    // 3rd parameter is pedestal. It is not calculated. This method
    // relies on input parameters for pedestals and gain ratio. Return
    // zero.
    //
    // 4th parameter is jitter which is a bad choice to call Tmax. It is
    // calculated by this method (in 25 nsec clock units)
    //
    // GF subtract 5 so that jitter==0 corresponds to synchronous hit
    //
    //
    // 5th parameter is chi2. It is possible to calculate chi2 for
    // Tmax. It is possible to calculate chi2 for Amax. However, these
    // values are not very useful without TmaxErr and AmaxErr. This
    // method can return one value for chi2 but there are 4 different
    // parameters that have useful information about the quality of Amax
    // ans Tmax. For now we can return TmaxErr. The quality of Tmax and
    // Amax can be judged from the magnitude of TmaxErr

    return EcalUncalibratedRecHit(dataFrame.id(),
                      calculatedRechit_.amplitudeMax,
                                      pedestal_,
                      calculatedRechit_.timeMax - 5,
                      calculatedRechit_.timeError);
}

Member Data Documentation

template<class C>

std::vector< double > EcalUncalibRecHitRatioMethodAlgo< C >::amplitudeErrors_

protected

Definition at line 62 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>

std::vector< double > EcalUncalibRecHitRatioMethodAlgo< C >::amplitudes_

protected

Definition at line 61 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>

double EcalUncalibRecHitRatioMethodAlgo< C >::ampMaxError_

protected

Definition at line 69 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>

CalculatedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::calculatedRechit_

protected

Definition at line 71 of file EcalUncalibRecHitRatioMethodAlgo.h.

Referenced by EcalUncalibRecHitRatioMethodAlgo< EBDataFrame >::getCalculatedRecHit().

template<class C>

int EcalUncalibRecHitRatioMethodAlgo< C >::num_

protected

Definition at line 68 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>

double EcalUncalibRecHitRatioMethodAlgo< C >::pedestal_

protected

Definition at line 67 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>

std::vector< Ratio > EcalUncalibRecHitRatioMethodAlgo< C >::ratios_

protected

Definition at line 63 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>

std::vector< Tmax > EcalUncalibRecHitRatioMethodAlgo< C >::times_

protected

Definition at line 64 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>

std::vector< Tmax > EcalUncalibRecHitRatioMethodAlgo< C >::timesAB_

protected

Definition at line 65 of file EcalUncalibRecHitRatioMethodAlgo.h.