hit::CCHitFinderAlg Class Reference

Hit finder algorithm designed to work with Cluster Crawler. More...

#include "CCHitFinderAlg.h"

Classes
struct	FitStats_t
struct	HitChannelInfo_t
	exchange data about the originating wire More...

Public Member Functions
	CCHitFinderAlg (fhicl::ParameterSet const &pset)
void	RunCCHitFinder (std::vector< recob::Wire > const &Wires)
std::vector< recob::Hit > &&	YieldHits ()
	Returns (and loses) the collection of reconstructed hits. More...
template<typename Stream >
void	PrintStats (Stream &out) const
	Print the fit statistics. More...

Public Attributes
std::vector< recob::Hit >	allhits

Private Member Functions
void	FitNG (unsigned short nGaus, unsigned short npt, float *ticks, float *signl)
void	MakeCrudeHit (unsigned short npt, float *ticks, float *signl)
void	StoreHits (unsigned short TStart, unsigned short npt, HitChannelInfo_t info, float adcsum)
void	StudyHits (unsigned short flag, unsigned short npt=0, float *signl=nullptr, unsigned short tstart=0)

Static Private Member Functions
static bool	FastGaussianFit (unsigned short npt, float const *ticks, float const *signl, std::array< double, 3 > &params, std::array< double, 3 > &paramerrors, float &chidof)
	Performs a "fast" fit. More...

Private Attributes
std::vector< float >	fMinPeak
std::vector< float >	fMinRMS
unsigned short	fMaxBumps
unsigned short	fMaxXtraHits
float	fChiSplit
	Estimated noise error on the Signal. More...
std::vector< float >	fChiNorms
std::vector< float >	fTimeOffsets
std::vector< float >	fChgNorms
raw::ChannelID_t	theChannel
unsigned short	theWireNum
unsigned short	thePlane
float	chinorm
bool	fUseChannelFilter
geo::WireReadoutGeom const *	wireReadoutGeom
std::vector< double >	par
std::vector< double >	parerr
std::vector< double >	parmin
std::vector< double >	parmax
float	chidof
int	dof
std::vector< unsigned short >	bumps
bool	fStudyHits
std::vector< short >	fUWireRange
std::vector< short >	fUTickRange
std::vector< short >	fVWireRange
std::vector< short >	fVTickRange
std::vector< short >	fWWireRange
std::vector< short >	fWTickRange
std::vector< int >	bumpCnt
std::vector< int >	RATCnt
std::vector< float >	bumpChi
std::vector< float >	bumpRMS
std::vector< int >	hitCnt
std::vector< float >	hitRMS
std::vector< float >	loWire
std::vector< float >	loTime
std::vector< float >	hiWire
std::vector< float >	hiTime
bool	SelRAT
bool	fUseFastFit
	whether to attempt using a fast fit on single gauss. More...
GausFitCache	FitCache
	a set of functions ready to be used More...
FitStats_t	FinalFitStats
	counts of the good fits More...
FitStats_t	TriedFitStats
	counts of the tried fits More...

Static Private Attributes
static constexpr float	Sqrt2Pi = 2.5066
static constexpr float	SqrtPi = 1.7725
static constexpr unsigned int	MaxGaussians = 20

Detailed Description

Hit finder algorithm designed to work with Cluster Crawler.

This algorithm used to store hits in a proprietary CCHit data structure. It has now been changed to use recob::Hit class directly. It is possible to translate the former into the latter, with one exception, as follows:

// this is the original CCHit definition
struct CCHit {
  float Charge;            // recob::Hit::Integral()
  float ChargeErr;         // recob::Hit::SigmaIntegral()
  float Amplitude;         // recob::Hit::PeakAmplitude()
  float AmplitudeErr;      // recob::Hit::SigmaPeakAmplitude()
  float Time;              // recob::Hit::PeakTime()
  float TimeErr;           // recob::Hit::SigmaPeakTime()
  float RMS;               // recob::Hit::RMS()
  float RMSErr;            // dropped
  float ChiDOF;            // recob::Hit::GoodnessOfFit()
  int   DOF;               // recob::Hit::DegreesOfFreedom()
  float ADCSum;            // recob::Hit::SummedADC()
  unsigned short WireNum;  // recob::Hit::WireID().Wire
  unsigned short numHits;  // recob::Hit::Multiplicity()
  unsigned int LoHitID;    // see below
  float LoTime;            // recob::Hit::StartTick()
  float HiTime;            // recob::Hit::EndTick()
  short InClus;            // dropped; see below
  geo::WireID WirID;       // recob::Hit::WireID()
  recob::Wire const* Wire; // dropped; see below
};

The uncertainty on RMS has been dropped for good.

The LoHitID member used to mean the index of the first hit in the "hit train" (that is the set of hits extracted from the same region of interest). That is a concept that is not portable. If your hit list is still the original one as produced by this algorithm, or if at least the hits from the same train are stored sorted and contiguously, for a hit with index iHit, the equivalent value of LoHitID is iHit - hit.LocalIndex().

There is no pointer to the wire any more in recob::Hit. The wire can be obtained through associations, that are typically produced by the art module that runs CCHitFinderAlg (e.g. CCHitFinder). The channel ID is also directly available as recob::Hit::Channel().

Definition at line 78 of file CCHitFinderAlg.h.

Constructor & Destructor Documentation

hit::CCHitFinderAlg::CCHitFinderAlg ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 43 of file CCHitFinderAlg.cxx.

