TotemTimingRecHitProducerAlgorithm.cc

Go to the documentation of this file.

/****************************************************************************
 *
 * This is a part of CTPPS offline software.
 * Authors:
 *   Laurent Forthomme (laurent.forthomme@cern.ch)
 *   Nicola Minafra
 *
 ****************************************************************************/

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "RecoCTPPS/TotemRPLocal/interface/TotemTimingRecHitProducerAlgorithm.h"

#include <numeric>

//----------------------------------------------------------------------------------------------------

TotemTimingRecHitProducerAlgorithm::TotemTimingRecHitProducerAlgorithm(const edm::ParameterSet& iConfig)
    : mergeTimePeaks_(iConfig.getParameter<bool>("mergeTimePeaks")),
      baselinePoints_(iConfig.getParameter<int>("baselinePoints")),
      saturationLimit_(iConfig.getParameter<double>("saturationLimit")),
      cfdFraction_(iConfig.getParameter<double>("cfdFraction")),
      smoothingPoints_(iConfig.getParameter<int>("smoothingPoints")),
      lowPassFrequency_(iConfig.getParameter<double>("lowPassFrequency")),
      hysteresis_(iConfig.getParameter<double>("hysteresis")) {}

//----------------------------------------------------------------------------------------------------

void TotemTimingRecHitProducerAlgorithm::setCalibration(const PPSTimingCalibration& calib) {
  sampicConversions_.reset(new TotemTimingConversions(mergeTimePeaks_, calib));
}

//----------------------------------------------------------------------------------------------------

void TotemTimingRecHitProducerAlgorithm::build(const CTPPSGeometry& geom,
                                               const edm::DetSetVector<TotemTimingDigi>& input,
                                               edm::DetSetVector<TotemTimingRecHit>& output) {
  for (const auto& vec : input) {
    const TotemTimingDetId detid(vec.detId());

    float x_pos = 0.f, y_pos = 0.f, z_pos = 0.f;
    float x_width = 0.f, y_width = 0.f, z_width = 0.f;

    // retrieve the geometry element associated to this DetID ( if present )
    const DetGeomDesc* det = geom.sensorNoThrow(detid);

    if (det) {
      x_pos = det->translation().x();
      y_pos = det->translation().y();
      z_pos = det->parentZPosition();  // retrieve the plane position;

      x_width = 2.0 * det->params()[0],  // parameters stand for half the size
          y_width = 2.0 * det->params()[1], z_width = 2.0 * det->params()[2];
    } else
      throw cms::Exception("TotemTimingRecHitProducerAlgorithm") << "Failed to retrieve a sensor for " << detid;

    if (!sampicConversions_)
      throw cms::Exception("TotemTimingRecHitProducerAlgorithm") << "No timing conversion retrieved.";

    edm::DetSet<TotemTimingRecHit>& rec_hits = output.find_or_insert(detid);

    for (const auto& digi : vec) {
      const float triggerCellTimeInstant(sampicConversions_->triggerTime(digi));
      const float timePrecision(sampicConversions_->timePrecision(digi));
      const std::vector<float> time(sampicConversions_->timeSamples(digi));
      std::vector<float> data(sampicConversions_->voltSamples(digi));

      auto max_it = std::max_element(data.begin(), data.end());

      RegressionResults baselineRegression = simplifiedLinearRegression(time, data, 0, baselinePoints_);

      // remove baseline
      std::vector<float> dataCorrected(data.size());
      for (unsigned int i = 0; i < data.size(); ++i)
        dataCorrected[i] = data[i] - (baselineRegression.q + baselineRegression.m * time[i]);
      auto max_corrected_it = std::max_element(dataCorrected.begin(), dataCorrected.end());

      float t = TotemTimingRecHit::NO_T_AVAILABLE;
      if (*max_it < saturationLimit_)
        t = constantFractionDiscriminator(time, dataCorrected);

      mode_ = TotemTimingRecHit::CFD;

      rec_hits.push_back(TotemTimingRecHit(x_pos,
                                           x_width,
                                           y_pos,
                                           y_width,
                                           z_pos,
                                           z_width,  // spatial information
                                           t,
                                           triggerCellTimeInstant,
                                           timePrecision,
                                           *max_corrected_it,
                                           baselineRegression.rms,
                                           mode_));
    }
  }
}

