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 79 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 232 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 93 of file CCHitFinderAlg.h.

93 { return std::move(allhits); }

hit::CCHitFinderAlg::allhits

std::vector< recob::Hit > allhits

Definition: CCHitFinderAlg.h:83

Member Data Documentation

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

Definition at line 83 of file CCHitFinderAlg.h.

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

private

Definition at line 166 of file CCHitFinderAlg.h.

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

private

Definition at line 164 of file CCHitFinderAlg.h.

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

private

Definition at line 167 of file CCHitFinderAlg.h.

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

private

Definition at line 137 of file CCHitFinderAlg.h.

float hit::CCHitFinderAlg::chidof

private

Definition at line 135 of file CCHitFinderAlg.h.

float hit::CCHitFinderAlg::chinorm

private

Definition at line 116 of file CCHitFinderAlg.h.

int hit::CCHitFinderAlg::dof

private

Definition at line 136 of file CCHitFinderAlg.h.

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

private

Definition at line 110 of file CCHitFinderAlg.h.

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

private

Definition at line 108 of file CCHitFinderAlg.h.

float hit::CCHitFinderAlg::fChiSplit

private

Estimated noise error on the Signal.

Definition at line 105 of file CCHitFinderAlg.h.

FitStats_t hit::CCHitFinderAlg::FinalFitStats

private

counts of the good fits

Definition at line 194 of file CCHitFinderAlg.h.

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

private

a set of functions ready to be used

Definition at line 179 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::fMaxBumps

private

Definition at line 103 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::fMaxXtraHits

private

Definition at line 104 of file CCHitFinderAlg.h.

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

private

Definition at line 101 of file CCHitFinderAlg.h.

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

private

Definition at line 102 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::fStudyHits

private

Definition at line 158 of file CCHitFinderAlg.h.

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

private

Definition at line 109 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::fUseChannelFilter

private

Definition at line 121 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::fUseFastFit

private

whether to attempt using a fast fit on single gauss.

Definition at line 177 of file CCHitFinderAlg.h.

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

private

Definition at line 159 of file CCHitFinderAlg.h.

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

private

Definition at line 159 of file CCHitFinderAlg.h.

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

private

Definition at line 160 of file CCHitFinderAlg.h.

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

private

Definition at line 160 of file CCHitFinderAlg.h.

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

private

Definition at line 161 of file CCHitFinderAlg.h.

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

private

Definition at line 161 of file CCHitFinderAlg.h.

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

private

Definition at line 125 of file CCHitFinderAlg.h.

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

private

Definition at line 168 of file CCHitFinderAlg.h.

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

private

Definition at line 174 of file CCHitFinderAlg.h.

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

private

Definition at line 169 of file CCHitFinderAlg.h.

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

private

Definition at line 173 of file CCHitFinderAlg.h.

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

private

Definition at line 172 of file CCHitFinderAlg.h.

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

private

Definition at line 171 of file CCHitFinderAlg.h.

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

staticprivate

Definition at line 221 of file CCHitFinderAlg.h.

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

private

Definition at line 131 of file CCHitFinderAlg.h.

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

private

Definition at line 132 of file CCHitFinderAlg.h.

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

private

Definition at line 134 of file CCHitFinderAlg.h.

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

private

Definition at line 133 of file CCHitFinderAlg.h.

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

private

Definition at line 165 of file CCHitFinderAlg.h.

bool hit::CCHitFinderAlg::SelRAT

private

Definition at line 175 of file CCHitFinderAlg.h.

constexpr float hit::CCHitFinderAlg::Sqrt2Pi = 2.5066

staticprivate

Definition at line 118 of file CCHitFinderAlg.h.

constexpr float hit::CCHitFinderAlg::SqrtPi = 1.7725

staticprivate

Definition at line 119 of file CCHitFinderAlg.h.

raw::ChannelID_t hit::CCHitFinderAlg::theChannel

private

Definition at line 112 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::thePlane

private

Definition at line 114 of file CCHitFinderAlg.h.

unsigned short hit::CCHitFinderAlg::theWireNum

private

Definition at line 113 of file CCHitFinderAlg.h.

FitStats_t hit::CCHitFinderAlg::TriedFitStats

private

counts of the tried fits

Definition at line 195 of file CCHitFinderAlg.h.

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

Classes

Public Member Functions

Public Attributes

Private Member Functions

Static Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation