LArSoft  v10_04_05
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 
81  };
82 
83  //----------------------------------------------------------------------
84  // Constructor.
86  {
87  fPlane = pset.get<size_t>("Plane", 0);
88  fMinDeltaTicks = pset.get<int>("MinDeltaTicks", 0);
89  fMaxDeltaTicks = pset.get<int>("MaxDeltaTicks", 30);
90  fMinDeltaPeaks = pset.get<float>("MinDeltaPeaks", 0.025);
91  fMinHitHeight = pset.get<float>("MinHitHeight", 2.0);
92  fNumInterveningTicks = pset.get<size_t>("NumInterveningTicks", 6);
93  fOutputHistograms = pset.get<bool>("OutputHistograms", false);
94 
95  // Recover the baseline tool
97  art::make_tool<reco_tool::IWaveformTool>(pset.get<fhicl::ParameterSet>("WaveformAlgs"));
98 
99  // If asked, define the global histograms
100  if (fOutputHistograms) {
101  // Access ART's TFileService, which will handle creating and writing
102  // histograms and n-tuples for us.
104 
105  fHistDirectory = tfs.get();
106 
107  // Make a directory for these histograms
108  art::TFileDirectory dir = fHistDirectory->mkdir(Form("HitPlane_%1zu", fPlane));
109 
111  dir.make<TH1F>(Form("DStopStart_%1zu", fPlane), ";Delta Stop/Start;", 200, 0., 200.);
113  dir.make<TH1F>(Form("DMaxTMinT_%1zu", fPlane), ";Delta Max/Min Tick;", 200, 0., 200.);
115  dir.make<TH1F>(Form("DMaxDMinD_%1zu", fPlane), ";Delta Max/Min Deriv;", 200, 0., 200.);
116  }
117 
118  return;
119  }
120 
122  const recob::Wire::RegionsOfInterest_t::datarange_t& dataRange,
123  const size_t roiStartTick,
124  const size_t channel,
125  const size_t eventCount,
126  HitCandidateVec& hitCandidateVec) const
127  {
128  // In this case we want to find hit candidates based on the derivative of of the input waveform
129  // We get this from our waveform algs too...
130  Waveform rawDerivativeVec;
131  Waveform derivativeVec;
132 
133  // Recover the actual waveform
134  const Waveform& waveform = dataRange.data();
135 
136  fWaveformTool->firstDerivative(waveform, rawDerivativeVec);
137  fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
138 
139  std::vector<geo::WireID> wids = fWireReadoutGeom->ChannelToWire(channel);
140  size_t plane = wids[0].Plane;
141  size_t cryo = wids[0].Cryostat;
142  size_t tpc = wids[0].TPC;
143  size_t wire = wids[0].Wire;
144 
145  // Just make sure the input candidate hit vector has been cleared
146  hitCandidateVec.clear();
147 
148  // Now find the hits
149  findHitCandidates(derivativeVec.begin(),
150  derivativeVec.end(),
151  roiStartTick,
154  hitCandidateVec);
155 
156  if (hitCandidateVec.empty()) {
157  if (plane == 0) {
158  std::cout << "** C/T/P: " << cryo << "/" << tpc << "/" << plane << ", wire: " << wire
159  << " has not hits with input size: " << waveform.size() << std::endl;
160  }
161  }
162 
163  // Reset the hit height from the input waveform
164  for (auto& hitCandidate : hitCandidateVec) {
165  size_t centerIdx = hitCandidate.hitCenter;
166 
167  hitCandidate.hitHeight = waveform.at(centerIdx);
168  }
169 
170  // Keep track of histograms if requested
171  if (fOutputHistograms) {
172 
173  size_t channelCnt = fChannelCntMap[channel]++;
174 
175  // Make a directory for these histograms
176  art::TFileDirectory dir = fHistDirectory->mkdir(
177  Form("HitPlane_%1zu/ev%04zu/c%1zut%1zuwire_%05zu", plane, eventCount, cryo, tpc, wire));
178 
179  size_t waveformSize = waveform.size();
180  int waveStart = roiStartTick;
181  int waveStop = waveStart + waveformSize;
182 
183  TProfile* waveHist = dir.make<TProfile>(
184  Form("HWfm_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
185  "Waveform",
186  waveformSize,
187  waveStart,
188  waveStop,
189  -500.,
190  500.);
191  TProfile* derivHist = dir.make<TProfile>(
192  Form("HDer_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
193  "Derivative",
194  waveformSize,
195  waveStart,
196  waveStop,
197  -500.,
198  500.);
199  TProfile* candHitHist = dir.make<TProfile>(
200  Form("HCan_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
201  "Cand Hits",
202  waveformSize,
203  waveStart,
204  waveStop,
205  -500.,
206  500.);
207  TProfile* maxDerivHist = dir.make<TProfile>(
208  Form("HMax_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
209  "Maxima",
210  waveformSize,
211  waveStart,
212  waveStop,
213  -500.,
214  500.);
215 
216  // Fill wave/derivative
217  for (size_t idx = 0; idx < waveform.size(); idx++) {
218  waveHist->Fill(roiStartTick + idx, waveform.at(idx));
219  derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
220  }
221 
222  // Fill hits
223  for (const auto& hitCandidate : hitCandidateVec) {
224  candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
225  maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
226  maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
227 
228  fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
229  fDMaxTickMinTickHist->Fill(hitCandidate.minTick - hitCandidate.maxTick, 1.);
230  fDMaxDerivMinDerivHist->Fill(hitCandidate.maxDerivative - hitCandidate.minDerivative, 1.);
231  }
232  }
233 
234  return;
235  }
236 
238  Waveform::const_iterator stopItr,
239  const size_t roiStartTick,
240  int dTicksThreshold,
241  float dPeakThreshold,
242  HitCandidateVec& hitCandidateVec) const
243  {
244  // Search for candidate hits...
245  // The idea will be to find the largest deviation in the input derivative waveform as the starting point. Depending
246  // on if a maximum or minimum, we search forward or backward to find the minimum or maximum that our extremum
247  // corresponds to.
248  std::pair<Waveform::const_iterator, Waveform::const_iterator> minMaxPair =
249  std::minmax_element(startItr, stopItr);
250 
251  Waveform::const_iterator maxItr = minMaxPair.second;
252  Waveform::const_iterator minItr = minMaxPair.first;
253 
254  // Use the larger of the two as the starting point and recover the nearest max or min
255  if (std::fabs(*maxItr) > std::fabs(*minItr))
256  minItr = findNearestMin(maxItr, stopItr);
257  else
258  maxItr = findNearestMax(minItr, startItr);
259 
260  int deltaTicks = std::distance(maxItr, minItr);
261  float range = *maxItr - *minItr;
262 
263  // At some point small rolling oscillations on the waveform need to be ignored...
264  if (deltaTicks >= dTicksThreshold && range > dPeakThreshold) {
265  // Need to back up to find zero crossing, this will be the starting point of our
266  // candidate hit but also the endpoint of the pre sub-waveform we'll search next
267  Waveform::const_iterator newEndItr = findStartTick(maxItr, startItr);
268 
269  int startTick = std::distance(startItr, newEndItr);
270 
271  // Now need to go forward to again get close to zero, this will then be the end point
272  // of our candidate hit and the starting point for the post sub-waveform to search
273  Waveform::const_iterator newStartItr = findStopTick(minItr, stopItr);
274 
275  int stopTick = std::distance(startItr, newStartItr);
276 
277  // Find hits in the section of the waveform leading up to this candidate hit
278  if (startTick > dTicksThreshold) {
279  // Special handling for merged hits
280  if (*(newEndItr - 1) > 0.) {
281  dTicksThreshold = 2;
282  dPeakThreshold = 0.;
283  }
284  else {
285  dTicksThreshold = fMinDeltaTicks;
286  dPeakThreshold = fMinDeltaPeaks;
287  }
288 
290  startItr, newEndItr + 1, roiStartTick, dTicksThreshold, dPeakThreshold, hitCandidateVec);
291  }
292 
293  // Create a new hit candidate and store away
294  HitCandidate hitCandidate;
295 
296  Waveform::const_iterator peakItr =
297  std::min_element(maxItr, minItr, [](const auto& left, const auto& right) {
298  return std::fabs(left) < std::fabs(right);
299  });
300 
301  // Check balance
302  if (2 * std::distance(peakItr, minItr) < std::distance(maxItr, peakItr))
303  peakItr--;
304  else if (2 * std::distance(maxItr, peakItr) < std::distance(peakItr, minItr))
305  peakItr++;
306 
307  hitCandidate.startTick = roiStartTick + startTick;
308  hitCandidate.stopTick = roiStartTick + stopTick;
309  hitCandidate.maxTick = roiStartTick + std::distance(startItr, maxItr);
310  hitCandidate.minTick = roiStartTick + std::distance(startItr, minItr);
311  hitCandidate.maxDerivative = maxItr != stopItr ? *maxItr : 0.;
312  hitCandidate.minDerivative = minItr != stopItr ? *minItr : 0.;
313  hitCandidate.hitCenter = roiStartTick + std::distance(startItr, peakItr) + 0.5;
314  hitCandidate.hitSigma = 0.5 * float(hitCandidate.minTick - hitCandidate.maxTick);
315  hitCandidate.hitHeight =
316  hitCandidate.hitSigma * (hitCandidate.maxDerivative - hitCandidate.minDerivative) / 1.2130;
317 
318  hitCandidateVec.push_back(hitCandidate);
319 
320  // Finally, search the section of the waveform following this candidate for more hits
321  if (std::distance(newStartItr, stopItr) > dTicksThreshold) {
322  // Special handling for merged hits
323  if (*(newStartItr + 1) < 0.) {
324  dTicksThreshold = 2;
325  dPeakThreshold = 0.;
326  }
327  else {
328  dTicksThreshold = fMinDeltaTicks;
329  dPeakThreshold = fMinDeltaPeaks;
330  }
331 
332  findHitCandidates(newStartItr,
333  stopItr,
334  roiStartTick + stopTick,
335  dTicksThreshold,
336  dPeakThreshold,
337  hitCandidateVec);
338  }
339  }
340 
341  return;
342  }
343 
344  void CandHitDerivative::MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&,
345  const HitCandidateVec& hitCandidateVec,
346  MergeHitCandidateVec& mergedHitsVec) const
347  {
348  // If nothing on the input end then nothing to do
349  if (hitCandidateVec.empty()) return;
350 
351  // The idea is to group hits that "touch" so they can be part of common fit, those that
352  // don't "touch" are fit independently. So here we build the output vector to achieve that
353  // Get a container for the hits...
354  HitCandidateVec groupedHitVec;
355 
356  // Initialize the end of the last hit which we'll set to the first input hit's stop
357  size_t lastStopTick = hitCandidateVec.front().stopTick;
358 
359  // Step through the input hit candidates and group them by proximity
360  for (const auto& hitCandidate : hitCandidateVec) {
361  // Small pulse height hits should not be considered?
362  if (hitCandidate.hitHeight > fMinHitHeight) {
363  // Check condition that we have a new grouping
364  if (hitCandidate.startTick > lastStopTick + fNumInterveningTicks &&
365  !groupedHitVec.empty()) {
366  mergedHitsVec.emplace_back(groupedHitVec);
367 
368  groupedHitVec.clear();
369  }
370 
371  // Add the current hit to the current group
372  groupedHitVec.emplace_back(hitCandidate);
373 
374  lastStopTick = hitCandidate.stopTick;
375  }
376  }
377 
378  // Check end condition
379  if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
380 
381  return;
382  }
383 
386  Waveform::const_iterator stopItr) const
387  {
388  // reset the min iterator and search forward to find the nearest minimum
389  Waveform::const_iterator lastItr = maxItr;
390 
391  // The strategy is simple...
392  // We are at a maximum so we search forward until we find the lowest negative point
393  while ((lastItr + 1) != stopItr) {
394  if (*(lastItr + 1) > *lastItr) break;
395 
396  lastItr++;
397  }
398 
399  // The minimum will be the last iterator value...
400  return lastItr;
401  }
402 
405  Waveform::const_iterator startItr) const
406  {
407  // Set the internal loop variable...
408  Waveform::const_iterator lastItr = minItr;
409 
410  // One extra condition to watch for here, make sure we can actually back up!
411  if (std::distance(startItr, minItr) > 0) {
412  // Similar to searching for a maximum, we loop backward over ticks looking for the waveform to start decreasing
413  while ((lastItr - 1) != startItr) {
414  if (*(lastItr - 1) < *lastItr) break;
415 
416  lastItr--;
417  }
418  }
419 
420  return lastItr;
421  }
422 
425  Waveform::const_iterator startItr) const
426  {
427  Waveform::const_iterator lastItr = maxItr;
428 
429  // If we can't back up then there is nothing to do
430  if (std::distance(startItr, lastItr) > 0) {
431  // In theory, we are starting at a maximum and want to find the "start" of the candidate peak
432  // Ideally we would look to search backward to the point where the (derivative) waveform crosses zero again.
433  // However, the complication is that we need to watch for the case where two peaks are merged together and
434  // we might run through another peak before crossing zero...
435  // So... loop until we hit the startItr...
436  Waveform::const_iterator loopItr = lastItr - 1;
437 
438  while (loopItr != startItr) {
439  // Ideal world case, we cross zero... but we might encounter a minimum... or an inflection point
440  if (*loopItr < 0. || !(*loopItr < *lastItr)) break;
441 
442  lastItr = loopItr--;
443  }
444  }
445  else
446  lastItr = startItr;
447 
448  return lastItr;
449  }
450 
453  Waveform::const_iterator stopItr) const
454  {
455  Waveform::const_iterator lastItr = minItr;
456 
457  // If we can't go forward then there is really nothing to do
458  if (std::distance(minItr, stopItr) > 1) {
459  // Pretty much the same strategy as for finding the start tick...
460  Waveform::const_iterator loopItr = lastItr + 1;
461 
462  while (loopItr != stopItr) {
463  // Ideal case that we have crossed zero coming from a minimum... but watch for a maximum as well
464  if (*loopItr > 0. || !(*loopItr > *lastItr)) break;
465 
466  lastItr = loopItr++;
467  }
468  }
469 
470  return lastItr;
471  }
472 
474 }
const geo::WireReadoutGeom * fWireReadoutGeom
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
void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
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
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
intermediate_table::const_iterator const_iterator
Interface for a class providing readout channel mapping to geometry.
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::unique_ptr< reco_tool::IWaveformTool > fWaveformTool
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
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
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