Public Member Functions | Private Member Functions | Private Attributes | List of all members
hit::GausHitFinder Class Reference
Inheritance diagram for hit::GausHitFinder:
art::SharedProducer art::detail::Producer art::detail::SharedModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 GausHitFinder (fhicl::ParameterSet const &pset, art::ProcessingFrame const &)
 
- Public Member Functions inherited from art::SharedProducer
 SharedProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 SharedProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 
- Public Member Functions inherited from art::detail::SharedModule
 SharedModule ()
 
 SharedModule (std::string const &moduleLabel)
 
hep::concurrency::SerialTaskQueueChain * serialTaskQueueChain () const
 
std::set< std::string > const & sharedResources () const
 
void createQueues (SharedResources const &resources)
 
template<BranchType , typename... T>
void serialize (T const &...resources)
 
template<BranchType , typename... T>
void serializeExternal (T const &...resources)
 

Private Member Functions

void produce (art::Event &evt, art::ProcessingFrame const &) override
 
void beginJob (art::ProcessingFrame const &) override
 
std::vector< double > FillOutHitParameterVector (const std::vector< double > &input)
 

Private Attributes

const bool fFilterHits
 
const bool fFillHists
 
const std::string fCalDataModuleLabel
 
const std::string fAllHitsInstanceName
 
const std::vector< int > fLongMaxHitsVec
 Maximum number hits on a really long pulse train. More...
 
const std::vector< int > fLongPulseWidthVec
 Sets width of hits used to describe long pulses. More...
 
const size_t fMaxMultiHit
 maximum hits for multi fit More...
 
const int fAreaMethod
 Type of area calculation. More...
 
const std::vector< double > fAreaNormsVec
 factors for converting area to same units as peak height More...
 
const double fChi2NDF
 
const std::vector< float > fPulseHeightCuts
 maximum Chisquared / NDF allowed for a hit to be saved More...
 
const std::vector< float > fPulseWidthCuts
 
const std::vector< float > fPulseRatioCuts
 
std::atomic< size_t > fEventCount {0}
 
std::vector< std::unique_ptr< reco_tool::ICandidateHitFinder > > fHitFinderToolVec
 For finding candidate hits. More...
 
std::unique_ptr< reco_tool::IPeakFitterfPeakFitterTool
 Perform fit to candidate peaks. More...
 
std::unique_ptr< HitFilterAlgfHitFilterAlg
 algorithm used to filter out noise hits More...
 
TH1F * fFirstChi2
 
TH1F * fChi2
 

Additional Inherited Members

- Public Types inherited from art::SharedProducer
using ModuleType = SharedProducer
 
using WorkerType = WorkerT< SharedProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
- Protected Member Functions inherited from art::detail::SharedModule
template<BranchType BT = InEvent, typename... T>
void serialize (T const &...)
 
template<BranchType BT = InEvent, typename... T>
void serializeExternal (T const &...)
 
template<BranchType BT = InEvent>
void async ()
 

Detailed Description

Definition at line 61 of file GausHitFinder_module.cc.

Constructor & Destructor Documentation

hit::GausHitFinder::GausHitFinder ( fhicl::ParameterSet const &  pset,
art::ProcessingFrame const &   
)
explicit

Definition at line 108 of file GausHitFinder_module.cc.

