LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
lar_cluster3d::SnippetHit3DBuilder Class Reference

SnippetHit3DBuilder class definiton. More...

Inheritance diagram for lar_cluster3d::SnippetHit3DBuilder:
lar_cluster3d::IHit3DBuilder

Public Types

enum  TimeValues { COLLECTARTHITS = 0, BUILDTHREEDHITS = 1, BUILDNEWHITS = 2, NUMTIMEVALUES }
 enumerate the possible values for time checking if monitoring timing More...
 
using RecobHitToPtrMap = std::unordered_map< const recob::Hit *, art::Ptr< recob::Hit >>
 Defines a structure mapping art representation to internal. More...
 

Public Member Functions

 SnippetHit3DBuilder (fhicl::ParameterSet const &pset)
 Constructor. More...
 
void produces (art::ProducesCollector &) override
 Each algorithm may have different objects it wants "produced" so use this to let the top level producer module "know" what it is outputting. More...
 
void Hit3DBuilder (art::Event &, reco::HitPairList &, RecobHitToPtrMap &) override
 Given a set of recob hits, run DBscan to form 3D clusters. More...
 
float getTimeToExecute (IHit3DBuilder::TimeValues index) const override
 If monitoring, recover the time to execute a particular function. More...
 

Private Types

using HitMatchTriplet = std::tuple< const reco::ClusterHit2D *, const reco::ClusterHit2D *, const reco::ClusterHit3D >
 This builds a list of candidate hit pairs from lists of hits on two planes. More...
 
using HitMatchTripletVec = std::vector< HitMatchTriplet >
 
using HitMatchTripletVecMap = std::map< geo::WireID, HitMatchTripletVec >
 
using ChannelStatusVec = std::vector< size_t >
 define data structure for keeping track of channel status More...
 
using ChannelStatusByPlaneVec = std::vector< ChannelStatusVec >
 

Private Member Functions

void CollectArtHits (const art::Event &evt) const
 Extract the ART hits and the ART hit-particle relationships. More...
 
void BuildHit3D (reco::HitPairList &hitPairList) const
 Given the ClusterHit2D objects, build the HitPairMap. More...
 
void CreateNewRecobHitCollection (art::Event &, reco::HitPairList &, std::vector< recob::Hit > &, RecobHitToPtrMap &)
 Create a new 2D hit collection from hits associated to 3D space points. More...
 
void makeWireAssns (const art::Event &, art::Assns< recob::Wire, recob::Hit > &, RecobHitToPtrMap &) const
 Create recob::Wire to recob::Hit associations. More...
 
void makeRawDigitAssns (const art::Event &, art::Assns< raw::RawDigit, recob::Hit > &, RecobHitToPtrMap &) const
 Create raw::RawDigit to recob::Hit associations. More...
 
size_t BuildHitPairMap (PlaneToSnippetHitMap &planeToHitVectorMap, reco::HitPairList &hitPairList) const
 Given the ClusterHit2D objects, build the HitPairMap. More...
 
size_t BuildHitPairMapByTPC (PlaneSnippetHitMapItrPairVec &planeSnippetHitMapItrPairVec, reco::HitPairList &hitPairList) const
 Given the ClusterHit2D objects, build the HitPairMap. More...
 
int findGoodHitPairs (SnippetHitMap::iterator &, SnippetHitMap::iterator &, SnippetHitMap::iterator &, HitMatchTripletVecMap &) const
 
void findGoodTriplets (HitMatchTripletVecMap &, HitMatchTripletVecMap &, reco::HitPairList &) const
 This algorithm takes lists of hit pairs and finds good triplets. More...
 
bool makeHitPair (reco::ClusterHit3D &pairOut, const reco::ClusterHit2D *hit1, const reco::ClusterHit2D *hit2, float hitWidthSclFctr=1., size_t hitPairCntr=0) const
 Make a HitPair object by checking two hits. More...
 
bool makeHitTriplet (reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pairIn, const reco::ClusterHit2D *hit2) const
 Make a 3D HitPair object by checking two hits. More...
 
bool makeDeadChannelPair (reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pair, size_t maxStatus, size_t minStatus) const
 Make a 3D HitPair object from a valid pair and a dead channel in the missing plane. More...
 
bool WireIDsIntersect (const geo::WireID &, const geo::WireID &, geo::WireIDIntersection &) const
 function to detemine if two wires "intersect" (in the 2D sense) More...
 
float closestApproach (const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, float &, float &) const
 function to compute the distance of closest approach and the arc length to the points of closest approach More...
 
const reco::ClusterHit2DFindBestMatchingHit (const Hit2DSet &hit2DSet, const reco::ClusterHit3D &pair, float pairDeltaTimeLimits) const
 A utility routine for finding a 2D hit closest in time to the given pair. More...
 
int FindNumberInRange (const Hit2DSet &hit2DSet, const reco::ClusterHit3D &pair, float range) const
 A utility routine for returning the number of 2D hits from the list in a given range. More...
 
geo::WireID NearestWireID (const Eigen::Vector3f &position, const geo::WireID &wireID) const
 Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range. More...
 
float DistanceFromPointToHitWire (const Eigen::Vector3f &position, const geo::WireID &wireID) const
 Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range. More...
 
void BuildChannelStatusVec () const
 Create the internal channel status vector (assume will eventually be event-by-event) More...
 
float chargeIntegral (float, float, float, int, int) const
 Perform charge integration between limits. More...
 
void clear ()
 clear the tuple vectors before processing next event More...
 

Private Attributes

std::vector< art::InputTagm_hitFinderTagVec
 Data members to follow. More...
 
float m_hitWidthSclFctr
 
float m_deltaPeakTimeSig
 
float m_rangeNumSig
 
float m_LongHitStretchFctr
 
float m_pulseHeightFrac
 
float m_PHLowSelection
 
std::vector< int > m_invalidTPCVec
 
float m_wirePitchScaleFactor
 Scaling factor to determine max distance allowed between candidate pairs. More...
 
float m_maxHit3DChiSquare
 Provide ability to select hits based on "chi square". More...
 
bool m_outputHistograms
 Take the time to create and fill some histograms for diagnostics. More...
 
bool m_enableMonitoring
 
float m_wirePitch [3]
 
std::vector< float > m_timeVector
 
float m_zPosOffset
 
TTree * m_tupleTree
 output analysis tree More...
 
std::vector< float > m_deltaTimeVec
 
std::vector< float > m_chiSquare3DVec
 
std::vector< float > m_maxPullVec
 
std::vector< float > m_overlapFractionVec
 
std::vector< float > m_overlapRangeVec
 
std::vector< float > m_maxDeltaPeakVec
 
std::vector< float > m_maxSideVecVec
 
std::vector< float > m_pairWireDistVec
 
std::vector< float > m_smallChargeDiffVec
 
std::vector< int > m_smallIndexVec
 
std::vector< float > m_qualityMetricVec
 
std::vector< float > m_spacePointChargeVec
 
std::vector< float > m_hitAsymmetryVec
 
Hit2DList m_clusterHit2DMasterList
 
PlaneToSnippetHitMap m_planeToSnippetHitMap
 
PlaneToWireToHitSetMap m_planeToWireToHitSetMap
 
ChannelStatusByPlaneVec m_channelStatus
 
size_t m_numBadChannels
 
bool m_weHaveAllBeenHereBefore = false
 
const geo::Geometrym_geometry
 
const geo::WireReadoutGeomm_wireReadoutGeom
 
const lariov::ChannelStatusProvider * m_channelFilter
 

Detailed Description

SnippetHit3DBuilder class definiton.

Definition at line 81 of file SnippetHit3DBuilder_tool.cc.

Member Typedef Documentation

using lar_cluster3d::SnippetHit3DBuilder::ChannelStatusVec = std::vector<size_t>
private

define data structure for keeping track of channel status

Definition at line 255 of file SnippetHit3DBuilder_tool.cc.

This builds a list of candidate hit pairs from lists of hits on two planes.

Definition at line 163 of file SnippetHit3DBuilder_tool.cc.

using lar_cluster3d::IHit3DBuilder::RecobHitToPtrMap = std::unordered_map<const recob::Hit*, art::Ptr<recob::Hit>>
inherited

Defines a structure mapping art representation to internal.

Definition at line 43 of file IHit3DBuilder.h.

Member Enumeration Documentation

enumerate the possible values for time checking if monitoring timing

Enumerator
COLLECTARTHITS 
BUILDTHREEDHITS 
BUILDNEWHITS 
NUMTIMEVALUES 

Definition at line 56 of file IHit3DBuilder.h.

Constructor & Destructor Documentation

lar_cluster3d::SnippetHit3DBuilder::SnippetHit3DBuilder ( fhicl::ParameterSet const &  pset)
explicit

Constructor.

Parameters
pset

Definition at line 317 of file SnippetHit3DBuilder_tool.cc.

References clear(), Get, art::ServiceHandle< T, SCOPE >::get(), fhicl::ParameterSet::get(), m_chiSquare3DVec, m_deltaPeakTimeSig, m_deltaTimeVec, m_enableMonitoring, m_geometry, m_hitAsymmetryVec, m_hitFinderTagVec, m_hitWidthSclFctr, m_invalidTPCVec, m_LongHitStretchFctr, m_maxDeltaPeakVec, m_maxHit3DChiSquare, m_maxPullVec, m_maxSideVecVec, m_outputHistograms, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_PHLowSelection, m_pulseHeightFrac, m_qualityMetricVec, m_rangeNumSig, m_smallChargeDiffVec, m_smallIndexVec, m_spacePointChargeVec, m_tupleTree, m_wirePitch, m_wirePitchScaleFactor, m_wireReadoutGeom, m_zPosOffset, and geo::WireReadoutGeom::Plane().

319  {
320  m_hitFinderTagVec = pset.get<std::vector<art::InputTag>>(
321  "HitFinderTagVec", std::vector<art::InputTag>() = {"gaushit"});
322  m_enableMonitoring = pset.get<bool>("EnableMonitoring", true);
323  m_hitWidthSclFctr = pset.get<float>("HitWidthScaleFactor", 6.);
324  m_rangeNumSig = pset.get<float>("RangeNumSigma", 3.);
325  m_LongHitStretchFctr = pset.get<float>("LongHitsStretchFactor", 1.5);
326  m_pulseHeightFrac = pset.get<float>("PulseHeightFraction", 0.5);
327  m_PHLowSelection = pset.get<float>("PHLowSelection", 20.);
328  m_deltaPeakTimeSig = pset.get<float>("DeltaPeakTimeSig", 1.7);
329  m_zPosOffset = pset.get<float>("ZPosOffset", 0.0);
330  m_invalidTPCVec = pset.get<std::vector<int>>("InvalidTPCVec", std::vector<int>());
331  m_wirePitchScaleFactor = pset.get<float>("WirePitchScaleFactor", 1.9);
332  m_maxHit3DChiSquare = pset.get<float>("MaxHitChiSquare", 6.0);
333  m_outputHistograms = pset.get<bool>("OutputHistograms", false);
334 
337 
338  // Returns the wire pitch per plane assuming they will be the same for all TPCs
339  constexpr geo::TPCID tpcid{0, 0};
340  m_wirePitch[0] = m_wireReadoutGeom->Plane(geo::PlaneID{tpcid, 0}).WirePitch();
341  m_wirePitch[1] = m_wireReadoutGeom->Plane(geo::PlaneID{tpcid, 1}).WirePitch();
342  m_wirePitch[2] = m_wireReadoutGeom->Plane(geo::PlaneID{tpcid, 2}).WirePitch();
343 
344  // Access ART's TFileService, which will handle creating and writing
345  // histograms and n-tuples for us.
346  if (m_outputHistograms) {
348 
349  m_tupleTree = tfs->make<TTree>("Hit3DBuilderTree", "Tree by SnippetHit3DBuilder");
350 
351  clear();
352 
353  m_tupleTree->Branch("DeltaTime2D", "std::vector<float>", &m_deltaTimeVec);
354  m_tupleTree->Branch("ChiSquare3D", "std::vector<float>", &m_chiSquare3DVec);
355  m_tupleTree->Branch("MaxPullValue", "std::vector<float>", &m_maxPullVec);
356  m_tupleTree->Branch("OverlapFraction", "std::vector<float>", &m_overlapFractionVec);
357  m_tupleTree->Branch("OverlapRange", "std::vector<float>", &m_overlapRangeVec);
358  m_tupleTree->Branch("MaxDeltaPeak", "std::vector<float>", &m_maxDeltaPeakVec);
359  m_tupleTree->Branch("MaxSideVec", "std::vector<float>", &m_maxSideVecVec);
360  m_tupleTree->Branch("PairWireDistVec", "std::vector<float>", &m_pairWireDistVec);
361  m_tupleTree->Branch("SmallChargeDiff", "std::vector<float>", &m_smallChargeDiffVec);
362  m_tupleTree->Branch("SmallChargeIdx", "std::vector<int>", &m_smallIndexVec);
363  m_tupleTree->Branch("QualityMetric", "std::vector<float>", &m_qualityMetricVec);
364  m_tupleTree->Branch("SPCharge", "std::vector<float>", &m_spacePointChargeVec);
365  m_tupleTree->Branch("HitAsymmetry", "std::vector<float>", &m_hitAsymmetryVec);
366  }
367  }
TTree * m_tupleTree
output analysis tree
void clear()
clear the tuple vectors before processing next event
T * get() const
Definition: ServiceHandle.h:69
const lariov::ChannelStatusProvider * m_channelFilter
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
float m_wirePitchScaleFactor
Scaling factor to determine max distance allowed between candidate pairs.
float m_maxHit3DChiSquare
Provide ability to select hits based on "chi square".
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
const geo::WireReadoutGeom * m_wireReadoutGeom
std::vector< art::InputTag > m_hitFinderTagVec
Data members to follow.
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
bool m_outputHistograms
Take the time to create and fill some histograms for diagnostics.