References art::errors::Configuration, fChiNorms, fChiSplit, FinalFitStats, fMaxBumps, fMaxXtraHits, fMinPeak, fMinRMS, fStudyHits, fUseChannelFilter, fUseFastFit, fUTickRange, fUWireRange, fVTickRange, fVWireRange, fWTickRange, fWWireRange, MaxGaussians, MF_LOG_WARNING, hit::CCHitFinderAlg::FitStats_t::Reset(), and TriedFitStats.

    : FitCache{"GausFitCache_CCHitFinderAlg"}
  {
    if (pset.has_key("MinSigInd"))
      throw art::Exception(art::errors::Configuration)
        << "CCHitFinderAlg: Using no-longer-valid fcl input: MinSigInd, MinSigCol, etc";

    fMinPeak = pset.get<std::vector<float>>("MinPeak");
    fMinRMS = pset.get<std::vector<float>>("MinRMS");
    fMaxBumps = pset.get<unsigned short>("MaxBumps");
    fMaxXtraHits = pset.get<unsigned short>("MaxXtraHits");
    fChiSplit = pset.get<float>("ChiSplit");
    fChiNorms = pset.get<std::vector<float>>("ChiNorms");
    fUseFastFit = pset.get<bool>("UseFastFit", false);
    fUseChannelFilter = pset.get<bool>("UseChannelFilter", true);
    fStudyHits = pset.get<bool>("StudyHits", false);
    // The following variables are only used in StudyHits mode
    fUWireRange = pset.get<std::vector<short>>("UWireRange");
    fUTickRange = pset.get<std::vector<short>>("UTickRange");
    fVWireRange = pset.get<std::vector<short>>("VWireRange");
    fVTickRange = pset.get<std::vector<short>>("VTickRange");
    fWWireRange = pset.get<std::vector<short>>("WWireRange");
    fWTickRange = pset.get<std::vector<short>>("WTickRange");

    if (fMinPeak.size() != fMinRMS.size()) {
      mf::LogError("CCTF") << "MinPeak size != MinRMS size";
      return;
    }

    if (fMaxBumps > MaxGaussians) {
      // MF_LOG_WARNING will point the user to this line of code.
      // Any value of MaxGaussians can be used.
      // That value is defined in the header file.
      MF_LOG_WARNING("CCHitFinderAlg")
        << "CCHitFinder algorithm is currently hard-coded to support at most " << MaxGaussians
        << " bumps per region of interest, but " << fMaxBumps << " have been requested.\n"
        << "We are forcing the parameter to " << MaxGaussians
        << ". If this is not acceptable, increase CCHitFinderAlg::MaxGaussians"
        << " value and recompile.";
      fMaxBumps = MaxGaussians;
    } // if too many gaussians

    FinalFitStats.Reset(MaxGaussians);
    TriedFitStats.Reset(MaxGaussians);

    // sanity check for StudyHits mode
    if (fStudyHits) {
      if (fUWireRange.size() != 2 || fUTickRange.size() != 2 || fVWireRange.size() != 2 ||
          fVTickRange.size() != 2 || fWWireRange.size() != 2 || fWTickRange.size() != 2) {
        mf::LogError("CCHF") << "Invalid vector size for StudyHits. Must be 2";
      }
    } // fStudyHits
  }

Member Function Documentation

bool hit::CCHitFinderAlg::FastGaussianFit ( unsigned short  npt, float const *  ticks, float const *  signl, std::array< double, 3 > &  params, std::array< double, 3 > &  paramerrors, float &  chidof  )

staticprivate

Performs a "fast" fit.

Parameters

npt	number of points to be fitted
ticks	tick coordinates
signl	signal amplitude
params	an array where the fit parameters will be stored
paramerrors	an array where the fit parameter errors will be stored
chidof	a variable where to store chi^2 over degrees of freedom

Returns

whether the fit was successful or not

Note that the fit will bail out and rteurn false if any of the input signal amplitudes is zero or negative.

Also note that currently the chi^2 is not the one from comparing the Gaussian to the signal, but from comparing a fitted parabola to the logarithm of the signal.

Definition at line 234 of file CCHitFinderAlg.cxx.

References lar::util::GaussianFit< T >::add(), lar::util::GaussianFit< T >::ChiSquare(), lar::util::GaussianFit< T >::FillResults(), MF_LOG_DEBUG, lar::util::GaussianFit< T >::NDF(), and par.

Referenced by FitNG().

  {
    // parameters: amplitude, mean, sigma

    lar::util::GaussianFit<double> fitter; // probably "double" is overdoing

    // apply a time shift so that the center of the time interval is 0
    const float time_shift = (ticks[npt - 1] + ticks[0]) / 2.F;

    // fill the input data, no uncertainty
    for (size_t i = 0; i < npt; ++i) {
      if (signl[i] <= 0) {
        MF_LOG_DEBUG("CCHitFinderAlg")
          << "Non-positive charge encountered. Backing up to ROOT fit.";
        return false;
      }
      // we could freely add a Poisson uncertainty (as third parameter)
      fitter.add(ticks[i] - time_shift, signl[i]);
    } // for

    // we might have found that we don't want the fast fit after all...
    if (!fitter.FillResults(params, paramerrors)) {
      // something went wrong...
      MF_LOG_DEBUG("CCHitFinderAlg") << "Fast Gaussian fit failed.";
      return false;
    }

    // note that this is not the full chi^2, but it is the chi^2 of the
    // parabolic fit underlying the Gaussian one
    const double chi2 = fitter.ChiSquare();
    chidof = chi2 / fitter.NDF();

    // remove the time shift
    params[1] += time_shift; // mean

    // GP: inflate the uncertainties on the fit parameters according to chi2/NDF
    // (not sure if this is in any way justified)
    if (chidof > 1.)
      for (double& par : paramerrors)
        par *= std::sqrt(chidof);

    return true;
  } // FastGaussianFit()

void hit::CCHitFinderAlg::FitNG ( unsigned short  nGaus, unsigned short  npt, float *  ticks, float *  signl  )

private

Definition at line 284 of file CCHitFinderAlg.cxx.