109  : SharedProducer{pset}
110  , fFilterHits(pset.get<bool>("FilterHits", false))
111  , fFillHists(pset.get<bool>("FillHists", false))
112  , fCalDataModuleLabel(pset.get<std::string>("CalDataModuleLabel"))
113  , fAllHitsInstanceName(pset.get<std::string>("AllHitsInstanceName", ""))
114  , fLongMaxHitsVec(pset.get<std::vector<int>>("LongMaxHits", std::vector<int>() = {25, 25, 25}))
116  pset.get<std::vector<int>>("LongPulseWidth", std::vector<int>() = {16, 16, 16}))
117  , fMaxMultiHit(pset.get<int>("MaxMultiHit"))
118  , fAreaMethod(pset.get<int>("AreaMethod"))
119  , fAreaNormsVec(FillOutHitParameterVector(pset.get<std::vector<double>>("AreaNorms")))
120  , fChi2NDF(pset.get<double>("Chi2NDF"))
122  pset.get<std::vector<float>>("PulseHeightCuts", std::vector<float>() = {3.0, 3.0, 3.0}))
123  , fPulseWidthCuts(
124  pset.get<std::vector<float>>("PulseWidthCuts", std::vector<float>() = {2.0, 1.5, 1.0}))
125  , fPulseRatioCuts(
126  pset.get<std::vector<float>>("PulseRatioCuts", std::vector<float>() = {0.35, 0.40, 0.20}))
127  {
128  if (fFillHists && art::Globals::instance()->nthreads() > 1u) {
130  << "Cannot fill histograms when multiple threads configured, please set fFillHists to "
131  "false or change number of threads to 1\n";
132  }
133  async<art::InEvent>();
134  if (fFilterHits) {
135  fHitFilterAlg = std::make_unique<HitFilterAlg>(pset.get<fhicl::ParameterSet>("HitFilterAlg"));
136  }
137 
138  // recover the tool to do the candidate hit finding
139  // Recover the vector of fhicl parameters for the ROI tools
140  const fhicl::ParameterSet& hitFinderTools = pset.get<fhicl::ParameterSet>("HitFinderToolVec");
141 
142  fHitFinderToolVec.resize(hitFinderTools.get_pset_names().size());
143 
144  for (const std::string& hitFinderTool : hitFinderTools.get_pset_names()) {
145  const fhicl::ParameterSet& hitFinderToolParamSet =
146  hitFinderTools.get<fhicl::ParameterSet>(hitFinderTool);
147  size_t planeIdx = hitFinderToolParamSet.get<size_t>("Plane");
148 
149  fHitFinderToolVec.at(planeIdx) =
150  art::make_tool<reco_tool::ICandidateHitFinder>(hitFinderToolParamSet);
151  }
152 
153  // Recover the peak fitting tool
155  art::make_tool<reco_tool::IPeakFitter>(pset.get<fhicl::ParameterSet>("PeakFitter"));
156 
157  // let HitCollectionCreator declare that we are going to produce
158  // hits and associations with wires and raw digits
159  // We want the option to output two hit collections, one filtered
160  // and one with all hits. The key to doing this will be a non-null
161  // instance name for the second collection
162  // (with no particular product label)
164  producesCollector(), fAllHitsInstanceName, true, false); //fMakeRawDigitAssns);
165 
166  // and now the filtered hits...
167  if (fAllHitsInstanceName != "")
169  producesCollector(), "", true, false); //fMakeRawDigitAssns);
170 
171  return;
172  } // GausHitFinder::GausHitFinder()
std::string string
Definition: nybbler.cc:12
std::unique_ptr< HitFilterAlg > fHitFilterAlg
algorithm used to filter out noise hits
const std::vector< int > fLongPulseWidthVec
Sets width of hits used to describe long pulses.
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:248
std::unique_ptr< reco_tool::IPeakFitter > fPeakFitterTool
Perform fit to candidate peaks.
const int fAreaMethod
Type of area calculation.
const std::string fCalDataModuleLabel
const std::vector< float > fPulseWidthCuts
std::vector< std::string > get_pset_names() const
SharedProducer(fhicl::ParameterSet const &pset)
std::vector< double > FillOutHitParameterVector(const std::vector< double > &input)
T get(std::string const &key) const
Definition: ParameterSet.h:271
const size_t fMaxMultiHit
maximum hits for multi fit
const std::vector< int > fLongMaxHitsVec
Maximum number hits on a really long pulse train.
ProducesCollector & producesCollector() noexcept
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
const std::vector< float > fPulseRatioCuts
const std::vector< float > fPulseHeightCuts
maximum Chisquared / NDF allowed for a hit to be saved
const std::string fAllHitsInstanceName
const std::vector< double > fAreaNormsVec
factors for converting area to same units as peak height
static Globals * instance()
Definition: Globals.cc:17
std::vector< std::unique_ptr< reco_tool::ICandidateHitFinder > > fHitFinderToolVec
For finding candidate hits.

Member Function Documentation

void hit::GausHitFinder::beginJob ( art::ProcessingFrame const &  )
overrideprivatevirtual