Member Function Documentation

void lar_cluster3d::SnippetHit3DBuilder::BuildChannelStatusVec ( ) const
private

Create the internal channel status vector (assume will eventually be event-by-event)

Definition at line 397 of file SnippetHit3DBuilder_tool.cc.

References geo::WireReadoutGeom::ChannelToWire(), m_channelFilter, m_channelStatus, m_numBadChannels, m_wireReadoutGeom, geo::WireReadoutGeom::Nchannels(), geo::WireReadoutGeom::Nplanes(), geo::WireReadoutGeom::Nwires(), geo::PlaneID::Plane, and geo::WireID::Wire.

Referenced by BuildHit3D().

398  {
399  // This is called each event, clear out the previous version and start over
400  m_channelStatus.clear();
401 
402  m_numBadChannels = 0;
404 
405  // Loop through views/planes to set the wire length vectors
406  constexpr geo::TPCID tpcid{0, 0};
407  for (unsigned int idx = 0; idx < m_channelStatus.size(); idx++) {
408  m_channelStatus[idx] =
410  }
411 
412  // Loop through the channels and mark those that are "bad"
413  for (size_t channel = 0; channel < m_wireReadoutGeom->Nchannels(); channel++) {
414  if (m_channelFilter->IsGood(channel)) continue;
415 
416  std::vector<geo::WireID> wireIDVec = m_wireReadoutGeom->ChannelToWire(channel);
417  geo::WireID wireID = wireIDVec[0];
418  lariov::ChannelStatusProvider::Status_t chanStat = m_channelFilter->Status(channel);
419 
420  m_channelStatus[wireID.Plane][wireID.Wire] = chanStat;
422  }
423  }
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
const lariov::ChannelStatusProvider * m_channelFilter
virtual unsigned int Nchannels() const =0
Returns the total number of channels present (not necessarily contiguous)
std::vector< size_t > ChannelStatusVec
define data structure for keeping track of channel status
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
const geo::WireReadoutGeom * m_wireReadoutGeom
unsigned int Nplanes(TPCID const &tpcid=details::tpc_zero) const
Returns the total number of planes in the specified TPC.
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
void lar_cluster3d::SnippetHit3DBuilder::BuildHit3D ( reco::HitPairList hitPairList) const
private

Given the ClusterHit2D objects, build the HitPairMap.

Driver for processing input 2D hits, transforming to 3D hits and building lists of associated 3D hits (candidate 3D clusters)

Definition at line 496 of file SnippetHit3DBuilder_tool.cc.

References BuildChannelStatusVec(), BuildHitPairMap(), lar_cluster3d::IHit3DBuilder::BUILDTHREEDHITS, m_enableMonitoring, m_planeToSnippetHitMap, and m_timeVector.

Referenced by Hit3DBuilder().

497  {
502  cet::cpu_timer theClockMakeHits;
503 
504  if (m_enableMonitoring) theClockMakeHits.start();
505 
506  // The first task is to take the lists of input 2D hits (a map of view to sorted lists of 2D hits)
507  // and then to build a list of 3D hits to be used in downstream processing
509 
510  size_t numHitPairs = BuildHitPairMap(m_planeToSnippetHitMap, hitPairList);
511 
512  if (m_enableMonitoring) {
513  theClockMakeHits.stop();
514 
515  m_timeVector[BUILDTHREEDHITS] = theClockMakeHits.accumulated_real_time();
516  }
517 
518  mf::LogDebug("Cluster3D") << ">>>>> 3D hit building done, found " << numHitPairs << " 3D Hits"
519  << std::endl;
520  }
void BuildChannelStatusVec() const
Create the internal channel status vector (assume will eventually be event-by-event) ...
size_t BuildHitPairMap(PlaneToSnippetHitMap &planeToHitVectorMap, reco::HitPairList &hitPairList) const
Given the ClusterHit2D objects, build the HitPairMap.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
size_t lar_cluster3d::SnippetHit3DBuilder::BuildHitPairMap ( PlaneToSnippetHitMap planeToHitVectorMap,
reco::HitPairList hitPairList 
) const
private

Given the ClusterHit2D objects, build the HitPairMap.

Given input 2D hits, build out the lists of possible 3D hits

   The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
   However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
   be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
   and then look for the match in the V plane. In the event we don't find the match in the V plane then we
   will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.

Definition at line 544 of file SnippetHit3DBuilder_tool.cc.

References BuildHitPairMapByTPC(), m_geometry, geo::GeometryCore::Ncryostats(), geo::GeometryCore::NTPC(), and lar_cluster3d::SetPairStartTimeOrder().

Referenced by BuildHit3D().

546  {
556  size_t totalNumHits(0);
557  size_t hitPairCntr(0);
558 
559  size_t nTriplets(0);
560  size_t nDeadChanHits(0);
561 
562  // Set up to loop over cryostats and tpcs...
563  for (size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++) {
564  for (size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++) {
566  planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 0));
568  planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 1));
570  planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 2));
571 
572  size_t nPlanesWithHits =
573  (mapItr0 != planeToSnippetHitMap.end() && !mapItr0->second.empty() ? 1 : 0) +
574  (mapItr1 != planeToSnippetHitMap.end() && !mapItr1->second.empty() ? 1 : 0) +
575  (mapItr2 != planeToSnippetHitMap.end() && !mapItr2->second.empty() ? 1 : 0);
576 
577  if (nPlanesWithHits < 2) continue;
578 
579  SnippetHitMap& snippetHitMap0 = mapItr0->second;
580  SnippetHitMap& snippetHitMap1 = mapItr1->second;
581  SnippetHitMap& snippetHitMap2 = mapItr2->second;
582 
583  PlaneSnippetHitMapItrPairVec hitItrVec = {
584  SnippetHitMapItrPair(snippetHitMap0.begin(), snippetHitMap0.end()),
585  SnippetHitMapItrPair(snippetHitMap1.begin(), snippetHitMap1.end()),
586  SnippetHitMapItrPair(snippetHitMap2.begin(), snippetHitMap2.end())};
587 
588  totalNumHits += BuildHitPairMapByTPC(hitItrVec, hitPairList);
589  }
590  }
591 
592  // Return the hit pair list but sorted by z and y positions (faster traversal in next steps)
593  hitPairList.sort(SetPairStartTimeOrder);
594 
595  // Where are we?
596  mf::LogDebug("Cluster3D") << "Total number hits: " << totalNumHits << std::endl;
597  mf::LogDebug("Cluster3D") << "Created a total of " << hitPairList.size()
598  << " hit pairs, counted: " << hitPairCntr << std::endl;
599  mf::LogDebug("Cluster3D") << "-- Triplets: " << nTriplets
600  << ", dead channel pairs: " << nDeadChanHits << std::endl;
601 
602  return hitPairList.size();
603  }
intermediate_table::iterator iterator
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:416
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
size_t BuildHitPairMapByTPC(PlaneSnippetHitMapItrPairVec &planeSnippetHitMapItrPairVec, reco::HitPairList &hitPairList) const
Given the ClusterHit2D objects, build the HitPairMap.
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:303
std::map< HitStartEndPair, HitVector > SnippetHitMap
std::vector< SnippetHitMapItrPair > PlaneSnippetHitMapItrPairVec
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
std::pair< SnippetHitMap::iterator, SnippetHitMap::iterator > SnippetHitMapItrPair
bool SetPairStartTimeOrder(const reco::ClusterHit3D &left, const reco::ClusterHit3D &right)
size_t lar_cluster3d::SnippetHit3DBuilder::BuildHitPairMapByTPC ( PlaneSnippetHitMapItrPairVec planeSnippetHitMapItrPairVec,
reco::HitPairList hitPairList 
) const
private

Given the ClusterHit2D objects, build the HitPairMap.

Given input 2D hits, build out the lists of possible 3D hits

   The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
   However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
   be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
   and then look for the match in the V plane. In the event we don't find the match in the V plane then we
   will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.

Definition at line 605 of file SnippetHit3DBuilder_tool.cc.

References findGoodHitPairs(), and findGoodTriplets().

Referenced by BuildHitPairMap().

608  {
619  // Define functions to set start/end iterators in the loop below
620  auto SetStartIterator =
621  [](SnippetHitMap::iterator startItr, SnippetHitMap::iterator endItr, float startTime) {
622  while (startItr != endItr) {
623  if (startItr->first.second < startTime)
624  startItr++;
625  else
626  break;
627  }
628  return startItr;
629  };
630 
631  auto SetEndIterator =
632  [](SnippetHitMap::iterator lastItr, SnippetHitMap::iterator endItr, float endTime) {
633  while (lastItr != endItr) {
634  if (lastItr->first.first < endTime)
635  lastItr++;
636  else
637  break;
638  }
639  return lastItr;
640  };
641 
642  //*********************************************************************************
643  // Basically, we try to loop until done...
644  while (1) {
645  // Sort so that the earliest hit time will be the first element, etc.
646  std::sort(snippetHitMapItrVec.begin(), snippetHitMapItrVec.end(), SetStartTimeOrder());
647 
648  // Make sure there are still hits on at least
649  int nPlanesWithHits(0);
650 
651  for (auto& pair : snippetHitMapItrVec)
652  if (pair.first != pair.second) nPlanesWithHits++;
653 
654  if (nPlanesWithHits < 2) break;
655 
656  // This loop iteration's snippet iterator
657  SnippetHitMap::iterator firstSnippetItr = snippetHitMapItrVec.front().first;
658 
659  // Set iterators to insure we'll be in the overlap ranges
660  SnippetHitMap::iterator snippetHitMapItr1Start = SetStartIterator(
661  snippetHitMapItrVec[1].first, snippetHitMapItrVec[1].second, firstSnippetItr->first.first);
662  SnippetHitMap::iterator snippetHitMapItr1End = SetEndIterator(
663  snippetHitMapItr1Start, snippetHitMapItrVec[1].second, firstSnippetItr->first.second);
664  SnippetHitMap::iterator snippetHitMapItr2Start = SetStartIterator(
665  snippetHitMapItrVec[2].first, snippetHitMapItrVec[2].second, firstSnippetItr->first.first);
666  SnippetHitMap::iterator snippetHitMapItr2End = SetEndIterator(
667  snippetHitMapItr2Start, snippetHitMapItrVec[2].second, firstSnippetItr->first.second);
668 
669  // Since we'll use these many times in the internal loops, pre make the pairs for the second set of hits
670  HitMatchTripletVecMap pair12Map;
671  HitMatchTripletVecMap pair13Map;
672 
673  size_t n12Pairs =
674  findGoodHitPairs(firstSnippetItr, snippetHitMapItr1Start, snippetHitMapItr1End, pair12Map);
675  size_t n13Pairs =
676  findGoodHitPairs(firstSnippetItr, snippetHitMapItr2Start, snippetHitMapItr2End, pair13Map);
677 
678  if (n12Pairs > n13Pairs)
679  findGoodTriplets(pair12Map, pair13Map, hitPairList);
680  else
681  findGoodTriplets(pair13Map, pair12Map, hitPairList);
682 
683  snippetHitMapItrVec.front().first++;
684  }
685 
686  return hitPairList.size();
687  }
intermediate_table::iterator iterator
std::map< geo::WireID, HitMatchTripletVec > HitMatchTripletVecMap
void findGoodTriplets(HitMatchTripletVecMap &, HitMatchTripletVecMap &, reco::HitPairList &) const
This algorithm takes lists of hit pairs and finds good triplets.
int findGoodHitPairs(SnippetHitMap::iterator &, SnippetHitMap::iterator &, SnippetHitMap::iterator &, HitMatchTripletVecMap &) const
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
float lar_cluster3d::SnippetHit3DBuilder::chargeIntegral ( float  peakMean,
float  peakAmp,
float  peakSigma,
int  low,
int  hi 
) const
private

Perform charge integration between limits.

Definition at line 1329 of file SnippetHit3DBuilder_tool.cc.

Referenced by makeHitTriplet().

1334  {
1335  float integral(0);
1336 
1337  for (int sigPos = low; sigPos < hi; sigPos++) {
1338  float arg = (float(sigPos) - peakMean + 0.5) / peakSigma;
1339  integral += peakAmp * std::exp(-0.5 * arg * arg);
1340  }
1341 
1342  return integral;
1343  }
void lar_cluster3d::SnippetHit3DBuilder::clear ( )
private

clear the tuple vectors before processing next event

Definition at line 380 of file SnippetHit3DBuilder_tool.cc.

References m_chiSquare3DVec, m_deltaTimeVec, m_hitAsymmetryVec, m_maxDeltaPeakVec, m_maxPullVec, m_maxSideVecVec, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_qualityMetricVec, m_smallChargeDiffVec, m_smallIndexVec, and m_spacePointChargeVec.

Referenced by Hit3DBuilder(), and SnippetHit3DBuilder().

381  {
382  m_deltaTimeVec.clear();
383  m_chiSquare3DVec.clear();
384  m_maxPullVec.clear();
385  m_overlapFractionVec.clear();
386  m_overlapRangeVec.clear();
387  m_maxDeltaPeakVec.clear();
388  m_maxSideVecVec.clear();
389  m_pairWireDistVec.clear();
390  m_smallChargeDiffVec.clear();
391  m_smallIndexVec.clear();
392  m_qualityMetricVec.clear();
393  m_spacePointChargeVec.clear();
394  m_hitAsymmetryVec.clear();
395  }
float lar_cluster3d::SnippetHit3DBuilder::closestApproach ( const Eigen::Vector3f &  P0,
const Eigen::Vector3f &  u0,
const Eigen::Vector3f &  P1,
const Eigen::Vector3f &  u1,
float &  arcLen0,
float &  arcLen1 
) const
private

function to compute the distance of closest approach and the arc length to the points of closest approach

Definition at line 1304 of file SnippetHit3DBuilder_tool.cc.

References d, den, and e.

Referenced by WireIDsIntersect().

1310  {
1311  // Technique is to compute the arclength to each point of closest approach
1312  Eigen::Vector3f w0 = P0 - P1;
1313  float a(1.);
1314  float b(u0.dot(u1));
1315  float c(1.);
1316  float d(u0.dot(w0));
1317  float e(u1.dot(w0));
1318  float den(a * c - b * b);
1319 
1320  arcLen0 = (b * e - c * d) / den;
1321  arcLen1 = (a * e - b * d) / den;
1322 
1323  Eigen::Vector3f poca0 = P0 + arcLen0 * u0;
1324  Eigen::Vector3f poca1 = P1 + arcLen1 * u1;
1325 
1326  return (poca0 - poca1).norm();
1327  }
Float_t den
Definition: plot.C:35
Float_t d
Definition: plot.C:235
Float_t e
Definition: plot.C:35
void lar_cluster3d::SnippetHit3DBuilder::CollectArtHits ( const art::Event evt) const
private

Extract the ART hits and the ART hit-particle relationships.

Parameters
evt- the ART event

Recover the 2D hits from art and fill out the local data structures for the 3D clustering

Definition at line 1557 of file SnippetHit3DBuilder_tool.cc.

References geo::WireReadoutGeom::ChannelToWire(), lar_cluster3d::IHit3DBuilder::COLLECTARTHITS, geo::CryostatID::Cryostat, art::ProductRetriever::getByLabel(), art::Handle< T >::isValid(), m_clusterHit2DMasterList, m_enableMonitoring, m_geometry, m_hitFinderTagVec, m_invalidTPCVec, m_planeToSnippetHitMap, m_planeToWireToHitSetMap, m_timeVector, m_weHaveAllBeenHereBefore, m_wireReadoutGeom, geo::GeometryCore::Ncryostats(), geo::GeometryCore::NTPC(), geo::PlaneID::Plane, geo::TPCID::TPC, and detinfo::trigger_offset().

Referenced by Hit3DBuilder().

1558  {
1563  // Start by getting a vector of valid, non empty hit collections to make sure we really have something to do here...
1564  // Here is a container for the hits...
1565  std::vector<const recob::Hit*> recobHitVec;
1566 
1567  // Loop through the list of input sources
1568  for (const auto& inputTag : m_hitFinderTagVec) {
1569  art::Handle<std::vector<recob::Hit>> recobHitHandle;
1570  evt.getByLabel(inputTag, recobHitHandle);
1571 
1572  if (!recobHitHandle.isValid() || recobHitHandle->size() == 0) continue;
1573 
1574  recobHitVec.reserve(recobHitVec.size() + recobHitHandle->size());
1575 
1576  for (const auto& hit : *recobHitHandle)
1577  recobHitVec.push_back(&hit);
1578  }
1579 
1580  // If the vector is empty there is nothing to do
1581  if (recobHitVec.empty()) return;
1582 
1583  cet::cpu_timer theClockMakeHits;
1584 
1585  if (m_enableMonitoring) theClockMakeHits.start();
1586 
1587  // We'll want to correct the hit times for the plane offsets
1588  // (note this is already taken care of when converting to position)
1589  std::map<geo::PlaneID, double> planeOffsetMap;
1590 
1591  // Need the detector properties which needs the clocks
1592  auto const clock_data =
1594  auto const det_prop =
1596 
1597  // Try to output a formatted string
1598  std::string debugMessage("");
1599 
1600  // Initialize the plane to hit vector map
1601  for (size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++) {
1602  for (size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++) {
1603  m_planeToSnippetHitMap[geo::PlaneID(cryoIdx, tpcIdx, 0)] = SnippetHitMap();
1604  m_planeToSnippetHitMap[geo::PlaneID(cryoIdx, tpcIdx, 1)] = SnippetHitMap();
1605  m_planeToSnippetHitMap[geo::PlaneID(cryoIdx, tpcIdx, 2)] = SnippetHitMap();
1606 
1607  // What we want here are the relative offsets between the planes
1608  // Note that plane 0 is assumed the "first" plane and is the reference
1609  planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 0)] = 0.;
1610  planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 1)] =
1611  det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 1)) -
1612  det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0));
1613  planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 2)] =
1614  det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 2)) -
1615  det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0));
1616 
1617  // Should we provide output?
1619  std::ostringstream outputString;
1620 
1621  outputString << "***> plane 0 offset: "
1622  << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 0)]
1623  << ", plane 1: " << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 1)]
1624  << ", plane 2: " << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 2)]
1625  << "\n";
1626  debugMessage += outputString.str();
1627  outputString << " Det prop plane 0: "
1628  << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0))
1629  << ", plane 1: "
1630  << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 1))
1631  << ", plane 2: "
1632  << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 2))
1633  << ", Trig: " << trigger_offset(clock_data) << "\n";
1634  debugMessage += outputString.str();
1635  }
1636  }
1637  }
1638 
1640  mf::LogDebug("Cluster3D") << debugMessage << std::endl;
1641 
1642  std::cout << debugMessage << std::endl;
1643 
1645  }
1646 
1647  // Cycle through the recob hits to build ClusterHit2D objects and insert
1648  // them into the map
1649  for (const auto& recobHit : recobHitVec) {
1650  // Reject hits with negative charge, these are misreconstructed
1651  if (recobHit->Integral() < 0.) continue;
1652 
1653  // For some detectors we can have multiple wire ID's associated to a given channel.
1654  // So we recover the list of these wire IDs
1655  const std::vector<geo::WireID>& wireIDs =
1656  m_wireReadoutGeom->ChannelToWire(recobHit->Channel());
1657 
1658  // Start/End ticks to identify the snippet
1659  HitStartEndPair hitStartEndPair(recobHit->StartTick(), recobHit->EndTick());
1660 
1661  // Can this really happen?
1662  if (hitStartEndPair.second <= hitStartEndPair.first) {
1663  std::cout << "Yes, found a hit with end time less than start time: "
1664  << hitStartEndPair.first << "/" << hitStartEndPair.second
1665  << ", mult: " << recobHit->Multiplicity() << std::endl;
1666  continue;
1667  }
1668 
1669  // And then loop over all possible to build out our maps
1670  for (auto wireID : wireIDs) {
1671  // Check if this is an invalid TPC
1672  // (for example, in protoDUNE there are logical TPC's which see no signal)
1673  if (std::find(m_invalidTPCVec.begin(), m_invalidTPCVec.end(), wireID.TPC) !=
1674  m_invalidTPCVec.end())
1675  continue;
1676 
1677  // Note that a plane ID will define cryostat, TPC and plane
1678  const geo::PlaneID& planeID = wireID.planeID();
1679 
1680  double hitPeakTime(recobHit->PeakTime() - planeOffsetMap[planeID]);
1681  double xPosition(det_prop.ConvertTicksToX(
1682  recobHit->PeakTime(), planeID.Plane, planeID.TPC, planeID.Cryostat));
1683 
1684  m_clusterHit2DMasterList.emplace_back(0, 0., 0., xPosition, hitPeakTime, wireID, recobHit);
1685 
1686  m_planeToSnippetHitMap[planeID][hitStartEndPair].emplace_back(
1687  &m_clusterHit2DMasterList.back());
1688  m_planeToWireToHitSetMap[planeID][wireID.Wire].insert(&m_clusterHit2DMasterList.back());
1689  }
1690  }
1691 
1692  if (m_enableMonitoring) {
1693  theClockMakeHits.stop();
1694 
1695  m_timeVector[COLLECTARTHITS] = theClockMakeHits.accumulated_real_time();
1696  }
1697 
1698  mf::LogDebug("Cluster3D") << ">>>>> Number of ART hits: " << m_clusterHit2DMasterList.size()
1699  << std::endl;
1700  }
std::pair< raw::TDCtick_t, raw::TDCtick_t > HitStartEndPair
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:416
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:303
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
bool isValid() const noexcept
Definition: Handle.h:203
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
const geo::WireReadoutGeom * m_wireReadoutGeom
std::map< HitStartEndPair, HitVector > SnippetHitMap
Detector simulation of raw signals on wires.
std::vector< art::InputTag > m_hitFinderTagVec
Data members to follow.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
int trigger_offset(DetectorClocksData const &data)
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
void lar_cluster3d::SnippetHit3DBuilder::CreateNewRecobHitCollection ( art::Event event,
reco::HitPairList hitPairList,
std::vector< recob::Hit > &  hitPtrVec,
RecobHitToPtrMap recobHitToPtrMap 
)
private

Create a new 2D hit collection from hits associated to 3D space points.

Definition at line 1704 of file SnippetHit3DBuilder_tool.cc.

References lar_cluster3d::IHit3DBuilder::BUILDNEWHITS, reco::ClusterHit2D::getHit(), m_clusterHit2DMasterList, m_enableMonitoring, and m_timeVector.

Referenced by Hit3DBuilder().