//----------------------------------------------------------------------------------------------------

TotemTimingRecHitProducerAlgorithm::RegressionResults TotemTimingRecHitProducerAlgorithm::simplifiedLinearRegression(
    const std::vector<float>& time,
    const std::vector<float>& data,
    const unsigned int start_at,
    const unsigned int points) const {
  RegressionResults results;
  if (time.size() != data.size())
    return results;
  if (start_at > time.size())
    return results;
  unsigned int stop_at = std::min((unsigned int)time.size(), start_at + points);
  unsigned int realPoints = stop_at - start_at;
  auto t_begin = std::next(time.begin(), start_at);
  auto t_end = std::next(time.begin(), stop_at);
  auto d_begin = std::next(data.begin(), start_at);
  auto d_end = std::next(data.begin(), stop_at);

  const float sx = std::accumulate(t_begin, t_end, 0.);
  const float sxx = std::inner_product(t_begin, t_end, t_begin, 0.);  // sum(t**2)
  const float sy = std::accumulate(d_begin, d_end, 0.);

  float sxy = 0.;
  for (unsigned int i = 0; i < realPoints; ++i)
    sxy += time[i] * data[i];

  // y = mx + q
  results.m = (sxy * realPoints - sx * sy) / (sxx * realPoints - sx * sx);
  results.q = sy / realPoints - results.m * sx / realPoints;

  float correctedSyy = 0.;
  for (unsigned int i = 0; i < realPoints; ++i)
    correctedSyy += pow(data[i] - (results.m * time[i] + results.q), 2);
  results.rms = sqrt(correctedSyy / realPoints);

  return results;
}

//----------------------------------------------------------------------------------------------------

int TotemTimingRecHitProducerAlgorithm::fastDiscriminator(const std::vector<float>& data, float threshold) const {
  int threholdCrossingIndex = -1;
  bool above = false;
  bool lockForHysteresis = false;

  for (unsigned int i = 0; i < data.size(); ++i) {
    // Look for first edge
    if (!above && !lockForHysteresis && data[i] > threshold) {
      threholdCrossingIndex = i;
      above = true;
      lockForHysteresis = true;
    }
    if (above && lockForHysteresis)  // NOTE: not else if!, "above" can be set in
                                     // the previous if
    {
      // Lock until above threshold_+hysteresis
      if (lockForHysteresis && data[i] > threshold + hysteresis_)
        lockForHysteresis = false;
      // Ignore noise peaks
      if (lockForHysteresis && data[i] < threshold) {
        above = false;
        lockForHysteresis = false;
        threholdCrossingIndex = -1;  // assigned because of noise
      }
    }
  }

  return threholdCrossingIndex;
}

float TotemTimingRecHitProducerAlgorithm::constantFractionDiscriminator(const std::vector<float>& time,
                                                                        const std::vector<float>& data) {
  std::vector<float> dataProcessed(data);
  if (lowPassFrequency_ != 0)  // Smoothing
    for (unsigned short i = 0; i < data.size(); ++i)
      for (unsigned short j = -smoothingPoints_ / 2; j <= +smoothingPoints_ / 2; ++j)
        if ((i + j) >= 0 && (i + j) < data.size() && j != 0) {
          float x = SINC_COEFFICIENT * lowPassFrequency_ * j;
          dataProcessed[i] += data[i + j] * std::sin(x) / x;
        }
  auto max_it = std::max_element(dataProcessed.begin(), dataProcessed.end());
  float max = *max_it;

  float threshold = cfdFraction_ * max;
  int indexOfThresholdCrossing = fastDiscriminator(dataProcessed, threshold);

  //--- necessary sanity checks before proceeding with the extrapolation
  return (indexOfThresholdCrossing >= 1 && indexOfThresholdCrossing >= baselinePoints_ &&
          indexOfThresholdCrossing < (int)time.size())
             ? (time[indexOfThresholdCrossing - 1] - time[indexOfThresholdCrossing]) /
                       (dataProcessed[indexOfThresholdCrossing - 1] - dataProcessed[indexOfThresholdCrossing]) *
                       (threshold - dataProcessed[indexOfThresholdCrossing]) +
                   time[indexOfThresholdCrossing]
             : (float)TotemTimingRecHit::NO_T_AVAILABLE;
}