Reimplemented from art::SharedProducer.

Definition at line 200 of file GausHitFinder_module.cc.

201  {
202  // get access to the TFile service
204 
205  // ======================================
206  // === Hit Information for Histograms ===
207  if (fFillHists) {
208  fFirstChi2 = tfs->make<TH1F>("fFirstChi2", "#chi^{2}", 10000, 0, 5000);
209  fChi2 = tfs->make<TH1F>("fChi2", "#chi^{2}", 10000, 0, 5000);
210  }
211  }
std::vector< double > hit::GausHitFinder::FillOutHitParameterVector ( const std::vector< double > &  input)
private

Definition at line 177 of file GausHitFinder_module.cc.

178  {
179  if (input.size() == 0)
180  throw std::runtime_error(
181  "GausHitFinder::FillOutHitParameterVector ERROR! Input config vector has zero size.");
182 
183  std::vector<double> output;
185  const unsigned int N_PLANES = geom->Nplanes();
186 
187  if (input.size() == 1)
188  output.resize(N_PLANES, input[0]);
189  else if (input.size() == N_PLANES)
190  output = input;
191  else
192  throw std::runtime_error("GausHitFinder::FillOutHitParameterVector ERROR! Input config "
193  "vector size !=1 and !=N_PLANES.");
194  return output;
195  }
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
void hit::GausHitFinder::produce ( art::Event evt,
art::ProcessingFrame const &   
)
overrideprivatevirtual

Implements art::SharedProducer.

Definition at line 218 of file GausHitFinder_module.cc.