1708  {
1709  // Set up the timing
1710  cet::cpu_timer theClockBuildNewHits;
1711 
1712  if (m_enableMonitoring) theClockBuildNewHits.start();
1713 
1714  // We want to build a unique list of hits from the input 3D points which, we know, reuse 2D points frequently
1715  // At the same time we need to create a new 2D hit with the "correct" WireID and replace the old 2D hit in the
1716  // ClusterHit2D object with this new one... while keeping track of the use of the old ones. My head is spinning...
1717  // Declare a set which will allow us to keep track of those CusterHit2D objects we have seen already
1718  std::set<const reco::ClusterHit2D*> visitedHit2DSet;
1719 
1720  // Use this handy art utility to make art::Ptr objects to the new recob::Hits for use in the output phase
1721  art::PtrMaker<recob::Hit> ptrMaker(event);
1722 
1723  // Reserve enough memory to replace every recob::Hit we have considered (this is upper limit)
1724  hitPtrVec.reserve(m_clusterHit2DMasterList.size());
1725 
1726  // Scheme is to loop through all 3D hits, then through each associated ClusterHit2D object
1727  for (reco::ClusterHit3D& hit3D : hitPairList) {
1728  reco::ClusterHit2DVec& hit2DVec = hit3D.getHits();
1729 
1730  // The loop is over the index so we can recover the correct WireID to associate to the new hit when made
1731  for (size_t idx = 0; idx < hit3D.getHits().size(); idx++) {
1732  const reco::ClusterHit2D* hit2D = hit2DVec[idx];
1733 
1734  // Have we seen this 2D hit already?
1735  if (visitedHit2DSet.find(hit2D) == visitedHit2DSet.end()) {
1736  visitedHit2DSet.insert(hit2D);
1737 
1738  // Create and save the new recob::Hit with the correct WireID
1739  hitPtrVec.emplace_back(
1740  recob::HitCreator(*hit2D->getHit(), hit3D.getWireIDs()[idx]).copy());
1741 
1742  // Recover a pointer to it...
1743  recob::Hit* newHit = &hitPtrVec.back();
1744 
1745  // Create a mapping from this hit to an art Ptr representing it
1746  recobHitToPtrMap[newHit] = ptrMaker(hitPtrVec.size() - 1);
1747 
1748  // And set the pointer to this hit in the ClusterHit2D object
1749  const_cast<reco::ClusterHit2D*>(hit2D)->setHit(newHit);
1750  }
1751  }
1752  }
1753 
1754  size_t numNewHits = hitPtrVec.size();
1755 
1756  if (m_enableMonitoring) {
1757  theClockBuildNewHits.stop();
1758 
1759  m_timeVector[BUILDNEWHITS] = theClockBuildNewHits.accumulated_real_time();
1760  }
1761 
1762  mf::LogDebug("Cluster3D") << ">>>>> New output recob::Hit size: " << numNewHits << " (vs "
1763  << m_clusterHit2DMasterList.size() << " input)" << std::endl;
1764 
1765  return;
1766  }
const recob::Hit * getHit() const
Definition: Cluster3D.h:77
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:87
std::vector< const reco::ClusterHit2D * > ClusterHit2DVec
Definition: Cluster3D.h:90
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
float lar_cluster3d::SnippetHit3DBuilder::DistanceFromPointToHitWire ( const Eigen::Vector3f &  position,
const geo::WireID wireID 
) const
private

Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range.

Definition at line 1501 of file SnippetHit3DBuilder_tool.cc.

References util::abs(), geo::WireGeo::Direction(), geo::vect::dot(), geo::WireGeo::GetCenter(), geo::WireGeo::HalfL(), m_wireReadoutGeom, and geo::WireReadoutGeom::Wire().

Referenced by makeHitTriplet().