References util::abs(), hit::CCHitFinderAlg::FitStats_t::AddFast(), hit::CCHitFinderAlg::FitStats_t::AddMultiGaus(), bumps, chidof, chinorm, dof, FastGaussianFit(), FinalFitStats, fMinPeak, fMinRMS, fUseFastFit, par, parerr, thePlane, and TriedFitStats.

Referenced by RunCCHitFinder().

  {
    // Fit the signal to n Gaussians

    dof = npt - 3 * nGaus;

    chidof = 9999.;

    if (dof < 3) return;
    if (bumps.size() == 0) return;

    // load the fit into a temp vector
    std::vector<double> partmp;
    std::vector<double> partmperr;

    //
    // if it is possible, we try first with the quick single Gaussian fit
    //
    TriedFitStats.AddMultiGaus(nGaus);

    bool bNeedROOTfit = (nGaus > 1) || !fUseFastFit;
    if (!bNeedROOTfit) {
      // so, we need only one puny Gaussian;
      std::array<double, 3> params, paramerrors;

      TriedFitStats.AddFast();

      if (FastGaussianFit(npt, ticks, signl, params, paramerrors, chidof)) {
        // success? copy the results in the proper structures
        partmp.resize(3);
        std::copy(params.begin(), params.end(), partmp.begin());
        partmperr.resize(3);
        std::copy(paramerrors.begin(), paramerrors.end(), partmperr.begin());
      }
      else
        bNeedROOTfit = true; // if we fail, let's schedule ROOT to back us up

      if (!bNeedROOTfit) FinalFitStats.AddFast();

    } // if we don't need ROOT to fit

    if (bNeedROOTfit) {
      // we may land here either because the simple Gaussian fit did not work (either
      // failed, or we chose not to trust it) or because the fit is multi-Gaussian

      // define the fit string to pass to TF1

      std::stringstream numConv;
      std::string eqn = "gaus";
      if (nGaus > 1) eqn = "gaus(0)";
      for (unsigned short ii = 3; ii < nGaus * 3; ii += 3) {
        eqn.append(" + gaus(");
        numConv.str("");
        numConv << ii;
        eqn.append(numConv.str());
        eqn.append(")");
      }

      auto Gn = std::make_unique<TF1>("gn", eqn.c_str());
      TGraph* fitn = new TGraph(npt, ticks, signl);

      // put in the bump parameters. Assume that nGaus >= bumps.size()
      for (unsigned short ii = 0; ii < bumps.size(); ++ii) {
        unsigned short index = ii * 3;
        unsigned short bumptime = bumps[ii];
        double amp = signl[bumptime];
        Gn->SetParameter(index, amp);
        Gn->SetParLimits(index, 0., 9999.);
        Gn->SetParameter(index + 1, (double)bumptime);
        Gn->SetParLimits(index + 1, 0, (double)npt);
        Gn->SetParameter(index + 2, (double)fMinRMS[thePlane]);
        Gn->SetParLimits(index + 2, 1., 3 * (double)fMinRMS[thePlane]);
      } // ii bumps

      // search for other bumps that may be hidden by the already found ones
      for (unsigned short ii = bumps.size(); ii < nGaus; ++ii) {
        // bump height must exceed fMinPeak
        float big = fMinPeak[thePlane];
        unsigned short imbig = 0;
        for (unsigned short jj = 0; jj < npt; ++jj) {
          float diff = signl[jj] - Gn->Eval((Double_t)jj, 0, 0, 0);
          if (diff > big) {
            big = diff;
            imbig = jj;
          }
        } // jj
        if (imbig > 0) {
          // set the parameters for the bump
          unsigned short index = ii * 3;
          Gn->SetParameter(index, (double)big);
          Gn->SetParLimits(index, 0., 9999.);
          Gn->SetParameter(index + 1, (double)imbig);
          Gn->SetParLimits(index + 1, 0, (double)npt);
          Gn->SetParameter(index + 2, (double)fMinRMS[thePlane]);
          Gn->SetParLimits(index + 2, 1., 5 * (double)fMinRMS[thePlane]);
        } // imbig > 0
      }   // ii

      // W = set weights to 1, N = no drawing or storing, Q = quiet
      // B = bounded parameters
      fitn->Fit(&*Gn, "WNQB");

      for (unsigned short ipar = 0; ipar < 3 * nGaus; ++ipar) {
        partmp.push_back(Gn->GetParameter(ipar));
        partmperr.push_back(Gn->GetParError(ipar));
      }
      chidof = Gn->GetChisquare() / (dof * chinorm);

      delete fitn;
      //  delete Gn;

    } // if ROOT fit

    // Sort by increasing time if necessary
    if (nGaus > 1) {
      std::vector<std::pair<unsigned short, unsigned short>> times;
      // fill the sort vector
      for (unsigned short ii = 0; ii < nGaus; ++ii) {
        unsigned short index = ii * 3;
        times.push_back(std::make_pair(partmp[index + 1], ii));
      } // ii
      std::sort(times.begin(), times.end());
      // see if re-arranging is necessary
      bool sortem = false;
      for (unsigned short ii = 0; ii < nGaus; ++ii) {
        if (times[ii].second != ii) {
          sortem = true;
          break;
        }
      } // ii
      if (sortem) {
        // temp temp vectors for putting things in the right time order
        std::vector<double> partmpt;
        std::vector<double> partmperrt;
        for (unsigned short ii = 0; ii < nGaus; ++ii) {
          unsigned short index = times[ii].second * 3;
          partmpt.push_back(partmp[index]);
          partmpt.push_back(partmp[index + 1]);
          partmpt.push_back(partmp[index + 2]);
          partmperrt.push_back(partmperr[index]);
          partmperrt.push_back(partmperr[index + 1]);
          partmperrt.push_back(partmperr[index + 2]);
        } // ii
        partmp = partmpt;
        partmperr = partmperrt;
      } // sortem
    }   // nGaus > 1

    // ensure that the fit is reasonable
    bool fitok = true;
    for (unsigned short ii = 0; ii < nGaus; ++ii) {
      unsigned short index = ii * 3;
      // ensure that the fitted time is within the signal bounds
      short fittime = partmp[index + 1];
      if (fittime < 0 || fittime > npt - 1) {
        fitok = false;
        break;
      }
      // ensure that the signal peak is large enough
      if (partmp[index] < fMinPeak[thePlane]) {
        fitok = false;
        break;
      }
      // ensure that the RMS is large enough but not too large
      float rms = partmp[index + 2];
      if (rms < 0.5 * fMinRMS[thePlane] || rms > 5 * fMinRMS[thePlane]) {
        fitok = false;
        break;
      }
      // ensure that the hits are not too similar in time (< 2 ticks)
      for (unsigned short jj = 0; jj < nGaus; ++jj) {
        if (jj == ii) continue;
        unsigned short jndex = jj * 3;
        float timediff = std::abs(partmp[jndex + 1] - partmp[index + 1]);
        if (timediff < 2.) {
          fitok = false;
          break;
        }
      }
      if (!fitok) break;
    }

    if (fitok) {
      par = partmp;
      parerr = partmperr;
    }
    else {
      chidof = 9999.;
      dof = -1;
    }
  } // FitNG