219  {
220  unsigned int count = fEventCount.fetch_add(1);
221  //==================================================================================================
222 
223  TH1::AddDirectory(kFALSE);
224 
225  // Instantiate and Reset a stop watch
226  //TStopwatch StopWatch;
227  //StopWatch.Reset();
228 
229  // ################################
230  // ### Calling Geometry service ###
231  // ################################
233 
234  // ###############################################
235  // ### Making a ptr vector to put on the event ###
236  // ###############################################
237  // this contains the hit collection
238  // and its associations to wires and raw digits
239  recob::HitCollectionCreator allHitCol(evt, fAllHitsInstanceName, true, false);
240 
241  // Handle the filtered hits collection...
242  recob::HitCollectionCreator hcol(evt, "", true, false);
243  recob::HitCollectionCreator* filteredHitCol = 0;
244 
245  if (fFilterHits) filteredHitCol = &hcol;
246 
247  //store in a thread safe way
248  struct hitstruct {
249  recob::Hit hit_tbb;
250  art::Ptr<recob::Wire> wire_tbb;
251  };
252 
253  tbb::concurrent_vector<hitstruct> hitstruct_vec;
254  tbb::concurrent_vector<hitstruct> filthitstruct_vec;
255 
256  // if (fAllHitsInstanceName != "") filteredHitCol = &hcol;
257 
258  // ##########################################
259  // ### Reading in the Wire List object(s) ###
260  // ##########################################
262  evt.getByLabel(fCalDataModuleLabel, wireVecHandle);
263 
264  //#################################################
265  //### Set the charge determination method ###
266  //### Default is to compute the normalized area ###
267  //#################################################
268  std::function<double(double, double, double, double, int, int)> chargeFunc =
269  [](double peakMean, double peakAmp, double peakWidth, double areaNorm, int low, int hi) {
270  return std::sqrt(2 * TMath::Pi()) * peakAmp * peakWidth / areaNorm;
271  };
272 
273  //##############################################
274  //### Alternative is to integrate over pulse ###
275  //##############################################
276  if (fAreaMethod == 0)
277  chargeFunc =
278  [](double peakMean, double peakAmp, double peakWidth, double areaNorm, int low, int hi) {
279  double charge(0);
280  for (int sigPos = low; sigPos < hi; sigPos++)
281  charge += peakAmp * TMath::Gaus(sigPos, peakMean, peakWidth);
282  return charge;
283  };
284 
285  //##############################
286  //### Looping over the wires ###
287  //##############################
288  //for(size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++)
289  //{
290  tbb::parallel_for(
291  static_cast<std::size_t>(0),
292  wireVecHandle->size(),
293  [&](size_t& wireIter) {
294  // ####################################
295  // ### Getting this particular wire ###
296  // ####################################
297  art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
298 
299  // --- Setting Channel Number and Signal type ---
300 
301  raw::ChannelID_t channel = wire->Channel();
302 
303  // get the WireID for this hit
304  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
305  // for now, just take the first option returned from ChannelToWire
306  geo::WireID wid = wids[0];
307  // We need to know the plane to look up parameters
308  geo::PlaneID::PlaneID_t plane = wid.Plane;
309 
310  // ----------------------------------------------------------
311  // -- Setting the appropriate signal widths and thresholds --
312  // -- for the right plane. --
313  // ----------------------------------------------------------
314 
315  // #################################################
316  // ### Set up to loop over ROI's for this wire ###
317  // #################################################
318  const recob::Wire::RegionsOfInterest_t& signalROI = wire->SignalROI();
319 
320  for (const auto& range : signalROI.get_ranges()) {
321  // ROI start time
322  raw::TDCtick_t roiFirstBinTick = range.begin_index();
323 
324  // ###########################################################
325  // ### Scan the waveform and find candidate peaks + merge ###
326  // ###########################################################
327 
330 
331  fHitFinderToolVec.at(plane)->findHitCandidates(range, 0, channel, count, hitCandidateVec);
332  fHitFinderToolVec.at(plane)->MergeHitCandidates(
333  range, hitCandidateVec, mergedCandidateHitVec);
334 
335  // #######################################################
336  // ### Lets loop over the pulses we found on this wire ###
337  // #######################################################
338 
339  for (auto& mergedCands : mergedCandidateHitVec) {
340  int startT = mergedCands.front().startTick;
341  int endT = mergedCands.back().stopTick;
342 
343  // ### Putting in a protection in case things went wrong ###
344  // ### In the end, this primarily catches the case where ###
345  // ### a fake pulse is at the start of the ROI ###
346  if (endT - startT < 5) continue;
347 
348  // #######################################################
349  // ### Clearing the parameter vector for the new pulse ###
350  // #######################################################
351 
352  // === Setting the number of Gaussians to try ===
353  int nGausForFit = mergedCands.size();
354 
355  // ##################################################
356  // ### Calling the function for fitting Gaussians ###
357  // ##################################################
358  double chi2PerNDF(0.);
359  int NDF(1);
360  /*stand alone
361  reco_tool::IPeakFitter::PeakParamsVec peakParamsVec(nGausForFit);
362  */
364 
365  // #######################################################
366  // ### If # requested Gaussians is too large then punt ###
367  // #######################################################
368  if (mergedCands.size() <= fMaxMultiHit) {
369  fPeakFitterTool->findPeakParameters(
370  range.data(), mergedCands, peakParamsVec, chi2PerNDF, NDF);
371 
372  // If the chi2 is infinite then there is a real problem so we bail
373  if (!(chi2PerNDF < std::numeric_limits<double>::infinity())) {
374  chi2PerNDF = 2. * fChi2NDF;
375  NDF = 2;
376  }
377 
378  if (fFillHists) fFirstChi2->Fill(chi2PerNDF);
379  }
380 
381  // #######################################################
382  // ### If too large then force alternate solution ###
383  // ### - Make n hits from pulse train where n will ###
384  // ### depend on the fhicl parameter fLongPulseWidth ###
385  // ### Also do this if chi^2 is too large ###
386  // #######################################################
387  if (mergedCands.size() > fMaxMultiHit || nGausForFit * chi2PerNDF > fChi2NDF) {
388  int longPulseWidth = fLongPulseWidthVec.at(plane);
389  int nHitsThisPulse = (endT - startT) / longPulseWidth;
390 
391  if (nHitsThisPulse > fLongMaxHitsVec.at(plane)) {
392  nHitsThisPulse = fLongMaxHitsVec.at(plane);
393  longPulseWidth = (endT - startT) / nHitsThisPulse;
394  }
395 
396  if (nHitsThisPulse * longPulseWidth < endT - startT) nHitsThisPulse++;
397 
398  int firstTick = startT;
399  int lastTick = std::min(firstTick + longPulseWidth, endT);
400 
401  peakParamsVec.clear();
402  nGausForFit = nHitsThisPulse;
403  NDF = 1.;
404  chi2PerNDF = chi2PerNDF > fChi2NDF ? chi2PerNDF : -1.;
405 
406  for (int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) {
407  // This hit parameters
408  double sumADC =
409  std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.);
410  double peakSigma = (lastTick - firstTick) / 3.; // Set the width...
411  double peakAmp = 0.3989 * sumADC / peakSigma; // Use gaussian formulation
412  double peakMean = (firstTick + lastTick) / 2.;
413 
414  // Store hit params
416 
417  peakParams.peakCenter = peakMean;
418  peakParams.peakCenterError = 0.1 * peakMean;
419  peakParams.peakSigma = peakSigma;
420  peakParams.peakSigmaError = 0.1 * peakSigma;
421  peakParams.peakAmplitude = peakAmp;
422  peakParams.peakAmplitudeError = 0.1 * peakAmp;
423 
424  peakParamsVec.push_back(peakParams);
425 
426  // set for next loop
427  firstTick = lastTick;
428  lastTick = std::min(lastTick + longPulseWidth, endT);
429  }
430  }
431 
432  // #######################################################
433  // ### Loop through returned peaks and make recob hits ###
434  // #######################################################
435 
436  int numHits(0);
437 
438  // Make a container for what will be the filtered collection
439  std::vector<recob::Hit> filteredHitVec;
440 
441  for (const auto& peakParams : peakParamsVec) {
442  // Extract values for this hit
443  float peakAmp = peakParams.peakAmplitude;
444  float peakMean = peakParams.peakCenter;
445  float peakWidth = peakParams.peakSigma;
446 
447  // Place one bit of protection here
448  if (std::isnan(peakAmp)) {
449  std::cout << "**** hit peak amplitude is a nan! Channel: " << channel
450  << ", start tick: " << startT << std::endl;
451  continue;
452  }
453 
454  // Extract errors
455  float peakAmpErr = peakParams.peakAmplitudeError;
456  float peakMeanErr = peakParams.peakCenterError;
457  float peakWidthErr = peakParams.peakSigmaError;
458 
459  // ### Charge ###
460  float charge =
461  chargeFunc(peakMean, peakAmp, peakWidth, fAreaNormsVec[plane], startT, endT);
462  ;
463  float chargeErr =
464  std::sqrt(TMath::Pi()) * (peakAmpErr * peakWidthErr + peakWidthErr * peakAmpErr);
465 
466  // ### limits for getting sums
467  std::vector<float>::const_iterator sumStartItr = range.begin() + startT;
468  std::vector<float>::const_iterator sumEndItr = range.begin() + endT;
469 
470  // ### Sum of ADC counts
471  double sumADC = std::accumulate(sumStartItr, sumEndItr, 0.);
472 
473  // ok, now create the hit
474  recob::HitCreator hitcreator(
475  *wire, // wire reference
476  wid, // wire ID
477  startT + roiFirstBinTick, // start_tick TODO check
478  endT + roiFirstBinTick, // end_tick TODO check
479  peakWidth, // rms
480  peakMean + roiFirstBinTick, // peak_time
481  peakMeanErr, // sigma_peak_time
482  peakAmp, // peak_amplitude
483  peakAmpErr, // sigma_peak_amplitude
484  charge, // hit_integral
485  chargeErr, // hit_sigma_integral
486  sumADC, // summedADC FIXME
487  nGausForFit, // multiplicity
488  numHits, // local_index TODO check that the order is correct
489  chi2PerNDF, // goodness_of_fit
490  NDF // dof
491  );
492 
493  if (filteredHitCol) filteredHitVec.push_back(hitcreator.copy());
494 
495  const recob::Hit hit(hitcreator.move());
496 
497  // This loop will store ALL hits
498  hitstruct tmp{std::move(hit), wire};
499  hitstruct_vec.push_back(std::move(tmp));
500 
501  numHits++;
502  } // <---End loop over gaussians
503 
504  // Should we filter hits?
505  if (filteredHitCol && !filteredHitVec.empty()) {
506  // #######################################################################
507  // Is all this sorting really necessary? Would it be faster to just loop
508  // through the hits and perform simple cuts on amplitude and width on a
509  // hit-by-hit basis, either here in the module (using fPulseHeightCuts and
510  // fPulseWidthCuts) or in HitFilterAlg?
511  // #######################################################################
512 
513  // Sort in ascending peak height
514  std::sort(filteredHitVec.begin(),
515  filteredHitVec.end(),
516  [](const auto& left, const auto& right) {
517  return left.PeakAmplitude() > right.PeakAmplitude();
518  });
519 
520  // Reject if the first hit fails the PH/wid cuts
521  if (filteredHitVec.front().PeakAmplitude() < fPulseHeightCuts.at(plane) ||
522  filteredHitVec.front().RMS() < fPulseWidthCuts.at(plane))
523  filteredHitVec.clear();
524 
525  // Now check other hits in the snippet
526  if (filteredHitVec.size() > 1) {
527  // The largest pulse height will now be at the front...
528  float largestPH = filteredHitVec.front().PeakAmplitude();
529 
530  // Find where the pulse heights drop below threshold
531  float threshold(fPulseRatioCuts.at(plane));
532 
533  std::vector<recob::Hit>::iterator smallHitItr = std::find_if(
534  filteredHitVec.begin(),
535  filteredHitVec.end(),
536  [largestPH, threshold](const auto& hit) {
537  return hit.PeakAmplitude() < 8. && hit.PeakAmplitude() / largestPH < threshold;
538  });
539 
540  // Shrink to fit
541  if (smallHitItr != filteredHitVec.end())
542  filteredHitVec.resize(std::distance(filteredHitVec.begin(), smallHitItr));
543 
544  // Resort in time order
545  std::sort(filteredHitVec.begin(),
546  filteredHitVec.end(),
547  [](const auto& left, const auto& right) {
548  return left.PeakTime() < right.PeakTime();
549  });
550  }
551 
552  // Copy the hits we want to keep to the filtered hit collection
553  for (const auto& filteredHit : filteredHitVec)
554  if (!fHitFilterAlg || fHitFilterAlg->IsGoodHit(filteredHit)) {
555  hitstruct tmp{std::move(filteredHit), wire};
556  filthitstruct_vec.push_back(std::move(tmp));
557  }
558 
559  if (fFillHists) fChi2->Fill(chi2PerNDF);
560  }
561  } //<---End loop over merged candidate hits
562 
563  } //<---End looping over ROI's
564  } //<---End looping over all the wires
565  ); //end tbb parallel for
566 
567  for (size_t i = 0; i < hitstruct_vec.size(); i++) {
568  allHitCol.emplace_back(hitstruct_vec[i].hit_tbb, hitstruct_vec[i].wire_tbb);
569  }
570 
571  for (size_t j = 0; j < filthitstruct_vec.size(); j++) {
572  filteredHitCol->emplace_back(filthitstruct_vec[j].hit_tbb, filthitstruct_vec[j].wire_tbb);
573  }
574 
575  //==================================================================================================
576  // End of the event -- move the hit collection and the associations into the event
577 
578  if (filteredHitCol) {
579 
580  // If we filtered hits but no instance name was
581  // specified for the "all hits" collection, then
582  // only save the filtered hits to the event
583  if (fAllHitsInstanceName == "") {
584  filteredHitCol->put_into(evt);
585 
586  // otherwise, save both
587  }
588  else {
589  filteredHitCol->put_into(evt);
590  allHitCol.put_into(evt);
591  }
592  }
593  else {
594  allHitCol.put_into(evt);
595  }
596 
597  // Keep track of events processed
598  //fEventCount++;
599 
600  } // End of produce()
intermediate_table::iterator iterator
std::unique_ptr< HitFilterAlg > fHitFilterAlg
algorithm used to filter out noise hits
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
intermediate_table::const_iterator const_iterator
uint8_t channel
Definition: CRTFragment.hh:201
const std::vector< int > fLongPulseWidthVec
Sets width of hits used to describe long pulses.
std::unique_ptr< reco_tool::IPeakFitter > fPeakFitterTool
Perform fit to candidate peaks.
const int fAreaMethod
Type of area calculation.
const std::string fCalDataModuleLabel
const std::vector< float > fPulseWidthCuts
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:83
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
A class handling a collection of hits and its associations.
Definition: HitCreator.h:508
std::atomic< size_t > fEventCount
def move(depos, offset)
Definition: depos.py:107
const size_t fMaxMultiHit
maximum hits for multi fit
const std::vector< int > fLongMaxHitsVec
Maximum number hits on a really long pulse train.
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
Definition: HitCreator.cxx:290
string tmp
Definition: languages.py:63
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:652
Detector simulation of raw signals on wires.
std::vector< PeakFitParams_t > PeakParamsVec
Definition: IPeakFitter.h:37
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
const std::vector< float > fPulseRatioCuts
const std::vector< float > fPulseHeightCuts
maximum Chisquared / NDF allowed for a hit to be saved
const std::string fAllHitsInstanceName
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
const std::vector< double > fAreaNormsVec
factors for converting area to same units as peak height
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
std::vector< HitCandidateVec > MergeHitCandidateVec
QTextStream & endl(QTextStream &s)
std::vector< std::unique_ptr< reco_tool::ICandidateHitFinder > > fHitFinderToolVec
For finding candidate hits.
std::vector< HitCandidate > HitCandidateVec

