hit::CCHitFinderAlg Class Reference

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

#include <CCHitFinderAlg.h>

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

Public Member Functions
	CCHitFinderAlg (fhicl::ParameterSet const &pset)
virtual	~CCHitFinderAlg ()=default
virtual void	reconfigure (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 *ticks=0, float *signl=0, 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
art::ServiceHandle< geo::Geometry const >	geom
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...
std::unique_ptr< 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

)

Definition at line 45 of file CCHitFinderAlg.cxx.

                                                             :
    FitCache(new
      GausFitCache // run-time, on demand TFormula cache
    //  CompiledGausFitCache<MaxGaussians> // precompiled Gaussian set
    //  CompiledTruncatedGausFitCache<MaxGaussians, 4> // precompiled truncated Gaussian set
    //  CompiledTruncatedGausFitCache<MaxGaussians, 5> // precompiled truncated Gaussian set
      ("GausFitCache_CCHitFinderAlg")
      )
  {
    this->reconfigure(pset);
  }

virtual hit::CCHitFinderAlg::~CCHitFinderAlg ( )

virtualdefault

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 276 of file CCHitFinderAlg.cxx.

    {
    // 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 325 of file CCHitFinderAlg.cxx.

  {
    // 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());
      /*
      TF1* Gn = FitCache->Get(nGaus);
      */
      TGraph *fitn = new TGraph(npt, ticks, signl);
  /*
    if(prt) mf::LogVerbatim("CCHitFinder")
      <<"FitNG nGaus "<<nGaus<<" nBumps "<<bumps.size();
  */
      // 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]);
  /*
    if(prt) mf::LogVerbatim("CCHitFinder")<<"Bump params "<<ii<<" "<<(short)amp
      <<" "<<(int)bumptime<<" "<<(int)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) {
  /*
    if(prt) mf::LogVerbatim("CCHitFinder")<<"Found bump "<<ii<<" "<<(short)big
      <<" "<<imbig;
  */
          // 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
/*
  if(prt) {
    mf::LogVerbatim("CCHitFinder")<<"Fit "<<nGaus<<" chi "<<chidof
      <<" npars "<<partmp.size();
    mf::LogVerbatim("CCHitFinder")<<"pars    errs ";
    for(unsigned short ii = 0; ii < partmp.size(); ++ii) {
      mf::LogVerbatim("CCHitFinder")<<ii<<" "<<partmp[ii]<<" "
        <<partmperr[ii];
    }
  }
*/
    // 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;
//      if(prt) mf::LogVerbatim("CCHitFinder")<<"Bad fit parameters";
    }

    return;
  } // FitNG

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

private

Definition at line 545 of file CCHitFinderAlg.cxx.

  {
    // 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();
/*
  if(prt) mf::LogVerbatim("CCHitFinder")<<"Crude hit Amp "<<(int)amp<<" mean "
    <<(int)mean<<" rms "<<rms;
*/
    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);
/*
  if(prt) mf::LogVerbatim("CCHitFinder")<<" errors Amp "<<amperr<<" mean "
    <<meanerr<<" rms "<<rmserr;
*/
    chidof = 9999.;
    dof = -1;
  } // MakeCrudeHit

template<typename Stream >

void hit::CCHitFinderAlg::PrintStats

(

Stream &

out

)

const

Print the fit statistics.

Definition at line 231 of file CCHitFinderAlg.h.

                                                    {

  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::reconfigure ( fhicl::ParameterSet const &  pset )

virtual

Definition at line 57 of file CCHitFinderAlg.cxx.

  {
    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";
        return;
      }
    } // fStudyHits

  }

void hit::CCHitFinderAlg::RunCCHitFinder

(

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

Wires

)

Definition at line 122 of file CCHitFinderAlg.cxx.

                                                                       {

    allhits.clear();

    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;

//    prt = false;
    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;
/*
      geo::SigType_t SigType = geom->SignalType(theChannel);
      minSig = 0.;
      minRMS = 0.;
      if(SigType == geo::kInduction){
        minSig = fMinSigInd;
        minRMS = fMinRMSInd;
      }//<-- End if Induction Plane
      else if(SigType == geo::kCollection){
        minSig = fMinSigCol;
        minRMS  = fMinRMSCol;
      }//<-- End if Collection Plane
*/

      std::vector<geo::WireID> wids = geom->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], *geom);

      // 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];

      // edit this line to debug hit fitting on a particular plane/wire
//      prt = (thePlane == 1 && theWireNum == 839);
      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);
//  if(prt) mf::LogVerbatim("CCHitFinder")<<"signl "<<ii<<" "<<signl[npt];
              ++npt;
            }
            // decide if this RAT should be studied
            if(fStudyHits) StudyHits(1, npt, ticks, 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, ticks, 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 589 of file CCHitFinderAlg.cxx.

    {
    // 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,   // summedADC
        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
        );
/*
  if(prt) {
    mf::LogVerbatim("CCHitFinder")<<"W:T "<<allhits.back().WireID().Wire
      <<":"<<(short)allhits.back().PeakTime()
      <<" Chg "<<(short)allhits.back().Integral()
      <<" RMS "<<allhits.back().RMS()
      <<" lo ID "<<allhits.back().LocalIndex()
      <<" numHits "<<allhits.back().Multiplicity()
      <<" loTime "<<allhits.back().StartTick()<<" hiTime "<<allhits.back().EndTick()
      <<" chidof "<<allhits.back().GoodnessOfFit()
      << " DOF " << allhits.back().DegreesOfFreedom();
  }
*/
    } // hit
  } // StoreHits

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

private

Definition at line 660 of file CCHitFinderAlg.cxx.

                                                         {
    // 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 92 of file CCHitFinderAlg.h.

{ return std::move(allhits); }

Member Data Documentation

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

Definition at line 82 of file CCHitFinderAlg.h.

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

private

Definition at line 165 of file CCHitFinderAlg.h.

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

private

Definition at line 163 of file CCHitFinderAlg.h.

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

private

Definition at line 166 of file CCHitFinderAlg.h.

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

private

Definition at line 136 of file CCHitFinderAlg.h.

float hit::CCHitFinderAlg::chidof

private

Definition at line 134 of file CCHitFinderAlg.h.

float hit::CCHitFinderAlg::chinorm

private

Definition at line 115 of file CCHitFinderAlg.h.

int hit::CCHitFinderAlg::dof

private

Definition at line 135 of file CCHitFinderAlg.h.

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

private

Definition at line 109 of file CCHitFinderAlg.h.

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

private

Definition at line 107 of file CCHitFinderAlg.h.

float hit::CCHitFinderAlg::fChiSplit

private

Estimated noise error on the Signal.

Definition at line 104 of file CCHitFinderAlg.h.

FitStats_t hit::CCHitFinderAlg::FinalFitStats

private

counts of the good fits

Definition at line 193 of file CCHitFinderAlg.h.

std::unique_ptr<GausFitCache> hit::CCHitFinderAlg::FitCache

private

a set of functions ready to be used

Definition at line 178 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::fMaxBumps

private

Definition at line 102 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::fMaxXtraHits

private

Definition at line 103 of file CCHitFinderAlg.h.

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

private

Definition at line 100 of file CCHitFinderAlg.h.

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

private

Definition at line 101 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::fStudyHits

private

Definition at line 157 of file CCHitFinderAlg.h.

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

private

Definition at line 108 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::fUseChannelFilter

private

Definition at line 120 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::fUseFastFit

private

whether to attempt using a fast fit on single gauss.

Definition at line 176 of file CCHitFinderAlg.h.

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

private

Definition at line 158 of file CCHitFinderAlg.h.

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

private

Definition at line 158 of file CCHitFinderAlg.h.

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

private

Definition at line 159 of file CCHitFinderAlg.h.

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

private

Definition at line 159 of file CCHitFinderAlg.h.

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

private

Definition at line 160 of file CCHitFinderAlg.h.

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

private

Definition at line 160 of file CCHitFinderAlg.h.

art::ServiceHandle<geo::Geometry const> hit::CCHitFinderAlg::geom

private

Definition at line 124 of file CCHitFinderAlg.h.

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

private

Definition at line 167 of file CCHitFinderAlg.h.

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

private

Definition at line 173 of file CCHitFinderAlg.h.

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

private

Definition at line 168 of file CCHitFinderAlg.h.

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

private

Definition at line 172 of file CCHitFinderAlg.h.

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

private

Definition at line 171 of file CCHitFinderAlg.h.

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

private

Definition at line 170 of file CCHitFinderAlg.h.

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

staticprivate

Definition at line 220 of file CCHitFinderAlg.h.

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

private

Definition at line 130 of file CCHitFinderAlg.h.

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

private

Definition at line 131 of file CCHitFinderAlg.h.

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

private

Definition at line 133 of file CCHitFinderAlg.h.

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

private

Definition at line 132 of file CCHitFinderAlg.h.

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

private

Definition at line 164 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::SelRAT

private

Definition at line 174 of file CCHitFinderAlg.h.

constexpr float hit::CCHitFinderAlg::Sqrt2Pi = 2.5066

staticprivate

Definition at line 117 of file CCHitFinderAlg.h.

constexpr float hit::CCHitFinderAlg::SqrtPi = 1.7725

staticprivate

Definition at line 118 of file CCHitFinderAlg.h.

raw::ChannelID_t hit::CCHitFinderAlg::theChannel

private

Definition at line 111 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::thePlane

private

Definition at line 113 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::theWireNum

private

Definition at line 112 of file CCHitFinderAlg.h.

FitStats_t hit::CCHitFinderAlg::TriedFitStats

private

counts of the tried fits

Definition at line 194 of file CCHitFinderAlg.h.

---

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

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