void hit::CCHitFinderAlg::MakeCrudeHit ( unsigned short  npt, float *  ticks, float *  signl  )

private

Definition at line 477 of file CCHitFinderAlg.cxx.

References chidof, dof, pmtana::mean(), par, parerr, and Sqrt2Pi.

Referenced by RunCCHitFinder().

  {
    // make a single crude hit if fitting failed
    float sumS = 0.;
    float sumST = 0.;
    for (unsigned short ii = 0; ii < npt; ++ii) {
      sumS += signl[ii];
      sumST += signl[ii] * ticks[ii];
    }
    float mean = sumST / sumS;
    float rms = 0.;
    for (unsigned short ii = 0; ii < npt; ++ii) {
      float arg = ticks[ii] - mean;
      rms += signl[ii] * arg * arg;
    }
    rms = std::sqrt(rms / sumS);
    float amp = sumS / (Sqrt2Pi * rms);
    par.clear();
    par.push_back(amp);
    par.push_back(mean);
    par.push_back(rms);
    // need to do the errors better
    parerr.clear();
    float amperr = npt;
    float meanerr = std::sqrt(1 / sumS);
    float rmserr = 0.2 * rms;
    parerr.push_back(amperr);
    parerr.push_back(meanerr);
    parerr.push_back(rmserr);
    chidof = 9999.;
    dof = -1;
  } // MakeCrudeHit

template<typename Stream >

void hit::CCHitFinderAlg::PrintStats

(

Stream &

out

)

const

Print the fit statistics.

Definition at line 212 of file CCHitFinderAlg.h.

References hit::CCHitFinderAlg::FitStats_t::FastFits, FinalFitStats, fUseFastFit, MaxGaussians, hit::CCHitFinderAlg::FitStats_t::MultiGausFits, and TriedFitStats.

Referenced by hit::HitFinder::endJob(), and YieldHits().

{
  out << "CCHitFinderAlg fit statistics:";
  if (fUseFastFit) {
    out << "\n  fast 1-Gaussian fits: " << FinalFitStats.FastFits << " succeeded ("
        << TriedFitStats.FastFits << " tried)";
  }
  else
    out << "\n  fast 1-Gaussian fits: disabled";

  for (unsigned int nGaus = 1; nGaus < MaxGaussians; ++nGaus) {
    if (TriedFitStats.MultiGausFits[nGaus - 1] == 0) continue;
    out << "\n  " << nGaus << "-Gaussian fits: " << FinalFitStats.MultiGausFits[nGaus - 1]
        << " accepted (" << TriedFitStats.MultiGausFits[nGaus - 1] << " tried)";
  } // for nGaus
  if (TriedFitStats.MultiGausFits.back() > 0) {
    out << "\n  " << FinalFitStats.MultiGausFits.size()
        << "-Gaussian fits or higher: " << FinalFitStats.MultiGausFits.back() << " accepted ("
        << TriedFitStats.MultiGausFits.back() << " tried)";
  }
  out << std::endl;
} // CCHitFinderAlg::FitStats_t::Print()

void hit::CCHitFinderAlg::RunCCHitFinder

(

std::vector< recob::Wire > const &

Wires

)

Definition at line 98 of file CCHitFinderAlg.cxx.

References hit::CCHitFinderAlg::FitStats_t::AddMultiGaus(), allhits, bumps, recob::Wire::Channel(), geo::WireReadoutGeom::ChannelToWire(), chidof, chinorm, dof, fChiNorms, fChiSplit, FinalFitStats, FitNG(), fMaxBumps, fMaxXtraHits, fMinPeak, fMinRMS, fStudyHits, MakeCrudeHit(), SelRAT, recob::Wire::Signal(), geo::WireReadoutGeom::SignalType(), StoreHits(), StudyHits(), theChannel, thePlane, theWireNum, and wireReadoutGeom.

