All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SnippetHit3DBuilderICARUS_tool.cc
Go to the documentation of this file.
1 /**
2  * @file SnippetHit3DBuilderICARUS_tool.cc
3  *
4  * @brief This tool provides "standard" 3D hits built (by this tool) from 2D hits
5  *
6  */
7 
8 // Framework Includes
9 #include "art/Framework/Core/ProducesCollector.h"
10 #include "art/Framework/Principal/Event.h"
11 #include "art/Framework/Principal/Handle.h"
12 #include "art/Framework/Services/Registry/ServiceHandle.h"
13 #include "art/Persistency/Common/PtrMaker.h"
14 #include "art/Utilities/ToolMacros.h"
15 #include "art_root_io/TFileDirectory.h"
16 #include "art_root_io/TFileService.h"
17 #include "canvas/Persistency/Common/Assns.h"
18 #include "canvas/Persistency/Common/FindOneP.h"
19 #include "canvas/Persistency/Common/Ptr.h"
20 #include "canvas/Utilities/InputTag.h"
21 #include "cetlib/cpu_timer.h"
22 #include "fhiclcpp/ParameterSet.h"
23 #include "messagefacility/MessageLogger/MessageLogger.h"
24 
25 // LArSoft includes
26 #include "larcore/Geometry/Geometry.h"
27 #include "lardata/ArtDataHelper/HitCreator.h"
28 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
29 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
30 #include "lardataobj/RecoBase/Hit.h"
31 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
32 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
33 #include "larreco/RecoAlg/Cluster3DAlgs/IHit3DBuilder.h"
34 
35 // Eigen
36 #include <Eigen/Core>
37 
38 // std includes
39 #include <string>
40 #include <iostream>
41 #include <memory>
42 #include <numeric> // std::accumulate
43 
44 // Ack!
45 #include "TH1F.h"
46 #include "TTree.h"
47 
48 //------------------------------------------------------------------------------------------------------------------------------------------
49 // implementation follows
50 
51 namespace lar_cluster3d {
52 
53 /**
54  * @brief What follows are several highly useful typedefs which we
55  * want to expose to the outside world
56  */
57 
58 // forward declaration to define an ordering function for our hit set
59 struct Hit2DSetCompare
60 {
61  bool operator() (const reco::ClusterHit2D*, const reco::ClusterHit2D*) const;
62 };
63 
64 using HitVector = std::vector<const reco::ClusterHit2D*>;
65 using HitStartEndPair = std::pair<raw::TDCtick_t,raw::TDCtick_t>;
66 using SnippetHitMap = std::map<HitStartEndPair,HitVector>;
67 using PlaneToSnippetHitMap = std::map<geo::PlaneID, SnippetHitMap>;
68 using TPCToPlaneToSnippetHitMap = std::map<geo::TPCID, PlaneToSnippetHitMap>;
69 using Hit2DList = std::list<reco::ClusterHit2D>;
70 using Hit2DSet = std::set<const reco::ClusterHit2D*, Hit2DSetCompare>;
71 using WireToHitSetMap = std::map<unsigned int, Hit2DSet>;
72 using PlaneToWireToHitSetMap = std::map<geo::PlaneID, WireToHitSetMap>;
73 using TPCToPlaneToWireToHitSetMap = std::map<geo::TPCID, PlaneToWireToHitSetMap>;
74 using HitVectorMap = std::map<size_t, HitVector>;
75 using SnippetHitMapItrPair = std::pair<SnippetHitMap::iterator,SnippetHitMap::iterator>;
76 using PlaneSnippetHitMapItrPairVec = std::vector<SnippetHitMapItrPair>;
77 
78 /**
79  * @brief SnippetHit3DBuilderICARUS class definiton
80  */
81 class SnippetHit3DBuilderICARUS : virtual public IHit3DBuilder
82 {
83 public:
84  /**
85  * @brief Constructor
86  *
87  * @param pset
88  */
89  explicit SnippetHit3DBuilderICARUS(fhicl::ParameterSet const &pset);
90 
91  /**
92  * @brief Each algorithm may have different objects it wants "produced" so use this to
93  * let the top level producer module "know" what it is outputting
94  */
95  virtual void produces(art::ProducesCollector&) override;
96 
97  virtual void configure(const fhicl::ParameterSet&) override;
98 
99  /**
100  * @brief Given a set of recob hits, run DBscan to form 3D clusters
101  *
102  * @param hitPairList The input list of 3D hits to run clustering on
103  * @param clusterParametersList A list of cluster objects (parameters from associated hits)
104  */
105  virtual void Hit3DBuilder(art::Event&, reco::HitPairList&, RecobHitToPtrMap&) override;
106 
107  /**
108  * @brief If monitoring, recover the time to execute a particular function
109  */
110  virtual float getTimeToExecute(IHit3DBuilder::TimeValues index) const override {return m_timeVector[index];}
111 
112 private:
113 
114  /**
115  * @brief Extract the ART hits and the ART hit-particle relationships
116  *
117  * @param evt - the ART event
118  */
119  void CollectArtHits(const art::Event& evt) const;
120 
121  /**
122  * @brief Given the ClusterHit2D objects, build the HitPairMap
123  */
124  void BuildHit3D(reco::HitPairList& hitPairList) const;
125 
126  /**
127  * @brief Create a new 2D hit collection from hits associated to 3D space points
128  */
129  void CreateNewRecobHitCollection(art::Event&, reco::HitPairList&, std::vector<recob::Hit>&, RecobHitToPtrMap&);
130 
131  /**
132  * @brief Create recob::Wire to recob::Hit associations
133  */
134  void makeWireAssns(const art::Event&, art::Assns<recob::Wire, recob::Hit>&, RecobHitToPtrMap&) const;
135 
136  /**
137  * @brief Create raw::RawDigit to recob::Hit associations
138  */
139  void makeRawDigitAssns(const art::Event&, art::Assns<raw::RawDigit, recob::Hit>&, RecobHitToPtrMap&) const;
140 
141  /**
142  * @brief Given the ClusterHit2D objects, build the HitPairMap
143  */
144  size_t BuildHitPairMap(PlaneToSnippetHitMap& planeToHitVectorMap, reco::HitPairList& hitPairList) const;
145 
146  /**
147  * @brief Given the ClusterHit2D objects, build the HitPairMap
148  */
149  size_t BuildHitPairMapByTPC(PlaneSnippetHitMapItrPairVec& planeSnippetHitMapItrPairVec, reco::HitPairList& hitPairList) const;
150 
151  /**
152  * @brief This builds a list of candidate hit pairs from lists of hits on two planes
153  */
154  using HitMatchTriplet = std::tuple<const reco::ClusterHit2D*,const reco::ClusterHit2D*,const reco::ClusterHit3D>;
155  using HitMatchTripletVec = std::vector<HitMatchTriplet>;
156  using HitMatchTripletVecMap = std::map<geo::WireID,HitMatchTripletVec>;
157 
158  int findGoodHitPairs(SnippetHitMap::iterator&, SnippetHitMap::iterator&, SnippetHitMap::iterator&, HitMatchTripletVecMap&) const;
159 
160  /**
161  * @brief This algorithm takes lists of hit pairs and finds good triplets
162  */
163  void findGoodTriplets(HitMatchTripletVecMap&, HitMatchTripletVecMap&, reco::HitPairList&, bool = false) const;
164 
165  /**
166  * @brief This will look at storing pair "orphans" where the 2D hits are otherwise unused
167  */
168 
169  int saveOrphanPairs(HitMatchTripletVecMap&, reco::HitPairList&) const;
170 
171  /**
172  * @brief Make a HitPair object by checking two hits
173  */
174  bool makeHitPair(reco::ClusterHit3D& pairOut,
175  const reco::ClusterHit2D* hit1,
176  const reco::ClusterHit2D* hit2,
177  float hitWidthSclFctr = 1.,
178  size_t hitPairCntr = 0) const;
179 
180  /**
181  * @brief Make a 3D HitPair object by checking two hits
182  */
183  bool makeHitTriplet(reco::ClusterHit3D& pairOut,
184  const reco::ClusterHit3D& pairIn,
185  const reco::ClusterHit2D* hit2) const;
186 
187  /**
188  * @brief Make a 3D HitPair object from a valid pair and a dead channel in the missing plane
189  */
190  bool makeDeadChannelPair(reco::ClusterHit3D& pairOut, const reco::ClusterHit3D& pair, size_t maxStatus = 4, size_t minStatus = 0, float minOverlap=0.2) const;
191 
192  /**
193  * @brief function to detemine if two wires "intersect" (in the 2D sense)
194  */
195 
196  bool WireIDsIntersect(const geo::WireID&, const geo::WireID&, geo::WireIDIntersection&) const;
197 
198  /**
199  * @brief function to compute the distance of closest approach and the arc length to the points of closest approach
200  */
201  float closestApproach(const Eigen::Vector3f&, const Eigen::Vector3f&, const Eigen::Vector3f&, const Eigen::Vector3f&, float&, float&) const;
202 
203  /**
204  * @brief A utility routine for finding a 2D hit closest in time to the given pair
205  */
206  const reco::ClusterHit2D* FindBestMatchingHit(const Hit2DSet& hit2DSet, const reco::ClusterHit3D& pair, float pairDeltaTimeLimits) const;
207 
208  /**
209  * @brief A utility routine for returning the number of 2D hits from the list in a given range
210  */
211  int FindNumberInRange(const Hit2DSet& hit2DSet, const reco::ClusterHit3D& pair, float range) const;
212 
213  /**
214  * @brief Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range
215  */
216  geo::WireID NearestWireID(const Eigen::Vector3f& position, const geo::WireID& wireID) const;
217 
218  /**
219  * @brief Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range
220  */
221  float DistanceFromPointToHitWire(const Eigen::Vector3f& position, const geo::WireID& wireID) const;
222 
223  /**
224  * @brief Create the internal channel status vector (assume will eventually be event-by-event)
225  */
226  void BuildChannelStatusVec(PlaneToWireToHitSetMap& planeToWiretoHitSetMap) const;
227 
228  /**
229  * @brief Perform charge integration between limits
230  */
231  float chargeIntegral(float,float,float,float,int,int) const;
232 
233  /**
234  * @brief define data structure for keeping track of channel status
235  */
236  using ChannelStatusVec = std::vector<size_t>;
237  using ChannelStatusByPlaneVec = std::vector<ChannelStatusVec>;
238 
239  /**
240  * @brief clear the tuple vectors before processing next event
241  */
242  void clear();
243 
244  /**
245  * @brief Data members to follow
246  */
247  using TickCorrectionArray = std::vector<std::vector<std::vector<float>>>;
248 
249  std::vector<art::InputTag> m_hitFinderTagVec;
250  float m_hitWidthSclFctr;
251  float m_deltaPeakTimeSig;
252  float m_rangeNumSig;
253  float m_LongHitStretchFctr;
254  float m_pulseHeightFrac;
255  float m_PHLowSelection;
256  float m_minPHFor2HitPoints; ///< Set a minimum pulse height for 2 hit space point candidates
257  std::vector<int> m_invalidTPCVec;
258  float m_wirePitchScaleFactor; ///< Scaling factor to determine max distance allowed between candidate pairs
259  float m_maxHit3DChiSquare; ///< Provide ability to select hits based on "chi square"
260  bool m_saveMythicalPoints; ///< Should we save valid 2 hit space points?
261  float m_maxMythicalChiSquare; ///< Selection cut on mythical points
262  bool m_useT0Offsets; ///< If true then we will use the LArSoft interplane offsets
263  bool m_outputHistograms; ///< Take the time to create and fill some histograms for diagnostics
264 
265  bool m_enableMonitoring; ///<
266  float m_wirePitch[3];
267  mutable std::vector<float> m_timeVector; ///<
268 
269  float m_zPosOffset;
270 
271  using PlaneToT0OffsetMap = std::map<geo::PlaneID,float>;
272 
273  PlaneToT0OffsetMap m_PlaneToT0OffsetMap;
274 
275  // Define some basic histograms
276  TTree* m_tupleTree; ///< output analysis tree
277 
278  mutable std::vector<float> m_deltaPeakTimePlane0Vec;
279  mutable std::vector<float> m_deltaPeakSigmaPlane0Vec;
280  mutable std::vector<float> m_deltaPeakTimePlane1Vec;
281  mutable std::vector<float> m_deltaPeakSigmaPlane1Vec;
282  mutable std::vector<float> m_deltaPeakTimePlane2Vec;
283  mutable std::vector<float> m_deltaPeakSigmaPlane2Vec;
284 
285  mutable std::vector<float> m_deltaTimeVec;
286  mutable std::vector<float> m_deltaTime0Vec;
287  mutable std::vector<float> m_deltaSigma0Vec;
288  mutable std::vector<float> m_deltaTime1Vec;
289  mutable std::vector<float> m_deltaSigma1Vec;
290  mutable std::vector<float> m_deltaTime2Vec;
291  mutable std::vector<float> m_deltaSigma2Vec;
292  mutable std::vector<float> m_chiSquare3DVec;
293  mutable std::vector<float> m_maxPullVec;
294  mutable std::vector<float> m_overlapFractionVec;
295  mutable std::vector<float> m_overlapRangeVec;
296  mutable std::vector<float> m_maxDeltaPeakVec;
297  mutable std::vector<float> m_maxSideVecVec;
298  mutable std::vector<float> m_pairWireDistVec;
299  mutable std::vector<float> m_smallChargeDiffVec;
300  mutable std::vector<int> m_smallIndexVec;
301  mutable std::vector<float> m_qualityMetricVec;
302  mutable std::vector<float> m_spacePointChargeVec;
303  mutable std::vector<float> m_hitAsymmetryVec;
304  mutable std::vector<float> m_2hit1stPHVec;
305  mutable std::vector<float> m_2hit2ndPHVec;
306  mutable std::vector<float> m_2hitDeltaPHVec;
307  mutable std::vector<float> m_2hitSumPHVec;
308 
309  // Get instances of the primary data structures needed
310  mutable Hit2DList m_clusterHit2DMasterList;
311  mutable PlaneToSnippetHitMap m_planeToSnippetHitMap;
312  mutable PlaneToWireToHitSetMap m_planeToWireToHitSetMap;
313 
314 
315  mutable ChannelStatusByPlaneVec m_channelStatus;
316  mutable size_t m_numBadChannels;
317 
318  mutable bool m_weHaveAllBeenHereBefore = false;
319 
320  const geo::Geometry* m_geometry; //< pointer to the Geometry service
321  const lariov::ChannelStatusProvider* m_channelFilter;
322 };
323 
324 SnippetHit3DBuilderICARUS::SnippetHit3DBuilderICARUS(fhicl::ParameterSet const &pset) :
325  m_channelFilter(&art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider())
326 
327 {
328  this->configure(pset);
329 }
330 
331 //------------------------------------------------------------------------------------------------------------------------------------------
332 
333 void SnippetHit3DBuilderICARUS::produces(art::ProducesCollector& collector)
334 {
335  collector.produces< std::vector<recob::Hit>>();
336  collector.produces< art::Assns<recob::Wire, recob::Hit>>();
337  collector.produces< art::Assns<raw::RawDigit, recob::Hit>>();
338 }
339 
340 //------------------------------------------------------------------------------------------------------------------------------------------
341 
342 void SnippetHit3DBuilderICARUS::configure(fhicl::ParameterSet const &pset)
343 {
344  m_hitFinderTagVec = pset.get<std::vector<art::InputTag>>("HitFinderTagVec", {"gaushit"});
345  m_enableMonitoring = pset.get<bool >("EnableMonitoring", true);
346  m_hitWidthSclFctr = pset.get<float >("HitWidthScaleFactor", 6. );
347  m_rangeNumSig = pset.get<float >("RangeNumSigma", 3. );
348  m_LongHitStretchFctr = pset.get<float >("LongHitsStretchFactor", 1.5 );
349  m_pulseHeightFrac = pset.get<float >("PulseHeightFraction", 0.5 );
350  m_PHLowSelection = pset.get<float >("PHLowSelection", 20. );
351  m_minPHFor2HitPoints = pset.get<float >("MinPHFor2HitPoints", 15. );
352  m_deltaPeakTimeSig = pset.get<float >("DeltaPeakTimeSig", 1.7 );
353  m_zPosOffset = pset.get<float >("ZPosOffset", 0.0 );
354  m_invalidTPCVec = pset.get<std::vector<int> >("InvalidTPCVec", std::vector<int>());
355  m_wirePitchScaleFactor = pset.get<float >("WirePitchScaleFactor", 1.9 );
356  m_maxHit3DChiSquare = pset.get<float >("MaxHitChiSquare", 6.0 );
357  m_saveMythicalPoints = pset.get<bool >("SaveMythicalPoints", true);
358  m_maxMythicalChiSquare = pset.get<float >("MaxMythicalChiSquare", 10.);
359  m_useT0Offsets = pset.get<bool >("UseT0Offsets", true);
360  m_outputHistograms = pset.get<bool >("OutputHistograms", false);
361 
362  m_geometry = art::ServiceHandle<geo::Geometry const>{}.get();
363 
364  // Returns the wire pitch per plane assuming they will be the same for all TPCs
365  m_wirePitch[0] = m_geometry->WirePitch(0);
366  m_wirePitch[1] = m_geometry->WirePitch(1);
367  m_wirePitch[2] = m_geometry->WirePitch(2);
368 
369  // Access ART's TFileService, which will handle creating and writing
370  // histograms and n-tuples for us.
371  if (m_outputHistograms)
372  {
373  art::ServiceHandle<art::TFileService> tfs;
374 
375  m_tupleTree = tfs->make<TTree>("Hit3DBuilderTree", "Tree by SnippetHit3DBuilderICARUS");
376 
377  clear();
378 
379  m_tupleTree->Branch("DeltaPeakTimePair0", "std::vector<float>", &m_deltaPeakTimePlane0Vec);
380  m_tupleTree->Branch("DeltaPeakSigmaPair0", "std::vector<float>", &m_deltaPeakSigmaPlane0Vec);
381  m_tupleTree->Branch("DeltaPeakTimePair1", "std::vector<float>", &m_deltaPeakTimePlane1Vec);
382  m_tupleTree->Branch("DeltaPeakSigmaPair1", "std::vector<float>", &m_deltaPeakSigmaPlane1Vec);
383  m_tupleTree->Branch("DeltaPeakTimePair2", "std::vector<float>", &m_deltaPeakTimePlane2Vec);
384  m_tupleTree->Branch("DeltaPeakSigmaPair2", "std::vector<float>", &m_deltaPeakSigmaPlane2Vec);
385 
386  m_tupleTree->Branch("DeltaTime2D", "std::vector<float>", &m_deltaTimeVec);
387  m_tupleTree->Branch("DeltaTime2D0", "std::vector<float>", &m_deltaTime0Vec);
388  m_tupleTree->Branch("DeltaSigma2D0", "std::vector<float>", &m_deltaSigma0Vec);
389  m_tupleTree->Branch("DeltaTime2D1", "std::vector<float>", &m_deltaTime1Vec);
390  m_tupleTree->Branch("DeltaSigma2D1", "std::vector<float>", &m_deltaSigma1Vec);
391  m_tupleTree->Branch("DeltaTime2D2", "std::vector<float>", &m_deltaTime2Vec);
392  m_tupleTree->Branch("DeltaSigma2D2", "std::vector<float>", &m_deltaSigma2Vec);
393  m_tupleTree->Branch("ChiSquare3D", "std::vector<float>", &m_chiSquare3DVec);
394  m_tupleTree->Branch("MaxPullValue", "std::vector<float>", &m_maxPullVec);
395  m_tupleTree->Branch("OverlapFraction", "std::vector<float>", &m_overlapFractionVec);
396  m_tupleTree->Branch("OverlapRange", "std::vector<float>", &m_overlapRangeVec);
397  m_tupleTree->Branch("MaxDeltaPeak", "std::vector<float>", &m_maxDeltaPeakVec);
398  m_tupleTree->Branch("MaxSideVec", "std::vector<float>", &m_maxSideVecVec);
399  m_tupleTree->Branch("PairWireDistVec", "std::vector<float>", &m_pairWireDistVec);
400  m_tupleTree->Branch("SmallChargeDiff", "std::vector<float>", &m_smallChargeDiffVec);
401  m_tupleTree->Branch("SmallChargeIdx", "std::vector<int>", &m_smallIndexVec);
402  m_tupleTree->Branch("QualityMetric", "std::vector<float>", &m_qualityMetricVec);
403  m_tupleTree->Branch("SPCharge", "std::vector<float>", &m_spacePointChargeVec);
404  m_tupleTree->Branch("HitAsymmetry", "std::vector<float>", &m_hitAsymmetryVec);
405 
406  m_tupleTree->Branch("2hit1stPH", "std::vector<float>", &m_2hit1stPHVec);
407  m_tupleTree->Branch("2hit2ndPH", "std::vector<float>", &m_2hit2ndPHVec);
408  m_tupleTree->Branch("2hitDeltaPH", "std::vector<float>", &m_2hitDeltaPHVec);
409  m_tupleTree->Branch("2hitSumPH", "std::vector<float>", &m_2hitSumPHVec);
410 
411  }
412 
413  return;
414 }
415 
416 void SnippetHit3DBuilderICARUS::clear()
417 {
418  m_deltaPeakTimePlane0Vec.clear();
419  m_deltaPeakSigmaPlane0Vec.clear();
420  m_deltaPeakTimePlane1Vec.clear();
421  m_deltaPeakSigmaPlane1Vec.clear();
422  m_deltaPeakTimePlane1Vec.clear();
423  m_deltaPeakSigmaPlane1Vec.clear();
424 
425  m_deltaTimeVec.clear();
426  m_deltaTime0Vec.clear();
427  m_deltaSigma0Vec.clear();
428  m_deltaTime1Vec.clear();
429  m_deltaSigma1Vec.clear();
430  m_deltaTime2Vec.clear();
431  m_deltaSigma2Vec.clear();
432  m_chiSquare3DVec.clear();
433  m_maxPullVec.clear();
434  m_overlapFractionVec.clear();
435  m_overlapRangeVec.clear();
436  m_maxDeltaPeakVec.clear();
437  m_maxSideVecVec.clear();
438  m_pairWireDistVec.clear();
439  m_smallChargeDiffVec.clear();
440  m_smallIndexVec.clear();
441  m_qualityMetricVec.clear();
442  m_spacePointChargeVec.clear();
443  m_hitAsymmetryVec.clear();
444 
445  m_2hit1stPHVec.clear();
446  m_2hit2ndPHVec.clear();
447  m_2hitDeltaPHVec.clear();
448  m_2hitSumPHVec.clear();
449 
450  return;
451 }
452 
453 void SnippetHit3DBuilderICARUS::BuildChannelStatusVec(PlaneToWireToHitSetMap& planeToWireToHitSetMap) const
454 {
455  // This is called each event, clear out the previous version and start over
456  m_channelStatus.clear();
457 
458  m_numBadChannels = 0;
459  m_channelStatus.resize(m_geometry->Nplanes());
460 
461  // Loop through views/planes to set the wire length vectors
462  for(size_t idx = 0; idx < m_channelStatus.size(); idx++)
463  {
464  m_channelStatus[idx] = ChannelStatusVec(m_geometry->Nwires(idx), 5);
465  }
466 
467  // Loop through the channels and mark those that are "bad"
468  for(size_t channel = 0; channel < m_geometry->Nchannels(); channel++)
469  {
470  if( !m_channelFilter->IsGood(channel))
471  {
472  std::vector<geo::WireID> wireIDVec = m_geometry->ChannelToWire(channel);
473  geo::WireID wireID = wireIDVec[0];
474  lariov::ChannelStatusProvider::Status_t chanStat = m_channelFilter->Status(channel);
475 
476  m_channelStatus[wireID.Plane][wireID.Wire] = chanStat;
477  m_numBadChannels++;
478  }
479  }
480 
481  return;
482 }
483 
484 
485 bool SetPeakHitPairIteratorOrder(const reco::HitPairList::iterator& left, const reco::HitPairList::iterator& right)
486 {
487  return (*left).getAvePeakTime() < (*right).getAvePeakTime();
488 }
489 
490 struct HitPairClusterOrder
491 {
492  bool operator()(const reco::HitPairClusterMap::iterator& left, const reco::HitPairClusterMap::iterator& right)
493  {
494  // Watch out for the case where two clusters can have the same number of hits!
495  if (left->second.size() == right->second.size())
496  return left->first < right->first;
497 
498  return left->second.size() > right->second.size();
499  }
500 };
501 
502 void SnippetHit3DBuilderICARUS::Hit3DBuilder(art::Event& evt, reco::HitPairList& hitPairList, RecobHitToPtrMap& clusterHitToArtPtrMap)
503 {
504  // Clear the internal data structures
505  m_clusterHit2DMasterList.clear();
506  m_planeToSnippetHitMap.clear();
507  m_planeToWireToHitSetMap.clear();
508 
509  // Do the one time initialization of the tick offsets.
510  if (m_PlaneToT0OffsetMap.empty())
511  {
512  // Need the detector properties which needs the clocks
513  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
514  auto const det_prop = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
515 
516  // Initialize the plane to hit vector map
517  for(size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++)
518  {
519  for(size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++)
520  {
521  for(size_t planeIdx = 0; planeIdx < m_geometry->Nplanes(); planeIdx++)
522  {
523  geo::PlaneID planeID(cryoIdx,tpcIdx,planeIdx);
524 
525  if (m_useT0Offsets) m_PlaneToT0OffsetMap[planeID] = det_prop.GetXTicksOffset(planeID) - det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,0));
526  else m_PlaneToT0OffsetMap[planeID] = 0.;
527  }
528  }
529  }
530  }
531 
532  m_timeVector.resize(NUMTIMEVALUES, 0.);
533 
534  // Get a hit refiner
535  std::unique_ptr<std::vector<recob::Hit>> outputHitPtrVec(new std::vector<recob::Hit>);
536 
537  // Recover the 2D hits and then organize them into data structures which will be used in the
538  // DBscan algorithm for building the 3D clusters
539  this->CollectArtHits(evt);
540 
541  // If there are no hits in our view/wire data structure then do not proceed with the full analysis
542  if (!m_planeToWireToHitSetMap.empty())
543  {
544  // Call the algorithm that builds 3D hits
545  this->BuildHit3D(hitPairList);
546 
547  // If we built 3D points then attempt to output a new hit list as well
548  if (!hitPairList.empty())
549  CreateNewRecobHitCollection(evt, hitPairList, *outputHitPtrVec, clusterHitToArtPtrMap);
550  }
551 
552  // Set up to make the associations (if desired)
553  /// Associations with wires.
554  std::unique_ptr<art::Assns<recob::Wire, recob::Hit>> wireAssns(new art::Assns<recob::Wire, recob::Hit>);
555 
556  makeWireAssns(evt, *wireAssns, clusterHitToArtPtrMap);
557 
558  /// Associations with raw digits.
559  std::unique_ptr<art::Assns<raw::RawDigit, recob::Hit>> rawDigitAssns(new art::Assns<raw::RawDigit, recob::Hit>);
560 
561  makeRawDigitAssns(evt, *rawDigitAssns, clusterHitToArtPtrMap);
562 
563  // Move everything into the event
564  evt.put(std::move(outputHitPtrVec));
565  evt.put(std::move(wireAssns));
566  evt.put(std::move(rawDigitAssns));
567 
568  // Handle tree output too
569  if (m_outputHistograms)
570  {
571  m_tupleTree->Fill();
572 
573  clear();
574  }
575 
576  return;
577 }
578 
579 void SnippetHit3DBuilderICARUS::BuildHit3D(reco::HitPairList& hitPairList) const
580 {
581  /**
582  * @brief Driver for processing input 2D hits, transforming to 3D hits and building lists
583  * of associated 3D hits (candidate 3D clusters)
584  */
585  cet::cpu_timer theClockMakeHits;
586 
587  if (m_enableMonitoring) theClockMakeHits.start();
588 
589  // The first task is to take the lists of input 2D hits (a map of view to sorted lists of 2D hits)
590  // and then to build a list of 3D hits to be used in downstream processing
591  BuildChannelStatusVec(m_planeToWireToHitSetMap);
592 
593  size_t numHitPairs = BuildHitPairMap(m_planeToSnippetHitMap, hitPairList);
594 
595  if (m_enableMonitoring)
596  {
597  theClockMakeHits.stop();
598 
599  m_timeVector[BUILDTHREEDHITS] = theClockMakeHits.accumulated_real_time();
600  }
601 
602  mf::LogDebug("Cluster3D") << ">>>>> 3D hit building done, found " << numHitPairs << " 3D Hits" << std::endl;
603 
604  return;
605 }
606 
607 //------------------------------------------------------------------------------------------------------------------------------------------
608 class SetStartTimeOrder
609 {
610 public:
611  SetStartTimeOrder() {}
612 
613  bool operator()(const SnippetHitMapItrPair& left, const SnippetHitMapItrPair& right) const
614  {
615  // Special case handling, there is nothing to compare for the left or right
616  if (left.first == left.second) return false;
617  if (right.first == right.second) return true;
618 
619  // de-referencing a bunch of pairs here...
620  return left.first->first.first < right.first->first.first;
621  }
622 
623 private:
624 };
625 
626 bool SetPairStartTimeOrder(const reco::ClusterHit3D& left, const reco::ClusterHit3D& right)
627 {
628  // Sort by "modified start time" of pulse
629  return left.getAvePeakTime() - left.getSigmaPeakTime() < right.getAvePeakTime() - right.getSigmaPeakTime();
630 }
631 
632 
633 //------------------------------------------------------------------------------------------------------------------------------------------
634 
635 size_t SnippetHit3DBuilderICARUS::BuildHitPairMap(PlaneToSnippetHitMap& planeToSnippetHitMap, reco::HitPairList& hitPairList) const
636 {
637  /**
638  * @brief Given input 2D hits, build out the lists of possible 3D hits
639  *
640  * The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
641  * However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
642  * be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
643  * 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
644  * will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.
645  */
646  size_t totalNumHits(0);
647  size_t hitPairCntr(0);
648 
649  size_t nTriplets(0);
650  size_t nDeadChanHits(0);
651 
652  // Set up to loop over cryostats and tpcs...
653  for(size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++)
654  {
655  for(size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++)
656  {
657  //************************************
658  // Kludge
659 // if (!(cryoIdx == 1 && tpcIdx == 0)) continue;
660 
661  PlaneToSnippetHitMap::iterator mapItr0 = planeToSnippetHitMap.find(geo::PlaneID(cryoIdx,tpcIdx,0));
662  PlaneToSnippetHitMap::iterator mapItr1 = planeToSnippetHitMap.find(geo::PlaneID(cryoIdx,tpcIdx,1));
663  PlaneToSnippetHitMap::iterator mapItr2 = planeToSnippetHitMap.find(geo::PlaneID(cryoIdx,tpcIdx,2));
664 
665  size_t nPlanesWithHits = (mapItr0 != planeToSnippetHitMap.end() && !mapItr0->second.empty() ? 1 : 0)
666  + (mapItr1 != planeToSnippetHitMap.end() && !mapItr1->second.empty() ? 1 : 0)
667  + (mapItr2 != planeToSnippetHitMap.end() && !mapItr2->second.empty() ? 1 : 0);
668 
669  if (nPlanesWithHits < 2) continue;
670 
671  SnippetHitMap& snippetHitMap0 = mapItr0->second;
672  SnippetHitMap& snippetHitMap1 = mapItr1->second;
673  SnippetHitMap& snippetHitMap2 = mapItr2->second;
674 
675  PlaneSnippetHitMapItrPairVec hitItrVec = {SnippetHitMapItrPair(snippetHitMap0.begin(),snippetHitMap0.end()),
676  SnippetHitMapItrPair(snippetHitMap1.begin(),snippetHitMap1.end()),
677  SnippetHitMapItrPair(snippetHitMap2.begin(),snippetHitMap2.end())};
678 
679  totalNumHits += BuildHitPairMapByTPC(hitItrVec, hitPairList);
680  }
681  }
682 
683  // Return the hit pair list but sorted by z and y positions (faster traversal in next steps)
684  hitPairList.sort(SetPairStartTimeOrder);
685 
686  // Where are we?
687  mf::LogDebug("Cluster3D") << "Total number hits: " << totalNumHits << std::endl;
688  mf::LogDebug("Cluster3D") << "Created a total of " << hitPairList.size() << " hit pairs, counted: " << hitPairCntr << std::endl;
689  mf::LogDebug("Cluster3D") << "-- Triplets: " << nTriplets << ", dead channel pairs: " << nDeadChanHits << std::endl;
690 
691  return hitPairList.size();
692 }
693 
694 size_t SnippetHit3DBuilderICARUS::BuildHitPairMapByTPC(PlaneSnippetHitMapItrPairVec& snippetHitMapItrVec, reco::HitPairList& hitPairList) const
695 {
696  /**
697  * @brief Given input 2D hits, build out the lists of possible 3D hits
698  *
699  * The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
700  * However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
701  * be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
702  * 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
703  * will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.
704  */
705 
706  // Define functions to set start/end iterators in the loop below
707  auto SetStartIterator = [](SnippetHitMap::iterator startItr, SnippetHitMap::iterator endItr, float startTime)
708  {
709  while(startItr != endItr)
710  {
711  if (startItr->first.second < startTime) startItr++;
712  else break;
713  }
714  return startItr;
715  };
716 
717  auto SetEndIterator = [](SnippetHitMap::iterator lastItr, SnippetHitMap::iterator endItr, float endTime)
718  {
719  while(lastItr != endItr)
720  {
721  if (lastItr->first.first < endTime) lastItr++;
722  else break;
723  }
724  return lastItr;
725  };
726 
727  size_t nTriplets(0);
728  size_t nDeadChanHits(0);
729  size_t nOrphanPairs(0);
730 
731  //*********************************************************************************
732  // Basically, we try to loop until done...
733  while(1)
734  {
735  // Sort so that the earliest hit time will be the first element, etc.
736  std::sort(snippetHitMapItrVec.begin(),snippetHitMapItrVec.end(),SetStartTimeOrder());
737 
738  // Make sure there are still hits on at least
739  int nPlanesWithHits(0);
740 
741  for(auto& pair : snippetHitMapItrVec)
742  if (pair.first != pair.second) nPlanesWithHits++;
743 
744  if (nPlanesWithHits < 2) break;
745 
746  // End condition: no more hit snippets
747 // if (snippetHitMapItrVec.front().first == snippetHitMapItrVec.front().second) break;
748 
749  // This loop iteration's snippet iterator
750  SnippetHitMap::iterator firstSnippetItr = snippetHitMapItrVec.front().first;
751 
752  // Set iterators to insure we'll be in the overlap ranges
753  SnippetHitMap::iterator snippetHitMapItr1Start = SetStartIterator(snippetHitMapItrVec[1].first, snippetHitMapItrVec[1].second, firstSnippetItr->first.first);
754  SnippetHitMap::iterator snippetHitMapItr1End = SetEndIterator( snippetHitMapItr1Start, snippetHitMapItrVec[1].second, firstSnippetItr->first.second);
755  SnippetHitMap::iterator snippetHitMapItr2Start = SetStartIterator(snippetHitMapItrVec[2].first, snippetHitMapItrVec[2].second, firstSnippetItr->first.first);
756  SnippetHitMap::iterator snippetHitMapItr2End = SetEndIterator( snippetHitMapItr2Start, snippetHitMapItrVec[2].second, firstSnippetItr->first.second);
757 
758  // Since we'll use these many times in the internal loops, pre make the pairs for the second set of hits
759  size_t curHitListSize(hitPairList.size());
760  HitMatchTripletVecMap pair12Map;
761  HitMatchTripletVecMap pair13Map;
762 
763  size_t n12Pairs = findGoodHitPairs(firstSnippetItr, snippetHitMapItr1Start, snippetHitMapItr1End, pair12Map);
764  size_t n13Pairs = findGoodHitPairs(firstSnippetItr, snippetHitMapItr2Start, snippetHitMapItr2End, pair13Map);
765 
766  nDeadChanHits += hitPairList.size() - curHitListSize;
767  curHitListSize = hitPairList.size();
768 
769  if (n12Pairs > n13Pairs) findGoodTriplets(pair12Map, pair13Map, hitPairList);
770  else findGoodTriplets(pair13Map, pair12Map, hitPairList);
771 
772  if (m_saveMythicalPoints)
773  {
774  nOrphanPairs += saveOrphanPairs(pair12Map, hitPairList);
775  nOrphanPairs += saveOrphanPairs(pair13Map, hitPairList);
776  }
777 
778  nTriplets += hitPairList.size() - curHitListSize;
779 
780  snippetHitMapItrVec.front().first++;
781  }
782 
783  mf::LogDebug("Cluster3D") << "--> Created " << nTriplets << " triplets of which " << nOrphanPairs << " are orphans" << std::endl;
784 
785  return hitPairList.size();
786 }
787 
788 int SnippetHit3DBuilderICARUS::findGoodHitPairs(SnippetHitMap::iterator& firstSnippetItr,
789  SnippetHitMap::iterator& startItr,
790  SnippetHitMap::iterator& endItr,
791  HitMatchTripletVecMap& hitMatchMap) const
792 {
793  int numPairs(0);
794 
795  HitVector::iterator firstMaxItr = std::max_element(firstSnippetItr->second.begin(),firstSnippetItr->second.end(),[](const auto& left, const auto& right){return left->getHit()->PeakAmplitude() < right->getHit()->PeakAmplitude();});
796  float firstPHCut = firstMaxItr != firstSnippetItr->second.end() ? m_pulseHeightFrac * (*firstMaxItr)->getHit()->PeakAmplitude() : 4096.;
797 
798  // Loop through the hits on the first snippet
799  for(const auto& hit1 : firstSnippetItr->second)
800  {
801  // Let's focus on the largest hit in the chain
802  if (hit1->getHit()->DegreesOfFreedom() > 1 && hit1->getHit()->PeakAmplitude() < firstPHCut && hit1->getHit()->PeakAmplitude() < m_PHLowSelection) continue;
803 
804  // Inside loop iterator
805  SnippetHitMap::iterator secondHitItr = startItr;
806 
807  // Loop through the input secon hits and make pairs
808  while(secondHitItr != endItr)
809  {
810  HitVector::iterator secondMaxItr = std::max_element(secondHitItr->second.begin(),secondHitItr->second.end(),[](const auto& left, const auto& right){return left->getHit()->PeakAmplitude() < right->getHit()->PeakAmplitude();});
811  float secondPHCut = secondMaxItr != secondHitItr->second.end() ? m_pulseHeightFrac * (*secondMaxItr)->getHit()->PeakAmplitude() : 0.;
812 
813  for(const auto& hit2 : secondHitItr->second)
814  {
815  // Again, focus on the large hits
816  if (hit2->getHit()->DegreesOfFreedom() > 1 && hit2->getHit()->PeakAmplitude() < secondPHCut && hit2->getHit()->PeakAmplitude() < m_PHLowSelection) continue;
817 
818  reco::ClusterHit3D pair;
819 
820  // pair returned with a negative ave time is signal of failure
821  if (!makeHitPair(pair, hit1, hit2, m_hitWidthSclFctr)) continue;
822 
823  std::vector<const recob::Hit*> recobHitVec = {nullptr,nullptr,nullptr};
824 
825  recobHitVec[hit1->WireID().Plane] = hit1->getHit();
826  recobHitVec[hit2->WireID().Plane] = hit2->getHit();
827 
828  geo::WireID wireID = hit2->WireID();
829 
830  hitMatchMap[wireID].emplace_back(hit1,hit2,pair);
831 
832  numPairs++;
833  }
834 
835  secondHitItr++;
836  }
837  }
838 
839  return numPairs;
840 }
841 
842 void SnippetHit3DBuilderICARUS::findGoodTriplets(HitMatchTripletVecMap& pair12Map, HitMatchTripletVecMap& pair13Map, reco::HitPairList& hitPairList, bool tagged) const
843 {
844  // Build triplets from the two lists of hit pairs
845  if (!pair12Map.empty())
846  {
847  // temporary container for dead channel hits
848  std::vector<reco::ClusterHit3D> tempDeadChanVec;
849  reco::ClusterHit3D deadChanPair;
850 
851  // Keep track of which third plane hits have been used
852  std::map<const reco::ClusterHit3D*,bool> usedPairMap;
853 
854  // Initial population of this map with the pair13Map hits
855  for(const auto& pair13 : pair13Map)
856  {
857  for(const auto& hit2Dhit3DPair : pair13.second) usedPairMap[&std::get<2>(hit2Dhit3DPair)] = false;
858  }
859 
860  // The outer loop is over all hit pairs made from the first two plane combinations
861  for(const auto& pair12 : pair12Map)
862  {
863  if (pair12.second.empty()) continue;
864 
865  // This loop is over hit pairs that share the same first two plane wires but may have different
866  // hit times on those wires
867  for(const auto& hit2Dhit3DPair12 : pair12.second)
868  {
869  const reco::ClusterHit3D& pair1 = std::get<2>(hit2Dhit3DPair12);
870 
871  // populate the map with initial value
872  usedPairMap[&pair1] = false;
873 
874  // The simplest approach here is to loop over all possibilities and let the triplet builder weed out the weak candidates
875  for(const auto& pair13 : pair13Map)
876  {
877  if (pair13.second.empty()) continue;
878 
879  for(const auto& hit2Dhit3DPair13 : pair13.second)
880  {
881  // Protect against double counting
882  if (std::get<0>(hit2Dhit3DPair12) != std::get<0>(hit2Dhit3DPair13)) continue;
883 
884  const reco::ClusterHit2D* hit2 = std::get<1>(hit2Dhit3DPair13);
885  const reco::ClusterHit3D& pair2 = std::get<2>(hit2Dhit3DPair13);
886 
887  // If success try for the triplet
888  reco::ClusterHit3D triplet;
889 
890  if (makeHitTriplet(triplet, pair1, hit2))
891  {
892  triplet.setID(hitPairList.size());
893  hitPairList.emplace_back(triplet);
894  usedPairMap[&pair1] = true;
895  usedPairMap[&pair2] = true;
896  }
897  }
898  }
899  }
900  }
901 
902  // One more loop through the other pairs to check for sick channels
903  if (m_numBadChannels > 0)
904  {
905  for(const auto& pairMapPair : usedPairMap)
906  {
907  if (pairMapPair.second) continue;
908 
909  const reco::ClusterHit3D* pair = pairMapPair.first;
910 
911  // Here we look to see if we failed to make a triplet because the partner wire was dead/noisy/sick
912  if (makeDeadChannelPair(deadChanPair, *pair, 4, 0, 0.)) tempDeadChanVec.emplace_back(deadChanPair);
913  }
914 
915  // Handle the dead wire triplets
916  if(!tempDeadChanVec.empty())
917  {
918  // If we have many then see if we can trim down a bit by keeping those with time significance
919  if (tempDeadChanVec.size() > 1)
920  {
921  // Sort by "significance" of agreement
922  std::sort(tempDeadChanVec.begin(),tempDeadChanVec.end(),[](const auto& left, const auto& right){return left.getDeltaPeakTime()/left.getSigmaPeakTime() < right.getDeltaPeakTime()/right.getSigmaPeakTime();});
923 
924  // What is the range of "significance" from first to last?
925  float firstSig = tempDeadChanVec.front().getDeltaPeakTime() / tempDeadChanVec.front().getSigmaPeakTime();
926  float lastSig = tempDeadChanVec.back().getDeltaPeakTime() / tempDeadChanVec.back().getSigmaPeakTime();
927  float sigRange = lastSig - firstSig;
928 
929  if (lastSig > 0.5 * m_deltaPeakTimeSig && sigRange > 0.5)
930  {
931  // Declare a maximum of 1.5 * the average of the first and last pairs...
932  float maxSignificance = std::max(0.75 * (firstSig + lastSig),1.0);
933 
934  std::vector<reco::ClusterHit3D>::iterator firstBadElem = std::find_if(tempDeadChanVec.begin(),tempDeadChanVec.end(),[&maxSignificance](const auto& pair){return pair.getDeltaPeakTime()/pair.getSigmaPeakTime() > maxSignificance;});
935 
936  // But only keep the best 10?
937  if (std::distance(tempDeadChanVec.begin(),firstBadElem) > 20) firstBadElem = tempDeadChanVec.begin() + 20;
938  // Keep at least one hit...
939  else if (firstBadElem == tempDeadChanVec.begin()) firstBadElem++;
940 
941  tempDeadChanVec.resize(std::distance(tempDeadChanVec.begin(),firstBadElem));
942  }
943  }
944 
945  for(auto& pair : tempDeadChanVec)
946  {
947  pair.setID(hitPairList.size());
948  hitPairList.emplace_back(pair);
949  }
950  }
951  }
952  }
953 
954  return;
955 }
956 
957 int SnippetHit3DBuilderICARUS::saveOrphanPairs(HitMatchTripletVecMap& pairMap, reco::HitPairList& hitPairList) const
958 {
959  int curTripletCount = hitPairList.size();
960 
961  // Build triplets from the two lists of hit pairs
962  if (!pairMap.empty())
963  {
964  // Initial population of this map with the pair13Map hits
965  for(const auto& pair : pairMap)
966  {
967  if (pair.second.empty()) continue;
968 
969  // This loop is over hit pairs that share the same first two plane wires but may have different
970  // hit times on those wires
971  for(const auto& hit2Dhit3DPair : pair.second)
972  {
973  const reco::ClusterHit3D& hit3D = std::get<2>(hit2Dhit3DPair);
974 
975  // No point considering a 3D hit that has been used to make a space point already
976  if (hit3D.getStatusBits() & reco::ClusterHit3D::MADESPACEPOINT) continue;
977 
978  const reco::ClusterHit2D* hit1 = std::get<0>(hit2Dhit3DPair);
979  const reco::ClusterHit2D* hit2 = std::get<1>(hit2Dhit3DPair);
980 
981  if (m_outputHistograms)
982  {
983  m_2hit1stPHVec.emplace_back(hit1->getHit()->PeakAmplitude());
984  m_2hit2ndPHVec.emplace_back(hit2->getHit()->PeakAmplitude());
985  m_2hitDeltaPHVec.emplace_back(hit2->getHit()->PeakAmplitude() - hit1->getHit()->PeakAmplitude());
986  m_2hitSumPHVec.emplace_back(hit2->getHit()->PeakAmplitude() + hit1->getHit()->PeakAmplitude());
987  }
988 
989  // If both hits already appear in a triplet then there is no gain here so reject
990  if (hit1->getHit()->PeakAmplitude() < m_minPHFor2HitPoints || hit2->getHit()->PeakAmplitude() < m_minPHFor2HitPoints) continue;
991 
992  // Require that one of the hits is on the collection plane
993  if (hit1->WireID().Plane == 2 || hit2->WireID().Plane == 2)
994  {
995  // Allow cut on the quality of the space point
996  if (hit3D.getHitChiSquare() < m_maxMythicalChiSquare)
997  {
998  // Add to the list
999  hitPairList.emplace_back(hit3D);
1000  hitPairList.back().setID(hitPairList.size()-1);
1001  }
1002  }
1003  }
1004  }
1005  }
1006 
1007  return hitPairList.size() - curTripletCount;
1008 }
1009 
1010 bool SnippetHit3DBuilderICARUS::makeHitPair(reco::ClusterHit3D& hitPair,
1011  const reco::ClusterHit2D* hit1,
1012  const reco::ClusterHit2D* hit2,
1013  float hitWidthSclFctr,
1014  size_t hitPairCntr) const
1015 {
1016  // Assume failure
1017  bool result(false);
1018 
1019  // Start by checking time consistency since this is fastest
1020  // Wires intersect so now we can check the timing
1021  float hit1Peak = hit1->getTimeTicks();
1022  float hit1Sigma = hit1->getHit()->RMS();
1023 
1024  float hit2Peak = hit2->getTimeTicks();
1025  float hit2Sigma = hit2->getHit()->RMS();
1026 
1027  // "Long hits" are an issue... so we deal with these differently
1028  int hit1NDF = hit1->getHit()->DegreesOfFreedom();
1029  int hit2NDF = hit2->getHit()->DegreesOfFreedom();
1030 
1031  // Basically, allow the range to extend to the nearest end of the snippet
1032  if (hit1NDF < 2) hit1Sigma *= m_LongHitStretchFctr; // This sets the range to the width of the pulse
1033  if (hit2NDF < 2) hit2Sigma *= m_LongHitStretchFctr;
1034 
1035  // The "hit sigma" is the gaussian fit sigma of the hit, we need to expand a bit to allow hit overlap efficiency
1036  float hit1Width = hitWidthSclFctr * hit1Sigma;
1037  float hit2Width = hitWidthSclFctr * hit2Sigma;
1038 
1039  // Coarse check hit times are "in range"
1040  if (fabs(hit1Peak - hit2Peak) <= (hit1Width + hit2Width))
1041  {
1042  // Check to see that hit peak times are consistent with each other
1043  float hit1SigSq = std::pow(hit1Sigma,2);
1044  float hit2SigSq = std::pow(hit2Sigma,2);
1045  float deltaPeakTime = std::fabs(hit1Peak - hit2Peak);
1046  float sigmaPeakTime = std::sqrt(hit1SigSq + hit2SigSq);
1047 
1048  if (m_outputHistograms)
1049  {
1050  // brute force... sigh...
1051  int plane1 = hit1->WireID().Plane;
1052  int plane2 = hit2->WireID().Plane;
1053  int planeIdx = (plane1 + plane2) - 1; // should be 0 for 0-1, 1 for 0-2 and 2 for 1-2
1054 
1055  float deltaTicks = hit2Peak - hit1Peak;
1056 
1057  if (plane1 > plane2) deltaTicks = -deltaTicks;
1058 
1059  if (planeIdx == 0)
1060  {
1061  m_deltaPeakTimePlane0Vec.emplace_back(deltaTicks);
1062  m_deltaPeakSigmaPlane0Vec.emplace_back(sigmaPeakTime);
1063  }
1064  else if (planeIdx == 1)
1065  {
1066  m_deltaPeakTimePlane1Vec.emplace_back(deltaTicks);
1067  m_deltaPeakSigmaPlane1Vec.emplace_back(sigmaPeakTime);
1068  }
1069  else
1070  {
1071  m_deltaPeakTimePlane2Vec.emplace_back(deltaTicks);
1072  m_deltaPeakSigmaPlane2Vec.emplace_back(sigmaPeakTime);
1073  }
1074  }
1075 
1076  // delta peak time consistency check here
1077  if (deltaPeakTime < m_deltaPeakTimeSig * sigmaPeakTime) // 2 sigma consistency? (do this way to avoid divide)
1078  {
1079  // We assume in this routine that we are looking at hits in different views
1080  // The first mission is to check that the wires intersect
1081  const geo::WireID& hit1WireID = hit1->WireID();
1082  const geo::WireID& hit2WireID = hit2->WireID();
1083 
1084  geo::WireIDIntersection widIntersect;
1085 
1086  if (WireIDsIntersect(hit1WireID, hit2WireID, widIntersect))
1087  {
1088  float oneOverWghts = hit1SigSq * hit2SigSq / (hit1SigSq + hit2SigSq);
1089  float avePeakTime = (hit1Peak / hit1SigSq + hit2Peak / hit2SigSq) * oneOverWghts;
1090  float averageCharge = 0.5 * (hit1->getHit()->Integral() + hit2->getHit()->Integral());
1091  float hitChiSquare = std::pow((hit1Peak - avePeakTime),2) / hit1SigSq
1092  + std::pow((hit2Peak - avePeakTime),2) / hit2SigSq;
1093 
1094  float xPositionHit1(hit1->getXPosition());
1095  float xPositionHit2(hit2->getXPosition());
1096  float xPosition = (xPositionHit1 / hit1SigSq + xPositionHit2 / hit2SigSq) * hit1SigSq * hit2SigSq / (hit1SigSq + hit2SigSq);
1097 
1098  Eigen::Vector3f position(xPosition, float(widIntersect.y), float(widIntersect.z)-m_zPosOffset);
1099 
1100  // If to here then we need to sort out the hit pair code telling what views are used
1101  unsigned statusBits = 1 << hit1->WireID().Plane | 1 << hit2->WireID().Plane;
1102 
1103  // handle status bits for the 2D hits
1104  if (hit1->getStatusBits() & reco::ClusterHit2D::USEDINPAIR) hit1->setStatusBit(reco::ClusterHit2D::SHAREDINPAIR);
1105  if (hit2->getStatusBits() & reco::ClusterHit2D::USEDINPAIR) hit2->setStatusBit(reco::ClusterHit2D::SHAREDINPAIR);
1106 
1107  hit1->setStatusBit(reco::ClusterHit2D::USEDINPAIR);
1108  hit2->setStatusBit(reco::ClusterHit2D::USEDINPAIR);
1109 
1110  reco::ClusterHit2DVec hitVector;
1111 
1112  hitVector.resize(3, NULL);
1113 
1114  hitVector[hit1->WireID().Plane] = hit1;
1115  hitVector[hit2->WireID().Plane] = hit2;
1116 
1117  unsigned int cryostatIdx = hit1->WireID().Cryostat;
1118  unsigned int tpcIdx = hit1->WireID().TPC;
1119 
1120  // Initialize the wireIdVec
1121  std::vector<geo::WireID> wireIDVec = {geo::WireID(cryostatIdx,tpcIdx,0,0),
1122  geo::WireID(cryostatIdx,tpcIdx,1,0),
1123  geo::WireID(cryostatIdx,tpcIdx,2,0)};
1124 
1125  wireIDVec[hit1->WireID().Plane] = hit1->WireID();
1126  wireIDVec[hit2->WireID().Plane] = hit2->WireID();
1127 
1128  // For compiling at the moment
1129  std::vector<float> hitDelTSigVec = {0.,0.,0.};
1130 
1131  hitDelTSigVec[hit1->WireID().Plane] = deltaPeakTime / sigmaPeakTime;
1132  hitDelTSigVec[hit2->WireID().Plane] = deltaPeakTime / sigmaPeakTime;
1133 
1134  // Create the 3D cluster hit
1135  hitPair.initialize(hitPairCntr,
1136  statusBits,
1137  position,
1138  averageCharge,
1139  avePeakTime,
1140  deltaPeakTime,
1141  sigmaPeakTime,
1142  hitChiSquare,
1143  0.,
1144  0.,
1145  0.,
1146  0.,
1147  hitVector,
1148  hitDelTSigVec,
1149  wireIDVec);
1150 
1151  result = true;
1152  }
1153  }
1154  }
1155 
1156  // Send it back
1157  return result;
1158 }
1159 
1160 
1161 bool SnippetHit3DBuilderICARUS::makeHitTriplet(reco::ClusterHit3D& hitTriplet,
1162  const reco::ClusterHit3D& pair,
1163  const reco::ClusterHit2D* hit) const
1164 {
1165  // Assume failure
1166  bool result(false);
1167 
1168  // We are going to force the wire pitch here, some time in the future we need to fix
1169  static const double wirePitch = 0.5 * m_wirePitchScaleFactor * *std::max_element(m_wirePitch,m_wirePitch+3);
1170 
1171  // Recover hit info
1172  float hitTimeTicks = hit->getTimeTicks();
1173  float hitSigma = hit->getHit()->RMS();
1174 
1175  // Special case check
1176  if (hitSigma > 2. * hit->getHit()->PeakAmplitude()) hitSigma = 2. * hit->getHit()->PeakAmplitude();
1177 
1178  // Let's do a quick consistency check on the input hits to make sure we are in range...
1179  // Require the W hit to be "in range" with the UV Pair
1180  if (fabs(hitTimeTicks - pair.getAvePeakTime()) < m_hitWidthSclFctr * (pair.getSigmaPeakTime() + hitSigma))
1181  {
1182  // Check the distance from the point to the wire the hit is on
1183  float hitWireDist = DistanceFromPointToHitWire(pair.getPosition(), hit->WireID());
1184 
1185  if (m_outputHistograms) m_maxSideVecVec.push_back(hitWireDist);
1186 
1187  // Reject hits that are not within range
1188  if (hitWireDist < wirePitch)
1189  {
1190  if (m_outputHistograms) m_pairWireDistVec.push_back(hitWireDist);
1191 
1192  // Let's mark that this pair had a valid 3 wire combination
1193  pair.setStatusBit(reco::ClusterHit3D::MADESPACEPOINT);
1194 
1195  // Use the existing code to see the U and W hits are willing to pair with the V hit
1196  reco::ClusterHit3D pair0h;
1197  reco::ClusterHit3D pair1h;
1198 
1199  // Recover all the hits involved
1200  const reco::ClusterHit2DVec& pairHitVec = pair.getHits();
1201  const reco::ClusterHit2D* hit0 = pairHitVec[0];
1202  const reco::ClusterHit2D* hit1 = pairHitVec[1];
1203 
1204  if (!hit0) hit0 = pairHitVec[2];
1205  else if (!hit1) hit1 = pairHitVec[2];
1206 
1207  // If good pairs made here then we can try to make a triplet
1208  if (makeHitPair(pair0h, hit0, hit, m_hitWidthSclFctr) && makeHitPair(pair1h, hit1, hit, m_hitWidthSclFctr))
1209  {
1210  // Get a copy of the input hit vector (note the order is by plane - by definition)
1211  reco::ClusterHit2DVec hitVector = pair.getHits();
1212 
1213  // include the new hit
1214  hitVector[hit->WireID().Plane] = hit;
1215 
1216  // Set up to get average peak time, hitChiSquare, etc.
1217  unsigned int statusBits(0x7);
1218  float avePeakTime(0.);
1219  float weightSum(0.);
1220  float xPosition(0.);
1221 
1222  // And get the wire IDs
1223  std::vector<geo::WireID> wireIDVec = {geo::WireID(), geo::WireID(), geo::WireID()};
1224 
1225  // First loop through the hits to get WireIDs and calculate the averages
1226  for(size_t planeIdx = 0; planeIdx < 3; planeIdx++)
1227  {
1228  const reco::ClusterHit2D* hit2D = hitVector[planeIdx];
1229 
1230  wireIDVec[planeIdx] = hit2D->WireID();
1231 
1232  float hitRMS = hit2D->getHit()->RMS();
1233  float peakTime = hit2D->getTimeTicks();
1234 
1235  // Basically, allow the range to extend to the nearest end of the snippet
1236  if (hit2D->getHit()->DegreesOfFreedom() < 2) hitRMS *= m_LongHitStretchFctr;
1237  //hitRMS = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),float(hit2D->getHit()->EndTick())-hit2D->getTimeTicks());
1238 
1239  float weight = 1. / (hitRMS * hitRMS);
1240 
1241  avePeakTime += peakTime * weight;
1242  xPosition += hit2D->getXPosition() * weight;
1243  weightSum += weight;
1244  }
1245 
1246  avePeakTime /= weightSum;
1247  xPosition /= weightSum;
1248 
1249  Eigen::Vector2f pair0hYZVec(pair0h.getPosition()[1],pair0h.getPosition()[2]);
1250  Eigen::Vector2f pair1hYZVec(pair1h.getPosition()[1],pair1h.getPosition()[2]);
1251  Eigen::Vector2f pairYZVec(pair.getPosition()[1],pair.getPosition()[2]);
1252  Eigen::Vector3f position(xPosition,
1253  float((pairYZVec[0] + pair0hYZVec[0] + pair1hYZVec[0]) / 3.),
1254  float((pairYZVec[1] + pair0hYZVec[1] + pair1hYZVec[1]) / 3.));
1255 
1256  // Armed with the average peak time, now get hitChiSquare and the sig vec
1257  float hitChiSquare(0.);
1258  float sigmaPeakTime(std::sqrt(1./weightSum));
1259  std::vector<float> hitDelTSigVec;
1260 
1261  for(const auto& hit2D : hitVector)
1262  {
1263  float hitRMS = hit2D->getHit()->RMS();
1264 
1265  // Basically, allow the range to extend to the nearest end of the snippet
1266  if (hit2D->getHit()->DegreesOfFreedom() < 2) hitRMS *= m_LongHitStretchFctr;
1267  //hitRMS = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),float(hit2D->getHit()->EndTick())-hit2D->getTimeTicks());
1268 
1269  float combRMS = std::sqrt(hitRMS*hitRMS - sigmaPeakTime*sigmaPeakTime);
1270  float peakTime = hit2D->getTimeTicks();
1271  float deltaTime = peakTime - avePeakTime;
1272  float hitSig = deltaTime / combRMS;
1273 
1274  hitChiSquare += hitSig * hitSig;
1275 
1276  hitDelTSigVec.emplace_back(std::fabs(hitSig));
1277  }
1278 
1279  if (m_outputHistograms) m_chiSquare3DVec.push_back(hitChiSquare);
1280 
1281  int lowMinIndex(std::numeric_limits<int>::max());
1282  int lowMaxIndex(std::numeric_limits<int>::min());
1283  int hiMinIndex(std::numeric_limits<int>::max());
1284  int hiMaxIndex(std::numeric_limits<int>::min());
1285 
1286  // First task is to get the min/max values for the common overlap region
1287  for(const auto& hit2D : hitVector)
1288  {
1289  float range = m_rangeNumSig * hit2D->getHit()->RMS();
1290 
1291  // Basically, allow the range to extend to the nearest end of the snippet
1292  if (hit2D->getHit()->DegreesOfFreedom() < 2) range *= m_LongHitStretchFctr;
1293  //range = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),float(hit2D->getHit()->EndTick())-hit2D->getTimeTicks());
1294 
1295  int hitStart = hit2D->getHit()->PeakTime() - range - 0.5;
1296  int hitStop = hit2D->getHit()->PeakTime() + range + 0.5;
1297 
1298  lowMinIndex = std::min(hitStart, lowMinIndex);
1299  lowMaxIndex = std::max(hitStart, lowMaxIndex);
1300  hiMinIndex = std::min(hitStop + 1, hiMinIndex);
1301  hiMaxIndex = std::max(hitStop + 1, hiMaxIndex);
1302  }
1303 
1304  // Keep only "good" hits...
1305  if (hitChiSquare < m_maxHit3DChiSquare && hiMinIndex > lowMaxIndex)
1306  {
1307  // One more pass through hits to get charge
1308  std::vector<float> chargeVec;
1309 
1310  for(const auto& hit2D : hitVector)
1311  chargeVec.push_back(chargeIntegral(hit2D->getHit()->PeakTime(),hit2D->getHit()->PeakAmplitude(),hit2D->getHit()->RMS(),1.,lowMaxIndex,hiMinIndex));
1312 
1313  float totalCharge = std::accumulate(chargeVec.begin(),chargeVec.end(),0.) / float(chargeVec.size());
1314  float overlapRange = float(hiMinIndex - lowMaxIndex);
1315  float overlapFraction = overlapRange / float(hiMaxIndex - lowMinIndex);
1316 
1317  // Set up to compute the charge asymmetry
1318  std::vector<float> smallestChargeDiffVec;
1319  std::vector<float> chargeAveVec;
1320  float smallestDiff(std::numeric_limits<float>::max());
1321  float maxDeltaPeak(0.);
1322  size_t chargeIndex(0);
1323 
1324  for(size_t idx = 0; idx < 3; idx++)
1325  {
1326  size_t leftIdx = (idx + 2) % 3;
1327  size_t rightIdx = (idx + 1) % 3;
1328 
1329  smallestChargeDiffVec.push_back(std::abs(chargeVec[leftIdx] - chargeVec[rightIdx]));
1330  chargeAveVec.push_back(float(0.5 * (chargeVec[leftIdx] + chargeVec[rightIdx])));
1331 
1332  if (smallestChargeDiffVec.back() < smallestDiff)
1333  {
1334  smallestDiff = smallestChargeDiffVec.back();
1335  chargeIndex = idx;
1336  }
1337 
1338  // Take opportunity to look at peak time diff
1339  float deltaPeakTime = hitVector[rightIdx]->getTimeTicks() - hitVector[leftIdx]->getTimeTicks();
1340 
1341  if (std::abs(deltaPeakTime) > maxDeltaPeak) maxDeltaPeak = std::abs(deltaPeakTime);
1342 
1343  if (m_outputHistograms)
1344  {
1345  int deltaTimeIdx = (hitVector[leftIdx]->WireID().Plane + hitVector[rightIdx]->WireID().Plane) - 1;
1346  float combRMS = std::sqrt(std::pow(hitVector[leftIdx]->getHit()->RMS(),2) + std::pow(hitVector[rightIdx]->getHit()->RMS(),2));
1347 
1348  m_deltaTimeVec.push_back(deltaPeakTime);
1349 
1350  // Want to get the sign of the difference correct...
1351  if (deltaTimeIdx == 0) // This is planes 1 and 0
1352  {
1353  m_deltaTime0Vec.emplace_back(float(hitVector[1]->getTimeTicks() - hitVector[0]->getTimeTicks()));
1354  m_deltaSigma0Vec.emplace_back(combRMS);
1355  }
1356  else if (deltaTimeIdx == 1) // This is planes 0 and 2
1357  {
1358  m_deltaTime1Vec.emplace_back(float(hitVector[2]->getTimeTicks() - hitVector[0]->getTimeTicks()));
1359  m_deltaSigma1Vec.emplace_back(combRMS);
1360  }
1361  else // must be planes 1 and 2
1362  {
1363  m_deltaTime2Vec.emplace_back(float(hitVector[2]->getTimeTicks() - hitVector[1]->getTimeTicks()));
1364  m_deltaSigma2Vec.emplace_back(combRMS);
1365  }
1366  }
1367  }
1368 
1369  float chargeAsymmetry = (chargeAveVec[chargeIndex] - chargeVec[chargeIndex]) / (chargeAveVec[chargeIndex] + chargeVec[chargeIndex]);
1370 
1371  // If this is true there has to be a negative charge that snuck in somehow
1372  if (chargeAsymmetry < -1. || chargeAsymmetry > 1.)
1373  {
1374  const geo::WireID& hitWireID = hitVector[chargeIndex]->WireID();
1375 
1376  std::cout << "============> Charge asymmetry out of range: " << chargeAsymmetry << " <============" << std::endl;
1377  std::cout << " hit C: " << hitWireID.Cryostat << ", TPC: " << hitWireID.TPC << ", Plane: " << hitWireID.Plane << ", Wire: " << hitWireID.Wire << std::endl;
1378  std::cout << " charge: " << chargeVec[0] << ", " << chargeVec[1] << ", " << chargeVec[2] << std::endl;
1379  std::cout << " index: " << chargeIndex << ", smallest diff: " << smallestDiff << std::endl;
1380  return result;
1381  }
1382 
1383  // Usurping "deltaPeakTime" to be the maximum pull
1384  float deltaPeakTime = *std::max_element(hitDelTSigVec.begin(),hitDelTSigVec.end());
1385 
1386  if (m_outputHistograms)
1387  {
1388  m_smallChargeDiffVec.push_back(smallestDiff);
1389  m_smallIndexVec.push_back(chargeIndex);
1390  m_maxPullVec.push_back(deltaPeakTime);
1391  m_qualityMetricVec.push_back(hitChiSquare);
1392  m_spacePointChargeVec.push_back(totalCharge);
1393  m_overlapFractionVec.push_back(overlapFraction);
1394  m_overlapRangeVec.push_back(overlapRange);
1395  m_maxDeltaPeakVec.push_back(maxDeltaPeak);
1396  m_hitAsymmetryVec.push_back(chargeAsymmetry);
1397  }
1398 
1399  // Try to weed out cases where overlap doesn't match peak separation
1400  if (maxDeltaPeak > 3. * overlapRange) return result;
1401 
1402  // Create the 3D cluster hit
1403  hitTriplet.initialize(0,
1404  statusBits,
1405  position,
1406  totalCharge,
1407  avePeakTime,
1408  deltaPeakTime,
1409  sigmaPeakTime,
1410  hitChiSquare,
1411  overlapFraction,
1412  chargeAsymmetry,
1413  0.,
1414  0.,
1415  hitVector,
1416  hitDelTSigVec,
1417  wireIDVec);
1418 
1419  // Since we are keeping the triplet, mark the hits as used
1420  for(const auto& hit2D : hitVector)
1421  {
1422  if (hit2D->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET) hit2D->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
1423 
1424  hit2D->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
1425  }
1426 
1427  // Mark the input pair
1428  pair.setStatusBit(reco::ClusterHit3D::MADESPACEPOINT);
1429 
1430  result = true;
1431  }
1432 // else std::cout << "-Rejecting triple with chiSquare: " << hitChiSquare << " and hiMinIndex: " << hiMinIndex << ", loMaxIndex: " << lowMaxIndex << std::endl;
1433  }
1434  }
1435  }
1436 // else std::cout << "-MakeTriplet hit cut, delta: " << hitTimeTicks - pair.getAvePeakTime() << ", min scale fctr: " <<m_hitWidthSclFctr << ", pair sig: " << pair.getSigmaPeakTime() << ", hitSigma: " << hitSigma << std::endl;
1437 
1438  // return success/fail
1439  return result;
1440 }
1441 
1442 bool SnippetHit3DBuilderICARUS::WireIDsIntersect(const geo::WireID& wireID0, const geo::WireID& wireID1, geo::WireIDIntersection& widIntersection) const
1443 {
1444  bool success(false);
1445 
1446  // Do quick check that things are in the same logical TPC
1447  if (wireID0.Cryostat != wireID1.Cryostat || wireID0.TPC != wireID1.TPC || wireID0.Plane == wireID1.Plane) return success;
1448 
1449  // Recover wire geometry information for each wire
1450  const geo::WireGeo& wireGeo0 = m_geometry->WireIDToWireGeo(wireID0);
1451  const geo::WireGeo& wireGeo1 = m_geometry->WireIDToWireGeo(wireID1);
1452 
1453  // Get wire position and direction for first wire
1454  double wirePosArr[3] = {0.,0.,0.};
1455  wireGeo0.GetCenter(wirePosArr);
1456 
1457  Eigen::Vector3f wirePos0(wirePosArr[0],wirePosArr[1],wirePosArr[2]);
1458  Eigen::Vector3f wireDir0(wireGeo0.Direction().X(),wireGeo0.Direction().Y(),wireGeo0.Direction().Z());
1459 
1460  //*********************************
1461  // Kludge
1462 // if (wireID0.Plane > 0) wireDir0[2] = -wireDir0[2];
1463 
1464  // And now the second one
1465  wireGeo1.GetCenter(wirePosArr);
1466 
1467  Eigen::Vector3f wirePos1(wirePosArr[0],wirePosArr[1],wirePosArr[2]);
1468  Eigen::Vector3f wireDir1(wireGeo1.Direction().X(),wireGeo1.Direction().Y(),wireGeo1.Direction().Z());
1469 
1470  //**********************************
1471  // Kludge
1472 // if (wireID1.Plane > 0) wireDir1[2] = -wireDir1[2];
1473 
1474  // Get the distance of closest approach
1475  float arcLen0;
1476  float arcLen1;
1477 
1478  if (closestApproach(wirePos0, wireDir0, wirePos1, wireDir1, arcLen0, arcLen1))
1479  {
1480  // Now check that arc lengths are within range
1481  if (std::abs(arcLen0) < wireGeo0.HalfL() && std::abs(arcLen1) < wireGeo1.HalfL())
1482  {
1483  Eigen::Vector3f poca0 = wirePos0 + arcLen0 * wireDir0;
1484 
1485  widIntersection.y = poca0[1];
1486  widIntersection.z = poca0[2];
1487 
1488  success = true;
1489  }
1490  }
1491 
1492  return success;
1493 }
1494 
1495 float SnippetHit3DBuilderICARUS::closestApproach(const Eigen::Vector3f& P0,
1496  const Eigen::Vector3f& u0,
1497  const Eigen::Vector3f& P1,
1498  const Eigen::Vector3f& u1,
1499  float& arcLen0,
1500  float& arcLen1) const
1501 {
1502  // Technique is to compute the arclength to each point of closest approach
1503  Eigen::Vector3f w0 = P0 - P1;
1504  float a(1.);
1505  float b(u0.dot(u1));
1506  float c(1.);
1507  float d(u0.dot(w0));
1508  float e(u1.dot(w0));
1509  float den(a * c - b * b);
1510 
1511  arcLen0 = (b * e - c * d) / den;
1512  arcLen1 = (a * e - b * d) / den;
1513 
1514  Eigen::Vector3f poca0 = P0 + arcLen0 * u0;
1515  Eigen::Vector3f poca1 = P1 + arcLen1 * u1;
1516 
1517  return (poca0 - poca1).norm();
1518 }
1519 
1520 float SnippetHit3DBuilderICARUS::chargeIntegral(float peakMean,
1521  float peakAmp,
1522  float peakSigma,
1523  float areaNorm,
1524  int low,
1525  int hi) const
1526 {
1527  float integral(0);
1528 
1529  for(int sigPos = low; sigPos < hi; sigPos++)
1530  {
1531  float arg = (float(sigPos) - peakMean + 0.5) / peakSigma;
1532  integral += peakAmp * std::exp(-0.5 * arg * arg);
1533  }
1534 
1535  return integral;
1536 }
1537 
1538 bool SnippetHit3DBuilderICARUS::makeDeadChannelPair(reco::ClusterHit3D& pairOut,
1539  const reco::ClusterHit3D& pair,
1540  size_t maxChanStatus,
1541  size_t minChanStatus,
1542  float minOverlap) const
1543 {
1544  // Assume failure (most common result)
1545  bool result(false);
1546 
1547  const reco::ClusterHit2D* hit0 = pair.getHits()[0];
1548  const reco::ClusterHit2D* hit1 = pair.getHits()[1];
1549 
1550  size_t missPlane(2);
1551 
1552  // u plane hit is missing
1553  if (!hit0)
1554  {
1555  hit0 = pair.getHits()[2];
1556  missPlane = 0;
1557  }
1558  // v plane hit is missing
1559  else if (!hit1)
1560  {
1561  hit1 = pair.getHits()[2];
1562  missPlane = 1;
1563  }
1564 
1565  // Which plane is missing?
1566  geo::WireID wireID0 = hit0->WireID();
1567  geo::WireID wireID1 = hit1->WireID();
1568 
1569  // Ok, recover the wireID expected in the third plane...
1570  geo::WireID wireIn(wireID0.Cryostat,wireID0.TPC,missPlane,0);
1571  geo::WireID wireID = NearestWireID(pair.getPosition(), wireIn);
1572 
1573  // There can be a round off issue so check the next wire as well
1574  bool wireStatus = m_channelStatus[wireID.Plane][wireID.Wire] < maxChanStatus && m_channelStatus[wireID.Plane][wireID.Wire] >= minChanStatus;
1575  bool wireOneStatus = m_channelStatus[wireID.Plane][wireID.Wire+1] < maxChanStatus && m_channelStatus[wireID.Plane][wireID.Wire+1] >= minChanStatus;
1576 
1577  // Make sure they are of at least the minimum status
1578  if(wireStatus || wireOneStatus)
1579  {
1580  // Sort out which is the wire we're dealing with
1581  if (!wireStatus) wireID.Wire += 1;
1582 
1583  // Want to refine position since we "know" the missing wire
1584  geo::WireIDIntersection widIntersect0;
1585 
1586  if (m_geometry->WireIDsIntersect(wireID0, wireID, widIntersect0))
1587  {
1588  geo::WireIDIntersection widIntersect1;
1589 
1590  if (m_geometry->WireIDsIntersect(wireID1, wireID, widIntersect1))
1591  {
1592  Eigen::Vector3f newPosition(pair.getPosition()[0],pair.getPosition()[1],pair.getPosition()[2]);
1593 
1594  newPosition[1] = (newPosition[1] + widIntersect0.y + widIntersect1.y) / 3.;
1595  newPosition[2] = (newPosition[2] + widIntersect0.z + widIntersect1.z - 2. * m_zPosOffset) / 3.;
1596 
1597  pairOut = pair;
1598  pairOut.setWireID(wireID);
1599  pairOut.setPosition(newPosition);
1600 
1601  if (hit0->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET) hit0->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
1602  if (hit1->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET) hit1->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
1603 
1604  hit0->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
1605  hit1->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
1606 
1607  result = true;
1608  }
1609  }
1610  }
1611 
1612  return result;
1613 }
1614 
1615 const reco::ClusterHit2D* SnippetHit3DBuilderICARUS::FindBestMatchingHit(const Hit2DSet& hit2DSet, const reco::ClusterHit3D& pair, float pairDeltaTimeLimits) const
1616 {
1617  static const float minCharge(0.);
1618 
1619  const reco::ClusterHit2D* bestVHit(0);
1620 
1621  float pairAvePeakTime(pair.getAvePeakTime());
1622 
1623  // Idea is to loop through the input set of hits and look for the best combination
1624  for (const auto& hit2D : hit2DSet)
1625  {
1626  if (hit2D->getHit()->Integral() < minCharge) continue;
1627 
1628  float hitVPeakTime(hit2D->getTimeTicks());
1629  float deltaPeakTime(pairAvePeakTime-hitVPeakTime);
1630 
1631  if (deltaPeakTime > pairDeltaTimeLimits) continue;
1632 
1633  if (deltaPeakTime < -pairDeltaTimeLimits) break;
1634 
1635  pairDeltaTimeLimits = fabs(deltaPeakTime);
1636  bestVHit = hit2D;
1637  }
1638 
1639  return bestVHit;
1640 }
1641 
1642 int SnippetHit3DBuilderICARUS::FindNumberInRange(const Hit2DSet& hit2DSet, const reco::ClusterHit3D& pair, float range) const
1643 {
1644  static const float minCharge(0.);
1645 
1646  int numberInRange(0);
1647  float pairAvePeakTime(pair.getAvePeakTime());
1648 
1649  // Idea is to loop through the input set of hits and look for the best combination
1650  for (const auto& hit2D : hit2DSet)
1651  {
1652  if (hit2D->getHit()->Integral() < minCharge) continue;
1653 
1654  float hitVPeakTime(hit2D->getTimeTicks());
1655  float deltaPeakTime(pairAvePeakTime-hitVPeakTime);
1656 
1657  if (deltaPeakTime > range) continue;
1658 
1659  if (deltaPeakTime < -range) break;
1660 
1661  numberInRange++;
1662  }
1663 
1664  return numberInRange;
1665 }
1666 
1667 geo::WireID SnippetHit3DBuilderICARUS::NearestWireID(const Eigen::Vector3f& position, const geo::WireID& wireIDIn) const
1668 {
1669  geo::WireID wireID = wireIDIn;
1670 
1671  // Embed the call to the geometry's services nearest wire id method in a try-catch block
1672  try
1673  {
1674  // Switch from NearestWireID to this method to avoid the roundoff error issues...
1675  double distanceToWire = m_geometry->Plane(wireIDIn).WireCoordinate(position.data());
1676 
1677  wireID.Wire = int(distanceToWire);
1678  }
1679  catch(std::exception& exc)
1680  {
1681  // This can happen, almost always because the coordinates are **just** out of range
1682  mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - " << exc.what() << std::endl;
1683 
1684  // Assume extremum for wire number depending on z coordinate
1685  if (position[2] < 0.5 * m_geometry->DetLength()) wireID.Wire = 0;
1686  else wireID.Wire = m_geometry->Nwires(wireIDIn.Plane) - 1;
1687  }
1688 
1689  return wireID;
1690 }
1691 
1692 float SnippetHit3DBuilderICARUS::DistanceFromPointToHitWire(const Eigen::Vector3f& position, const geo::WireID& wireIDIn) const
1693 {
1694  float distance = std::numeric_limits<float>::max();
1695 
1696  // Embed the call to the geometry's services nearest wire id method in a try-catch block
1697  try
1698  {
1699  // Recover wire geometry information for each wire
1700  const geo::WireGeo& wireGeo = m_geometry->WireIDToWireGeo(wireIDIn);
1701 
1702  // Get wire position and direction for first wire
1703  double wirePosArr[3] = {0.,0.,0.};
1704  wireGeo.GetCenter(wirePosArr);
1705 
1706  Eigen::Vector3f wirePos(wirePosArr[0],wirePosArr[1],wirePosArr[2]);
1707  Eigen::Vector3f wireDir(wireGeo.Direction().X(),wireGeo.Direction().Y(),wireGeo.Direction().Z());
1708 
1709  //*********************************
1710  // Kludge
1711 // if (wireIDIn.Plane > 0) wireDir[2] = -wireDir[2];
1712 
1713 
1714  // Want the hit position to have same x value as wire coordinates
1715  Eigen::Vector3f hitPosition(wirePos[0],position[1],position[2]);
1716 
1717  // Get arc length to doca
1718  double arcLen = (hitPosition - wirePos).dot(wireDir);
1719 
1720  // Make sure arclen is in range
1721  if (abs(arcLen) < wireGeo.HalfL())
1722  {
1723  Eigen::Vector3f docaVec = hitPosition - (wirePos + arcLen * wireDir);
1724 
1725  distance = docaVec.norm();
1726  }
1727  }
1728  catch(std::exception& exc)
1729  {
1730  // This can happen, almost always because the coordinates are **just** out of range
1731  mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - " << exc.what() << std::endl;
1732 
1733  // Assume extremum for wire number depending on z coordinate
1734  distance = 0.;
1735  }
1736 
1737  return distance;
1738 }
1739 
1740 //------------------------------------------------------------------------------------------------------------------------------------------
1741 bool SetHitTimeOrder(const reco::ClusterHit2D* left, const reco::ClusterHit2D* right)
1742 {
1743  // Sort by "modified start time" of pulse
1744  return left->getHit()->PeakTime() < right->getHit()->PeakTime();
1745 }
1746 
1747 bool Hit2DSetCompare::operator() (const reco::ClusterHit2D* left, const reco::ClusterHit2D* right) const
1748 {
1749  return left->getHit()->PeakTime() < right->getHit()->PeakTime();
1750 }
1751 
1752 //------------------------------------------------------------------------------------------------------------------------------------------
1753 void SnippetHit3DBuilderICARUS::CollectArtHits(const art::Event& evt) const
1754 {
1755  /**
1756  * @brief Recover the 2D hits from art and fill out the local data structures for the 3D clustering
1757  */
1758 
1759  // Start by getting a vector of valid, non empty hit collections to make sure we really have something to do here...
1760  // Here is a container for the hits...
1761  std::vector<const recob::Hit*> recobHitVec;
1762 
1763  // Loop through the list of input sources
1764  for(const auto& inputTag : m_hitFinderTagVec)
1765  {
1766  art::Handle< std::vector<recob::Hit> > recobHitHandle;
1767  evt.getByLabel(inputTag, recobHitHandle);
1768 
1769  if (!recobHitHandle.isValid() || recobHitHandle->size() == 0) continue;
1770 
1771  recobHitVec.reserve(recobHitVec.size() + recobHitHandle->size());
1772 
1773  for(const auto& hit : *recobHitHandle) recobHitVec.push_back(&hit);
1774  }
1775 
1776  // If the vector is empty there is nothing to do
1777  if (recobHitVec.empty()) return;
1778 
1779  cet::cpu_timer theClockMakeHits;
1780 
1781  if (m_enableMonitoring) theClockMakeHits.start();
1782 
1783  // Need the detector properties which needs the clocks
1784  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
1785  auto const det_prop = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
1786 
1787  // Try to output a formatted string
1788  std::string debugMessage("");
1789 
1790  // Keep track of x position limits
1791  std::map<geo::PlaneID,double> planeIDToPositionMap;
1792 
1793  // Initialize the plane to hit vector map
1794  for(size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++)
1795  {
1796  for(size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++)
1797  {
1798  m_planeToSnippetHitMap[geo::PlaneID(cryoIdx,tpcIdx,0)] = SnippetHitMap();
1799  m_planeToSnippetHitMap[geo::PlaneID(cryoIdx,tpcIdx,1)] = SnippetHitMap();
1800  m_planeToSnippetHitMap[geo::PlaneID(cryoIdx,tpcIdx,2)] = SnippetHitMap();
1801 
1802  // Should we provide output?
1803  if (!m_weHaveAllBeenHereBefore)
1804  {
1805  std::ostringstream outputString;
1806 
1807  outputString << "***> plane 0 offset: " << m_PlaneToT0OffsetMap.find(geo::PlaneID(cryoIdx,tpcIdx,0))->second
1808  << ", plane 1: " << m_PlaneToT0OffsetMap.find(geo::PlaneID(cryoIdx,tpcIdx,1))->second
1809  << ", plane 2: " << m_PlaneToT0OffsetMap.find(geo::PlaneID(cryoIdx,tpcIdx,2))->second << "\n";
1810  outputString << " Det prop plane 0: " << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,0)) << ", plane 1: " << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,1)) << ", plane 2: " << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,2)) << ", Trig: " << trigger_offset(clock_data) << "\n";
1811  debugMessage += outputString.str() + "\n";
1812  }
1813 
1814  double xPosition(det_prop.ConvertTicksToX(0., 2, tpcIdx, cryoIdx));
1815 
1816  planeIDToPositionMap[geo::PlaneID(cryoIdx,tpcIdx,2)] = xPosition;
1817  }
1818  }
1819 
1820  if (!m_weHaveAllBeenHereBefore)
1821  {
1822  for(const auto& planeToPositionPair : planeIDToPositionMap)
1823  {
1824  std::ostringstream outputString;
1825 
1826  outputString << "***> Plane " << planeToPositionPair.first << " has time=0 position: " << planeToPositionPair.second << "\n";
1827 
1828  debugMessage += outputString.str();
1829  }
1830 
1831  mf::LogDebug("Cluster3D") << debugMessage << std::endl;
1832 
1833  m_weHaveAllBeenHereBefore = true;
1834  }
1835 
1836  // Cycle through the recob hits to build ClusterHit2D objects and insert
1837  // them into the map
1838  for (const auto& recobHit : recobHitVec)
1839  {
1840  // Reject hits with negative charge, these are misreconstructed
1841  if (recobHit->Integral() < 0.) continue;
1842 
1843  // For some detectors we can have multiple wire ID's associated to a given channel.
1844  // So we recover the list of these wire IDs
1845  const std::vector<geo::WireID>& wireIDs = m_geometry->ChannelToWire(recobHit->Channel());
1846 
1847  // Start/End ticks to identify the snippet
1848  HitStartEndPair hitStartEndPair(recobHit->StartTick(),recobHit->EndTick());
1849 
1850  // Can this really happen?
1851  if (hitStartEndPair.second <= hitStartEndPair.first)
1852  {
1853  std::cout << "Yes, found a hit with end time less than start time: " << hitStartEndPair.first << "/" << hitStartEndPair.second << ", mult: " << recobHit->Multiplicity() << std::endl;
1854  continue;
1855  }
1856 
1857  // And then loop over all possible to build out our maps
1858  //for(const auto& wireID : wireIDs)
1859  for(auto wireID : wireIDs)
1860  {
1861  // Check if this is an invalid TPC
1862  // (for example, in protoDUNE there are logical TPC's which see no signal)
1863  if (std::find(m_invalidTPCVec.begin(),m_invalidTPCVec.end(),wireID.TPC) != m_invalidTPCVec.end()) continue;
1864 
1865  // Note that a plane ID will define cryostat, TPC and plane
1866  const geo::PlaneID& planeID = wireID.planeID();
1867 
1868  double hitPeakTime(recobHit->PeakTime() - m_PlaneToT0OffsetMap.find(planeID)->second);
1869  double xPosition(det_prop.ConvertTicksToX(recobHit->PeakTime(), planeID.Plane, planeID.TPC, planeID.Cryostat));
1870 
1871  m_clusterHit2DMasterList.emplace_back(0, 0., 0., xPosition, hitPeakTime, wireID, recobHit);
1872 
1873  m_planeToSnippetHitMap[planeID][hitStartEndPair].emplace_back(&m_clusterHit2DMasterList.back());
1874  m_planeToWireToHitSetMap[planeID][wireID.Wire].insert(&m_clusterHit2DMasterList.back());
1875  }
1876  }
1877 
1878  // Make a loop through to sort the recover hits in time order
1879 // for(auto& hitVectorMap : m_planeToSnippetHitMap)
1880 // std::sort(hitVectorMap.second.begin(), hitVectorMap.second.end(), SetHitTimeOrder);
1881 
1882  if (m_enableMonitoring)
1883  {
1884  theClockMakeHits.stop();
1885 
1886  m_timeVector[COLLECTARTHITS] = theClockMakeHits.accumulated_real_time();
1887  }
1888 
1889  mf::LogDebug("Cluster3D") << ">>>>> Number of ART hits: " << m_clusterHit2DMasterList.size() << std::endl;
1890 }
1891 
1892 //------------------------------------------------------------------------------------------------------------------------------------------
1893 
1894 void SnippetHit3DBuilderICARUS::CreateNewRecobHitCollection(art::Event& event,
1895  reco::HitPairList& hitPairList,
1896  std::vector<recob::Hit>& hitPtrVec,
1897  RecobHitToPtrMap& recobHitToPtrMap)
1898 {
1899  // Set up the timing
1900  cet::cpu_timer theClockBuildNewHits;
1901 
1902  if (m_enableMonitoring) theClockBuildNewHits.start();
1903 
1904  // We want to build a unique list of hits from the input 3D points which, we know, reuse 2D points frequently
1905  // At the same time we need to create a new 2D hit with the "correct" WireID and replace the old 2D hit in the
1906  // ClusterHit2D object with this new one... while keeping track of the use of the old ones. My head is spinning...
1907  // Declare a set which will allow us to keep track of those CusterHit2D objects we have seen already
1908  std::set<const reco::ClusterHit2D*> visitedHit2DSet;
1909 
1910  // Use this handy art utility to make art::Ptr objects to the new recob::Hits for use in the output phase
1911  art::PtrMaker<recob::Hit> ptrMaker(event);
1912 
1913  // Reserve enough memory to replace every recob::Hit we have considered (this is upper limit)
1914  hitPtrVec.reserve(m_clusterHit2DMasterList.size());
1915 
1916  // Scheme is to loop through all 3D hits, then through each associated ClusterHit2D object
1917  for(reco::ClusterHit3D& hit3D : hitPairList)
1918  {
1919  reco::ClusterHit2DVec& hit2DVec = hit3D.getHits();
1920 
1921  // The loop is over the index so we can recover the correct WireID to associate to the new hit when made
1922  for(size_t idx = 0; idx < hit3D.getHits().size(); idx++)
1923  {
1924  const reco::ClusterHit2D* hit2D = hit2DVec[idx];
1925 
1926  if (!hit2D) continue;
1927 
1928  // Have we seen this 2D hit already?
1929  if (visitedHit2DSet.find(hit2D) == visitedHit2DSet.end())
1930  {
1931  visitedHit2DSet.insert(hit2D);
1932 
1933  // Create and save the new recob::Hit with the correct WireID
1934  hitPtrVec.emplace_back(recob::HitCreator(*hit2D->getHit(), hit3D.getWireIDs()[idx]).copy());
1935 
1936  // Recover a pointer to it...
1937  recob::Hit* newHit = &hitPtrVec.back();
1938 
1939  // Create a mapping from this hit to an art Ptr representing it
1940  recobHitToPtrMap[newHit] = ptrMaker(hitPtrVec.size()-1);
1941 
1942  // And set the pointer to this hit in the ClusterHit2D object
1943  const_cast<reco::ClusterHit2D*>(hit2D)->setHit(newHit);
1944  }
1945  }
1946  }
1947 
1948  size_t numNewHits = hitPtrVec.size();
1949 
1950  if (m_enableMonitoring)
1951  {
1952  theClockBuildNewHits.stop();
1953 
1954  m_timeVector[BUILDNEWHITS] = theClockBuildNewHits.accumulated_real_time();
1955  }
1956 
1957  mf::LogDebug("Cluster3D") << ">>>>> New output recob::Hit size: " << numNewHits << " (vs " << m_clusterHit2DMasterList.size() << " input)" << std::endl;
1958 
1959  return;
1960 }
1961 
1962 void SnippetHit3DBuilderICARUS::makeWireAssns(const art::Event& evt, art::Assns<recob::Wire, recob::Hit>& wireAssns, RecobHitToPtrMap& recobHitPtrMap) const
1963 {
1964  // Let's make sure the input associations container is empty
1965  wireAssns = art::Assns<recob::Wire, recob::Hit>();
1966 
1967  // First task is to recover all of the previous wire <--> hit associations and map them by channel number
1968  // Create the temporary container
1969  std::unordered_map<raw::ChannelID_t, art::Ptr<recob::Wire>> channelToWireMap;
1970 
1971  // Go through the list of input sources and fill out the map
1972  for(const auto& inputTag : m_hitFinderTagVec)
1973  {
1974  art::ValidHandle<std::vector<recob::Hit>> hitHandle = evt.getValidHandle<std::vector<recob::Hit>>(inputTag);
1975 
1976  art::FindOneP<recob::Wire> hitToWireAssns(hitHandle, evt, inputTag);
1977 
1978  if (hitToWireAssns.isValid())
1979  {
1980  for(size_t wireIdx = 0; wireIdx < hitToWireAssns.size(); wireIdx++)
1981  {
1982  art::Ptr<recob::Wire> wire = hitToWireAssns.at(wireIdx);
1983 
1984  channelToWireMap[wire->Channel()] = wire;
1985  }
1986  }
1987  }
1988 
1989  // Now fill the container
1990  for(const auto& hitPtrPair : recobHitPtrMap)
1991  {
1992  raw::ChannelID_t channel = hitPtrPair.first->Channel();
1993 
1994  std::unordered_map<raw::ChannelID_t, art::Ptr<recob::Wire>>::iterator chanWireItr = channelToWireMap.find(channel);
1995 
1996  if (!(chanWireItr != channelToWireMap.end()))
1997  {
1998  //mf::LogDebug("Cluster3D") << "** Did not find channel to wire match! Skipping..." << std::endl;
1999  continue;
2000  }
2001 
2002  wireAssns.addSingle(chanWireItr->second, hitPtrPair.second);
2003  }
2004 
2005  return;
2006 }
2007 
2008 void SnippetHit3DBuilderICARUS::makeRawDigitAssns(const art::Event& evt, art::Assns<raw::RawDigit, recob::Hit>& rawDigitAssns, RecobHitToPtrMap& recobHitPtrMap) const
2009 {
2010  // Let's make sure the input associations container is empty
2011  rawDigitAssns = art::Assns<raw::RawDigit, recob::Hit>();
2012 
2013  // First task is to recover all of the previous wire <--> hit associations and map them by channel number
2014  // Create the temporary container
2015  std::unordered_map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> channelToRawDigitMap;
2016 
2017  // Go through the list of input sources and fill out the map
2018  for(const auto& inputTag : m_hitFinderTagVec)
2019  {
2020  art::ValidHandle<std::vector<recob::Hit>> hitHandle = evt.getValidHandle<std::vector<recob::Hit>>(inputTag);
2021 
2022  art::FindOneP<raw::RawDigit> hitToRawDigitAssns(hitHandle, evt, inputTag);
2023 
2024  if (hitToRawDigitAssns.isValid())
2025  {
2026  for(size_t rawDigitIdx = 0; rawDigitIdx < hitToRawDigitAssns.size(); rawDigitIdx++)
2027  {
2028  art::Ptr<raw::RawDigit> rawDigit = hitToRawDigitAssns.at(rawDigitIdx);
2029 
2030  channelToRawDigitMap[rawDigit->Channel()] = rawDigit;
2031  }
2032  }
2033  }
2034 
2035  // Now fill the container
2036  for(const auto& hitPtrPair : recobHitPtrMap)
2037  {
2038  raw::ChannelID_t channel = hitPtrPair.first->Channel();
2039 
2040  std::unordered_map<raw::ChannelID_t, art::Ptr<raw::RawDigit>>::iterator chanRawDigitItr = channelToRawDigitMap.find(channel);
2041 
2042  if (chanRawDigitItr == channelToRawDigitMap.end())
2043  {
2044  //mf::LogDebug("Cluster3D") << "** Did not find channel to wire match! Skipping..." << std::endl;
2045  continue;
2046  }
2047 
2048  rawDigitAssns.addSingle(chanRawDigitItr->second, hitPtrPair.second);
2049  }
2050 
2051  return;
2052 }
2053 
2054 //------------------------------------------------------------------------------------------------------------------------------------------
2055 
2056 DEFINE_ART_CLASS_TOOL(SnippetHit3DBuilderICARUS)
2057 } // namespace lar_cluster3d
reco::ClusterHit3D::initialize
void initialize(size_t id, unsigned int statusBits, const Eigen::Vector3f &position, float totalCharge, float avePeakTime, float deltaPeakTime, float sigmaPeakTime, float hitChiSquare, float overlapFraction, float chargeAsymmetry, float docaToAxis, float arclenToPoca, const ClusterHit2DVec &hitVec, const std::vector< float > &hitDelTSigVec, const std::vector< geo::WireID > &wireIDVec)
Definition: Cluster3D.cxx:137
lar_cluster3d::SnippetHit3DBuilderICARUS::makeHitPair
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.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1010
geo::WireGeo
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
lar_cluster3d::SnippetHit3DBuilderICARUS::PlaneToT0OffsetMap
std::map< geo::PlaneID, float > PlaneToT0OffsetMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:271
lar_cluster3d::SnippetHit3DBuilderICARUS::makeHitTriplet
bool makeHitTriplet(reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pairIn, const reco::ClusterHit2D *hit2) const
Make a 3D HitPair object by checking two hits.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1161
lar_cluster3d::SnippetHit3DBuilderICARUS::m_PlaneToT0OffsetMap
PlaneToT0OffsetMap m_PlaneToT0OffsetMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:273
lar_cluster3d::SnippetHit3DBuilderICARUS::Hit3DBuilder
virtual void Hit3DBuilder(art::Event &, reco::HitPairList &, RecobHitToPtrMap &) override
Given a set of recob hits, run DBscan to form 3D clusters.
Definition: SnippetHit3DBuilderICARUS_tool.cc:502
reco::ClusterHit3D::setWireID
void setWireID(const geo::WireID &wid) const
Definition: Cluster3D.cxx:172
lar_cluster3d::SnippetHit3DBuilderICARUS::m_LongHitStretchFctr
float m_LongHitStretchFctr
Definition: SnippetHit3DBuilderICARUS_tool.cc:253
lar_cluster3d::SnippetHit3DBuilderICARUS::m_timeVector
std::vector< float > m_timeVector
Definition: SnippetHit3DBuilderICARUS_tool.cc:267
geo::WireIDIntersection::z
double z
z position of intersection
Definition: geo_types.h:805
reco::HitPairList
std::list< reco::ClusterHit3D > HitPairList
Definition: Cluster3D.h:338
lar_cluster3d::SnippetHit3DBuilderICARUS::SnippetHit3DBuilderICARUS
SnippetHit3DBuilderICARUS(fhicl::ParameterSet const &pset)
Constructor.
Definition: SnippetHit3DBuilderICARUS_tool.cc:324
second
second
Definition: detectorproperties_icarus.fcl:8
geo::vect::dot
constexpr auto dot(Vector const &a, Vector const &b)
Return cross product of two vectors.
Definition: geo_vectors_utils.h:1194
lar_cluster3d::SnippetHit3DBuilderICARUS::clear
void clear()
clear the tuple vectors before processing next event
Definition: SnippetHit3DBuilderICARUS_tool.cc:416
lar_cluster3d::SnippetHit3DBuilderICARUS::m_smallIndexVec
std::vector< int > m_smallIndexVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:300
geo::GeometryCore::Plane
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Definition: GeometryCore.h:3000
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaPeakTimePlane1Vec
std::vector< float > m_deltaPeakTimePlane1Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:280
reco::ClusterHit2D::getTimeTicks
float getTimeTicks() const
Definition: Cluster3D.h:75
lar_cluster3d::SnippetHit3DBuilderICARUS::CreateNewRecobHitCollection
void CreateNewRecobHitCollection(art::Event &, reco::HitPairList &, std::vector< recob::Hit > &, RecobHitToPtrMap &)
Create a new 2D hit collection from hits associated to 3D space points.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1894
lar_cluster3d::SnippetHit3DBuilderICARUS::m_hitAsymmetryVec
std::vector< float > m_hitAsymmetryVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:303
lar_cluster3d::SnippetHit3DBuilderICARUS::m_tupleTree
TTree * m_tupleTree
output analysis tree
Definition: SnippetHit3DBuilderICARUS_tool.cc:276
lar_cluster3d::SetStartTimeOrder
Definition: SnippetHit3DBuilderICARUS_tool.cc:608
lar_cluster3d::SnippetHit3DBuilderICARUS::m_2hit1stPHVec
std::vector< float > m_2hit1stPHVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:304
lar_cluster3d::SnippetHit3DBuilderICARUS::m_rangeNumSig
float m_rangeNumSig
Definition: SnippetHit3DBuilderICARUS_tool.cc:252
lar_cluster3d::SnippetHit3DBuilderICARUS::m_enableMonitoring
bool m_enableMonitoring
Definition: SnippetHit3DBuilderICARUS_tool.cc:265
lar_cluster3d::SnippetHit3DBuilderICARUS::m_overlapRangeVec
std::vector< float > m_overlapRangeVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:295
lar_cluster3d::SnippetHit3DBuilderICARUS::m_chiSquare3DVec
std::vector< float > m_chiSquare3DVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:292
reco::ClusterHit3D::getPosition
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:157
lar_cluster3d::SetHitTimeOrder
bool SetHitTimeOrder(const reco::ClusterHit2D *left, const reco::ClusterHit2D *right)
Definition: SnippetHit3DBuilderICARUS_tool.cc:1741
lar_cluster3d::SnippetHit3DBuilderICARUS::m_2hitSumPHVec
std::vector< float > m_2hitSumPHVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:307
recob::Hit::RMS
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:220
lar_cluster3d::Hit2DSetCompare
What follows are several highly useful typedefs which we want to expose to the outside world...
Definition: SnippetHit3DBuilderICARUS_tool.cc:59
Hit.h
Declaration of signal hit object.
lar_cluster3d::SnippetHit3DBuilderICARUS::getTimeToExecute
virtual float getTimeToExecute(IHit3DBuilder::TimeValues index) const override
If monitoring, recover the time to execute a particular function.
Definition: SnippetHit3DBuilderICARUS_tool.cc:110
right
walls no right
Definition: selectors.fcl:105
lar_cluster3d::SnippetHit3DBuilderICARUS::BuildHit3D
void BuildHit3D(reco::HitPairList &hitPairList) const
Given the ClusterHit2D objects, build the HitPairMap.
Definition: SnippetHit3DBuilderICARUS_tool.cc:579
lar_cluster3d::SnippetHit3DBuilderICARUS::m_maxDeltaPeakVec
std::vector< float > m_maxDeltaPeakVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:296
lar_cluster3d::SnippetHit3DBuilderICARUS::m_overlapFractionVec
std::vector< float > m_overlapFractionVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:294
lar_cluster3d::SnippetHit3DBuilderICARUS::makeWireAssns
void makeWireAssns(const art::Event &, art::Assns< recob::Wire, recob::Hit > &, RecobHitToPtrMap &) const
Create recob::Wire to recob::Hit associations.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1962
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
reco::ClusterHit2D::WireID
const geo::WireID & WireID() const
Definition: Cluster3D.h:76
lar_cluster3d::Hit2DSetCompare::operator()
bool operator()(const reco::ClusterHit2D *, const reco::ClusterHit2D *) const
Definition: SnippetHit3DBuilderICARUS_tool.cc:1747
lar_cluster3d::SnippetHit3DBuilderICARUS::BuildHitPairMap
size_t BuildHitPairMap(PlaneToSnippetHitMap &planeToHitVectorMap, reco::HitPairList &hitPairList) const
Given the ClusterHit2D objects, build the HitPairMap.
Definition: SnippetHit3DBuilderICARUS_tool.cc:635
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaTime2Vec
std::vector< float > m_deltaTime2Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:290
lar_cluster3d::SnippetHit3DBuilderICARUS::BuildChannelStatusVec
void BuildChannelStatusVec(PlaneToWireToHitSetMap &planeToWiretoHitSetMap) const
Create the internal channel status vector (assume will eventually be event-by-event) ...
Definition: SnippetHit3DBuilderICARUS_tool.cc:453
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaPeakSigmaPlane2Vec
std::vector< float > m_deltaPeakSigmaPlane2Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:283
recob::Hit::DegreesOfFreedom
int DegreesOfFreedom() const
Definition: Hit.h:229
geo::GeometryCore::ChannelToWire
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition: GeometryCore.cxx:1196
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
lar_cluster3d::SnippetHit3DBuilderICARUS::configure
virtual void configure(const fhicl::ParameterSet &) override
Interface for configuring the particular algorithm tool.
Definition: SnippetHit3DBuilderICARUS_tool.cc:342
lar_cluster3d::SnippetHit3DBuilderICARUS::m_zPosOffset
float m_zPosOffset
Definition: SnippetHit3DBuilderICARUS_tool.cc:269
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaSigma0Vec
std::vector< float > m_deltaSigma0Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:287
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
lar_cluster3d::SnippetHit3DBuilderICARUS::m_geometry
const geo::Geometry * m_geometry
Definition: SnippetHit3DBuilderICARUS_tool.cc:320
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
reco::ClusterHit3D::getSigmaPeakTime
float getSigmaPeakTime() const
Definition: Cluster3D.h:164
lar_cluster3d::SnippetHit3DBuilderICARUS::m_2hitDeltaPHVec
std::vector< float > m_2hitDeltaPHVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:306
lar_cluster3d::SnippetHit3DBuilderICARUS::m_weHaveAllBeenHereBefore
bool m_weHaveAllBeenHereBefore
Definition: SnippetHit3DBuilderICARUS_tool.cc:318
hit
process_name hit
Definition: cheaterreco.fcl:51
lar_cluster3d::PlaneToWireToHitSetMap
std::map< geo::PlaneID, WireToHitSetMap > PlaneToWireToHitSetMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:72
geo::GeometryCore::Ncryostats
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:2002
lar_cluster3d::SnippetHit3DBuilderICARUS::m_planeToSnippetHitMap
PlaneToSnippetHitMap m_planeToSnippetHitMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:311
lar_cluster3d::SnippetHit3DBuilderICARUS::findGoodHitPairs
int findGoodHitPairs(SnippetHitMap::iterator &, SnippetHitMap::iterator &, SnippetHitMap::iterator &, HitMatchTripletVecMap &) const
Definition: SnippetHit3DBuilderICARUS_tool.cc:788
lar_cluster3d::SnippetHit3DBuilderICARUS
SnippetHit3DBuilderICARUS class definiton.
Definition: SnippetHit3DBuilderICARUS_tool.cc:81
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaSigma2Vec
std::vector< float > m_deltaSigma2Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:291
lar_cluster3d::SnippetHit3DBuilderICARUS::m_minPHFor2HitPoints
float m_minPHFor2HitPoints
Set a minimum pulse height for 2 hit space point candidates.
Definition: SnippetHit3DBuilderICARUS_tool.cc:256
lar_cluster3d::SnippetHit3DBuilderICARUS::ChannelStatusByPlaneVec
std::vector< ChannelStatusVec > ChannelStatusByPlaneVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:237
reco::ClusterHit3D::getStatusBits
unsigned int getStatusBits() const
Definition: Cluster3D.h:156
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaPeakSigmaPlane1Vec
std::vector< float > m_deltaPeakSigmaPlane1Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:281
lar_cluster3d::IHit3DBuilder::NUMTIMEVALUES
Definition: IHit3DBuilder.h:78
lar_cluster3d::SnippetHit3DBuilderICARUS::m_numBadChannels
size_t m_numBadChannels
Definition: SnippetHit3DBuilderICARUS_tool.cc:316
geo::GeometryCore::Nwires
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
Definition: GeometryCore.h:3345
lar_cluster3d::SnippetHit3DBuilderICARUS::m_wirePitch
float m_wirePitch[3]
Definition: SnippetHit3DBuilderICARUS_tool.cc:266
lar_cluster3d::SetPeakHitPairIteratorOrder
bool SetPeakHitPairIteratorOrder(const reco::HitPairList::iterator &left, const reco::HitPairList::iterator &right)
Definition: SnippetHit3DBuilderICARUS_tool.cc:485
lar_cluster3d::SnippetHit3DBuilderICARUS::closestApproach
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...
Definition: SnippetHit3DBuilderICARUS_tool.cc:1495
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaTime1Vec
std::vector< float > m_deltaTime1Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:288
low
standard_dbscan3dalg useful for diagnostics hits not in a line will not be clustered on on only for track like only for track like on on the smaller the less shower like tracks low
Definition: clusteralgorithms.fcl:253
lar_cluster3d::SnippetHit3DBuilderICARUS::m_channelStatus
ChannelStatusByPlaneVec m_channelStatus
Definition: SnippetHit3DBuilderICARUS_tool.cc:315
reco::ClusterHit2D::getHit
const recob::Hit * getHit() const
Definition: Cluster3D.h:77
geo::GeometryCore::Nchannels
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
Definition: GeometryCore.cxx:208
recob::Hit::PeakAmplitude
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaSigma1Vec
std::vector< float > m_deltaSigma1Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:289
geo::fhicl::PlaneID
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
Definition: geo_types_fhicl.h:375
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
lar_cluster3d::PlaneSnippetHitMapItrPairVec
std::vector< SnippetHitMapItrPair > PlaneSnippetHitMapItrPairVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:76
lariov::ChannelStatusProvider::Status_t
unsigned short Status_t
type representing channel status
Definition: ChannelStatusProvider.h:53
lar_cluster3d::IHit3DBuilder::COLLECTARTHITS
Definition: IHit3DBuilder.h:75
lar_cluster3d::SnippetHit3DBuilderICARUS::m_hitFinderTagVec
std::vector< art::InputTag > m_hitFinderTagVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:249
lar_cluster3d::SnippetHit3DBuilderICARUS::m_pairWireDistVec
std::vector< float > m_pairWireDistVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:298
lar_cluster3d::SnippetHit3DBuilderICARUS::FindBestMatchingHit
const reco::ClusterHit2D * FindBestMatchingHit(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.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1615
lar_cluster3d::SnippetHit3DBuilderICARUS::HitMatchTripletVecMap
std::map< geo::WireID, HitMatchTripletVec > HitMatchTripletVecMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:156
reco::ClusterHit3D::getAvePeakTime
float getAvePeakTime() const
Definition: Cluster3D.h:162
lar_cluster3d::SnippetHit3DBuilderICARUS::makeRawDigitAssns
void makeRawDigitAssns(const art::Event &, art::Assns< raw::RawDigit, recob::Hit > &, RecobHitToPtrMap &) const
Create raw::RawDigit to recob::Hit associations.
Definition: SnippetHit3DBuilderICARUS_tool.cc:2008
lar_cluster3d::HitVectorMap
std::map< size_t, HitVector > HitVectorMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:74
recob::HitCreator
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:83
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
HitCreator.h
Helper functions to create a hit.
lar_cluster3d::SnippetHit3DBuilderICARUS::produces
virtual void produces(art::ProducesCollector &) override
Each algorithm may have different objects it wants &quot;produced&quot; so use this to let the top level produc...
Definition: SnippetHit3DBuilderICARUS_tool.cc:333
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaPeakTimePlane2Vec
std::vector< float > m_deltaPeakTimePlane2Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:282
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2854
lar_cluster3d::SnippetHit3DBuilderICARUS::m_clusterHit2DMasterList
Hit2DList m_clusterHit2DMasterList
Definition: SnippetHit3DBuilderICARUS_tool.cc:310
lar_cluster3d::Hit2DList
std::list< reco::ClusterHit2D > Hit2DList
Definition: SnippetHit3DBuilderICARUS_tool.cc:69
lar_cluster3d::SnippetHit3DBuilderICARUS::m_pulseHeightFrac
float m_pulseHeightFrac
Definition: SnippetHit3DBuilderICARUS_tool.cc:254
lar_cluster3d::SnippetHit3DBuilderICARUS::m_maxPullVec
std::vector< float > m_maxPullVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:293
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
lar_cluster3d::SnippetHit3DBuilderICARUS::m_qualityMetricVec
std::vector< float > m_qualityMetricVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:301
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
lar_cluster3d::SnippetHit3DBuilderICARUS::m_spacePointChargeVec
std::vector< float > m_spacePointChargeVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:302
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaPeakSigmaPlane0Vec
std::vector< float > m_deltaPeakSigmaPlane0Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:279
lar_cluster3d::SnippetHit3DBuilderICARUS::m_hitWidthSclFctr
float m_hitWidthSclFctr
Definition: SnippetHit3DBuilderICARUS_tool.cc:250
lar_cluster3d::SnippetHit3DBuilderICARUS::m_outputHistograms
bool m_outputHistograms
Take the time to create and fill some histograms for diagnostics.
Definition: SnippetHit3DBuilderICARUS_tool.cc:263
lar_cluster3d::SnippetHit3DBuilderICARUS::BuildHitPairMapByTPC
size_t BuildHitPairMapByTPC(PlaneSnippetHitMapItrPairVec &planeSnippetHitMapItrPairVec, reco::HitPairList &hitPairList) const
Given the ClusterHit2D objects, build the HitPairMap.
Definition: SnippetHit3DBuilderICARUS_tool.cc:694
lar_cluster3d::SnippetHit3DBuilderICARUS::ChannelStatusVec
std::vector< size_t > ChannelStatusVec
define data structure for keeping track of channel status
Definition: SnippetHit3DBuilderICARUS_tool.cc:236
lariov::ChannelStatusProvider::Status
virtual Status_t Status(raw::ChannelID_t channel) const
Returns a status integer with arbitrary meaning.
Definition: ChannelStatusProvider.h:89
geo::GeometryCore::DetLength
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
Definition: GeometryCore.cxx:785
lariov::ChannelStatusProvider::IsGood
virtual bool IsGood(raw::ChannelID_t channel) const
Returns whether the specified channel is physical and good.
Definition: ChannelStatusProvider.h:84
lar_cluster3d::SnippetHit3DBuilderICARUS::m_channelFilter
const lariov::ChannelStatusProvider * m_channelFilter
Definition: SnippetHit3DBuilderICARUS_tool.cc:321
lar_cluster3d::SnippetHit3DBuilderICARUS::m_useT0Offsets
bool m_useT0Offsets
If true then we will use the LArSoft interplane offsets.
Definition: SnippetHit3DBuilderICARUS_tool.cc:262
left
walls no left
Definition: selectors.fcl:105
lar_cluster3d::SnippetHit3DBuilderICARUS::m_wirePitchScaleFactor
float m_wirePitchScaleFactor
Scaling factor to determine max distance allowed between candidate pairs.
Definition: SnippetHit3DBuilderICARUS_tool.cc:258
lar_cluster3d::SnippetHit3DBuilderICARUS::WireIDsIntersect
bool WireIDsIntersect(const geo::WireID &, const geo::WireID &, geo::WireIDIntersection &) const
function to detemine if two wires &quot;intersect&quot; (in the 2D sense)
Definition: SnippetHit3DBuilderICARUS_tool.cc:1442
reco::ClusterHit2D::getXPosition
float getXPosition() const
Definition: Cluster3D.h:74
lar_cluster3d::IHit3DBuilder::TimeValues
TimeValues
enumerate the possible values for time checking if monitoring timing
Definition: IHit3DBuilder.h:75
lariov::ChannelStatusProvider
Class providing information about the quality of channels.
Definition: ChannelStatusProvider.h:49
lar_cluster3d::SnippetHit3DBuilderICARUS::CollectArtHits
void CollectArtHits(const art::Event &evt) const
Extract the ART hits and the ART hit-particle relationships.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1753
geo::Geometry
The geometry of one entire detector, as served by art.
Definition: Geometry.h:181
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
lar_cluster3d::SnippetHit3DBuilderICARUS::DistanceFromPointToHitWire
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...
Definition: SnippetHit3DBuilderICARUS_tool.cc:1692
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaPeakTimeSig
float m_deltaPeakTimeSig
Definition: SnippetHit3DBuilderICARUS_tool.cc:251
lar_cluster3d::IHit3DBuilder
IHit3DBuilder interface class definiton.
Definition: IHit3DBuilder.h:38
geo::WireGeo::Direction
Vector Direction() const
Definition: WireGeo.h:587
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
IHit3DBuilder.h
This provides an art tool interface definition for tools which construct 3D hits used in 3D clusterin...
reco::ClusterHit2DVec
std::vector< const reco::ClusterHit2D * > ClusterHit2DVec
Definition: Cluster3D.h:91
lar_cluster3d::Hit2DSet
std::set< const reco::ClusterHit2D *, Hit2DSetCompare > Hit2DSet
Definition: SnippetHit3DBuilderICARUS_tool.cc:70
geo::GeometryCore::NTPC
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:2311
lar_cluster3d::TPCToPlaneToWireToHitSetMap
std::map< geo::TPCID, PlaneToWireToHitSetMap > TPCToPlaneToWireToHitSetMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:73
lar_cluster3d::HitVector
std::vector< const reco::ClusterHit2D * > HitVector
Definition: SnippetHit3DBuilderICARUS_tool.cc:64
geo::WireGeo::HalfL
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:128
lar_cluster3d::SnippetHit3DBuilderICARUS::HitMatchTripletVec
std::vector< HitMatchTriplet > HitMatchTripletVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:155
geo::GeometryCore::WireIDsIntersect
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
Definition: GeometryCore.cxx:1512
lar_cluster3d::SnippetHit3DBuilderICARUS::m_maxSideVecVec
std::vector< float > m_maxSideVecVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:297
lar_cluster3d::HitPairClusterOrder::operator()
bool operator()(const reco::HitPairClusterMap::iterator &left, const reco::HitPairClusterMap::iterator &right)
Definition: SnippetHit3DBuilderICARUS_tool.cc:492
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
lar_cluster3d::SetStartTimeOrder::operator()
bool operator()(const SnippetHitMapItrPair &left, const SnippetHitMapItrPair &right) const
Definition: SnippetHit3DBuilderICARUS_tool.cc:613
reco::ClusterHit3D::setID
void setID(const size_t &id) const
Definition: Cluster3D.h:178
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaTimeVec
std::vector< float > m_deltaTimeVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:285
lar_cluster3d::PlaneToSnippetHitMap
std::map< geo::PlaneID, SnippetHitMap > PlaneToSnippetHitMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:67
reco::ClusterHit3D::MADESPACEPOINT
Hit has been made into Space Point.
Definition: Cluster3D.h:102
geo::WireIDIntersection::y
double y
y position of intersection
Definition: geo_types.h:804
lar_cluster3d::TPCToPlaneToSnippetHitMap
std::map< geo::TPCID, PlaneToSnippetHitMap > TPCToPlaneToSnippetHitMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:68
lar_cluster3d::SnippetHit3DBuilderICARUS::m_2hit2ndPHVec
std::vector< float > m_2hit2ndPHVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:305
lar_cluster3d::SnippetHit3DBuilderICARUS::m_maxHit3DChiSquare
float m_maxHit3DChiSquare
Provide ability to select hits based on &quot;chi square&quot;.
Definition: SnippetHit3DBuilderICARUS_tool.cc:259
lar_cluster3d::SnippetHit3DBuilderICARUS::NearestWireID
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...
Definition: SnippetHit3DBuilderICARUS_tool.cc:1667
lar_cluster3d::SnippetHit3DBuilderICARUS::findGoodTriplets
void findGoodTriplets(HitMatchTripletVecMap &, HitMatchTripletVecMap &, reco::HitPairList &, bool=false) const
This algorithm takes lists of hit pairs and finds good triplets.
Definition: SnippetHit3DBuilderICARUS_tool.cc:842
ChannelStatusProvider.h
Interface for experiment-specific channel quality info provider.
e
do i e
Definition: updateVersionList.sh:68
copy
T copy(T const &v)
Definition: RangeForWrapper_test.cc:98
lar_cluster3d::SnippetHitMapItrPair
std::pair< SnippetHitMap::iterator, SnippetHitMap::iterator > SnippetHitMapItrPair
Definition: SnippetHit3DBuilderICARUS_tool.cc:75
lar_cluster3d::HitPairClusterOrder
Definition: SnippetHit3DBuilderICARUS_tool.cc:490
detinfo::trigger_offset
int trigger_offset(DetectorClocksData const &data)
Definition: DetectorClocksData.h:632
geo::WireID::WireID
constexpr WireID()=default
Default constructor: an invalid TPC ID.
lar_cluster3d::HitStartEndPair
std::pair< raw::TDCtick_t, raw::TDCtick_t > HitStartEndPair
Definition: SnippetHit3DBuilderICARUS_tool.cc:65
lar_cluster3d::IHit3DBuilder::BUILDTHREEDHITS
Definition: IHit3DBuilder.h:76
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaPeakTimePlane0Vec
std::vector< float > m_deltaPeakTimePlane0Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:278
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
reco::ClusterHit3D::getHitChiSquare
float getHitChiSquare() const
Definition: Cluster3D.h:165
lar_cluster3d::SnippetHit3DBuilderICARUS::m_planeToWireToHitSetMap
PlaneToWireToHitSetMap m_planeToWireToHitSetMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:312
geo::WireGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
lar_cluster3d::SnippetHit3DBuilderICARUS::chargeIntegral
float chargeIntegral(float, float, float, float, int, int) const
Perform charge integration between limits.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1520
geo::PlaneGeo::WireCoordinate
double WireCoordinate(Point const &point) const
Returns the coordinate of the point on the plane, in wire units.
Definition: PlaneGeo.h:882
lar_cluster3d::SnippetHit3DBuilderICARUS::makeDeadChannelPair
bool makeDeadChannelPair(reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pair, size_t maxStatus=4, size_t minStatus=0, float minOverlap=0.2) const
Make a 3D HitPair object from a valid pair and a dead channel in the missing plane.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1538
lar_cluster3d::SnippetHitMap
std::map< HitStartEndPair, HitVector > SnippetHitMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:66
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
ChannelStatusService.h
Interface for experiment-specific service for channel quality info.
P0
process_name physics producers generator hPHist_pi physics producers generator P0
Definition: icaruscode/fcl/gen/single/prodsingle_electron_pi+_bnblike_forward.fcl:10
lar_cluster3d::SetStartTimeOrder::SetStartTimeOrder
SetStartTimeOrder()
Definition: SnippetHit3DBuilderICARUS_tool.cc:611
lar_cluster3d::SnippetHit3DBuilderICARUS::HitMatchTriplet
std::tuple< const reco::ClusterHit2D *, const reco::ClusterHit2D *, const reco::ClusterHit3D > HitMatchTriplet
This builds a list of candidate hit pairs from lists of hits on two planes.
Definition: SnippetHit3DBuilderICARUS_tool.cc:154
lar_cluster3d::IHit3DBuilder::BUILDNEWHITS
Definition: IHit3DBuilder.h:77
reco::ClusterHit3D::getHits
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:170
lar_cluster3d::IHit3DBuilder::RecobHitToPtrMap
std::unordered_map< const recob::Hit *, art::Ptr< recob::Hit >> RecobHitToPtrMap
Defines a structure mapping art representation to internal.
Definition: IHit3DBuilder.h:62
lar_cluster3d::SnippetHit3DBuilderICARUS::m_smallChargeDiffVec
std::vector< float > m_smallChargeDiffVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:299
lar_cluster3d::SnippetHit3DBuilderICARUS::m_saveMythicalPoints
bool m_saveMythicalPoints
Should we save valid 2 hit space points?
Definition: SnippetHit3DBuilderICARUS_tool.cc:260
reco::ClusterHit3D::setPosition
void setPosition(const Eigen::Vector3f &pos) const
Definition: Cluster3D.h:185
lar_cluster3d::SnippetHit3DBuilderICARUS::TickCorrectionArray
std::vector< std::vector< std::vector< float >>> TickCorrectionArray
Data members to follow.
Definition: SnippetHit3DBuilderICARUS_tool.cc:247
lar_cluster3d::SetPairStartTimeOrder
bool SetPairStartTimeOrder(const reco::ClusterHit3D &left, const reco::ClusterHit3D &right)
Definition: SnippetHit3DBuilderICARUS_tool.cc:626
lar_cluster3d::SnippetHit3DBuilderICARUS::FindNumberInRange
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.
Definition: SnippetHit3DBuilderICARUS_tool.cc:1642
lar_cluster3d::SnippetHit3DBuilderICARUS::m_invalidTPCVec
std::vector< int > m_invalidTPCVec
Definition: SnippetHit3DBuilderICARUS_tool.cc:257
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
lar_cluster3d::SnippetHit3DBuilderICARUS::m_PHLowSelection
float m_PHLowSelection
Definition: SnippetHit3DBuilderICARUS_tool.cc:255
lar_cluster3d::SnippetHit3DBuilderICARUS::m_maxMythicalChiSquare
float m_maxMythicalChiSquare
Selection cut on mythical points.
Definition: SnippetHit3DBuilderICARUS_tool.cc:261
lar_cluster3d::WireToHitSetMap
std::map< unsigned int, Hit2DSet > WireToHitSetMap
Definition: SnippetHit3DBuilderICARUS_tool.cc:71
lar_cluster3d::SnippetHit3DBuilderICARUS::m_deltaTime0Vec
std::vector< float > m_deltaTime0Vec
Definition: SnippetHit3DBuilderICARUS_tool.cc:286
geo::GeometryCore::WireIDToWireGeo
WireGeo const & WireIDToWireGeo(geo::WireID const &wireid) const
Definition: GeometryCore.h:3424
lar_cluster3d::SnippetHit3DBuilderICARUS::saveOrphanPairs
int saveOrphanPairs(HitMatchTripletVecMap &, reco::HitPairList &) const
This will look at storing pair &quot;orphans&quot; where the 2D hits are otherwise unused.
Definition: SnippetHit3DBuilderICARUS_tool.cc:957
reco::ClusterHit2D::getStatusBits
unsigned getStatusBits() const
Definition: Cluster3D.h:71
reco::ClusterHit3D::setStatusBit
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:179
reco::ClusterHit2D::setStatusBit
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:79
Generated by   doxygen 1.8.5