Member Data Documentation

const std::string hit::GausHitFinder::fAllHitsInstanceName
private

Definition at line 75 of file GausHitFinder_module.cc.

const int hit::GausHitFinder::fAreaMethod
private

Type of area calculation.

Definition at line 81 of file GausHitFinder_module.cc.

const std::vector<double> hit::GausHitFinder::fAreaNormsVec
private

factors for converting area to same units as peak height

Definition at line 83 of file GausHitFinder_module.cc.

const std::string hit::GausHitFinder::fCalDataModuleLabel
private

Definition at line 74 of file GausHitFinder_module.cc.

TH1F* hit::GausHitFinder::fChi2
private

Definition at line 102 of file GausHitFinder_module.cc.

const double hit::GausHitFinder::fChi2NDF
private

Definition at line 84 of file GausHitFinder_module.cc.

std::atomic<size_t> hit::GausHitFinder::fEventCount {0}
private

Definition at line 90 of file GausHitFinder_module.cc.

const bool hit::GausHitFinder::fFillHists
private

Definition at line 72 of file GausHitFinder_module.cc.

const bool hit::GausHitFinder::fFilterHits
private

Definition at line 71 of file GausHitFinder_module.cc.

TH1F* hit::GausHitFinder::fFirstChi2
private

Definition at line 101 of file GausHitFinder_module.cc.