Referenced by hit::HitFinder::produce(), and cluster::ClusterCrawler::produce().

  {
    allhits.clear();

    constexpr unsigned short maxticks = 1000;
    float* ticks = new float[maxticks];
    // define the ticks array used for fitting
    for (unsigned short ii = 0; ii < maxticks; ++ii) {
      ticks[ii] = ii;
    }
    float* signl = new float[maxticks];
    float adcsum = 0;
    // initialize the vectors for the hit study
    if (fStudyHits) StudyHits(0);
    bool first;

    lariov::ChannelStatusProvider const& channelStatus =
      art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider();

    for (size_t wireIter = 0; wireIter < Wires.size(); wireIter++) {
      recob::Wire const& theWire = Wires[wireIter];
      theChannel = theWire.Channel();
      // ignore bad channels
      if (channelStatus.IsBad(theChannel)) continue;

      std::vector<geo::WireID> wids = wireReadoutGeom->ChannelToWire(theChannel);
      thePlane = wids[0].Plane;
      if (thePlane > fMinPeak.size() - 1) {
        mf::LogError("CCHF") << "MinPeak vector too small for plane " << thePlane;
        return;
      }
      theWireNum = wids[0].Wire;
      HitChannelInfo_t WireInfo{&theWire, wids[0], wireReadoutGeom->SignalType(theWire.Channel())};

      // minimum number of time samples
      unsigned short minSamples = 2 * fMinRMS[thePlane];

      // factor used to normalize the chi/dof fits for each plane
      chinorm = fChiNorms[thePlane];

      std::vector<float> signal(theWire.Signal());

      unsigned short nabove = 0;
      unsigned short tstart = 0;
      unsigned short maxtime = signal.size() - 2;
      // find the min time when the signal is below threshold
      unsigned short mintime = 3;
      for (unsigned short time = 3; time < maxtime; ++time) {
        if (signal[time] < fMinPeak[thePlane]) {
          mintime = time;
          break;
        }
      }
      for (unsigned short time = mintime; time < maxtime; ++time) {
        if (signal[time] > fMinPeak[thePlane]) {
          if (nabove == 0) tstart = time;
          ++nabove;
        }
        else {
          // check for a wide enough signal above threshold
          if (nabove > minSamples) {
            // skip this wire if the RAT is too long
            if (nabove > maxticks)
              mf::LogError("CCHitFinder")
                << "Long RAT " << nabove << " " << maxticks << " No signal on wire " << theWireNum
                << " after time " << time;
            if (nabove > maxticks) break;
            unsigned short npt = 0;
            // look for bumps to inform the fit
            bumps.clear();
            adcsum = 0;
            for (unsigned short ii = tstart; ii < time; ++ii) {
              signl[npt] = signal[ii];
              adcsum += signl[npt];
              if (signal[ii] > signal[ii - 1] && signal[ii - 1] > signal[ii - 2] &&
                  signal[ii] > signal[ii + 1] && signal[ii + 1] > signal[ii + 2])
                bumps.push_back(npt);
              ++npt;
            }
            // decide if this RAT should be studied
            if (fStudyHits) StudyHits(1, npt, signl, tstart);
            // just make a crude hit if too many bumps
            if (bumps.size() > fMaxBumps) {
              MakeCrudeHit(npt, ticks, signl);
              StoreHits(tstart, npt, WireInfo, adcsum);
              nabove = 0;
              continue;
            }
            // start looking for hits with the found bumps
            unsigned short nHitsFit = bumps.size();
            unsigned short nfit = 0;
            chidof = 0.;
            dof = -1;
            bool HitStored = false;
            unsigned short nMaxFit = bumps.size() + fMaxXtraHits;
            // only used in StudyHits mode
            first = true;
            while (nHitsFit <= nMaxFit) {

              FitNG(nHitsFit, npt, ticks, signl);
              if (fStudyHits && first && SelRAT) {
                first = false;
                StudyHits(2, npt, signl, tstart);
              }
              // good chisq so store it
              if (chidof < fChiSplit) {
                StoreHits(tstart, npt, WireInfo, adcsum);
                HitStored = true;
                break;
              }
              // the previous fit was better, so revert to it and store it
              ++nHitsFit;
              ++nfit;
            } // nHitsFit < fMaxXtraHits
            if (!HitStored && npt < maxticks) {
              // failed all fitting. Make a crude hit
              MakeCrudeHit(npt, ticks, signl);
              StoreHits(tstart, npt, WireInfo, adcsum);
            }
            else if (nHitsFit > 0)
              FinalFitStats.AddMultiGaus(nHitsFit);
          } // nabove > minSamples
          nabove = 0;
        } // signal < fMinPeak
      }   // time
    }     // wireIter

    // print out
    if (fStudyHits) StudyHits(4);

    delete[] ticks;
    delete[] signl;

  } //RunCCHitFinder

void hit::CCHitFinderAlg::StoreHits ( unsigned short  TStart, unsigned short  npt, HitChannelInfo_t  info, float  adcsum  )

private

Definition at line 511 of file CCHitFinderAlg.cxx.

References allhits, recob::Wire::Channel(), chidof, dof, fStudyHits, hiTime, loTime, par, parerr, hit::CCHitFinderAlg::HitChannelInfo_t::sigType, Sqrt2Pi, SqrtPi, StudyHits(), recob::Wire::View(), hit::CCHitFinderAlg::HitChannelInfo_t::wire, and hit::CCHitFinderAlg::HitChannelInfo_t::wireID.

Referenced by RunCCHitFinder().

  {
    // store the hits in the struct
    size_t nhits = par.size() / 3;

    if (allhits.max_size() - allhits.size() < nhits) {
      mf::LogError("CCHitFinder") << "Too many hits: existing " << allhits.size() << " plus new "
                                  << nhits << " beyond the maximum " << allhits.max_size();
      return;
    }

    if (nhits == 0) return;

    // fill RMS for single hits
    if (fStudyHits) StudyHits(3);

    const float loTime = TStart;
    const float hiTime = TStart + npt;

    // Find sum of the areas of all Gaussians
    float gsum = 0.;
    for (size_t hit = 0; hit < nhits; ++hit) {
      const unsigned short index = 3 * hit;
      gsum += Sqrt2Pi * par[index] * par[index + 2];
    }
    for (size_t hit = 0; hit < nhits; ++hit) {
      const size_t index = 3 * hit;
      const float charge = Sqrt2Pi * par[index] * par[index + 2];
      const float charge_err =
        SqrtPi * (parerr[index] * par[index + 2] + par[index] * parerr[index + 2]);

      allhits.emplace_back(info.wire->Channel(),    // channel
                           loTime,                  // start_tick
                           hiTime,                  // end_tick
                           par[index + 1] + TStart, // peak_time
                           parerr[index + 1],       // sigma_peak_time
                           par[index + 2],          // rms
                           par[index],              // peak_amplitude
                           parerr[index],           // sigma_peak_amplitude
                           adcsum * charge / gsum,  // ROIsummedADC
                           adcsum * charge / gsum,  // HitsummedADC  NOT CORRECTLY FILLED
                           charge,                  // hit_integral
                           charge_err,              // hit_sigma_integral
                           nhits,                   // multiplicity
                           hit,                     // local_index
                           chidof,                  // goodness_of_fit
                           dof,                     // dof
                           info.wire->View(),       // view
                           info.sigType,            // signal_type
                           info.wireID              // wireID
      );
    } // hit
  }   // StoreHits

