LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CandHitDerivative_tool.cc
Go to the documentation of this file.
1 //note for MT: this implementation is not thread-safe
6 
10 
14 #include "art_root_io/TFileService.h"
15 #include "cetlib_except/exception.h"
17 
18 #include "TProfile.h"
19 
20 #include <cmath>
21 
22 namespace reco_tool {
23 
25  public:
26  explicit CandHitDerivative(const fhicl::ParameterSet& pset);
27 
28  void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&,
29  const size_t,
30  const size_t,
31  const size_t,
32  HitCandidateVec&) const override;
33 
34  void MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&,
35  const HitCandidateVec&,
36  MergeHitCandidateVec&) const override;
37 
38  private:
39  // Internal functions
42  const size_t,
43  int,
44  float,
45  HitCandidateVec&) const;
46 
47  // Finding the nearest maximum/minimum from current point
52 
53  // handle finding the "start" and "stop" of a candidate hit
57 
58  // some control variables
59  size_t fPlane; //< Identifies the plane this tool is meant to operate on
60  int fMinDeltaTicks; //< minimum ticks from max to min to consider
61  int fMaxDeltaTicks; //< maximum ticks from max to min to consider
62  float fMinDeltaPeaks; //< minimum maximum to minimum peak distance
63  float fMinHitHeight; //< Drop candidate hits with height less than this
64  size_t fNumInterveningTicks; //< Number ticks between candidate hits to merge
65  bool fOutputHistograms; //< If true will generate a very large file of hists!
66 
67  art::TFileDirectory* fHistDirectory;
68 
69  // Global histograms
73 
74  mutable std::map<size_t, int> fChannelCntMap;
75 
76  // Member variables from the fhicl file
77  std::unique_ptr<reco_tool::IWaveformTool> fWaveformTool;
78 
79  const geo::GeometryCore* fGeometry = lar::providerFrom<geo::Geometry>();
80  };
81 
82  //----------------------------------------------------------------------
83  // Constructor.
85  {
86  fPlane = pset.get<size_t>("Plane", 0);
87  fMinDeltaTicks = pset.get<int>("MinDeltaTicks", 0);
88  fMaxDeltaTicks = pset.get<int>("MaxDeltaTicks", 30);
89  fMinDeltaPeaks = pset.get<float>("MinDeltaPeaks", 0.025);
90  fMinHitHeight = pset.get<float>("MinHitHeight", 2.0);
91  fNumInterveningTicks = pset.get<size_t>("NumInterveningTicks", 6);
92  fOutputHistograms = pset.get<bool>("OutputHistograms", false);
93 
94  // Recover the baseline tool
96  art::make_tool<reco_tool::IWaveformTool>(pset.get<fhicl::ParameterSet>("WaveformAlgs"));
97 
98  // If asked, define the global histograms
99  if (fOutputHistograms) {
100  // Access ART's TFileService, which will handle creating and writing
101  // histograms and n-tuples for us.
103 
104  fHistDirectory = tfs.get();
105 
106  // Make a directory for these histograms
107  art::TFileDirectory dir = fHistDirectory->mkdir(Form("HitPlane_%1zu", fPlane));
108 
110  dir.make<TH1F>(Form("DStopStart_%1zu", fPlane), ";Delta Stop/Start;", 200, 0., 200.);
112  dir.make<TH1F>(Form("DMaxTMinT_%1zu", fPlane), ";Delta Max/Min Tick;", 200, 0., 200.);
114  dir.make<TH1F>(Form("DMaxDMinD_%1zu", fPlane), ";Delta Max/Min Deriv;", 200, 0., 200.);
115  }
116 
117  return;
118  }
119 
121  const recob::Wire::RegionsOfInterest_t::datarange_t& dataRange,
122  const size_t roiStartTick,
123  const size_t channel,
124  const size_t eventCount,
125  HitCandidateVec& hitCandidateVec) const
126  {
127  // In this case we want to find hit candidates based on the derivative of of the input waveform
128  // We get this from our waveform algs too...
129  Waveform rawDerivativeVec;
130  Waveform derivativeVec;
131 
132  // Recover the actual waveform
133  const Waveform& waveform = dataRange.data();
134 
135  fWaveformTool->firstDerivative(waveform, rawDerivativeVec);
136  fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
137 
138  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
139  size_t plane = wids[0].Plane;
140  size_t cryo = wids[0].Cryostat;
141  size_t tpc = wids[0].TPC;
142  size_t wire = wids[0].Wire;
143 
144  // Just make sure the input candidate hit vector has been cleared
145  hitCandidateVec.clear();
146 
147  // Now find the hits
148  findHitCandidates(derivativeVec.begin(),
149  derivativeVec.end(),
150  roiStartTick,
153  hitCandidateVec);
154 
155  if (hitCandidateVec.empty()) {
156  if (plane == 0) {
157  std::cout << "** C/T/P: " << cryo << "/" << tpc << "/" << plane << ", wire: " << wire
158  << " has not hits with input size: " << waveform.size() << std::endl;
159  }
160  }
161 
162  // Reset the hit height from the input waveform
163  for (auto& hitCandidate : hitCandidateVec) {
164  size_t centerIdx = hitCandidate.hitCenter;
165 
166  hitCandidate.hitHeight = waveform.at(centerIdx);
167  }
168 
169  // Keep track of histograms if requested
170  if (fOutputHistograms) {
171  // Recover the details...
172  // std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
173  // size_t plane = wids[0].Plane;
174  // size_t cryo = wids[0].Cryostat;
175  // size_t tpc = wids[0].TPC;
176  // size_t wire = wids[0].Wire;
177 
178  size_t channelCnt = fChannelCntMap[channel]++;
179 
180  // Make a directory for these histograms
181  art::TFileDirectory dir = fHistDirectory->mkdir(
182  Form("HitPlane_%1zu/ev%04zu/c%1zut%1zuwire_%05zu", plane, eventCount, cryo, tpc, wire));
183 
184  size_t waveformSize = waveform.size();
185  int waveStart = roiStartTick;
186  int waveStop = waveStart + waveformSize;
187 
188  TProfile* waveHist = dir.make<TProfile>(
189  Form("HWfm_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
190  "Waveform",
191  waveformSize,
192  waveStart,
193  waveStop,
194  -500.,
195  500.);
196  TProfile* derivHist = dir.make<TProfile>(
197  Form("HDer_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
198  "Derivative",
199  waveformSize,
200  waveStart,
201  waveStop,
202  -500.,
203  500.);
204  TProfile* candHitHist = dir.make<TProfile>(
205  Form("HCan_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
206  "Cand Hits",
207  waveformSize,
208  waveStart,
209  waveStop,
210  -500.,
211  500.);
212  TProfile* maxDerivHist = dir.make<TProfile>(
213  Form("HMax_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
214  "Maxima",
215  waveformSize,
216  waveStart,
217  waveStop,
218  -500.,
219  500.);
220 
221  // Fill wave/derivative
222  for (size_t idx = 0; idx < waveform.size(); idx++) {
223  waveHist->Fill(roiStartTick + idx, waveform.at(idx));
224  derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
225  }
226 
227  // Fill hits
228  for (const auto& hitCandidate : hitCandidateVec) {
229  candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
230  maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
231  maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
232 
233  fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
234  fDMaxTickMinTickHist->Fill(hitCandidate.minTick - hitCandidate.maxTick, 1.);
235  fDMaxDerivMinDerivHist->Fill(hitCandidate.maxDerivative - hitCandidate.minDerivative, 1.);
236  }
237  }
238 
239  return;
240  }
241 
243  Waveform::const_iterator stopItr,
244  const size_t roiStartTick,
245  int dTicksThreshold,
246  float dPeakThreshold,
247  HitCandidateVec& hitCandidateVec) const
248  {
249  // Search for candidate hits...
250  // The idea will be to find the largest deviation in the input derivative waveform as the starting point. Depending
251  // on if a maximum or minimum, we search forward or backward to find the minimum or maximum that our extremum
252  // corresponds to.
253  std::pair<Waveform::const_iterator, Waveform::const_iterator> minMaxPair =
254  std::minmax_element(startItr, stopItr);
255 
256  Waveform::const_iterator maxItr = minMaxPair.second;
257  Waveform::const_iterator minItr = minMaxPair.first;
258 
259  // Use the larger of the two as the starting point and recover the nearest max or min
260  if (std::fabs(*maxItr) > std::fabs(*minItr))
261  minItr = findNearestMin(maxItr, stopItr);
262  else
263  maxItr = findNearestMax(minItr, startItr);
264 
265  int deltaTicks = std::distance(maxItr, minItr);
266  float range = *maxItr - *minItr;
267 
268  // At some point small rolling oscillations on the waveform need to be ignored...
269  if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) {
270  // Need to back up to find zero crossing, this will be the starting point of our
271  // candidate hit but also the endpoint of the pre sub-waveform we'll search next
272  Waveform::const_iterator newEndItr = findStartTick(maxItr, startItr);
273 
274  int startTick = std::distance(startItr, newEndItr);
275 
276  // Now need to go forward to again get close to zero, this will then be the end point
277  // of our candidate hit and the starting point for the post sub-waveform to search
278  Waveform::const_iterator newStartItr = findStopTick(minItr, stopItr);
279 
280  int stopTick = std::distance(startItr, newStartItr);
281 
282  // Find hits in the section of the waveform leading up to this candidate hit
283  if (startTick > dTicksThreshold) {
284  // Special handling for merged hits
285  if (*(newEndItr - 1) > 0.) {
286  dTicksThreshold = 2;
287  dPeakThreshold = 0.;
288  }
289  else {
290  dTicksThreshold = fMinDeltaTicks;
291  dPeakThreshold = fMinDeltaPeaks;
292  }
293 
295  startItr, newEndItr + 1, roiStartTick, dTicksThreshold, dPeakThreshold, hitCandidateVec);
296  }
297 
298  // Create a new hit candidate and store away
299  HitCandidate hitCandidate;
300 
301  Waveform::const_iterator peakItr =
302  std::min_element(maxItr, minItr, [](const auto& left, const auto& right) {
303  return std::fabs(left) < std::fabs(right);
304  });
305 
306  // Check balance
307  if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
308  peakItr--;
309  else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
310  peakItr++;
311 
312  hitCandidate.startTick = roiStartTick + startTick;
313  hitCandidate.stopTick = roiStartTick + stopTick;
314  hitCandidate.maxTick = roiStartTick + std::distance(startItr, maxItr);
315  hitCandidate.minTick = roiStartTick + std::distance(startItr, minItr);
316  hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
317  hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
318  hitCandidate.hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
319  hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
320  hitCandidate.hitHeight =
321  hitCandidate.hitSigma * (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
322 
323  hitCandidateVec.push_back(hitCandidate);
324 
325  // Finally, search the section of the waveform following this candidate for more hits
326  if (std::distance(newStartItr, stopItr) > dTicksThreshold) {
327  // Special handling for merged hits
328  if (*(newStartItr + 1) < 0.) {
329  dTicksThreshold = 2;
330  dPeakThreshold = 0.;
331  }
332  else {
333  dTicksThreshold = fMinDeltaTicks;
334  dPeakThreshold = fMinDeltaPeaks;
335  }
336 
337  findHitCandidates(newStartItr,
338  stopItr,
339  roiStartTick + stopTick,
340  dTicksThreshold,
341  dPeakThreshold,
342  hitCandidateVec);
343  }
344  }
345 
346  return;
347  }
348 
349  void CandHitDerivative::MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&,
350  const HitCandidateVec& hitCandidateVec,
351  MergeHitCandidateVec& mergedHitsVec) const
352  {
353  // If nothing on the input end then nothing to do
354  if (hitCandidateVec.empty()) return;
355 
356  // The idea is to group hits that "touch" so they can be part of common fit, those that
357  // don't "touch" are fit independently. So here we build the output vector to achieve that
358  // Get a container for the hits...
359  HitCandidateVec groupedHitVec;
360 
361  // Initialize the end of the last hit which we'll set to the first input hit's stop
362  size_t lastStopTick = hitCandidateVec.front().stopTick;
363 
364  // Step through the input hit candidates and group them by proximity
365  for (const auto& hitCandidate : hitCandidateVec) {
366  // Small pulse height hits should not be considered?
367  if (hitCandidate.hitHeight > fMinHitHeight) {
368  // Check condition that we have a new grouping
369  if (hitCandidate.startTick > lastStopTick + fNumInterveningTicks &&
370  !groupedHitVec.empty()) {
371  mergedHitsVec.emplace_back(groupedHitVec);
372 
373  groupedHitVec.clear();
374  }
375 
376  // Add the current hit to the current group
377  groupedHitVec.emplace_back(hitCandidate);
378 
379  lastStopTick = hitCandidate.stopTick;
380  }
381  }
382 
383  // Check end condition
384  if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
385 
386  return;
387  }
388 
391  Waveform::const_iterator stopItr) const
392  {
393  // reset the min iterator and search forward to find the nearest minimum
394  Waveform::const_iterator lastItr = maxItr;
395 
396  // The strategy is simple...
397  // We are at a maximum so we search forward until we find the lowest negative point
398  while ((lastItr + 1) != stopItr) {
399  if (*(lastItr + 1) > *lastItr) break;
400 
401  lastItr++;
402  }
403 
404  // The minimum will be the last iterator value...
405  return lastItr;
406  }
407 
410  Waveform::const_iterator startItr) const
411  {
412  // Set the internal loop variable...
413  Waveform::const_iterator lastItr = minItr;
414 
415  // One extra condition to watch for here, make sure we can actually back up!
416  if (std::distance(startItr, minItr) > 0) {
417  // Similar to searching for a maximum, we loop backward over ticks looking for the waveform to start decreasing
418  while ((lastItr - 1) != startItr) {
419  if (*(lastItr - 1) < *lastItr) break;
420 
421  lastItr--;
422  }
423  }
424 
425  return lastItr;
426  }
427 
430  Waveform::const_iterator startItr) const
431  {
432  Waveform::const_iterator lastItr = maxItr;
433 
434  // If we can't back up then there is nothing to do
435  if (std::distance(startItr, lastItr) > 0) {
436  // In theory, we are starting at a maximum and want to find the "start" of the candidate peak
437  // Ideally we would look to search backward to the point where the (derivative) waveform crosses zero again.
438  // However, the complication is that we need to watch for the case where two peaks are merged together and
439  // we might run through another peak before crossing zero...
440  // So... loop until we hit the startItr...
441  Waveform::const_iterator loopItr = lastItr - 1;
442 
443  while (loopItr != startItr) {
444  // Ideal world case, we cross zero... but we might encounter a minimum... or an inflection point
445  if (*loopItr < 0. || !(*loopItr < *lastItr)) break;
446 
447  lastItr = loopItr--;
448  }
449  }
450  else
451  lastItr = startItr;
452 
453  return lastItr;
454  }
455 
458  Waveform::const_iterator stopItr) const
459  {
460  Waveform::const_iterator lastItr = minItr;
461 
462  // If we can't go forward then there is really nothing to do
463  if (std::distance(minItr, stopItr) > 1) {
464  // Pretty much the same strategy as for finding the start tick...
465  Waveform::const_iterator loopItr = lastItr + 1;
466 
467  while (loopItr != stopItr) {
468  // Ideal case that we have crossed zero coming from a minimum... but watch for a maximum as well
469  if (*loopItr > 0. || !(*loopItr > *lastItr)) break;
470 
471  lastItr = loopItr++;
472  }
473  }
474 
475  return lastItr;
476  }
477 
479 }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
Utilities related to art service access.
void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
T * get() const
Definition: ServiceHandle.h:69
This is the interface class for tools/algorithms that perform various operations on waveforms...
Waveform::const_iterator findNearestMax(Waveform::const_iterator, Waveform::const_iterator) const
const geo::GeometryCore * fGeometry
intermediate_table::const_iterator const_iterator
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::unique_ptr< reco_tool::IWaveformTool > fWaveformTool
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
std::map< size_t, int > fChannelCntMap
void MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const HitCandidateVec &, MergeHitCandidateVec &) const override
Waveform::const_iterator findStopTick(Waveform::const_iterator, Waveform::const_iterator) const
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
TDirectory * dir
Definition: macro.C:5
CandHitDerivative(const fhicl::ParameterSet &pset)
This provides an interface for tools which are tasked with finding candidate hits on input waveforms...
Waveform::const_iterator findStartTick(Waveform::const_iterator, Waveform::const_iterator) const
art::TFileDirectory * fHistDirectory
std::vector< HitCandidateVec > MergeHitCandidateVec
art framework interface to geometry description
Waveform::const_iterator findNearestMin(Waveform::const_iterator, Waveform::const_iterator) const
std::vector< HitCandidate > HitCandidateVec