1503  {
1504  float distance = std::numeric_limits<float>::max();
1505 
1506  // Embed the call to the geometry's services nearest wire id method in a try-catch block
1507  try {
1508  // Recover wire geometry information for each wire
1509  const geo::WireGeo& wireGeo = m_wireReadoutGeom->Wire(wireIDIn);
1510 
1511  // Get wire position and direction for first wire
1512  auto const wirePosArr = wireGeo.GetCenter();
1513 
1514  Eigen::Vector3f wirePos(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
1515  Eigen::Vector3f wireDir(
1516  wireGeo.Direction().X(), wireGeo.Direction().Y(), wireGeo.Direction().Z());
1517 
1518  // Want the hit position to have same x value as wire coordinates
1519  Eigen::Vector3f hitPosition(wirePos[0], position[1], position[2]);
1520 
1521  // Get arc length to doca
1522  double arcLen = (hitPosition - wirePos).dot(wireDir);
1523 
1524  // Make sure arclen is in range
1525  if (abs(arcLen) < wireGeo.HalfL()) {
1526  Eigen::Vector3f docaVec = hitPosition - (wirePos + arcLen * wireDir);
1527 
1528  distance = docaVec.norm();
1529  }
1530  }
1531  catch (std::exception& exc) {
1532  // This can happen, almost always because the coordinates are **just** out of range
1533  mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - "
1534  << exc.what() << std::endl;
1535 
1536  // Assume extremum for wire number depending on z coordinate
1537  distance = 0.;
1538  }
1539 
1540  return distance;
1541  }
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:112
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:219
constexpr auto abs(T v)
Returns the absolute value of the argument.
WireGeo const & Wire(WireID const &wireid) const
Returns the specified wire.
Vector_t Direction() const
Definition: WireGeo.h:287
constexpr auto dot(Vector const &a, OtherVector const &b)
Return cross product of two vectors.
const geo::WireReadoutGeom * m_wireReadoutGeom
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:170
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
const reco::ClusterHit2D * lar_cluster3d::SnippetHit3DBuilder::FindBestMatchingHit ( const Hit2DSet hit2DSet,
const reco::ClusterHit3D pair,
float  pairDeltaTimeLimits 
) const
private

A utility routine for finding a 2D hit closest in time to the given pair.

Definition at line 1418 of file SnippetHit3DBuilder_tool.cc.

References reco::ClusterHit3D::getAvePeakTime().

1422  {
1423  static const float minCharge(0.);
1424 
1425  const reco::ClusterHit2D* bestVHit(0);
1426 
1427  float pairAvePeakTime(pair.getAvePeakTime());
1428 
1429  // Idea is to loop through the input set of hits and look for the best combination
1430  for (const auto& hit2D : hit2DSet) {
1431  if (hit2D->getHit()->Integral() < minCharge) continue;
1432 
1433  float hitVPeakTime(hit2D->getTimeTicks());
1434  float deltaPeakTime(pairAvePeakTime - hitVPeakTime);
1435 
1436  if (deltaPeakTime > pairDeltaTimeLimits) continue;
1437 
1438  if (deltaPeakTime < -pairDeltaTimeLimits) break;
1439 
1440  pairDeltaTimeLimits = fabs(deltaPeakTime);
1441  bestVHit = hit2D;
1442  }
1443 
1444  return bestVHit;
1445  }
float getAvePeakTime() const
Definition: Cluster3D.h:160
int lar_cluster3d::SnippetHit3DBuilder::findGoodHitPairs ( SnippetHitMap::iterator firstSnippetItr,
SnippetHitMap::iterator startItr,
SnippetHitMap::iterator endItr,
HitMatchTripletVecMap hitMatchMap 
) const
private

Definition at line 689 of file SnippetHit3DBuilder_tool.cc.

References art::left(), m_hitWidthSclFctr, m_PHLowSelection, m_pulseHeightFrac, makeHitPair(), art::right(), and geo::WireID::WireID().

Referenced by BuildHitPairMapByTPC().

693  {
694  int numPairs(0);
695 
696  HitVector::iterator firstMaxItr =
697  std::max_element(firstSnippetItr->second.begin(),
698  firstSnippetItr->second.end(),
699  [](const auto& left, const auto& right) {
700  return left->getHit()->PeakAmplitude() < right->getHit()->PeakAmplitude();
701  });
702  float firstPHCut = firstMaxItr != firstSnippetItr->second.end() ?
703  m_pulseHeightFrac * (*firstMaxItr)->getHit()->PeakAmplitude() :
704  4096.;
705 
706  // Loop through the hits on the first snippet
707  for (const auto& hit1 : firstSnippetItr->second) {
708  // Let's focus on the largest hit in the chain
709  if (hit1->getHit()->DegreesOfFreedom() > 1 && hit1->getHit()->PeakAmplitude() < firstPHCut &&
710  hit1->getHit()->PeakAmplitude() < m_PHLowSelection)
711  continue;
712 
713  // Inside loop iterator
714  SnippetHitMap::iterator secondHitItr = startItr;
715 
716  // Loop through the input secon hits and make pairs
717  while (secondHitItr != endItr) {
718  HitVector::iterator secondMaxItr = std::max_element(
719  secondHitItr->second.begin(),
720  secondHitItr->second.end(),
721  [](const auto& left, const auto& right) {
722  return left->getHit()->PeakAmplitude() < right->getHit()->PeakAmplitude();
723  });
724  float secondPHCut = secondMaxItr != secondHitItr->second.end() ?
725  m_pulseHeightFrac * (*secondMaxItr)->getHit()->PeakAmplitude() :
726  0.;
727 
728  for (const auto& hit2 : secondHitItr->second) {
729  // Again, focus on the large hits
730  if (hit2->getHit()->DegreesOfFreedom() > 1 &&
731  hit2->getHit()->PeakAmplitude() < secondPHCut &&
732  hit2->getHit()->PeakAmplitude() < m_PHLowSelection)
733  continue;
734 
735  reco::ClusterHit3D pair;
736 
737  // pair returned with a negative ave time is signal of failure
738  if (!makeHitPair(pair, hit1, hit2, m_hitWidthSclFctr)) continue;
739 
740  std::vector<const recob::Hit*> recobHitVec = {nullptr, nullptr, nullptr};
741 
742  recobHitVec[hit1->WireID().Plane] = hit1->getHit();
743  recobHitVec[hit2->WireID().Plane] = hit2->getHit();
744 
745  geo::WireID wireID = hit2->WireID();
746 
747  hitMatchMap[wireID].emplace_back(hit1, hit2, pair);
748 
749  numPairs++;
750  }
751 
752  secondHitItr++;
753  }
754  }
755 
756  return numPairs;
757  }
intermediate_table::iterator iterator
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
bool makeHitPair(reco::ClusterHit3D &pairOut, const reco::ClusterHit2D *hit1, const reco::ClusterHit2D *hit2, float hitWidthSclFctr=1., size_t hitPairCntr=0) const
Make a HitPair object by checking two hits.
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
constexpr WireID()=default
Default constructor: an invalid TPC ID.
void lar_cluster3d::SnippetHit3DBuilder::findGoodTriplets ( HitMatchTripletVecMap pair12Map,
HitMatchTripletVecMap pair13Map,
reco::HitPairList hitPairList 
) const
private

This algorithm takes lists of hit pairs and finds good triplets.

Definition at line 759 of file SnippetHit3DBuilder_tool.cc.

References art::left(), m_deltaPeakTimeSig, m_numBadChannels, makeDeadChannelPair(), makeHitTriplet(), art::right(), and reco::ClusterHit3D::setID().

Referenced by BuildHitPairMapByTPC().

762  {
763  // Build triplets from the two lists of hit pairs
764  if (!pair12Map.empty()) {
765  // temporary container for dead channel hits
766  std::vector<reco::ClusterHit3D> tempDeadChanVec;
767  reco::ClusterHit3D deadChanPair;
768 
769  // Keep track of which third plane hits have been used
770  std::map<const reco::ClusterHit3D*, bool> usedPairMap;
771 
772  // Initial population of this map with the pair13Map hits
773  for (const auto& pair13 : pair13Map) {
774  for (const auto& hit2Dhit3DPair : pair13.second)
775  usedPairMap[&std::get<2>(hit2Dhit3DPair)] = false;
776  }
777 
778  // The outer loop is over all hit pairs made from the first two plane combinations
779  for (const auto& pair12 : pair12Map) {
780  if (pair12.second.empty()) continue;
781 
782  // This loop is over hit pairs that share the same first two plane wires but may have different
783  // hit times on those wires
784  for (const auto& hit2Dhit3DPair12 : pair12.second) {
785  const reco::ClusterHit3D& pair1 = std::get<2>(hit2Dhit3DPair12);
786 
787  // populate the map with initial value
788  usedPairMap[&pair1] = false;
789 
790  // The simplest approach here is to loop over all possibilities and let the triplet builder weed out the weak candidates
791  for (const auto& pair13 : pair13Map) {
792  if (pair13.second.empty()) continue;
793 
794  for (const auto& hit2Dhit3DPair13 : pair13.second) {
795  // Protect against double counting
796  if (std::get<0>(hit2Dhit3DPair12) != std::get<0>(hit2Dhit3DPair13)) continue;
797 
798  const reco::ClusterHit2D* hit2 = std::get<1>(hit2Dhit3DPair13);
799  const reco::ClusterHit3D& pair2 = std::get<2>(hit2Dhit3DPair13);
800 
801  // If success try for the triplet
802  reco::ClusterHit3D triplet;
803 
804  if (makeHitTriplet(triplet, pair1, hit2)) {
805  triplet.setID(hitPairList.size());
806  hitPairList.emplace_back(triplet);
807  usedPairMap[&pair1] = true;
808  usedPairMap[&pair2] = true;
809  }
810  }
811  }
812  }
813  }
814 
815  // One more loop through the other pairs to check for sick channels
816  if (m_numBadChannels > 0) {
817  for (const auto& pairMapPair : usedPairMap) {
818  if (pairMapPair.second) continue;
819 
820  const reco::ClusterHit3D* pair = pairMapPair.first;
821 
822  // Here we look to see if we failed to make a triplet because the partner wire was dead/noisy/sick
823  if (makeDeadChannelPair(deadChanPair, *pair, 4, 0))
824  tempDeadChanVec.emplace_back(deadChanPair);
825  }
826 
827  // Handle the dead wire triplets
828  if (!tempDeadChanVec.empty()) {
829  // If we have many then see if we can trim down a bit by keeping those with time significance
830  if (tempDeadChanVec.size() > 1) {
831  // Sort by "significance" of agreement
832  std::sort(tempDeadChanVec.begin(),
833  tempDeadChanVec.end(),
834  [](const auto& left, const auto& right) {
835  return left.getDeltaPeakTime() / left.getSigmaPeakTime() <
836  right.getDeltaPeakTime() / right.getSigmaPeakTime();
837  });
838 
839  // What is the range of "significance" from first to last?
840  float firstSig = tempDeadChanVec.front().getDeltaPeakTime() /
841  tempDeadChanVec.front().getSigmaPeakTime();
842  float lastSig =
843  tempDeadChanVec.back().getDeltaPeakTime() / tempDeadChanVec.back().getSigmaPeakTime();
844  float sigRange = lastSig - firstSig;
845 
846  if (lastSig > 0.5 * m_deltaPeakTimeSig && sigRange > 0.5) {
847  // Declare a maximum of 1.5 * the average of the first and last pairs...
848  float maxSignificance = std::max(0.75 * (firstSig + lastSig), 1.0);
849 
850  std::vector<reco::ClusterHit3D>::iterator firstBadElem = std::find_if(
851  tempDeadChanVec.begin(),
852  tempDeadChanVec.end(),
853  [&maxSignificance](const auto& pair) {
854  return pair.getDeltaPeakTime() / pair.getSigmaPeakTime() > maxSignificance;
855  });
856 
857  // But only keep the best 10?
858  if (std::distance(tempDeadChanVec.begin(), firstBadElem) > 20)
859  firstBadElem = tempDeadChanVec.begin() + 20;
860  // Keep at least one hit...
861  else if (firstBadElem == tempDeadChanVec.begin())
862  firstBadElem++;
863 
864  tempDeadChanVec.resize(std::distance(tempDeadChanVec.begin(), firstBadElem));
865  }
866  }
867 
868  for (auto& pair : tempDeadChanVec) {
869  pair.setID(hitPairList.size());
870  hitPairList.emplace_back(pair);
871  }
872  }
873  }
874  }
875 
876  return;
877  }
intermediate_table::iterator iterator
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
bool makeHitTriplet(reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pairIn, const reco::ClusterHit2D *hit2) const
Make a 3D HitPair object by checking two hits.
float getSigmaPeakTime() const
Definition: Cluster3D.h:162
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
bool makeDeadChannelPair(reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pair, size_t maxStatus, size_t minStatus) const
Make a 3D HitPair object from a valid pair and a dead channel in the missing plane.
void setID(const size_t &id) const
Definition: Cluster3D.h:176
float getDeltaPeakTime() const
Definition: Cluster3D.h:161
int lar_cluster3d::SnippetHit3DBuilder::FindNumberInRange ( const Hit2DSet hit2DSet,
const reco::ClusterHit3D pair,
float  range 
) const
private

A utility routine for returning the number of 2D hits from the list in a given range.

Definition at line 1447 of file SnippetHit3DBuilder_tool.cc.

References reco::ClusterHit3D::getAvePeakTime().

1450  {
1451  static const float minCharge(0.);
1452 
1453  int numberInRange(0);
1454  float pairAvePeakTime(pair.getAvePeakTime());
1455 
1456  // Idea is to loop through the input set of hits and look for the best combination
1457  for (const auto& hit2D : hit2DSet) {
1458  if (hit2D->getHit()->Integral() < minCharge) continue;
1459 
1460  float hitVPeakTime(hit2D->getTimeTicks());
1461  float deltaPeakTime(pairAvePeakTime - hitVPeakTime);
1462 
1463  if (deltaPeakTime > range) continue;
1464 
1465  if (deltaPeakTime < -range) break;
1466 
1467  numberInRange++;
1468  }
1469 
1470  return numberInRange;
1471  }
float getAvePeakTime() const
Definition: Cluster3D.h:160
float lar_cluster3d::SnippetHit3DBuilder::getTimeToExecute ( IHit3DBuilder::TimeValues  index) const
inlineoverridevirtual

If monitoring, recover the time to execute a particular function.

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 107 of file SnippetHit3DBuilder_tool.cc.

References tca::evt.

108  {
109  return m_timeVector[index];
110  }
void lar_cluster3d::SnippetHit3DBuilder::Hit3DBuilder ( art::Event evt,
reco::HitPairList hitPairList,
RecobHitToPtrMap clusterHitToArtPtrMap 
)
overridevirtual

Given a set of recob hits, run DBscan to form 3D clusters.

Parameters
hitPairListThe input list of 3D hits to run clustering on
clusterParametersListA list of cluster objects (parameters from associated hits)

Associations with wires.

Associations with raw digits.

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 442 of file SnippetHit3DBuilder_tool.cc.

References BuildHit3D(), clear(), CollectArtHits(), CreateNewRecobHitCollection(), m_clusterHit2DMasterList, m_outputHistograms, m_planeToSnippetHitMap, m_planeToWireToHitSetMap, m_timeVector, m_tupleTree, makeRawDigitAssns(), makeWireAssns(), lar_cluster3d::IHit3DBuilder::NUMTIMEVALUES, and art::Event::put().

445  {
446  // Clear the internal data structures
447  m_clusterHit2DMasterList.clear();
448  m_planeToSnippetHitMap.clear();
449  m_planeToWireToHitSetMap.clear();
450 
451  m_timeVector.resize(NUMTIMEVALUES, 0.);
452 
453  // Get a hit refiner
454  std::unique_ptr<std::vector<recob::Hit>> outputHitPtrVec(new std::vector<recob::Hit>);
455 
456  // Recover the 2D hits and then organize them into data structures which will be used in the
457  // DBscan algorithm for building the 3D clusters
458  CollectArtHits(evt);
459 
460  // If there are no hits in our view/wire data structure then do not proceed with the full analysis
461  if (!m_planeToWireToHitSetMap.empty()) {
462  // Call the algorithm that builds 3D hits
463  BuildHit3D(hitPairList);
464 
465  // If we built 3D points then attempt to output a new hit list as well
466  if (!hitPairList.empty())
467  CreateNewRecobHitCollection(evt, hitPairList, *outputHitPtrVec, clusterHitToArtPtrMap);
468  }
469 
470  // Set up to make the associations (if desired)
472  std::unique_ptr<art::Assns<recob::Wire, recob::Hit>> wireAssns(
474 
475  makeWireAssns(evt, *wireAssns, clusterHitToArtPtrMap);
476 
478  std::unique_ptr<art::Assns<raw::RawDigit, recob::Hit>> rawDigitAssns(
480 
481  makeRawDigitAssns(evt, *rawDigitAssns, clusterHitToArtPtrMap);
482 
483  // Move everything into the event
484  evt.put(std::move(outputHitPtrVec));
485  evt.put(std::move(wireAssns));
486  evt.put(std::move(rawDigitAssns));
487 
488  // Handle tree output too
489  if (m_outputHistograms) {
490  m_tupleTree->Fill();
491 
492  clear();
493  }
494  }
TTree * m_tupleTree
output analysis tree
void clear()
clear the tuple vectors before processing next event
void CreateNewRecobHitCollection(art::Event &, reco::HitPairList &, std::vector< recob::Hit > &, RecobHitToPtrMap &)
Create a new 2D hit collection from hits associated to 3D space points.
void BuildHit3D(reco::HitPairList &hitPairList) const
Given the ClusterHit2D objects, build the HitPairMap.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
void makeWireAssns(const art::Event &, art::Assns< recob::Wire, recob::Hit > &, RecobHitToPtrMap &) const
Create recob::Wire to recob::Hit associations.
void makeRawDigitAssns(const art::Event &, art::Assns< raw::RawDigit, recob::Hit > &, RecobHitToPtrMap &) const
Create raw::RawDigit to recob::Hit associations.
void CollectArtHits(const art::Event &evt) const
Extract the ART hits and the ART hit-particle relationships.
bool m_outputHistograms
Take the time to create and fill some histograms for diagnostics.
bool lar_cluster3d::SnippetHit3DBuilder::makeDeadChannelPair ( reco::ClusterHit3D pairOut,
const reco::ClusterHit3D pair,
size_t  maxStatus,
size_t  minStatus 
) const
private

Make a 3D HitPair object from a valid pair and a dead channel in the missing plane.

Definition at line 1345 of file SnippetHit3DBuilder_tool.cc.

References geo::CryostatID::Cryostat, reco::ClusterHit3D::getHits(), reco::ClusterHit3D::getPosition(), reco::ClusterHit2D::getStatusBits(), m_channelStatus, m_wireReadoutGeom, m_zPosOffset, NearestWireID(), geo::PlaneID::Plane, reco::ClusterHit3D::setPosition(), reco::ClusterHit2D::setStatusBit(), reco::ClusterHit3D::setWireID(), reco::ClusterHit2D::SHAREDINTRIPLET, geo::TPCID::TPC, reco::ClusterHit2D::USEDINTRIPLET, geo::WireID::Wire, reco::ClusterHit2D::WireID(), and geo::WireReadoutGeom::WireIDsIntersect().

Referenced by findGoodTriplets().

1349  {
1350  // Assume failure (most common result)
1351  bool result(false);
1352 
1353  const reco::ClusterHit2D* hit0 = pair.getHits()[0];
1354  const reco::ClusterHit2D* hit1 = pair.getHits()[1];
1355 
1356  size_t missPlane(2);
1357 
1358  // u plane hit is missing
1359  if (!hit0) {
1360  hit0 = pair.getHits()[2];
1361  missPlane = 0;
1362  }
1363  // v plane hit is missing
1364  else if (!hit1) {
1365  hit1 = pair.getHits()[2];
1366  missPlane = 1;
1367  }
1368 
1369  // Which plane is missing?
1370  geo::WireID wireID0 = hit0->WireID();
1371  geo::WireID wireID1 = hit1->WireID();
1372 
1373  // Ok, recover the wireID expected in the third plane...
1374  geo::WireID wireIn(wireID0.Cryostat, wireID0.TPC, missPlane, 0);
1375  geo::WireID wireID = NearestWireID(pair.getPosition(), wireIn);
1376 
1377  // There can be a round off issue so check the next wire as well
1378  bool wireStatus = m_channelStatus[wireID.Plane][wireID.Wire] < maxChanStatus &&
1379  m_channelStatus[wireID.Plane][wireID.Wire] >= minChanStatus;
1380  bool wireOneStatus = m_channelStatus[wireID.Plane][wireID.Wire + 1] < maxChanStatus &&
1381  m_channelStatus[wireID.Plane][wireID.Wire + 1] >= minChanStatus;
1382 
1383  // Make sure they are of at least the minimum status
1384  if (wireStatus || wireOneStatus) {
1385  // Sort out which is the wire we're dealing with
1386  if (!wireStatus) wireID.Wire += 1;
1387 
1388  // Want to refine position since we "know" the missing wire
1389  if (auto widIntersect0 = m_wireReadoutGeom->WireIDsIntersect(wireID0, wireID)) {
1390  if (auto widIntersect1 = m_wireReadoutGeom->WireIDsIntersect(wireID1, wireID)) {
1391  Eigen::Vector3f newPosition(
1392  pair.getPosition()[0], pair.getPosition()[1], pair.getPosition()[2]);
1393 
1394  newPosition[1] = (newPosition[1] + widIntersect0->y + widIntersect1->y) / 3.;
1395  newPosition[2] =
1396  (newPosition[2] + widIntersect0->z + widIntersect1->z - 2. * m_zPosOffset) / 3.;
1397 
1398  pairOut = pair;
1399  pairOut.setWireID(wireID);
1400  pairOut.setPosition(newPosition);
1401 
1406 
1409 
1410  result = true;
1411  }
1412  }
1413  }
1414 
1415  return result;
1416  }
void setWireID(const geo::WireID &wid) const
Definition: Cluster3D.cxx:191
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:155
const geo::WireID & WireID() const
Definition: Cluster3D.h:76
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
const geo::WireReadoutGeom * m_wireReadoutGeom
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
geo::WireID NearestWireID(const Eigen::Vector3f &position, const geo::WireID &wireID) const
Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range...
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:168
void setPosition(const Eigen::Vector3f &pos) const
Definition: Cluster3D.h:183
unsigned getStatusBits() const
Definition: Cluster3D.h:71
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:79
bool lar_cluster3d::SnippetHit3DBuilder::makeHitPair ( reco::ClusterHit3D pairOut,
const reco::ClusterHit2D hit1,
const reco::ClusterHit2D hit2,
float  hitWidthSclFctr = 1.,
size_t  hitPairCntr = 0 
) const
private

Make a HitPair object by checking two hits.

Definition at line 879 of file SnippetHit3DBuilder_tool.cc.

References geo::CryostatID::Cryostat, recob::Hit::DegreesOfFreedom(), reco::ClusterHit2D::getHit(), reco::ClusterHit2D::getStatusBits(), reco::ClusterHit2D::getTimeTicks(), reco::ClusterHit2D::getXPosition(), reco::ClusterHit3D::initialize(), recob::Hit::Integral(), m_deltaPeakTimeSig, m_LongHitStretchFctr, m_zPosOffset, geo::PlaneID::Plane, recob::Hit::RMS(), reco::ClusterHit2D::setStatusBit(), reco::ClusterHit2D::SHAREDINPAIR, geo::TPCID::TPC, reco::ClusterHit2D::USEDINPAIR, reco::ClusterHit2D::WireID(), WireIDsIntersect(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by findGoodHitPairs(), and makeHitTriplet().

884  {
885  // Assume failure
886  bool result(false);
887 
888  // Start by checking time consistency since this is fastest
889  // Wires intersect so now we can check the timing
890  float hit1Peak = hit1->getTimeTicks();
891  float hit1Sigma = hit1->getHit()->RMS();
892 
893  float hit2Peak = hit2->getTimeTicks();
894  float hit2Sigma = hit2->getHit()->RMS();
895 
896  // "Long hits" are an issue... so we deal with these differently
897  int hit1NDF = hit1->getHit()->DegreesOfFreedom();
898  int hit2NDF = hit2->getHit()->DegreesOfFreedom();
899 
900  // Basically, allow the range to extend to the nearest end of the snippet
901  if (hit1NDF < 2)
902  hit1Sigma *= m_LongHitStretchFctr; // This sets the range to the width of the pulse
903  if (hit2NDF < 2) hit2Sigma *= m_LongHitStretchFctr;
904 
905  // The "hit sigma" is the gaussian fit sigma of the hit, we need to expand a bit to allow hit overlap efficiency
906  float hit1Width = hitWidthSclFctr * hit1Sigma;
907  float hit2Width = hitWidthSclFctr * hit2Sigma;
908 
909  // Coarse check hit times are "in range"
910  if (fabs(hit1Peak - hit2Peak) <= (hit1Width + hit2Width)) {
911  // Check to see that hit peak times are consistent with each other
912  float hit1SigSq = hit1Sigma * hit1Sigma;
913  float hit2SigSq = hit2Sigma * hit2Sigma;
914  float deltaPeakTime = std::fabs(hit1Peak - hit2Peak);
915  float sigmaPeakTime = std::sqrt(hit1SigSq + hit2SigSq);
916 
917  // delta peak time consistency check here
918  if (deltaPeakTime <
919  m_deltaPeakTimeSig * sigmaPeakTime) // 2 sigma consistency? (do this way to avoid divide)
920  {
921  // We assume in this routine that we are looking at hits in different views
922  // The first mission is to check that the wires intersect
923  const geo::WireID& hit1WireID = hit1->WireID();
924  const geo::WireID& hit2WireID = hit2->WireID();
925 
926  geo::WireIDIntersection widIntersect;
927 
928  if (WireIDsIntersect(hit1WireID, hit2WireID, widIntersect)) {
929  float oneOverWghts = hit1SigSq * hit2SigSq / (hit1SigSq + hit2SigSq);
930  float avePeakTime = (hit1Peak / hit1SigSq + hit2Peak / hit2SigSq) * oneOverWghts;
931  float totalCharge = hit1->getHit()->Integral() + hit2->getHit()->Integral();
932  float hitChiSquare = std::pow((hit1Peak - avePeakTime), 2) / hit1SigSq +
933  std::pow((hit2Peak - avePeakTime), 2) / hit2SigSq;
934 
935  float xPositionHit1(hit1->getXPosition());
936  float xPositionHit2(hit2->getXPosition());
937  float xPosition = (xPositionHit1 / hit1SigSq + xPositionHit2 / hit2SigSq) * hit1SigSq *
938  hit2SigSq / (hit1SigSq + hit2SigSq);
939 
940  Eigen::Vector3f position(
941  xPosition, float(widIntersect.y), float(widIntersect.z) - m_zPosOffset);
942 
943  // If to here then we need to sort out the hit pair code telling what views are used
944  unsigned statusBits = 1 << hit1->WireID().Plane | 1 << hit2->WireID().Plane;
945 
946  // handle status bits for the 2D hits
951 
954 
955  reco::ClusterHit2DVec hitVector;
956 
957  hitVector.resize(3, NULL);
958 
959  hitVector[hit1->WireID().Plane] = hit1;
960  hitVector[hit2->WireID().Plane] = hit2;
961 
962  unsigned int cryostatIdx = hit1->WireID().Cryostat;
963  unsigned int tpcIdx = hit1->WireID().TPC;
964 
965  // Initialize the wireIdVec
966  std::vector<geo::WireID> wireIDVec = {geo::WireID(cryostatIdx, tpcIdx, 0, 0),
967  geo::WireID(cryostatIdx, tpcIdx, 1, 0),
968  geo::WireID(cryostatIdx, tpcIdx, 2, 0)};
969 
970  wireIDVec[hit1->WireID().Plane] = hit1->WireID();
971  wireIDVec[hit2->WireID().Plane] = hit2->WireID();
972 
973  // For compiling at the moment
974  std::vector<float> hitDelTSigVec = {0., 0., 0.};
975 
976  hitDelTSigVec[hit1->WireID().Plane] = deltaPeakTime / sigmaPeakTime;
977  hitDelTSigVec[hit2->WireID().Plane] = deltaPeakTime / sigmaPeakTime;
978 
979  // Create the 3D cluster hit
980  hitPair.initialize(hitPairCntr,
981  statusBits,
982  position,
983  totalCharge,
984  avePeakTime,
985  deltaPeakTime,
986  sigmaPeakTime,
987  hitChiSquare,
988  0.,
989  0.,
990  0.,
991  0.,
992  hitVector,
993  hitDelTSigVec,
994  wireIDVec);
995 
996  result = true;
997  }
998  }
999  }
1000 
1001  // Send it back
1002  return result;
1003  }
double z
z position of intersection
Definition: geo_types.h:584
float getTimeTicks() const
Definition: Cluster3D.h:75
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:234
const geo::WireID & WireID() const
Definition: Cluster3D.h:76
int DegreesOfFreedom() const
Initial tdc tick for hit.
Definition: Hit.h:274
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:254
const recob::Hit * getHit() const
Definition: Cluster3D.h:77
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
float getXPosition() const
Definition: Cluster3D.h:74
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
std::vector< const reco::ClusterHit2D * > ClusterHit2DVec
Definition: Cluster3D.h:90
bool WireIDsIntersect(const geo::WireID &, const geo::WireID &, geo::WireIDIntersection &) const
function to detemine if two wires "intersect" (in the 2D sense)
double y
y position of intersection
Definition: geo_types.h:583
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
unsigned getStatusBits() const
Definition: Cluster3D.h:71
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:79
bool lar_cluster3d::SnippetHit3DBuilder::makeHitTriplet ( reco::ClusterHit3D pairOut,
const reco::ClusterHit3D pairIn,
const reco::ClusterHit2D hit2 
) const
private

Make a 3D HitPair object by checking two hits.

Definition at line 1005 of file SnippetHit3DBuilder_tool.cc.

References util::abs(), chargeIntegral(), geo::CryostatID::Cryostat, recob::Hit::DegreesOfFreedom(), DistanceFromPointToHitWire(), reco::ClusterHit3D::getAvePeakTime(), reco::ClusterHit2D::getHit(), reco::ClusterHit3D::getHits(), reco::ClusterHit3D::getPosition(), reco::ClusterHit3D::getSigmaPeakTime(), reco::ClusterHit2D::getStatusBits(), reco::ClusterHit2D::getTimeTicks(), reco::ClusterHit2D::getXPosition(), reco::ClusterHit3D::initialize(), m_chiSquare3DVec, m_deltaTimeVec, m_hitAsymmetryVec, m_hitWidthSclFctr, m_LongHitStretchFctr, m_maxDeltaPeakVec, m_maxPullVec, m_maxSideVecVec, m_outputHistograms, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_qualityMetricVec, m_rangeNumSig, m_smallChargeDiffVec, m_smallIndexVec, m_spacePointChargeVec, m_wirePitch, m_wirePitchScaleFactor, makeHitPair(), recob::Hit::PeakAmplitude(), geo::PlaneID::Plane, recob::Hit::RMS(), reco::ClusterHit2D::setStatusBit(), reco::ClusterHit2D::SHAREDINTRIPLET, geo::TPCID::TPC, reco::ClusterHit2D::USEDINTRIPLET, weight, geo::WireID::Wire, reco::ClusterHit2D::WireID(), and geo::WireID::WireID().

Referenced by findGoodTriplets().

1008  {
1009  // Assume failure
1010  bool result(false);
1011 
1012  // We are going to force the wire pitch here, some time in the future we need to fix
1013  static const double wirePitch =
1014  0.5 * m_wirePitchScaleFactor * *std::max_element(m_wirePitch, m_wirePitch + 3);
1015 
1016  // Recover hit info
1017  float hitTimeTicks = hit->getTimeTicks();
1018  float hitSigma = hit->getHit()->RMS();
1019 
1020  // Special case check
1021  if (hitSigma > 2. * hit->getHit()->PeakAmplitude())
1022  hitSigma = 2. * hit->getHit()->PeakAmplitude();
1023 
1024  // Let's do a quick consistency check on the input hits to make sure we are in range...
1025  // Require the W hit to be "in range" with the UV Pair
1026  if (fabs(hitTimeTicks - pair.getAvePeakTime()) <
1027  m_hitWidthSclFctr * (pair.getSigmaPeakTime() + hitSigma)) {
1028  // Check the distance from the point to the wire the hit is on
1029  float hitWireDist = DistanceFromPointToHitWire(pair.getPosition(), hit->WireID());
1030 
1031  if (m_outputHistograms) m_maxSideVecVec.push_back(hitWireDist);
1032 
1033  // Reject hits that are not within range
1034  if (hitWireDist < wirePitch) {
1035  if (m_outputHistograms) m_pairWireDistVec.push_back(hitWireDist);
1036 
1037  // Use the existing code to see the U and W hits are willing to pair with the V hit
1038  reco::ClusterHit3D pair0h;
1039  reco::ClusterHit3D pair1h;
1040 
1041  // Recover all the hits involved
1042  const reco::ClusterHit2DVec& pairHitVec = pair.getHits();
1043  const reco::ClusterHit2D* hit0 = pairHitVec[0];
1044  const reco::ClusterHit2D* hit1 = pairHitVec[1];
1045 
1046  if (!hit0)
1047  hit0 = pairHitVec[2];
1048  else if (!hit1)
1049  hit1 = pairHitVec[2];
1050 
1051  // If good pairs made here then we can try to make a triplet
1052  if (makeHitPair(pair0h, hit0, hit, m_hitWidthSclFctr) &&
1053  makeHitPair(pair1h, hit1, hit, m_hitWidthSclFctr)) {
1054  // Get a copy of the input hit vector (note the order is by plane - by definition)
1055  reco::ClusterHit2DVec hitVector = pair.getHits();
1056 
1057  // include the new hit
1058  hitVector[hit->WireID().Plane] = hit;
1059 
1060  // Set up to get average peak time, hitChiSquare, etc.
1061  unsigned int statusBits(0x7);
1062  float avePeakTime(0.);
1063  float weightSum(0.);
1064  float xPosition(0.);
1065 
1066  // And get the wire IDs
1067  std::vector<geo::WireID> wireIDVec = {geo::WireID(), geo::WireID(), geo::WireID()};
1068 
1069  // First loop through the hits to get WireIDs and calculate the averages
1070  for (size_t planeIdx = 0; planeIdx < 3; planeIdx++) {
1071  const reco::ClusterHit2D* hit2D = hitVector[planeIdx];
1072 
1073  wireIDVec[planeIdx] = hit2D->WireID();
1074 
1077 
1079 
1080  float hitRMS = hit2D->getHit()->RMS();
1081  float peakTime = hit2D->getTimeTicks();
1082 
1083  // Basically, allow the range to extend to the nearest end of the snippet
1084  if (hit2D->getHit()->DegreesOfFreedom() < 2) hitRMS *= m_LongHitStretchFctr;
1085 
1086  float weight = 1. / (hitRMS * hitRMS);
1087 
1088  avePeakTime += peakTime * weight;
1089  xPosition += hit2D->getXPosition() * weight;
1090  weightSum += weight;
1091  }
1092 
1093  avePeakTime /= weightSum;
1094  xPosition /= weightSum;
1095 
1096  Eigen::Vector2f pair0hYZVec(pair0h.getPosition()[1], pair0h.getPosition()[2]);
1097  Eigen::Vector2f pair1hYZVec(pair1h.getPosition()[1], pair1h.getPosition()[2]);
1098  Eigen::Vector2f pairYZVec(pair.getPosition()[1], pair.getPosition()[2]);
1099  Eigen::Vector3f position(xPosition,
1100  float((pairYZVec[0] + pair0hYZVec[0] + pair1hYZVec[0]) / 3.),
1101  float((pairYZVec[1] + pair0hYZVec[1] + pair1hYZVec[1]) / 3.));
1102 
1103  // Armed with the average peak time, now get hitChiSquare and the sig vec
1104  float hitChiSquare(0.);
1105  float sigmaPeakTime(std::sqrt(1. / weightSum));
1106  std::vector<float> hitDelTSigVec;
1107 
1108  for (const auto& hit2D : hitVector) {
1109  float hitRMS = hit2D->getHit()->RMS();
1110 
1111  // Basically, allow the range to extend to the nearest end of the snippet
1112  if (hit2D->getHit()->DegreesOfFreedom() < 2) hitRMS *= m_LongHitStretchFctr;
1113 
1114  float combRMS = std::sqrt(hitRMS * hitRMS - sigmaPeakTime * sigmaPeakTime);
1115  float peakTime = hit2D->getTimeTicks();
1116  float deltaTime = peakTime - avePeakTime;
1117  float hitSig = deltaTime / combRMS;
1118 
1119  hitChiSquare += hitSig * hitSig;
1120 
1121  hitDelTSigVec.emplace_back(std::fabs(hitSig));
1122  }
1123 
1124  if (m_outputHistograms) m_chiSquare3DVec.push_back(hitChiSquare);
1125 
1126  int lowMinIndex(std::numeric_limits<int>::max());
1127  int lowMaxIndex(std::numeric_limits<int>::min());
1128  int hiMinIndex(std::numeric_limits<int>::max());
1129  int hiMaxIndex(std::numeric_limits<int>::min());
1130 
1131  // First task is to get the min/max values for the common overlap region
1132  for (const auto& hit2D : hitVector) {
1133  float range = m_rangeNumSig * hit2D->getHit()->RMS();
1134 
1135  // Basically, allow the range to extend to the nearest end of the snippet
1136  if (hit2D->getHit()->DegreesOfFreedom() < 2) range *= m_LongHitStretchFctr;
1137 
1138  int hitStart = hit2D->getHit()->PeakTime() - range - 0.5;
1139  int hitStop = hit2D->getHit()->PeakTime() + range + 0.5;
1140 
1141  lowMinIndex = std::min(hitStart, lowMinIndex);
1142  lowMaxIndex = std::max(hitStart, lowMaxIndex);
1143  hiMinIndex = std::min(hitStop + 1, hiMinIndex);
1144  hiMaxIndex = std::max(hitStop + 1, hiMaxIndex);
1145  }
1146 
1147  // Keep only "good" hits...
1148  if (hitChiSquare < m_maxHit3DChiSquare && hiMinIndex > lowMaxIndex) {
1149  // One more pass through hits to get charge
1150  std::vector<float> chargeVec;
1151 
1152  for (const auto& hit2D : hitVector)
1153  chargeVec.push_back(chargeIntegral(hit2D->getHit()->PeakTime(),
1154  hit2D->getHit()->PeakAmplitude(),
1155  hit2D->getHit()->RMS(),
1156  lowMaxIndex,
1157  hiMinIndex));
1158 
1159  float totalCharge =
1160  std::accumulate(chargeVec.begin(), chargeVec.end(), 0.) / float(chargeVec.size());
1161  float overlapRange = float(hiMinIndex - lowMaxIndex);
1162  float overlapFraction = overlapRange / float(hiMaxIndex - lowMinIndex);
1163 
1164  // Set up to compute the charge asymmetry
1165  std::vector<float> smallestChargeDiffVec;
1166  std::vector<float> chargeAveVec;
1167  float smallestDiff(std::numeric_limits<float>::max());
1168  float maxDeltaPeak(0.);
1169  size_t chargeIndex(0);
1170 
1171  for (size_t idx = 0; idx < 3; idx++) {
1172  size_t leftIdx = (idx + 2) % 3;
1173  size_t rightIdx = (idx + 1) % 3;
1174 
1175  smallestChargeDiffVec.push_back(std::abs(chargeVec[leftIdx] - chargeVec[rightIdx]));
1176  chargeAveVec.push_back(float(0.5 * (chargeVec[leftIdx] + chargeVec[rightIdx])));
1177 
1178  if (smallestChargeDiffVec.back() < smallestDiff) {
1179  smallestDiff = smallestChargeDiffVec.back();
1180  chargeIndex = idx;
1181  }
1182 
1183  // Take opportunity to look at peak time diff
1184  float deltaPeakTime =
1185  hitVector[leftIdx]->getTimeTicks() - hitVector[rightIdx]->getTimeTicks();
1186 
1187  if (std::abs(deltaPeakTime) > maxDeltaPeak) maxDeltaPeak = std::abs(deltaPeakTime);
1188 
1189  if (m_outputHistograms) m_deltaTimeVec.push_back(deltaPeakTime);
1190  }
1191 
1192  float chargeAsymmetry = (chargeAveVec[chargeIndex] - chargeVec[chargeIndex]) /
1193  (chargeAveVec[chargeIndex] + chargeVec[chargeIndex]);
1194 
1195  // If this is true there has to be a negative charge that snuck in somehow
1196  if (chargeAsymmetry < -1. || chargeAsymmetry > 1.) {
1197  const geo::WireID& hitWireID = hitVector[chargeIndex]->WireID();
1198 
1199  std::cout << "============> Charge asymmetry out of range: " << chargeAsymmetry
1200  << " <============" << std::endl;
1201  std::cout << " hit C: " << hitWireID.Cryostat << ", TPC: " << hitWireID.TPC
1202  << ", Plane: " << hitWireID.Plane << ", Wire: " << hitWireID.Wire
1203  << std::endl;
1204  std::cout << " charge: " << chargeVec[0] << ", " << chargeVec[1] << ", "
1205  << chargeVec[2] << std::endl;
1206  std::cout << " index: " << chargeIndex << ", smallest diff: " << smallestDiff
1207  << std::endl;
1208  return result;
1209  }
1210 
1211  // Usurping "deltaPeakTime" to be the maximum pull
1212  float deltaPeakTime = *std::max_element(hitDelTSigVec.begin(), hitDelTSigVec.end());
1213 
1214  if (m_outputHistograms) {
1215  m_smallChargeDiffVec.push_back(smallestDiff);
1216  m_smallIndexVec.push_back(chargeIndex);
1217  m_maxPullVec.push_back(deltaPeakTime);
1218  m_qualityMetricVec.push_back(hitChiSquare);
1219  m_spacePointChargeVec.push_back(totalCharge);
1220  m_overlapFractionVec.push_back(overlapFraction);
1221  m_overlapRangeVec.push_back(overlapRange);
1222  m_maxDeltaPeakVec.push_back(maxDeltaPeak);
1223  m_hitAsymmetryVec.push_back(chargeAsymmetry);
1224  }
1225 
1226  // Try to weed out cases where overlap doesn't match peak separation
1227  if (maxDeltaPeak > 3. * overlapRange) return result;
1228 
1229  // Create the 3D cluster hit
1230  hitTriplet.initialize(0,
1231  statusBits,
1232  position,
1233  totalCharge,
1234  avePeakTime,
1235  deltaPeakTime,
1236  sigmaPeakTime,
1237  hitChiSquare,
1238  overlapFraction,
1239  chargeAsymmetry,
1240  0.,
1241  0.,
1242  hitVector,
1243  hitDelTSigVec,
1244  wireIDVec);
1245 
1246  result = true;
1247  }
1248  }
1249  }
1250  }
1251 
1252  // return success/fail
1253  return result;
1254  }
float getTimeTicks() const
Definition: Cluster3D.h:75
bool makeHitPair(reco::ClusterHit3D &pairOut, const reco::ClusterHit2D *hit1, const reco::ClusterHit2D *hit2, float hitWidthSclFctr=1., size_t hitPairCntr=0) const
Make a HitPair object by checking two hits.
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:155
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:234
float chargeIntegral(float, float, float, int, int) const
Perform charge integration between limits.
const geo::WireID & WireID() const
Definition: Cluster3D.h:76
constexpr auto abs(T v)
Returns the absolute value of the argument.
int DegreesOfFreedom() const
Initial tdc tick for hit.
Definition: Hit.h:274
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
float DistanceFromPointToHitWire(const Eigen::Vector3f &position, const geo::WireID &wireID) const
Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range...
float m_wirePitchScaleFactor
Scaling factor to determine max distance allowed between candidate pairs.
const recob::Hit * getHit() const
Definition: Cluster3D.h:77
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:238
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
float getXPosition() const
Definition: Cluster3D.h:74
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
Detector simulation of raw signals on wires.
std::vector< const reco::ClusterHit2D * > ClusterHit2DVec
Definition: Cluster3D.h:90
double weight
Definition: plottest35.C:25
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:226
constexpr WireID()=default
Default constructor: an invalid TPC ID.
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
bool m_outputHistograms
Take the time to create and fill some histograms for diagnostics.
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:168
unsigned getStatusBits() const
Definition: Cluster3D.h:71
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:79
void lar_cluster3d::SnippetHit3DBuilder::makeRawDigitAssns ( const art::Event evt,
art::Assns< raw::RawDigit, recob::Hit > &  rawDigitAssns,
RecobHitToPtrMap recobHitPtrMap 
) const
private

Create raw::RawDigit to recob::Hit associations.

Definition at line 1808 of file SnippetHit3DBuilder_tool.cc.

References art::Assns< L, R, D >::addSingle(), raw::RawDigit::Channel(), DEFINE_ART_CLASS_TOOL, art::Assns< L, R, D >::end(), art::ProductRetriever::getValidHandle(), and m_hitFinderTagVec.

Referenced by Hit3DBuilder().

1811  {
1812  // Let's make sure the input associations container is empty
1813  rawDigitAssns = art::Assns<raw::RawDigit, recob::Hit>();
1814 
1815  // First task is to recover all of the previous wire <--> hit associations and map them by channel number
1816  // Create the temporary container
1817  std::unordered_map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> channelToRawDigitMap;
1818 
1819  // Go through the list of input sources and fill out the map
1820  for (const auto& inputTag : m_hitFinderTagVec) {
1822  evt.getValidHandle<std::vector<recob::Hit>>(inputTag);
1823 
1824  art::FindOneP<raw::RawDigit> hitToRawDigitAssns(hitHandle, evt, inputTag);
1825 
1826  if (hitToRawDigitAssns.isValid()) {
1827  for (size_t rawDigitIdx = 0; rawDigitIdx < hitToRawDigitAssns.size(); rawDigitIdx++) {
1828  art::Ptr<raw::RawDigit> rawDigit = hitToRawDigitAssns.at(rawDigitIdx);
1829 
1830  channelToRawDigitMap[rawDigit->Channel()] = rawDigit;
1831  }
1832  }
1833  }
1834 
1835  // Now fill the container
1836  for (const auto& hitPtrPair : recobHitPtrMap) {
1837  raw::ChannelID_t channel = hitPtrPair.first->Channel();
1838 
1839  std::unordered_map<raw::ChannelID_t, art::Ptr<raw::RawDigit>>::iterator chanRawDigitItr =
1840  channelToRawDigitMap.find(channel);
1841 
1842  if (chanRawDigitItr == channelToRawDigitMap.end()) { continue; }
1843 
1844  rawDigitAssns.addSingle(chanRawDigitItr->second, hitPtrPair.second);
1845  }
1846  }
ChannelID_t Channel() const
DAQ channel this raw data was read from.
Definition: RawDigit.h:213
std::vector< art::InputTag > m_hitFinderTagVec
Data members to follow.
const_iterator end() const
Definition: Assns.h:514
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
void addSingle(Ptr< left_t > const &left, Ptr< right_t > const &right, data_t const &data)
Definition: Assns.h:549
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
void lar_cluster3d::SnippetHit3DBuilder::makeWireAssns ( const art::Event evt,
art::Assns< recob::Wire, recob::Hit > &  wireAssns,
RecobHitToPtrMap recobHitPtrMap 
) const
private

Create recob::Wire to recob::Hit associations.

Definition at line 1768 of file SnippetHit3DBuilder_tool.cc.

References art::Assns< L, R, D >::addSingle(), recob::Wire::Channel(), art::Assns< L, R, D >::end(), art::ProductRetriever::getValidHandle(), and m_hitFinderTagVec.

Referenced by Hit3DBuilder().

1771  {
1772  // Let's make sure the input associations container is empty
1774 
1775  // First task is to recover all of the previous wire <--> hit associations and map them by channel number
1776  // Create the temporary container
1777  std::unordered_map<raw::ChannelID_t, art::Ptr<recob::Wire>> channelToWireMap;
1778 
1779  // Go through the list of input sources and fill out the map
1780  for (const auto& inputTag : m_hitFinderTagVec) {
1782  evt.getValidHandle<std::vector<recob::Hit>>(inputTag);
1783 
1784  art::FindOneP<recob::Wire> hitToWireAssns(hitHandle, evt, inputTag);
1785 
1786  if (hitToWireAssns.isValid()) {
1787  for (size_t wireIdx = 0; wireIdx < hitToWireAssns.size(); wireIdx++) {
1788  art::Ptr<recob::Wire> wire = hitToWireAssns.at(wireIdx);
1789 
1790  channelToWireMap[wire->Channel()] = wire;
1791  }
1792  }
1793  }
1794 
1795  // Now fill the container
1796  for (const auto& hitPtrPair : recobHitPtrMap) {
1797  raw::ChannelID_t channel = hitPtrPair.first->Channel();
1798 
1799  std::unordered_map<raw::ChannelID_t, art::Ptr<recob::Wire>>::iterator chanWireItr =
1800  channelToWireMap.find(channel);
1801 
1802  if (!(chanWireItr != channelToWireMap.end())) { continue; }
1803 
1804  wireAssns.addSingle(chanWireItr->second, hitPtrPair.second);
1805  }
1806  }
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
Definition: Wire.h:223
std::vector< art::InputTag > m_hitFinderTagVec
Data members to follow.
const_iterator end() const
Definition: Assns.h:514
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
void addSingle(Ptr< left_t > const &left, Ptr< right_t > const &right, data_t const &data)
Definition: Assns.h:549
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
Definition: fwd.h:26
geo::WireID lar_cluster3d::SnippetHit3DBuilder::NearestWireID ( const Eigen::Vector3f &  position,
const geo::WireID wireID 
) const
private

Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range.

Definition at line 1473 of file SnippetHit3DBuilder_tool.cc.

References geo::WireID::asPlaneID(), geo::TPCGeo::Length(), m_geometry, m_wireReadoutGeom, geo::WireReadoutGeom::Nwires(), geo::WireReadoutGeom::Plane(), geo::vect::toPoint(), geo::GeometryCore::TPC(), geo::WireID::Wire, and geo::PlaneGeo::WireCoordinate().

Referenced by makeDeadChannelPair().

1475  {
1476  geo::WireID wireID = wireIDIn;
1477 
1478  // Embed the call to the geometry's services nearest wire id method in a try-catch block
1479  try {
1480  // Switch from NearestWireID to this method to avoid the roundoff error issues...
1481  double distanceToWire =
1482  m_wireReadoutGeom->Plane(wireIDIn).WireCoordinate(geo::vect::toPoint(position.data()));
1483 
1484  wireID.Wire = int(distanceToWire);
1485  }
1486  catch (std::exception& exc) {
1487  // This can happen, almost always because the coordinates are **just** out of range
1488  mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - "
1489  << exc.what() << std::endl;
1490 
1491  // Assume extremum for wire number depending on z coordinate
1492  if (position[2] < 0.5 * m_geometry->TPC().Length())
1493  wireID.Wire = 0;
1494  else
1495  wireID.Wire = m_wireReadoutGeom->Nwires(wireIDIn.asPlaneID()) - 1;
1496  }
1497 
1498  return wireID;
1499  }
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
double WireCoordinate(Point_t const &point) const
Returns the coordinate of the point on the plane, in wire units.
Definition: PlaneGeo.h:704
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:110
Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
const geo::WireReadoutGeom * m_wireReadoutGeom
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:448
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void lar_cluster3d::SnippetHit3DBuilder::produces ( art::ProducesCollector collector)
overridevirtual

Each algorithm may have different objects it wants "produced" so use this to let the top level producer module "know" what it is outputting.

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 371 of file SnippetHit3DBuilder_tool.cc.

References art::ProducesCollector::produces().

372  {
373  collector.produces<std::vector<recob::Hit>>();
376  }
void produces(std::string const &instanceName={}, Persistable const persistable=Persistable::Yes)
bool lar_cluster3d::SnippetHit3DBuilder::WireIDsIntersect ( const geo::WireID wireID0,
const geo::WireID wireID1,
geo::WireIDIntersection widIntersection 
) const
private

function to detemine if two wires "intersect" (in the 2D sense)

Definition at line 1256 of file SnippetHit3DBuilder_tool.cc.

References util::abs(), closestApproach(), geo::CryostatID::Cryostat, geo::WireGeo::Direction(), geo::WireGeo::GetCenter(), geo::WireGeo::HalfL(), m_wireReadoutGeom, geo::PlaneID::Plane, geo::TPCID::TPC, geo::WireReadoutGeom::Wire(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by makeHitPair().

1259  {
1260  bool success(false);
1261 
1262  // Do quick check that things are in the same logical TPC
1263  if (wireID0.Cryostat != wireID1.Cryostat || wireID0.TPC != wireID1.TPC ||
1264  wireID0.Plane == wireID1.Plane)
1265  return success;
1266 
1267  // Recover wire geometry information for each wire
1268  const geo::WireGeo& wireGeo0 = m_wireReadoutGeom->Wire(wireID0);
1269  const geo::WireGeo& wireGeo1 = m_wireReadoutGeom->Wire(wireID1);
1270 
1271  // Get wire position and direction for first wire
1272  auto wirePosArr = wireGeo0.GetCenter();
1273 
1274  Eigen::Vector3f wirePos0(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
1275  Eigen::Vector3f wireDir0(
1276  wireGeo0.Direction().X(), wireGeo0.Direction().Y(), wireGeo0.Direction().Z());
1277 
1278  // And now the second one
1279  wirePosArr = wireGeo1.GetCenter();
1280 
1281  Eigen::Vector3f wirePos1(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
1282  Eigen::Vector3f wireDir1(
1283  wireGeo1.Direction().X(), wireGeo1.Direction().Y(), wireGeo1.Direction().Z());
1284 
1285  // Get the distance of closest approach
1286  float arcLen0;
1287  float arcLen1;
1288 
1289  if (closestApproach(wirePos0, wireDir0, wirePos1, wireDir1, arcLen0, arcLen1)) {
1290  // Now check that arc lengths are within range
1291  if (std::abs(arcLen0) < wireGeo0.HalfL() && std::abs(arcLen1) < wireGeo1.HalfL()) {
1292  Eigen::Vector3f poca0 = wirePos0 + arcLen0 * wireDir0;
1293 
1294  widIntersection.y = poca0[1];
1295  widIntersection.z = poca0[2];
1296 
1297  success = true;
1298  }
1299  }
1300 
1301  return success;
1302  }
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:112
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:219
double z
z position of intersection
Definition: geo_types.h:584
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
WireGeo const & Wire(WireID const &wireid) const
Returns the specified wire.
Vector_t Direction() const
Definition: WireGeo.h:287
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
const geo::WireReadoutGeom * m_wireReadoutGeom
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:170
double y
y position of intersection
Definition: geo_types.h:583
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
float closestApproach(const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, float &, float &) const
function to compute the distance of closest approach and the arc length to the points of closest appr...

Member Data Documentation

const lariov::ChannelStatusProvider* lar_cluster3d::SnippetHit3DBuilder::m_channelFilter
private

Definition at line 314 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec().

ChannelStatusByPlaneVec lar_cluster3d::SnippetHit3DBuilder::m_channelStatus
mutableprivate

Definition at line 307 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), and makeDeadChannelPair().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_chiSquare3DVec
mutableprivate

Definition at line 289 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

Hit2DList lar_cluster3d::SnippetHit3DBuilder::m_clusterHit2DMasterList
mutableprivate
float lar_cluster3d::SnippetHit3DBuilder::m_deltaPeakTimeSig
private

Definition at line 268 of file SnippetHit3DBuilder_tool.cc.

Referenced by findGoodTriplets(), makeHitPair(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_deltaTimeVec
mutableprivate

Definition at line 288 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

bool lar_cluster3d::SnippetHit3DBuilder::m_enableMonitoring
private
const geo::Geometry* lar_cluster3d::SnippetHit3DBuilder::m_geometry
private
std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_hitAsymmetryVec
mutableprivate

Definition at line 300 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<art::InputTag> lar_cluster3d::SnippetHit3DBuilder::m_hitFinderTagVec
private

Data members to follow.

Definition at line 266 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), makeRawDigitAssns(), makeWireAssns(), and SnippetHit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_hitWidthSclFctr
private
std::vector<int> lar_cluster3d::SnippetHit3DBuilder::m_invalidTPCVec
private

Definition at line 273 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), and SnippetHit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_LongHitStretchFctr
private

Definition at line 270 of file SnippetHit3DBuilder_tool.cc.

Referenced by makeHitPair(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_maxDeltaPeakVec
mutableprivate

Definition at line 293 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_maxHit3DChiSquare
private

Provide ability to select hits based on "chi square".

Definition at line 276 of file SnippetHit3DBuilder_tool.cc.

Referenced by SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_maxPullVec
mutableprivate

Definition at line 290 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_maxSideVecVec
mutableprivate

Definition at line 294 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

size_t lar_cluster3d::SnippetHit3DBuilder::m_numBadChannels
mutableprivate

Definition at line 308 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), and findGoodTriplets().

bool lar_cluster3d::SnippetHit3DBuilder::m_outputHistograms
private

Take the time to create and fill some histograms for diagnostics.

Definition at line 277 of file SnippetHit3DBuilder_tool.cc.

Referenced by Hit3DBuilder(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_overlapFractionVec
mutableprivate

Definition at line 291 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_overlapRangeVec
mutableprivate

Definition at line 292 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_pairWireDistVec
mutableprivate

Definition at line 295 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_PHLowSelection
private

Definition at line 272 of file SnippetHit3DBuilder_tool.cc.

Referenced by findGoodHitPairs(), and SnippetHit3DBuilder().

PlaneToSnippetHitMap lar_cluster3d::SnippetHit3DBuilder::m_planeToSnippetHitMap
mutableprivate

Definition at line 304 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildHit3D(), CollectArtHits(), and Hit3DBuilder().

PlaneToWireToHitSetMap lar_cluster3d::SnippetHit3DBuilder::m_planeToWireToHitSetMap
mutableprivate

Definition at line 305 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), and Hit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_pulseHeightFrac
private

Definition at line 271 of file SnippetHit3DBuilder_tool.cc.

Referenced by findGoodHitPairs(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_qualityMetricVec
mutableprivate

Definition at line 298 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_rangeNumSig
private

Definition at line 269 of file SnippetHit3DBuilder_tool.cc.

Referenced by makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_smallChargeDiffVec
mutableprivate

Definition at line 296 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<int> lar_cluster3d::SnippetHit3DBuilder::m_smallIndexVec
mutableprivate

Definition at line 297 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_spacePointChargeVec
mutableprivate

Definition at line 299 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), makeHitTriplet(), and SnippetHit3DBuilder().

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_timeVector
mutableprivate
TTree* lar_cluster3d::SnippetHit3DBuilder::m_tupleTree
private

output analysis tree

Definition at line 286 of file SnippetHit3DBuilder_tool.cc.

Referenced by Hit3DBuilder(), and SnippetHit3DBuilder().

bool lar_cluster3d::SnippetHit3DBuilder::m_weHaveAllBeenHereBefore = false
mutableprivate

Definition at line 310 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits().

float lar_cluster3d::SnippetHit3DBuilder::m_wirePitch[3]
private

Definition at line 280 of file SnippetHit3DBuilder_tool.cc.

Referenced by makeHitTriplet(), and SnippetHit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_wirePitchScaleFactor
private

Scaling factor to determine max distance allowed between candidate pairs.

Definition at line 275 of file SnippetHit3DBuilder_tool.cc.

Referenced by makeHitTriplet(), and SnippetHit3DBuilder().

const geo::WireReadoutGeom* lar_cluster3d::SnippetHit3DBuilder::m_wireReadoutGeom
private
float lar_cluster3d::SnippetHit3DBuilder::m_zPosOffset
private

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