void hit::CCHitFinderAlg::StudyHits ( unsigned short  flag, unsigned short  npt = 0, float *  signl = nullptr, unsigned short  tstart = 0  )

private

Definition at line 569 of file CCHitFinderAlg.cxx.

References util::abs(), bumpChi, bumpCnt, bumpRMS, bumps, chidof, fChiNorms, fMinPeak, fUTickRange, fUWireRange, fVTickRange, fVWireRange, fWTickRange, fWWireRange, hitCnt, hiTime, hitRMS, hiWire, loTime, loWire, par, RATCnt, art::right(), SelRAT, sum, thePlane, and theWireNum.

Referenced by RunCCHitFinder(), and StoreHits().

  {
    // study hits in user-selected ranges of wires and ticks in each plane. The user should identify
    // a shallow-angle isolated track, e.g. using the event display, to determine the wire/tick ranges.
    // One hit should be reconstructed on each wire when the hit finding fcl parameters are set correctly.
    // The intent of this study is to determine the correct fcl parameters. The flag variable determines
    // the operation performed.
    // flag = 0: Initialize study vectors
    // flag = 1: Set SelRat true if the Region Above Threshold resides within a wire/hit range
    // flag = 2: Find the maximum signal and calculate the RMS. Also find the low and high ticks of signals
    //           in the wire range to allow a later calculation of the track angle. This isn't strictly
    //           necessary for the study and presumes that the user has selected compatible regions in each plane.
    // flag = 3: Accumulate the RMS from the first Gaussian fit
    // flag = 4: Calculate recommended fcl parameters and print the results to the screen

    // init
    if (flag == 0) {
      for (unsigned short ipl = 0; ipl < 3; ++ipl) {
        // Average chisq of the first fit on a single bump in each plane
        bumpChi.push_back(0.);
        // Average RMS of the dump
        bumpRMS.push_back(0.);
        // The number of single bumps in each plane
        bumpCnt.push_back(0.);
        // number of RATs
        RATCnt.push_back(0);
        // The number of single hits found in each plane
        hitCnt.push_back(0.);
        // Average reconstructed hit RMS
        hitRMS.push_back(0.);
        // lo/hi wire/time
        loWire.push_back(9999.);
        loTime.push_back(0.);
        hiWire.push_back(-1.);
        hiTime.push_back(0.);
      } // ii
      return;
    } // flag == 0

    if (flag == 1) {
      SelRAT = false;
      if (thePlane == 0) {
        if (theWireNum > fUWireRange[0] && theWireNum < fUWireRange[1] && tstart > fUTickRange[0] &&
            tstart < fUTickRange[1]) {
          SelRAT = true;
          RATCnt[thePlane] += 1;
        }
        return;
      } // thePlane == 0
      if (thePlane == 1) {
        if (theWireNum > fVWireRange[0] && theWireNum < fVWireRange[1] && tstart > fVTickRange[0] &&
            tstart < fVTickRange[1]) {
          SelRAT = true;
          RATCnt[thePlane] += 1;
        }
        return;
      } // thePlane == 1
      if (thePlane == 2) {
        if (theWireNum > fWWireRange[0] && theWireNum < fWWireRange[1] && tstart > fWTickRange[0] &&
            tstart < fWTickRange[1]) {
          SelRAT = true;
          RATCnt[thePlane] += 1;
        }
        return;
      } // thePlane == 2
    }   // flag == 1

    if (flag == 2) {
      if (!SelRAT) return;
      // in this section we find the low/hi wire/time for a signal. This can be used to calculate
      // the slope dT/dW to study hit width, fraction of crude hits, etc vs dT/dW
      float big = 0.;
      float imbig = 0.;
      for (unsigned short ii = 0; ii < npt; ++ii) {
        if (signl[ii] > big) {
          big = signl[ii];
          imbig = ii;
        }
      } // ii
      // require a significant PH
      if (big > fMinPeak[0]) {
        // get the Lo info
        if (theWireNum < loWire[thePlane]) {
          loWire[thePlane] = theWireNum;
          loTime[thePlane] = tstart + imbig;
        }
        // get the Hi info
        if (theWireNum > hiWire[thePlane]) {
          hiWire[thePlane] = theWireNum;
          hiTime[thePlane] = tstart + imbig;
        }
      } // big > fMinPeak[0]
      if (bumps.size() == 1 && chidof < 9999.) {
        bumpCnt[thePlane] += bumps.size();
        bumpChi[thePlane] += chidof;
        // calculate the average bin
        float sumt = 0.;
        float sum = 0.;
        for (unsigned short ii = 0; ii < npt; ++ii) {
          sum += signl[ii];
          sumt += signl[ii] * ii;
        } // ii
        float aveb = sumt / sum;
        // now calculate the RMS
        sumt = 0.;
        for (unsigned short ii = 0; ii < npt; ++ii) {
          float dbin = (float)ii - aveb;
          sumt += signl[ii] * dbin * dbin;
        } // ii
        bumpRMS[thePlane] += std::sqrt(sumt / sum);
      } // bumps.size() == 1 && chidof < 9999.
      return;
    } // flag == 2

    // fill info for single hits
    if (flag == 3) {
      if (!SelRAT) return;
      if (par.size() == 3) {
        hitCnt[thePlane] += 1;
        hitRMS[thePlane] += par[2];
      }
      return;
    }

    if (flag == 4) {
      // The goal is to adjust the fcl inputs so that the number of single
      // hits found is ~equal to the number of single bumps found for shallow
      // angle tracks. The ChiNorm inputs should be adjusted so the average
      //  chisq/DOF is ~1 in each plane.
      std::cout << "Check lo and hi W/T for each plane" << std::endl;
      for (unsigned short ipl = 0; ipl < 3; ++ipl) {
        std::cout << ipl << " lo " << loWire[ipl] << " " << loTime[ipl] << " hi " << hiWire[ipl]
                  << " " << hiTime[ipl] << std::endl;
      }
      std::cout << " ipl nRAT bCnt   bChi   bRMS hCnt   hRMS  dT/dW New_ChiNorm" << std::endl;
      for (unsigned short ipl = 0; ipl < 3; ++ipl) {
        if (bumpCnt[ipl] > 0) {
          bumpChi[ipl] = bumpChi[ipl] / (float)bumpCnt[ipl];
          bumpRMS[ipl] = bumpRMS[ipl] / (float)bumpCnt[ipl];
          hitRMS[ipl] = hitRMS[ipl] / (float)hitCnt[ipl];
          // calculate the slope
          float dTdW = std::abs((hiTime[ipl] - loTime[ipl]) / (hiWire[ipl] - loWire[ipl]));
          std::cout << ipl << std::right << std::setw(5) << RATCnt[ipl] << std::setw(5)
                    << bumpCnt[ipl] << std::setw(7) << std::fixed << std::setprecision(2)
                    << bumpChi[ipl] << std::setw(7) << bumpRMS[ipl] << std::setw(7) << hitCnt[ipl]
                    << std::setw(7) << std::setprecision(1) << hitRMS[ipl] << std::setw(7) << dTdW
                    << std::setw(7) << std::setprecision(2) << bumpChi[ipl] * fChiNorms[ipl]
                    << std::endl;
        } //
      }   // ipl
      std::cout << "nRAT is the number of Regions Above Threshold (RAT) used in the study.\n";
      std::cout << "bCnt is the number of single bumps that were successfully fitted \n";
      std::cout << "bChi is the average chisq/DOF of the first fit\n";
      std::cout << "bRMS is the average calculated RMS of the bumps\n";
      std::cout << "hCnt is the number of RATs that have a single hit\n";
      std::cout << "hRMS is the average RMS from the Gaussian fit -> use this value for "
                   "fMinRMS[plane] in the fcl file\n";
      std::cout << "dTdW is the slope of the track\n";
      std::cout
        << "New_ChiNorm is the recommended values of ChiNorm that should be used in the fcl file\n";
      bumpChi.clear();
      bumpRMS.clear();
      bumpCnt.clear();
      RATCnt.clear();
      hitRMS.clear();
      hitCnt.clear();
      loWire.clear();
      loTime.clear();
      hiWire.clear();
      hiTime.clear();
    }
  } // StudyHits