std::unique_ptr<HitFilterAlg> hit::GausHitFinder::fHitFilterAlg
private

algorithm used to filter out noise hits

Definition at line 98 of file GausHitFinder_module.cc.

std::vector<std::unique_ptr<reco_tool::ICandidateHitFinder> > hit::GausHitFinder::fHitFinderToolVec
private

For finding candidate hits.

Definition at line 94 of file GausHitFinder_module.cc.

const std::vector<int> hit::GausHitFinder::fLongMaxHitsVec
private

Maximum number hits on a really long pulse train.

Definition at line 77 of file GausHitFinder_module.cc.

const std::vector<int> hit::GausHitFinder::fLongPulseWidthVec
private

Sets width of hits used to describe long pulses.

Definition at line 78 of file GausHitFinder_module.cc.

const size_t hit::GausHitFinder::fMaxMultiHit
private

maximum hits for multi fit

Definition at line 80 of file GausHitFinder_module.cc.

std::unique_ptr<reco_tool::IPeakFitter> hit::GausHitFinder::fPeakFitterTool
private

Perform fit to candidate peaks.

Definition at line 96 of file GausHitFinder_module.cc.

const std::vector<float> hit::GausHitFinder::fPulseHeightCuts
private

maximum Chisquared / NDF allowed for a hit to be saved

Definition at line 86 of file GausHitFinder_module.cc.

const std::vector<float> hit::GausHitFinder::fPulseRatioCuts
private

Definition at line 88 of file GausHitFinder_module.cc.

const std::vector<float> hit::GausHitFinder::fPulseWidthCuts
private

Definition at line 87 of file GausHitFinder_module.cc.


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