Public Member Functions | Private Member Functions | Private Attributes | List of all members
reco_tool::CandHitDerivative Class Reference
Inheritance diagram for reco_tool::CandHitDerivative:
reco_tool::ICandidateHitFinder

Public Member Functions

 CandHitDerivative (const fhicl::ParameterSet &pset)
 
void findHitCandidates (const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
 
void MergeHitCandidates (const recob::Wire::RegionsOfInterest_t::datarange_t &, const HitCandidateVec &, MergeHitCandidateVec &) const override
 

Private Member Functions

void findHitCandidates (Waveform::const_iterator, Waveform::const_iterator, const size_t, int, float, HitCandidateVec &) const
 
Waveform::const_iterator findNearestMax (Waveform::const_iterator, Waveform::const_iterator) const
 
Waveform::const_iterator findNearestMin (Waveform::const_iterator, Waveform::const_iterator) const
 
Waveform::const_iterator findStartTick (Waveform::const_iterator, Waveform::const_iterator) const
 
Waveform::const_iterator findStopTick (Waveform::const_iterator, Waveform::const_iterator) const
 
- Private Member Functions inherited from reco_tool::ICandidateHitFinder
virtual ~ICandidateHitFinder () noexcept=default
 

Private Attributes

size_t fPlane
 
int fMinDeltaTicks
 
int fMaxDeltaTicks
 
float fMinDeltaPeaks
 
float fMinHitHeight
 
size_t fNumInterveningTicks
 
bool fOutputHistograms
 
art::TFileDirectory * fHistDirectory
 
TH1F * fDStopStartHist
 
TH1F * fDMaxTickMinTickHist
 
TH1F * fDMaxDerivMinDerivHist
 
std::map< size_t, int > fChannelCntMap
 
std::unique_ptr< reco_tool::IWaveformToolfWaveformTool
 
const geo::GeometryCorefGeometry = lar::providerFrom<geo::Geometry>()
 

Additional Inherited Members

- Private Types inherited from reco_tool::ICandidateHitFinder
using HitCandidateVec = std::vector< HitCandidate >
 
using MergeHitCandidateVec = std::vector< HitCandidateVec >
 
using Waveform = std::vector< float >
 

Detailed Description

Definition at line 22 of file CandHitDerivative_tool.cc.

Constructor & Destructor Documentation

reco_tool::CandHitDerivative::CandHitDerivative ( const fhicl::ParameterSet pset)
explicit

Definition at line 82 of file CandHitDerivative_tool.cc.

83  {
84  fPlane = pset.get<size_t>("Plane", 0);
85  fMinDeltaTicks = pset.get<int>("MinDeltaTicks", 0);
86  fMaxDeltaTicks = pset.get<int>("MaxDeltaTicks", 30);
87  fMinDeltaPeaks = pset.get<float>("MinDeltaPeaks", 0.025);
88  fMinHitHeight = pset.get<float>("MinHitHeight", 2.0);
89  fNumInterveningTicks = pset.get<size_t>("NumInterveningTicks", 6);
90  fOutputHistograms = pset.get<bool>("OutputHistograms", false);
91 
92  // Recover the baseline tool
94  art::make_tool<reco_tool::IWaveformTool>(pset.get<fhicl::ParameterSet>("WaveformAlgs"));
95 
96  // If asked, define the global histograms
97  if (fOutputHistograms) {
98  // Access ART's TFileService, which will handle creating and writing
99  // histograms and n-tuples for us.
101 
102  fHistDirectory = tfs.get();
103 
104  // Make a directory for these histograms
105  art::TFileDirectory dir = fHistDirectory->mkdir(Form("HitPlane_%1zu", fPlane));
106 
108  dir.make<TH1F>(Form("DStopStart_%1zu", fPlane), ";Delta Stop/Start;", 200, 0., 200.);
110  dir.make<TH1F>(Form("DMaxTMinT_%1zu", fPlane), ";Delta Max/Min Tick;", 200, 0., 200.);
112  dir.make<TH1F>(Form("DMaxDMinD_%1zu", fPlane), ";Delta Max/Min Deriv;", 200, 0., 200.);
113  }
114 
115  return;
116  }
T * get() const
Definition: ServiceHandle.h:63
string dir
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::unique_ptr< reco_tool::IWaveformTool > fWaveformTool
art::TFileDirectory * fHistDirectory

Member Function Documentation

void reco_tool::CandHitDerivative::findHitCandidates ( const recob::Wire::RegionsOfInterest_t::datarange_t &  dataRange,
const size_t  roiStartTick,
const size_t  channel,
const size_t  eventCount,
HitCandidateVec hitCandidateVec 
) const
overridevirtual

Implements reco_tool::ICandidateHitFinder.

Definition at line 119 of file CandHitDerivative_tool.cc.

125  {
126  // In this case we want to find hit candidates based on the derivative of of the input waveform
127  // We get this from our waveform algs too...
128  Waveform rawDerivativeVec;
129  Waveform derivativeVec;
130 
131  // Recover the actual waveform
132  const Waveform& waveform = dataRange.data();
133 
134  fWaveformTool->firstDerivative(waveform, rawDerivativeVec);
135  fWaveformTool->triangleSmooth(rawDerivativeVec, derivativeVec);
136 
137  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
138  size_t plane = wids[0].Plane;
139  size_t cryo = wids[0].Cryostat;
140  size_t tpc = wids[0].TPC;
141  size_t wire = wids[0].Wire;
142 
143  // Just make sure the input candidate hit vector has been cleared
144  hitCandidateVec.clear();
145 
146  // Now find the hits
147  findHitCandidates(derivativeVec.begin(),
148  derivativeVec.end(),
149  roiStartTick,
152  hitCandidateVec);
153 
154  if (hitCandidateVec.empty()) {
155  if (plane == 0) {
156  std::cout << "** C/T/P: " << cryo << "/" << tpc << "/" << plane << ", wire: " << wire
157  << " has not hits with input size: " << waveform.size() << std::endl;
158  }
159  }
160 
161  // Reset the hit height from the input waveform
162  for (auto& hitCandidate : hitCandidateVec) {
163  size_t centerIdx = hitCandidate.hitCenter;
164 
165  hitCandidate.hitHeight = waveform.at(centerIdx);
166  }
167 
168  // Keep track of histograms if requested
169  if (fOutputHistograms) {
170  // Recover the details...
171  // std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
172  // size_t plane = wids[0].Plane;
173  // size_t cryo = wids[0].Cryostat;
174  // size_t tpc = wids[0].TPC;
175  // size_t wire = wids[0].Wire;
176 
177  size_t channelCnt = fChannelCntMap[channel]++;
178 
179  // Make a directory for these histograms
180  art::TFileDirectory dir = fHistDirectory->mkdir(
181  Form("HitPlane_%1zu/ev%04zu/c%1zut%1zuwire_%05zu", plane, eventCount, cryo, tpc, wire));
182 
183  size_t waveformSize = waveform.size();
184  int waveStart = roiStartTick;
185  int waveStop = waveStart + waveformSize;
186 
187  TProfile* waveHist = dir.make<TProfile>(
188  Form("HWfm_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
189  "Waveform",
190  waveformSize,
191  waveStart,
192  waveStop,
193  -500.,
194  500.);
195  TProfile* derivHist = dir.make<TProfile>(
196  Form("HDer_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
197  "Derivative",
198  waveformSize,
199  waveStart,
200  waveStop,
201  -500.,
202  500.);
203  TProfile* candHitHist = dir.make<TProfile>(
204  Form("HCan_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
205  "Cand Hits",
206  waveformSize,
207  waveStart,
208  waveStop,
209  -500.,
210  500.);
211  TProfile* maxDerivHist = dir.make<TProfile>(
212  Form("HMax_%03zu_ctw%01zu-%01zu-%01zu-%05zu", channelCnt, cryo, tpc, plane, wire),
213  "Maxima",
214  waveformSize,
215  waveStart,
216  waveStop,
217  -500.,
218  500.);
219 
220  // Fill wave/derivative
221  for (size_t idx = 0; idx < waveform.size(); idx++) {
222  waveHist->Fill(roiStartTick + idx, waveform.at(idx));
223  derivHist->Fill(roiStartTick + idx, derivativeVec.at(idx));
224  }
225 
226  // Fill hits
227  for (const auto& hitCandidate : hitCandidateVec) {
228  candHitHist->Fill(hitCandidate.hitCenter, hitCandidate.hitHeight);
229  maxDerivHist->Fill(hitCandidate.maxTick, hitCandidate.maxDerivative);
230  maxDerivHist->Fill(hitCandidate.minTick, hitCandidate.minDerivative);
231 
232  fDStopStartHist->Fill(hitCandidate.stopTick - hitCandidate.startTick, 1.);
233  fDMaxTickMinTickHist->Fill(hitCandidate.minTick - hitCandidate.maxTick, 1.);
234  fDMaxDerivMinDerivHist->Fill(hitCandidate.maxDerivative - hitCandidate.minDerivative, 1.);
235  }
236  }
237 
238  return;
239  }
void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
const geo::GeometryCore * fGeometry
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
uint8_t channel
Definition: CRTFragment.hh:201
string dir
std::unique_ptr< reco_tool::IWaveformTool > fWaveformTool
std::map< size_t, int > fChannelCntMap
art::TFileDirectory * fHistDirectory
QTextStream & endl(QTextStream &s)
void reco_tool::CandHitDerivative::findHitCandidates ( Waveform::const_iterator  startItr,
Waveform::const_iterator  stopItr,
const size_t  roiStartTick,
int  dTicksThreshold,
float  dPeakThreshold,
HitCandidateVec hitCandidateVec 
) const
private

Definition at line 242 of file CandHitDerivative_tool.cc.

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  }
void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t &, const size_t, const size_t, const size_t, HitCandidateVec &) const override
Waveform::const_iterator findNearestMax(Waveform::const_iterator, Waveform::const_iterator) const
intermediate_table::const_iterator const_iterator
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Waveform::const_iterator findStopTick(Waveform::const_iterator, Waveform::const_iterator) const
Waveform::const_iterator findStartTick(Waveform::const_iterator, Waveform::const_iterator) const
Waveform::const_iterator findNearestMin(Waveform::const_iterator, Waveform::const_iterator) const
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findNearestMax ( Waveform::const_iterator  minItr,
Waveform::const_iterator  startItr 
) const
private

Definition at line 411 of file CandHitDerivative_tool.cc.

413  {
414  // Set the internal loop variable...
415  Waveform::const_iterator lastItr = minItr;
416 
417  // One extra condition to watch for here, make sure we can actually back up!
418  if (std::distance(startItr, minItr) > 0) {
419  // Similar to searching for a maximum, we loop backward over ticks looking for the waveform to start decreasing
420  while ((lastItr - 1) != startItr) {
421  if (*(lastItr - 1) < *lastItr) break;
422 
423  lastItr--;
424  }
425  }
426 
427  return lastItr;
428  }
intermediate_table::const_iterator const_iterator
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findNearestMin ( Waveform::const_iterator  maxItr,
Waveform::const_iterator  stopItr 
) const
private

Definition at line 392 of file CandHitDerivative_tool.cc.

394  {
395  // reset the min iterator and search forward to find the nearest minimum
396  Waveform::const_iterator lastItr = maxItr;
397 
398  // The strategy is simple...
399  // We are at a maximum so we search forward until we find the lowest negative point
400  while ((lastItr + 1) != stopItr) {
401  if (*(lastItr + 1) > *lastItr) break;
402 
403  lastItr++;
404  }
405 
406  // The minimum will be the last iterator value...
407  return lastItr;
408  }
intermediate_table::const_iterator const_iterator
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findStartTick ( Waveform::const_iterator  maxItr,
Waveform::const_iterator  startItr 
) const
private

Definition at line 431 of file CandHitDerivative_tool.cc.

433  {
434  Waveform::const_iterator lastItr = maxItr;
435 
436  // If we can't back up then there is nothing to do
437  if (std::distance(startItr, lastItr) > 0) {
438  // In theory, we are starting at a maximum and want to find the "start" of the candidate peak
439  // Ideally we would look to search backward to the point where the (derivative) waveform crosses zero again.
440  // However, the complication is that we need to watch for the case where two peaks are merged together and
441  // we might run through another peak before crossing zero...
442  // So... loop until we hit the startItr...
443  Waveform::const_iterator loopItr = lastItr - 1;
444 
445  while (loopItr != startItr) {
446  // Ideal world case, we cross zero... but we might encounter a minimum... or an inflection point
447  if (*loopItr < 0. || !(*loopItr < *lastItr)) break;
448 
449  lastItr = loopItr--;
450  }
451  }
452  else
453  lastItr = startItr;
454 
455  return lastItr;
456  }
intermediate_table::const_iterator const_iterator
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
ICandidateHitFinder::Waveform::const_iterator reco_tool::CandHitDerivative::findStopTick ( Waveform::const_iterator  minItr,
Waveform::const_iterator  stopItr 
) const
private

Definition at line 459 of file CandHitDerivative_tool.cc.

461  {
462  Waveform::const_iterator lastItr = minItr;
463 
464  // If we can't go forward then there is really nothing to do
465  if (std::distance(minItr, stopItr) > 1) {
466  // Pretty much the same strategy as for finding the start tick...
467  Waveform::const_iterator loopItr = lastItr + 1;
468 
469  while (loopItr != stopItr) {
470  // Ideal case that we have crossed zero coming from a minimum... but watch for a maximum as well
471  if (*loopItr > 0. || !(*loopItr > *lastItr)) break;
472 
473  lastItr = loopItr++;
474  }
475  }
476 
477  return lastItr;
478  }
intermediate_table::const_iterator const_iterator
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
void reco_tool::CandHitDerivative::MergeHitCandidates ( const recob::Wire::RegionsOfInterest_t::datarange_t &  rangeData,
const HitCandidateVec hitCandidateVec,
MergeHitCandidateVec mergedHitsVec 
) const
overridevirtual

Implements reco_tool::ICandidateHitFinder.

Definition at line 350 of file CandHitDerivative_tool.cc.

354  {
355  // If nothing on the input end then nothing to do
356  if (hitCandidateVec.empty()) return;
357 
358  // The idea is to group hits that "touch" so they can be part of common fit, those that
359  // don't "touch" are fit independently. So here we build the output vector to achieve that
360  // Get a container for the hits...
361  HitCandidateVec groupedHitVec;
362 
363  // Initialize the end of the last hit which we'll set to the first input hit's stop
364  size_t lastStopTick = hitCandidateVec.front().stopTick;
365 
366  // Step through the input hit candidates and group them by proximity
367  for (const auto& hitCandidate : hitCandidateVec) {
368  // Small pulse height hits should not be considered?
369  if (hitCandidate.hitHeight > fMinHitHeight) {
370  // Check condition that we have a new grouping
371  if (hitCandidate.startTick > lastStopTick + fNumInterveningTicks &&
372  !groupedHitVec.empty()) {
373  mergedHitsVec.emplace_back(groupedHitVec);
374 
375  groupedHitVec.clear();
376  }
377 
378  // Add the current hit to the current group
379  groupedHitVec.emplace_back(hitCandidate);
380 
381  lastStopTick = hitCandidate.stopTick;
382  }
383  }
384 
385  // Check end condition
386  if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec);
387 
388  return;
389  }
std::vector< HitCandidate > HitCandidateVec

Member Data Documentation

std::map<size_t, int> reco_tool::CandHitDerivative::fChannelCntMap
mutableprivate

Definition at line 72 of file CandHitDerivative_tool.cc.

TH1F* reco_tool::CandHitDerivative::fDMaxDerivMinDerivHist
private

Definition at line 70 of file CandHitDerivative_tool.cc.

TH1F* reco_tool::CandHitDerivative::fDMaxTickMinTickHist
private

Definition at line 69 of file CandHitDerivative_tool.cc.

TH1F* reco_tool::CandHitDerivative::fDStopStartHist
private

Definition at line 68 of file CandHitDerivative_tool.cc.

const geo::GeometryCore* reco_tool::CandHitDerivative::fGeometry = lar::providerFrom<geo::Geometry>()
private

Definition at line 77 of file CandHitDerivative_tool.cc.

art::TFileDirectory* reco_tool::CandHitDerivative::fHistDirectory
private

Definition at line 65 of file CandHitDerivative_tool.cc.

int reco_tool::CandHitDerivative::fMaxDeltaTicks
private

Definition at line 59 of file CandHitDerivative_tool.cc.

float reco_tool::CandHitDerivative::fMinDeltaPeaks
private

Definition at line 60 of file CandHitDerivative_tool.cc.

int reco_tool::CandHitDerivative::fMinDeltaTicks
private

Definition at line 58 of file CandHitDerivative_tool.cc.

float reco_tool::CandHitDerivative::fMinHitHeight
private

Definition at line 61 of file CandHitDerivative_tool.cc.

size_t reco_tool::CandHitDerivative::fNumInterveningTicks
private

Definition at line 62 of file CandHitDerivative_tool.cc.

bool reco_tool::CandHitDerivative::fOutputHistograms
private

Definition at line 63 of file CandHitDerivative_tool.cc.

size_t reco_tool::CandHitDerivative::fPlane
private

Definition at line 57 of file CandHitDerivative_tool.cc.

std::unique_ptr<reco_tool::IWaveformTool> reco_tool::CandHitDerivative::fWaveformTool
private

Definition at line 75 of file CandHitDerivative_tool.cc.


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