std::vector<recob::Hit>&& hit::CCHitFinderAlg::YieldHits ( )

inline

Returns (and loses) the collection of reconstructed hits.

Definition at line 86 of file CCHitFinderAlg.h.

References PrintStats().

Referenced by hit::HitFinder::produce(), and cluster::ClusterCrawler::produce().

{ return std::move(allhits); }

Member Data Documentation

std::vector<recob::Hit> hit::CCHitFinderAlg::allhits

Definition at line 80 of file CCHitFinderAlg.h.

Referenced by RunCCHitFinder(), and StoreHits().

std::vector<float> hit::CCHitFinderAlg::bumpChi

private

Definition at line 151 of file CCHitFinderAlg.h.

Referenced by StudyHits().

std::vector<int> hit::CCHitFinderAlg::bumpCnt

private

Definition at line 149 of file CCHitFinderAlg.h.

Referenced by StudyHits().

std::vector<float> hit::CCHitFinderAlg::bumpRMS

private

Definition at line 152 of file CCHitFinderAlg.h.

Referenced by StudyHits().

std::vector<unsigned short> hit::CCHitFinderAlg::bumps

private

Definition at line 126 of file CCHitFinderAlg.h.

Referenced by FitNG(), RunCCHitFinder(), and StudyHits().

float hit::CCHitFinderAlg::chidof

private

Definition at line 124 of file CCHitFinderAlg.h.

Referenced by FitNG(), MakeCrudeHit(), RunCCHitFinder(), StoreHits(), and StudyHits().

float hit::CCHitFinderAlg::chinorm

private

Definition at line 108 of file CCHitFinderAlg.h.

Referenced by FitNG(), and RunCCHitFinder().

int hit::CCHitFinderAlg::dof

private

Definition at line 125 of file CCHitFinderAlg.h.

Referenced by FitNG(), MakeCrudeHit(), RunCCHitFinder(), and StoreHits().

std::vector<float> hit::CCHitFinderAlg::fChgNorms

private

Definition at line 102 of file CCHitFinderAlg.h.

std::vector<float> hit::CCHitFinderAlg::fChiNorms

private

Definition at line 100 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), RunCCHitFinder(), and StudyHits().

float hit::CCHitFinderAlg::fChiSplit

private

Estimated noise error on the Signal.

Definition at line 97 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and RunCCHitFinder().

FitStats_t hit::CCHitFinderAlg::FinalFitStats

private

counts of the good fits

