LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CRHitRemoval_module.cc
Go to the documentation of this file.
1 //
3 // Class: CRHitRemoval
4 // Module Type: producer
5 // File: CRHitRemoval_module.cc
6 //
7 // This module produces RecoBase/Hit objects after removing those
8 // deemed to be due to CR muons.
9 //
10 // Configuration parameters:
11 //
12 // CosmicProducerLabels - a list of cosmic ray producers which should be or'ed
13 // HitProducerLabel - the producer of the recob::Hit objects
14 // PFParticleProducerLabel - the producer of the recob::PFParticles to consider
15 // TrackProducerLabel - the producer of the recob::Track objects
16 // CosmicTagThresholds - a vector of thresholds to apply to label as cosmic
17 // EndTickPadding - # ticks to "pad" the end tick to account for possible
18 // uncertainty in drift velocity
19 //
20 // Created by Tracy Usher (usher@slac.stanford.edu) on September 18, 2014
21 //
23 
24 #include <algorithm>
25 
36 
47 
48 class CRHitRemoval : public art::EDProducer {
49 public:
50  // Copnstructors, destructor.
51  explicit CRHitRemoval(fhicl::ParameterSet const& pset);
52 
53  // Overrides.
54  virtual void produce(art::Event& e);
55  virtual void beginJob();
56  virtual void endJob();
57 
58 private:
59  // define vector for hits to make sure of uniform use
60  using HitPtrVector = std::vector<art::Ptr<recob::Hit>>;
61 
62  // Methods
63  void collectPFParticleHits(const recob::PFParticle* pfParticle,
64  const art::Handle<std::vector<recob::PFParticle>>& pfParticleHandle,
65  const art::FindManyP<recob::Cluster>& partToClusAssns,
66  const art::FindManyP<recob::Hit>& clusToHitAssns,
67  HitPtrVector& hitVec);
68 
72 
76 
77  void FilterHits(HitPtrVector& hits, HitPtrVector& used_hits);
78 
79  // Fcl parameters.
80  std::vector<std::string> fCosmicProducerLabels;
81  std::string fHitProducerLabel;
83  std::vector<std::string> fTrackProducerLabels;
84  std::vector<std::string> fAssnProducerLabels;
85 
86  std::vector<double> fCosmicTagThresholds;
87 
89 
94 
95  // Statistics.
96  int fNumEvent;
98 };
99 
101 
102 //----------------------------------------------------------------------------
109 CRHitRemoval::CRHitRemoval(fhicl::ParameterSet const& pset)
110  : EDProducer{pset}, fNumEvent(0), fNumCRRejects(0)
111 {
112  fCosmicProducerLabels = pset.get<std::vector<std::string>>("CosmicProducerLabels");
113  fHitProducerLabel = pset.get<std::string>("HitProducerLabel");
114  fPFParticleProducerLabel = pset.get<std::string>("PFParticleProducerLabel");
115  fTrackProducerLabels = pset.get<std::vector<std::string>>("TrackProducerLabels");
116  fAssnProducerLabels = pset.get<std::vector<std::string>>("AssnProducerLabels");
117  fCosmicTagThresholds = pset.get<std::vector<double>>("CosmicTagThresholds");
118  fEndTickPadding = pset.get<int>("EndTickPadding", 50);
119  fMaxOutOfTime = pset.get<int>("MaxOutOfTime", 4);
120 
121  // let HitCollectionCreator declare that we are going to produce
122  // hits and associations with wires and raw digits
123  // (with no particular product label)
125 
126  // Report.
127  mf::LogInfo("CRHitRemoval") << "CRHitRemoval configured\n";
128 }
129 
130 //----------------------------------------------------------------------------
133 {
134  auto const* geo = lar::providerFrom<geo::Geometry>();
135  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
136  auto const detp =
138 
139  float const samplingRate = sampling_rate(clock_data);
140  float const driftVelocity = detp.DriftVelocity(detp.Efield(), detp.Temperature()); // cm/us
141 
142  fDetectorWidthTicks = 2 * geo->DetHalfWidth() / (driftVelocity * samplingRate / 1000);
143  fMinTickDrift =
144  clock_data.Time2Tick(clock_data.TriggerTime()); // this is the hardware trigger time
146 }
147 
148 //----------------------------------------------------------------------------
162 {
163  ++fNumEvent;
164 
165  // Start by looking up the original hits
167  evt.getByLabel(fHitProducerLabel, hitHandle);
168 
169  // If there are no hits then there should be no output
170  if (!hitHandle.isValid()) return;
171 
172  // also get the associated wires
173  // Note that we no longer assume RawDigit associations and will not preserve them if they exist
174  // we assume they have been created by the same module as the hits
175  art::FindOneP<recob::Wire> ChannelHitWires(hitHandle, evt, fHitProducerLabel);
176 
177  HitPtrVector ChHits;
178  art::fill_ptr_vector(ChHits, hitHandle);
179 
180  // this object contains the hit collection
181  // and its associations to wires and raw digits:
182  recob::HitCollectionCreator hcol(evt, ChannelHitWires.isValid(), false);
183 
184  // Now recover thre remaining collections of objects in the event store that we need
185  // Recover the PFParticles that are ultimately associated to tracks
186  // This will be the key to recovering the hits
188  evt.getByLabel(fPFParticleProducerLabel, pfParticleHandle);
189 
190  // Without a valid collection of PFParticles we can't do the hit removal
191  if (!pfParticleHandle.isValid()) {
192  copyAllHits(ChHits, ChannelHitWires, hcol);
193 
194  // put the hit collection and associations into the event
195  hcol.put_into(evt);
196 
197  return;
198  }
199 
200  // Recover the clusters so we can do associations to the hits
201  // In theory the clusters come from the same producer as the PFParticles
203  evt.getByLabel(fPFParticleProducerLabel, clusterHandle);
204 
205  // If there are no clusters then something is really wrong
206  if (!clusterHandle.isValid()) {
207  copyAllHits(ChHits, ChannelHitWires, hcol);
208 
209  // put the hit collection and associations into the event
210  hcol.put_into(evt);
211 
212  return;
213  }
214 
215  // Build the list of cosmic tag producers, associations and thresholds
216  std::vector<art::Handle<std::vector<anab::CosmicTag>>> cosmicHandleVec;
217  std::vector<std::map<int, std::vector<art::Ptr<recob::Track>>>> cosmicToTrackVecMapVec;
218  std::vector<std::map<int, std::vector<art::Ptr<recob::PFParticle>>>> trackToPFParticleVecMapVec;
219  std::vector<double> thresholdVec;
220 
221  for (size_t idx = 0; idx != fCosmicProducerLabels.size(); idx++) {
222  std::string& handleLabel(fCosmicProducerLabels[idx]);
223 
225  evt.getByLabel(handleLabel, cosmicHandle);
226 
227  if (cosmicHandle.isValid()) {
228  // Look up the associations to tracks
229  art::FindManyP<recob::Track> cosmicTrackAssns(cosmicHandle, evt, handleLabel);
230 
231  // Get the handle to the tracks
233  evt.getByLabel(fTrackProducerLabels[idx], trackHandle);
234 
235  // This should be the case
236  if (trackHandle.isValid()) {
237  std::map<int, std::vector<art::Ptr<recob::Track>>> cosmicToTrackVecMap;
238  std::map<int, std::vector<art::Ptr<recob::PFParticle>>> trackToPFParticleVecMap;
239 
240  // Loop through the cosmic tags
241  for (size_t cosmicIdx = 0; cosmicIdx < cosmicHandle->size(); cosmicIdx++) {
242  art::Ptr<anab::CosmicTag> cosmicTag(cosmicHandle, cosmicIdx);
243 
244  std::vector<art::Ptr<recob::Track>> cosmicToTrackVec =
245  cosmicTrackAssns.at(cosmicTag.key());
246 
247  cosmicToTrackVecMap[cosmicTag.key()] = cosmicToTrackVec;
248 
249  art::FindManyP<recob::PFParticle> trackPFParticleAssns(
250  trackHandle, evt, fAssnProducerLabels[idx]);
251 
252  for (auto& track : cosmicToTrackVec) {
253  std::vector<art::Ptr<recob::PFParticle>> trackToPFParticleVec =
254  trackPFParticleAssns.at(track.key());
255 
256  for (auto& pfParticle : trackToPFParticleVec)
257  trackToPFParticleVecMap[track.key()].push_back(pfParticle);
258  }
259  }
260 
261  // Store these collections
262  cosmicHandleVec.emplace_back(cosmicHandle);
263  cosmicToTrackVecMapVec.emplace_back(cosmicToTrackVecMap);
264  trackToPFParticleVecMapVec.emplace_back(trackToPFParticleVecMap);
265  thresholdVec.push_back(fCosmicTagThresholds[idx]);
266  }
267  }
268  }
269 
270  // No cosmic tags then nothing to do here
271  if (cosmicHandleVec.empty()) {
272  copyAllHits(ChHits, ChannelHitWires, hcol);
273 
274  // put the hit collection and associations into the event
275  hcol.put_into(evt);
276 
277  return;
278  }
279 
280  // Now we walk up the remaining list of associations needed to go from CR tags to hits
281  // We'll need to go from tracks to PFParticles
282  // From PFParticles we go to clusters
283  art::FindManyP<recob::Cluster> clusterAssns(pfParticleHandle, evt, fPFParticleProducerLabel);
284 
285  // Likewise, recover the collection of associations to hits
286  art::FindManyP<recob::Hit> clusterHitAssns(clusterHandle, evt, fPFParticleProducerLabel);
287 
288  // No point double counting hits
289  std::set<const recob::PFParticle*> taggedSet;
290 
291  // Start the identification of hits to remove. The outer loop is over the various producers of
292  // the CosmicTag objects we're examininig
293  for (size_t idx = 0; idx != cosmicHandleVec.size(); idx++) {
294  // Obviously, dereference the handle and associations
295  const art::Handle<std::vector<anab::CosmicTag>>& cosmicHandle(cosmicHandleVec[idx]);
296  const std::map<int, std::vector<art::Ptr<recob::Track>>>& crTagTrackVec(
297  cosmicToTrackVecMapVec[idx]);
298  const std::map<int, std::vector<art::Ptr<recob::PFParticle>>>& trackToPFParticleVecMap(
299  trackToPFParticleVecMapVec[idx]);
300 
301  for (size_t crIdx = 0; crIdx != cosmicHandle->size(); crIdx++) {
302  art::Ptr<anab::CosmicTag> cosmicTag(cosmicHandle, crIdx);
303 
304  // If this was tagged as a CR muon then we have work to do!
305  if (cosmicTag->CosmicScore() > thresholdVec[idx]) {
306  // Recover the associated track
307  std::vector<art::Ptr<recob::Track>> trackVec(crTagTrackVec.at(cosmicTag.key()));
308 
309  // Loop over the tracks (almost always only 1)
310  for (const auto& track : trackVec) {
311  std::map<int, std::vector<art::Ptr<recob::PFParticle>>>::const_iterator
312  trackToPFParticleVecItr = trackToPFParticleVecMap.find(track.key());
313 
314  if (trackToPFParticleVecItr != trackToPFParticleVecMap.end()) {
315  // Loop through them
316  for (const auto& pfParticlePtr : trackToPFParticleVecItr->second) {
317  // Get bare pointer
318  const recob::PFParticle* pfParticle = pfParticlePtr.get();
319 
320  // A cosmic ray must be a primary (by fiat)
321  while (!pfParticle->IsPrimary())
322  pfParticle =
323  art::Ptr<recob::PFParticle>(pfParticleHandle, pfParticle->Parent()).get();
324 
325  // Add to our list of tagged PFParticles
326  taggedSet.insert(pfParticle);
327  }
328  }
329  }
330  }
331  }
332  }
333 
334  // If no PFParticles have been tagged then nothing to do
335  if (!taggedSet.empty()) {
336  // This may all seem backwards... but what you want to do is remove all the hits which are associated to tagged
337  // cosmic ray tracks/PFParticles and what you want to leave is all other hits. This includes hits on other PFParticles
338  // as well as unassociated (as yet) hits.
339  // One also needs to deal with a slight complication with hits in 2D which are shared between untagged PFParticles and
340  // tagged ones... we should leave these in.
341 
342  // All that is left is to go through the PFParticles which have not been tagged and collect up all their hits
343  // to output to our new collection
344  // Note that this SHOULD take care of the case of shared 2D hits automagically since if the PFParticle has not
345  // been tagged and it shares hits we'll pick those up here.
346  HitPtrVector taggedHits;
347  HitPtrVector untaggedHits;
348 
349  // Loop through the PFParticles and build out the list of hits on untagged PFParticle trees
350  for (const auto& pfParticle : *pfParticleHandle) {
351  // Start with only primaries
352  if (!pfParticle.IsPrimary()) continue;
353 
354  // Temporary container for these hits
355  HitPtrVector tempHits;
356 
357  // Find the hits associated to this untagged PFParticle
358  collectPFParticleHits(&pfParticle, pfParticleHandle, clusterAssns, clusterHitAssns, tempHits);
359 
360  // One more possible chance at identifying tagged hits...
361  // Check these hits to see if any lie outside time window
362  bool goodHits(true);
363 
364  if (taggedSet.find(&pfParticle) != taggedSet.end()) goodHits = false;
365 
366  if (goodHits) {
367  int nOutOfTime(0);
368 
369  for (const auto& hit : tempHits) {
370  // Check on out of time hits
371  if (hit->PeakTimeMinusRMS() < fMinTickDrift || hit->PeakTimePlusRMS() > fMaxTickDrift)
372  nOutOfTime++;
373 
374  if (nOutOfTime > fMaxOutOfTime) {
375  goodHits = false;
376  break;
377  }
378  }
379  }
380 
381  if (goodHits)
382  std::copy(tempHits.begin(), tempHits.end(), std::back_inserter(untaggedHits));
383  else
384  std::copy(tempHits.begin(), tempHits.end(), std::back_inserter(taggedHits));
385  }
386 
387  // First task - remove hits from the tagged hit collection that are inthe untagged hits (shared hits)
388  FilterHits(taggedHits, untaggedHits);
389 
390  // Now filter the tagged hits from the total hit collection
391  FilterHits(ChHits, taggedHits);
392  }
393 
394  // Copy our new hit collection to the output
395  copyInTimeHits(ChHits, ChannelHitWires, hcol);
396 
397  // put the hit collection and associations into the event
398  hcol.put_into(evt);
399 }
400 
401 //----------------------------------------------------------------------------
417  const recob::PFParticle* pfParticle,
418  const art::Handle<std::vector<recob::PFParticle>>& pfParticleHandle,
419  const art::FindManyP<recob::Cluster>& partToClusAssns,
420  const art::FindManyP<recob::Hit>& clusToHitAssns,
421  HitPtrVector& hitVec)
422 {
423  // Recover the clusters associated to the input PFParticle
424  std::vector<art::Ptr<recob::Cluster>> clusterVec = partToClusAssns.at(pfParticle->Self());
425 
426  int minTick(fMinTickDrift);
427  int maxTick(fMaxTickDrift);
428 
429  // Loop over the clusters and grab the associated hits
430  for (const auto& cluster : clusterVec) {
431  std::vector<art::Ptr<recob::Hit>> clusHitVec = clusToHitAssns.at(cluster.key());
432  hitVec.insert(hitVec.end(), clusHitVec.begin(), clusHitVec.end());
433 
434  int minHitTick = (*std::min_element(clusHitVec.begin(),
435  clusHitVec.end(),
436  [](const auto& hit, const auto& min) {
437  return hit->PeakTimeMinusRMS() < min->PeakTimeMinusRMS();
438  }))
439  ->PeakTimeMinusRMS();
440  int maxHitTick = (*std::max_element(clusHitVec.begin(),
441  clusHitVec.end(),
442  [](const auto& hit, const auto& max) {
443  return hit->PeakTimePlusRMS() > max->PeakTimePlusRMS();
444  }))
445  ->PeakTimePlusRMS();
446 
447  if (minHitTick < minTick) minTick = minHitTick;
448  if (maxHitTick < maxTick) maxTick = maxHitTick;
449  }
450 
451  // Loop over the daughters of this particle and remove their hits as well
452  for (const auto& daughterId : pfParticle->Daughters()) {
453  art::Ptr<recob::PFParticle> daughter(pfParticleHandle, daughterId);
454 
456  daughter.get(), pfParticleHandle, partToClusAssns, clusToHitAssns, hitVec);
457  }
458 
459  return;
460 }
461 
463  art::FindOneP<recob::Wire>& wireAssns,
464  recob::HitCollectionCreator& newHitCollection)
465 {
466  for (const auto& hitPtr : inputHits) {
467  art::Ptr<recob::Wire> wire = wireAssns.at(hitPtr.key());
468 
469  // just copy it
470  newHitCollection.emplace_back(*hitPtr, wire);
471  }
472 
473  return;
474 }
475 
477  art::FindOneP<recob::Wire>& wireAssns,
478  recob::HitCollectionCreator& newHitCollection)
479 {
480  for (const auto& hitPtr : inputHits) {
481  // Check on out of time hits
482  if (hitPtr->PeakTimeMinusRMS() < fMinTickDrift || hitPtr->PeakTimePlusRMS() > fMaxTickDrift)
483  continue;
484 
485  art::Ptr<recob::Wire> wire = wireAssns.at(hitPtr.key());
486 
487  // just copy it
488  newHitCollection.emplace_back(*hitPtr, wire);
489  }
490 
491  return;
492 }
493 
494 //----------------------------------------------------------------------------
495 // Filter a collection of hits (set difference).
496 // This function is copied from Track3DKalmanHit_module.cc
497 //
498 // Arguments:
499 //
500 // hits - Hit collection from which hits should be removed.
501 // used_hits - Hits to remove.
502 //
504 {
505  if (used_hits.size() > 0) {
506  // Make sure both hit collections are sorted.
507  std::stable_sort(hits.begin(), hits.end());
508  std::stable_sort(used_hits.begin(), used_hits.end());
509 
510  // Do set difference operation.
511  std::vector<art::Ptr<recob::Hit>>::iterator it = std::set_difference(
512  hits.begin(), hits.end(), used_hits.begin(), used_hits.end(), hits.begin());
513 
514  // Truncate hit collection.
515  hits.erase(it, hits.end());
516  }
517 }
518 
519 //----------------------------------------------------------------------------
522 {
523  double aveCRPerEvent = fNumEvent > 0 ? double(fNumCRRejects) / double(fNumEvent) : 0.;
524 
525  mf::LogInfo("CRHitRemoval") << "CRHitRemoval statistics:\n"
526  << " Number of events = " << fNumEvent << "\n"
527  << " Number of Cosmic Rays found = " << fNumCRRejects << ", "
528  << aveCRPerEvent << " average/event";
529 }
const std::vector< size_t > & Daughters() const
Returns the collection of daughter particles.
Definition: PFParticle.h:110
std::string fPFParticleProducerLabel
PFParticle producer.
Utilities related to art service access.
size_t Self() const
Returns the index of this particle.
Definition: PFParticle.h:88
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
int fEndTickPadding
Padding the end tick.
void copyAllHits(std::vector< art::Ptr< recob::Hit >> &, art::FindOneP< recob::Wire > &, recob::HitCollectionCreator &)
Cluster finding and building.
int fNumCRRejects
Number of tracks produced.
float & CosmicScore()
Definition: CosmicTag.h:54
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
bool isValid() const noexcept
Definition: Handle.h:203
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Helper functions to create a hit.
void hits()
Definition: readHits.C:15
A class handling a collection of hits and its associations.
Definition: HitCreator.h:481
std::vector< std::string > fAssnProducerLabels
Track to PFParticle assns producer.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
size_t Parent() const
Definition: PFParticle.h:92
parameter set interface
std::vector< art::Ptr< recob::Hit >> HitPtrVector
std::string fHitProducerLabel
The full collection of hits.
key_type key() const noexcept
Definition: Ptr.h:166
int fMaxTickDrift
Ending tick.
Provides recob::Track data product.
bool IsPrimary() const
Returns whether the particle is the root of the flow.
Definition: PFParticle.h:82
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:286
virtual void endJob()
End job method.
std::vector< double > fCosmicTagThresholds
Thresholds for tagging.
Declaration of cluster object.
Detector simulation of raw signals on wires.
ProducesCollector & producesCollector() noexcept
void collectPFParticleHits(const recob::PFParticle *pfParticle, const art::Handle< std::vector< recob::PFParticle >> &pfParticleHandle, const art::FindManyP< recob::Cluster > &partToClusAssns, const art::FindManyP< recob::Hit > &clusToHitAssns, HitPtrVector &hitVec)
virtual void produce(art::Event &e)
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
int fMaxOutOfTime
Max hits that can be out of time before rejecting.
CRHitRemoval(fhicl::ParameterSet const &pset)
void copyInTimeHits(std::vector< art::Ptr< recob::Hit >> &, art::FindOneP< recob::Wire > &, recob::HitCollectionCreator &)
int fDetectorWidthTicks
Effective drift time in ticks.
std::vector< std::string > fTrackProducerLabels
Track producer.
int fNumEvent
Number of events seen.
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Float_t e
Definition: plot.C:35
T const * get() const
Definition: Ptr.h:138
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Namespace collecting geometry-related classes utilities.
Float_t track
Definition: plot.C:35
int fMinTickDrift
Starting tick.
void FilterHits(HitPtrVector &hits, HitPtrVector &used_hits)
virtual void beginJob()
Begin job method.
art framework interface to geometry description
std::vector< std::string > fCosmicProducerLabels
List of cosmic tagger producers.