RFFHitFinderAlg.cxx

Go to the documentation of this file.

#include "fhiclcpp/ParameterSet.h"
#include "larcore/Geometry/Geometry.h"
#include "larreco/HitFinder/RFFHitFitter.h"

#include "RFFHitFinderAlg.h"

#include <numeric>

hit::RFFHitFinderAlg::RFFHitFinderAlg(fhicl::ParameterSet const& p)
{
  fMatchThresholdVec = p.get<std::vector<float>>("MeanMatchThreshold");
  fMergeMultiplicityVec = p.get<std::vector<unsigned int>>("MinMergeMultiplicity");
  fAmpThresholdVec = p.get<std::vector<float>>("AmplitudeThreshold", std::vector<float>(1, 0.0));
}

void hit::RFFHitFinderAlg::SetFitterParamsVectors(geo::Geometry const& geo)
{
  const unsigned int n_planes = geo.Nplanes();

  //If size zero, throw. If size one, assume same for all planes.
  //If size > 1 but < n_planes, throw. If size = n_plane, good.

  if (fMatchThresholdVec.size() == 0 || fMergeMultiplicityVec.size() == 0 ||
      fAmpThresholdVec.size() == 0)
    throw std::runtime_error("Error in RFFHitFinderAlg: Configured with zero planes.");

  if ((fMatchThresholdVec.size() > 1 && fMatchThresholdVec.size() < n_planes) ||
      (fMergeMultiplicityVec.size() > 1 && fMergeMultiplicityVec.size() < n_planes) ||
      (fAmpThresholdVec.size() > 1 && fAmpThresholdVec.size() < n_planes))
    throw std::runtime_error("Error in RFFHitFinderAlg: Configured with incorrect n_planes.");

  if (fMatchThresholdVec.size() == 1) fMatchThresholdVec.resize(n_planes, fMatchThresholdVec[0]);

  if (fMergeMultiplicityVec.size() == 1)
    fMergeMultiplicityVec.resize(n_planes, fMergeMultiplicityVec[0]);

  if (fAmpThresholdVec.size() == 1) fAmpThresholdVec.resize(n_planes, fAmpThresholdVec[0]);
}

void hit::RFFHitFinderAlg::SetFitterParams(unsigned int p)
{
  fFitter.SetFitterParams(fMatchThresholdVec[p], fMergeMultiplicityVec[p], fAmpThresholdVec[p]);
}

void hit::RFFHitFinderAlg::Run(std::vector<recob::Wire> const& wireVector,
                               std::vector<recob::Hit>& hitVector,
                               geo::Geometry const& geo)
{
  hitVector.reserve(wireVector.size());
  for (auto const& wire : wireVector) {
    geo::SigType_t const& sigtype = geo.SignalType(wire.Channel());
    geo::WireID const wireID = geo.ChannelToWire(wire.Channel()).at(0);

    SetFitterParams(wire.View());

    for (auto const& roi : wire.SignalROI().get_ranges()) {
      fFitter.RunFitter(roi.data());

      const float summedADCTotal = std::accumulate(roi.data().begin(), roi.data().end(), 0.0);
      const raw::TDCtick_t startTick = roi.begin_index();
      const raw::TDCtick_t endTick = roi.begin_index() + roi.size();

      EmplaceHit(hitVector, wire, summedADCTotal, startTick, endTick, sigtype, wireID);
    } //end loop over ROIs on wire

  } //end loop over wires
}

void hit::RFFHitFinderAlg::EmplaceHit(std::vector<recob::Hit>& hitVector,
                                      recob::Wire const& wire,
                                      float const& summedADCTotal,
                                      raw::TDCtick_t const& startTick,
                                      raw::TDCtick_t const& endTick,
                                      geo::SigType_t const& sigtype,
                                      geo::WireID const& wireID)
{

  float totalArea = 0.0;
  std::vector<float> areaVector(fFitter.NHits());
  std::vector<float> areaErrorVector(fFitter.NHits());
  std::vector<float> areaFracVector(fFitter.NHits());

  for (size_t ihit = 0; ihit < fFitter.NHits(); ihit++) {
    areaVector[ihit] = fFitter.AmplitudeVector()[ihit] * fFitter.SigmaVector()[ihit] * SQRT_TWO_PI;
    areaErrorVector[ihit] =
      SQRT_TWO_PI * std::sqrt(fFitter.AmplitudeVector()[ihit] * fFitter.SigmaErrorVector()[ihit] *
                                fFitter.AmplitudeVector()[ihit] * fFitter.SigmaErrorVector()[ihit] +
                              fFitter.AmplitudeErrorVector()[ihit] * fFitter.SigmaVector()[ihit] *
                                fFitter.AmplitudeErrorVector()[ihit] * fFitter.SigmaVector()[ihit]);
    totalArea += areaVector[ihit];
  }

  for (size_t ihit = 0; ihit < fFitter.NHits(); ihit++) {
    areaFracVector[ihit] = areaVector[ihit] / totalArea;

    hitVector.emplace_back(wire.Channel(),
                           startTick,
                           endTick,
                           fFitter.MeanVector()[ihit] + (float)startTick,
                           fFitter.MeanErrorVector()[ihit],
                           fFitter.SigmaVector()[ihit],
                           fFitter.AmplitudeVector()[ihit],
                           fFitter.AmplitudeErrorVector()[ihit],
                           summedADCTotal * areaFracVector[ihit],
                           areaVector[ihit],
                           areaErrorVector[ihit],
                           fFitter.NHits(),
                           ihit,
                           -999.,
                           -999,
                           wire.View(),
                           sigtype,
                           wireID);
  }
}