Definition at line 175 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), FitNG(), PrintStats(), and RunCCHitFinder().

GausFitCache hit::CCHitFinderAlg::FitCache

private

a set of functions ready to be used

Definition at line 164 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::fMaxBumps

private

Definition at line 95 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and RunCCHitFinder().

unsigned short hit::CCHitFinderAlg::fMaxXtraHits

private

Definition at line 96 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and RunCCHitFinder().

std::vector<float> hit::CCHitFinderAlg::fMinPeak

private

Definition at line 93 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), FitNG(), RunCCHitFinder(), and StudyHits().

std::vector<float> hit::CCHitFinderAlg::fMinRMS

private

Definition at line 94 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), FitNG(), and RunCCHitFinder().

bool hit::CCHitFinderAlg::fStudyHits

private

Definition at line 141 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), RunCCHitFinder(), and StoreHits().

std::vector<float> hit::CCHitFinderAlg::fTimeOffsets

private

Definition at line 101 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::fUseChannelFilter

private

Definition at line 112 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg().

bool hit::CCHitFinderAlg::fUseFastFit

private

whether to attempt using a fast fit on single gauss.

Definition at line 162 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), FitNG(), and PrintStats().

std::vector<short> hit::CCHitFinderAlg::fUTickRange

private

Definition at line 142 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and StudyHits().

std::vector<short> hit::CCHitFinderAlg::fUWireRange

private

Definition at line 142 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and StudyHits().

std::vector<short> hit::CCHitFinderAlg::fVTickRange

private

Definition at line 143 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and StudyHits().

std::vector<short> hit::CCHitFinderAlg::fVWireRange

private

Definition at line 143 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and StudyHits().

std::vector<short> hit::CCHitFinderAlg::fWTickRange

private

Definition at line 144 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and StudyHits().

std::vector<short> hit::CCHitFinderAlg::fWWireRange

private

Definition at line 144 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and StudyHits().

std::vector<int> hit::CCHitFinderAlg::hitCnt

private

Definition at line 153 of file CCHitFinderAlg.h.

Referenced by StudyHits().

std::vector<float> hit::CCHitFinderAlg::hiTime

private

Definition at line 159 of file CCHitFinderAlg.h.

Referenced by StoreHits(), and StudyHits().

std::vector<float> hit::CCHitFinderAlg::hitRMS

private

Definition at line 154 of file CCHitFinderAlg.h.

Referenced by StudyHits().

std::vector<float> hit::CCHitFinderAlg::hiWire

private

Definition at line 158 of file CCHitFinderAlg.h.

Referenced by StudyHits().

std::vector<float> hit::CCHitFinderAlg::loTime

private

Definition at line 157 of file CCHitFinderAlg.h.

Referenced by StoreHits(), and StudyHits().

std::vector<float> hit::CCHitFinderAlg::loWire

private

Definition at line 156 of file CCHitFinderAlg.h.

Referenced by StudyHits().

constexpr unsigned int hit::CCHitFinderAlg::MaxGaussians = 20

staticprivate

Definition at line 202 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), and PrintStats().

std::vector<double> hit::CCHitFinderAlg::par

private

Definition at line 120 of file CCHitFinderAlg.h.

Referenced by FastGaussianFit(), FitNG(), MakeCrudeHit(), StoreHits(), and StudyHits().

std::vector<double> hit::CCHitFinderAlg::parerr

private

Definition at line 121 of file CCHitFinderAlg.h.

Referenced by FitNG(), MakeCrudeHit(), and StoreHits().

std::vector<double> hit::CCHitFinderAlg::parmax

private

Definition at line 123 of file CCHitFinderAlg.h.

std::vector<double> hit::CCHitFinderAlg::parmin

private

Definition at line 122 of file CCHitFinderAlg.h.

std::vector<int> hit::CCHitFinderAlg::RATCnt

private

Definition at line 150 of file CCHitFinderAlg.h.

Referenced by StudyHits().

bool hit::CCHitFinderAlg::SelRAT

private

Definition at line 160 of file CCHitFinderAlg.h.

Referenced by RunCCHitFinder(), and StudyHits().

constexpr float hit::CCHitFinderAlg::Sqrt2Pi = 2.5066

staticprivate

Definition at line 109 of file CCHitFinderAlg.h.

Referenced by MakeCrudeHit(), and StoreHits().

constexpr float hit::CCHitFinderAlg::SqrtPi = 1.7725

staticprivate

Definition at line 110 of file CCHitFinderAlg.h.

Referenced by StoreHits().

raw::ChannelID_t hit::CCHitFinderAlg::theChannel

private

Definition at line 104 of file CCHitFinderAlg.h.

Referenced by RunCCHitFinder().

unsigned short hit::CCHitFinderAlg::thePlane

private

Definition at line 106 of file CCHitFinderAlg.h.

Referenced by FitNG(), RunCCHitFinder(), and StudyHits().

unsigned short hit::CCHitFinderAlg::theWireNum

private

Definition at line 105 of file CCHitFinderAlg.h.

Referenced by RunCCHitFinder(), and StudyHits().

FitStats_t hit::CCHitFinderAlg::TriedFitStats

private

counts of the tried fits

Definition at line 176 of file CCHitFinderAlg.h.

Referenced by CCHitFinderAlg(), FitNG(), and PrintStats().

geo::WireReadoutGeom const* hit::CCHitFinderAlg::wireReadoutGeom

private

Initial value:

{
      &art::ServiceHandle<geo::WireReadout const>()->Get()}

Definition at line 114 of file CCHitFinderAlg.h.

Referenced by RunCCHitFinder().

---

The documentation for this class was generated from the following files:

• larreco/v10_01_07/source/larreco/RecoAlg/CCHitFinderAlg.h
• larreco/v10_01_07/source/larreco/RecoAlg/CCHitFinderAlg.cxx