LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
FeatureVertexFinder_module.cc
Go to the documentation of this file.
1 //
3 // FeatureVertexFinder class
4 //
5 // jasaadi@fnal.gov
6 //
7 // This algorithm is designed to reconstruct the vertices using the
8 // 2D cluster information and CornerFinderAlg to find 2-d and 3-d verticies
9 //
10 // 06/05/14
11 // This is a major rewrite of the original code and will be factored to take
12 // advantage of improvements in LArSoft. The code now goes something like this:
13 //
14 // 1)Take up EndPoint2d's from Cluster Crawler and keep ones that make sense in 3d
15 // (this algorithm wants ClusterCrawler to have been run, but should totally work
16 // even if it hasn't)
17 //
18 // 2)Take up all EndPoint2d's from the image finding CornerFinder algorithm
19 // (this algortithm wants CornerFinder to have been run, but should totally work
20 // even if it hasn't)
21 //
22 // 3)Take up an additional clustering algorithm (dBcluster by default but can work
23 // with any clustering algorithm) and use slopes and endpoints to find intersections
24 //
25 // 4)Now merge together any duplicates and sort them by their weight (which needs its units still
26 // decided somehow)
27 //
28 // 5)Return a list of verticies based on what mode you run this algorithm in
29 // a) Primary mode: This mode will return only a single 3d vertex and an appropriate
30 // number of EndPoint2d's corresponding to the most likely neutrino vertex
31 // b) All mode: This mode returns all 3d/2d verticies that this algorithm has found but
32 // sorts the list by the most likely candidate
33 //
34 //
35 //
36 // This is Preliminary Work and needs modifications
37 //
38 // ////////////////////////////////////////////////////////////////////////
39 
40 // Basic C++ Includes
41 #include <cmath>
42 #include <iomanip>
43 #include <iostream>
44 #include <string>
45 
46 // Framework Includes
55 #include "fhiclcpp/ParameterSet.h"
57 
58 // LArSoft Includes
70 
71 // ROOT Includes
72 #include "TF1.h"
73 #include "TGraph.h"
74 
75 namespace vertex {
76 
78  public:
79  explicit FeatureVertexFinder(fhicl::ParameterSet const& pset);
80 
81  private:
82  void produce(art::Event& evt) override;
83 
84  // This function will take in and EndPoint2d from either cluster crawler
85  // or corner finder and only save points that make a 3d-candidate
87  std::vector<art::Ptr<recob::EndPoint2D>> const& EndPoints,
88  bool PlaneDet);
89 
90  // This function will take in 2d Clusters (by default dBCluster)
91  // and return 2d vertex candidates for sorting later
95 
96  // This function takes in 2d Vertex candidates (found by
97  // Find2dClusterVertexCandidates), does some basic merging and
98  // then finds 3d-candidates for use later
100  std::vector<double> const& Wire_2dvtx,
101  std::vector<double> const& Time_2dvtx,
102  std::vector<double> const& Plane_2dvtx);
103 
104  // This function merges and sorts all the 3d vertex candidates, it
105  // merges them if they are within 2.0 cm in X, Y, and Z
106  // simultaneously and then sorts the list by vertex strength and Z
107  // location (lowest Z is first on the list
108  void MergeAndSort3dVtxCandidate(std::vector<double> const& merge_vtxX,
109  std::vector<double> const& merge_vtxY,
110  std::vector<double> const& merge_vtxZ,
111  std::vector<double> const& merge_vtxStgth);
112 
113  std::string fClusterModuleLabel;
114  std::string fHitModuleLabel;
117  Double_t fRunningMode;
118 
119  bool GT2PlaneDetector = false;
120 
121  // Unsorted Raw list of 3d-Vertex candidates filled
122 
123  std::vector<double> candidate_x = {0.};
124  std::vector<double> candidate_y = {0.};
125  std::vector<double> candidate_z = {0.};
126  std::vector<double> candidate_strength = {0.};
127 
128  // Merged and sorted list of 3d-Vertex candidates filled
129 
130  std::vector<double> MergeSort3dVtx_xpos = {0.};
131  std::vector<double> MergeSort3dVtx_ypos = {0.};
132  std::vector<double> MergeSort3dVtx_zpos = {0.};
133  std::vector<double> MergeSort3dVtx_strength = {0.};
134 
135  // 2d Cluster based Variables
136 
137  std::vector<std::vector<int>> Cls; //<---- Index to clusters in each view
138  std::vector<double> dtdwstart = {0.}; //<----Slope (delta Time Tick vs delta Wire)
139  std::vector<double> dtdwstartError = {0.}; //<---Error on the slope
140  std::vector<double> Clu_Plane = {0.}; //<---Plane of the current cluster
141  std::vector<double> Clu_StartPos_Wire = {0.}; //<---Starting wire number of cluster
142  std::vector<double> Clu_StartPos_TimeTick = {0.}; //<---Starting TDC value of the cluster
143 
144  std::vector<double> Clu_EndPos_Wire = {0.}; //<---Ending wire number of cluster
145  std::vector<double> Clu_EndPos_TimeTick = {0.}; //<---Ending TDC value of the cluster
146 
147  std::vector<double> Clu_Slope = {0.}; //<---Calculated Slope of the cluster (TDC/Wire)
148  std::vector<double> Clu_Yintercept = {0.}; //<---Clusters Y Intercept using start positions
149  std::vector<double> Clu_Yintercept2 = {0.}; //<---Clusters Y Intercept using end positions
150  std::vector<double> Clu_Length = {0.}; //<---Calculated Length of the cluster
151 
152  // 2d Cluster based verticies
153 
154  std::vector<double> TwoDvtx_wire = {0.};
155  std::vector<double> TwoDvtx_time = {0.};
156  std::vector<double> TwoDvtx_plane = {0.};
157  };
158 
159 } //<---End namespace vertex
160 
161 namespace vertex {
162 
164  {
165  fCornerFinderModuleLabel = pset.get<std::string>("CornerFinderModuleLabel");
166  fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
167  fHitModuleLabel = pset.get<std::string>("HitModuleLabel");
168  fCCrawlerEndPoint2dModuleLabel = pset.get<std::string>("CCrawlerEndPoint2dModuleLabel");
169  fRunningMode = pset.get<double>("RunningMode");
170 
171  produces<std::vector<recob::Vertex>>();
172  produces<std::vector<recob::EndPoint2D>>();
173  produces<art::Assns<recob::EndPoint2D, recob::Hit>>();
174 
175  // === Don't think I'll actually have any of these associations ===
176  // === should consider removing in the future
177  produces<art::Assns<recob::Vertex, recob::Hit>>();
178  produces<art::Assns<recob::Vertex, recob::Shower>>();
179  produces<art::Assns<recob::Vertex, recob::Track>>();
180 
182  Cls.resize(geom->Nplanes(), std::vector<int>());
183  }
184 
185  // -----------------------------------------------------------------------------
186  // Produce
188  {
190  auto const detProp =
192 
193  // Figuring out if I have a 2 or 3 plane detector
194 
195  GT2PlaneDetector = false;
196 
197  for (auto const& tpc : geom->Iterate<geo::TPCGeo>()) {
198  if (tpc.Nplanes() > 2) {
199  GT2PlaneDetector = true;
200  std::cout << "yeah" << std::endl;
201  break;
202  }
203  }
204 
205  // These are the things I want to put on the event
206 
207  auto vcol = std::make_unique<std::vector<recob::Vertex>>(); // 3D vertex
208  auto epcol = std::make_unique<std::vector<recob::EndPoint2D>>(); // 2D vertex
209  auto assnep = std::make_unique<art::Assns<recob::EndPoint2D, recob::Hit>>();
210  auto assnsh = std::make_unique<art::Assns<recob::Vertex, recob::Shower>>();
211  auto assntr = std::make_unique<art::Assns<recob::Vertex, recob::Track>>();
212  auto assnh = std::make_unique<art::Assns<recob::Vertex, recob::Hit>>();
213 
214  // ClusterCrawler
215 
216  // Getting the EndPoint2d's for cluster crawler
217 
218  try {
219  art::Handle<std::vector<recob::EndPoint2D>> ccrawlerFinderHandle;
220  evt.getByLabel(fCCrawlerEndPoint2dModuleLabel, ccrawlerFinderHandle);
221  std::vector<art::Ptr<recob::EndPoint2D>> ccrawlerEndPoints;
222  art::fill_ptr_vector(ccrawlerEndPoints, ccrawlerFinderHandle);
223 
224  // Passing in the EndPoint2d's from Cluster Crawler
225  Get3dVertexCandidates(detProp, ccrawlerEndPoints, GT2PlaneDetector);
226  }
227  catch (...) {
228  mf::LogWarning("FeatureVertexFinder") << "Failed to get EndPoint2d's from Cluster Crawler";
229  }
230 
231  // CornerFinder EndPoint2d's
232 
233  // Getting the EndPoint2d's for Corner Finder
234  try {
235  art::Handle<std::vector<recob::EndPoint2D>> CornerFinderHandle;
236  evt.getByLabel(fCornerFinderModuleLabel, CornerFinderHandle);
237  std::vector<art::Ptr<recob::EndPoint2D>> cornerEndPoints;
238  art::fill_ptr_vector(cornerEndPoints, CornerFinderHandle);
239 
240  // Passing in the EndPoint2d's from Corner Finder
241  Get3dVertexCandidates(detProp, cornerEndPoints, GT2PlaneDetector);
242  }
243  catch (...) {
244  mf::LogWarning("FeatureVertexFinder") << "Failed to get EndPoint2d's from Corner Finder";
245  }
246 
247  // Making Cluster Slope EndPoint2d's
248 
249  // Retrieving the Cluster Module for the event
250  try {
251  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
252  evt.getByLabel(fClusterModuleLabel, clusterListHandle);
253 
254  // Finding hits associated with the clusters
255  art::FindManyP<recob::Hit> fmh(clusterListHandle, evt, fClusterModuleLabel);
256 
257  // Filling the Cluster Vector
259  for (unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
260  art::Ptr<recob::Cluster> clusterHolder(clusterListHandle, ii);
261  clusters.push_back(clusterHolder);
262  }
263 
264  // Passing in the clusters to find 2d Vertices
265  Find2dClusterVertexCandidates(detProp, clusters, fmh);
266 
267  // Finding 3d Candidates from 2d cluster vertex candidates
269  }
270  catch (...) {
271  mf::LogWarning("FeatureVertexFinder") << "Failed to get Cluster from default cluster module";
272  }
273 
274  // Merging and sorting 3d vertex candidates
276 
277  // Putting Vertices on the event
278  // ------------------------------------------------------------
279  // Now that I have a list of 3d vertex candidates I will return
280  // 3d/2d vertices based on which option the user has chosen
281  // fRunningMode == 0 (this returns a full list of all 3d/2d vertex
282  // candidates) fRunningMode == 1 (this returns only one vertex
283  // which is established as the most likely primary)
284 
285  // Returning all vertex candidates
286  if (fRunningMode == 0) {
287 
288  // Looping over Primary Verticies
289  for (size_t pri = 0; pri < MergeSort3dVtx_xpos.size(); pri++) {
290 
291  // Push each primary vertex onto the event
292 
293  double tempxyz[3] = {
295 
296  // Skipping a vertex that is zero
297 
298  if (tempxyz[0] == 0 && tempxyz[1] == 0 && tempxyz[2] == 0) { continue; }
299  recob::Vertex the3Dvertex(tempxyz, vcol->size());
300  vcol->push_back(the3Dvertex);
301 
302  // ---------------------------------------------------------------------
303  // --- Now go make the 2DEndPoints that correspond to each 3d vertex ---
304  // ---------------------------------------------------------------------
305 
306  for (auto const& planeID : geom->Iterate<geo::PlaneID>()) {
307  geo::Point_t const temp2dXYZ{
308  MergeSort3dVtx_xpos[pri], MergeSort3dVtx_ypos[pri], MergeSort3dVtx_zpos[pri]};
309  double temp2dStrength = MergeSort3dVtx_strength[pri];
310 
311  // Skipping a vertex that is zero
312 
313  if (temp2dXYZ.X() == 0 && temp2dXYZ.Y() == 0 && temp2dXYZ.Z() == 0) { continue; }
314 
315  // Converting the 3d vertex into 2d time ticks, wire, and
316  // channel
317 
318  double EndPoint2d_TimeTick = detProp.ConvertXToTicks(temp2dXYZ.X(), planeID);
319  int EndPoint2d_Wire = 0;
320  int EndPoint2d_Channel = 0;
321  // Putting in protection in case NearestWire Fails
322  try {
323  EndPoint2d_Wire = geom->NearestWireID(temp2dXYZ, planeID).Wire;
324  }
325  catch (...) {
326  mf::LogWarning("FeatureVertexFinder") << "2dWire failed";
327  continue;
328  }
329  // Putting in protection in case NearestChannel Fails
330  try {
331  EndPoint2d_Channel = geom->NearestChannel(temp2dXYZ, planeID);
332  }
333  catch (...) {
334  mf::LogWarning("FeatureVertexFinder") << "2dWire failed";
335  continue;
336  }
337 
338  // Making geo::WireID and getting the current View number
339  geo::View_t View = geom->View(EndPoint2d_Channel);
340  geo::WireID wireID(planeID, EndPoint2d_Wire);
341 
342  // Putting the 2d Vertex found on the event
343 
344  recob::EndPoint2D vertex(EndPoint2d_TimeTick, //<---TimeTick
345  wireID, //<---geo::WireID
346  temp2dStrength, //<---Vtx strength (JA: ?)
347  epcol->size(), //<---Vtx ID (JA: ?)
348  View, //<---Vtx View
349  1); //<---Total Charge (JA: Need to figure this one?)
350  epcol->push_back(vertex);
351  } //<---End Plane loop
352  } //<---End pri loop
353  } //<---End fRunningMode == 0
354 
355  // ================================================
356  // === Returning only primary vertex candidates ===
357  // ================================================
358  if (fRunningMode != 0) {
359  int position = 0;
360  int bail = 0;
361 
362  // Looping over Primary Verticies
363 
364  for (size_t pri = 0; pri < MergeSort3dVtx_xpos.size(); pri++) {
365 
366  // Push each primary vertex onto the event
367 
368  double tempxyz[3] = {
370 
371  // Skipping a vertex that is zero
372 
373  if (bail > 0) { continue; }
374  if (tempxyz[0] == 0 && tempxyz[1] == 0 && tempxyz[2] == 0) { continue; }
375  position = pri;
376  bail++;
377  recob::Vertex the3Dvertex(tempxyz, vcol->size());
378  vcol->push_back(the3Dvertex);
379  }
380 
381  // ---------------------------------------------------------------------
382  // --- Now go make the 2DEndPoints that correspond to each 3d vertex ---
383  // ---------------------------------------------------------------------
384 
385  // Looping over cryostats
386 
387  for (auto const& planeID : geom->Iterate<geo::PlaneID>()) {
388  geo::Point_t const temp2dXYZ{MergeSort3dVtx_xpos[position],
389  MergeSort3dVtx_ypos[position],
390  MergeSort3dVtx_zpos[position]};
391  double temp2dStrength = MergeSort3dVtx_strength[position];
392 
393  // Converting the 3d vertex into 2d time ticks, wire, and
394  // channel
395  double EndPoint2d_TimeTick = detProp.ConvertXToTicks(temp2dXYZ.X(), planeID);
396  int EndPoint2d_Wire = 0;
397  int EndPoint2d_Channel = 0;
398  // Putting in protection in case NearestWire Fails
399  try {
400  EndPoint2d_Wire = geom->NearestWireID(temp2dXYZ, planeID).Wire;
401  }
402  catch (...) {
403  mf::LogWarning("FeatureVertexFinder") << "2dWire failed";
404  continue;
405  }
406  // Putting in protection in case NearestChannel Fails
407  try {
408  EndPoint2d_Channel = geom->NearestChannel(temp2dXYZ, planeID);
409  }
410  catch (...) {
411  mf::LogWarning("FeatureVertexFinder") << "2dWire failed";
412  continue;
413  }
414 
415  // Making geo::WireID and getting the current View number
416  geo::View_t View = geom->View(EndPoint2d_Channel);
417  geo::WireID wireID(planeID, EndPoint2d_Wire);
418 
419  // Putting the 2d Vertex found on the event
420 
421  recob::EndPoint2D vertex(EndPoint2d_TimeTick, //<---TimeTick
422  wireID, //<---geo::WireID
423  temp2dStrength, //<---Vtx strength (JA: ?)
424  epcol->size(), //<---Vtx ID (JA: ?)
425  View, //<---Vtx View
426  1); //<---Total Charge (JA: Need to figure this one?)
427  epcol->push_back(vertex);
428  } //<---End Plane loop
429  } //<---End fRunningMode == 1
430 
431  mf::LogVerbatim("Summary") << std::setfill('-') << std::setw(175) << "-" << std::setfill(' ');
432  mf::LogVerbatim("Summary") << "FeatureVertexFinder Summary:";
433  for (size_t i = 0; i < epcol->size(); ++i)
434  mf::LogVerbatim("Summary") << epcol->at(i);
435  for (size_t i = 0; i < vcol->size(); ++i)
436  mf::LogVerbatim("Summary") << vcol->at(i);
437 
438  evt.put(std::move(epcol));
439  evt.put(std::move(vcol));
440  evt.put(std::move(assnep));
441  evt.put(std::move(assntr));
442  evt.put(std::move(assnsh));
443  evt.put(std::move(assnh));
444 
445  // Clearing vectors at the end of the event
446 
447  vcol.reset();
448  epcol.reset();
449  candidate_x.clear();
450  candidate_y.clear();
451  candidate_z.clear();
452  candidate_strength.clear();
453  MergeSort3dVtx_xpos.clear();
454  MergeSort3dVtx_ypos.clear();
455  MergeSort3dVtx_zpos.clear();
456  MergeSort3dVtx_strength.clear();
457  dtdwstart.clear();
458  dtdwstartError.clear();
459  Clu_Plane.clear();
460  Clu_StartPos_Wire.clear();
461  Clu_StartPos_TimeTick.clear();
462  Clu_EndPos_Wire.clear();
463  Clu_EndPos_TimeTick.clear();
464  Clu_Slope.clear();
465  Clu_Yintercept.clear();
466  Clu_Yintercept2.clear();
467  Clu_Length.clear();
468  TwoDvtx_wire.clear();
469  TwoDvtx_time.clear();
470  TwoDvtx_plane.clear();
471 
472  } //<--End FeatureVertexFinder::produce
473 
474  // -----------------------------------------------------------------------------
475  // Get 3d Vertex Candidates
476  // -----------------------------------------------------------------------------
478  detinfo::DetectorPropertiesData const& detProp,
479  std::vector<art::Ptr<recob::EndPoint2D>> const& EndPoints,
480  bool PlaneDet)
481  {
483 
484  double y = 0., z = 0.;
485  double yy = 0., zz = 0.;
486  double yy2 = 0., zz2 = 0.;
487  double yy3 = 0., zz3 = 0.;
488 
489  auto const numEndPoints = size(EndPoints);
490 
491  for (auto const& tpcid : geom->Iterate<geo::TPCID>()) {
492  for (size_t iendpt1 = 0; iendpt1 < numEndPoints; ++iendpt1) {
493  for (size_t iendpt2 = iendpt1 + 1; iendpt2 < numEndPoints; ++iendpt2) {
494 
495  auto const& endpt1 = *EndPoints[iendpt1];
496  auto const& endpt2 = *EndPoints[iendpt2];
497 
498  geo::PlaneID const endpt1_planeid{tpcid, endpt1.WireID().Plane};
499  geo::PlaneID const endpt2_planeid{tpcid, endpt2.WireID().Plane};
500 
501  // Check to make sure we are comparing features from different planes
502  if (endpt1_planeid.Plane == endpt2_planeid.Plane) continue;
503 
504  // Get the appropriate time offset for the two planes we are considering
505  float tempXFeature1 = detProp.ConvertTicksToX(endpt1.DriftTime(), endpt1_planeid);
506  float tempXFeature2 = detProp.ConvertTicksToX(endpt2.DriftTime(), endpt2_planeid);
507 
508  // Skip features that are not within 0.5 cm in projected X
509  if (std::abs(tempXFeature1 - tempXFeature2) >= 0.5) continue;
510 
511  // Checking to see if these features have intersecting channels
512  geo::WireID const endpt1_wireid{endpt1_planeid, endpt1.WireID().Wire};
513  geo::WireID const endpt2_wireid{endpt2_planeid, endpt2.WireID().Wire};
514 
515  if (!geom->ChannelsIntersect(geom->PlaneWireToChannel(endpt1_wireid),
516  geom->PlaneWireToChannel(endpt2_wireid),
517  yy,
518  zz))
519  continue;
520 
521  // Use this fill if we are in a detector with fewer than 3 plane (e.g. ArgoNeuT)
522 
523  if (!PlaneDet) {
524  candidate_x.push_back(tempXFeature1);
525  candidate_y.push_back(yy);
526  candidate_z.push_back(zz);
527  candidate_strength.push_back(endpt1.Strength() + endpt2.Strength());
528  continue;
529  } //<---End fill for 2 plane detector
530 
531  // Now need to check for a match across more than 2 planes
532 
533  // Looping over the rest of the list
534 
535  for (size_t iendpt3 = iendpt2 + 1; iendpt3 < numEndPoints; ++iendpt3) {
536  auto const& endpt3 = *EndPoints[iendpt3];
537 
538  geo::PlaneID const endpt3_planeid{tpcid, endpt3.WireID().Plane};
539 
540  // Check to make sure we are comparing features from different planes
541  // N.B. We do not need to check between endpt1 and endpt2 as that has been done above.
542 
543  if (endpt3_planeid.Plane == endpt2_planeid.Plane ||
544  endpt3_planeid.Plane == endpt1_planeid.Plane)
545  continue;
546 
547  // Check to see that the third feature is within 1.0 cm projected in X
548  geo::WireID const endpt3_wireid{endpt3_planeid, endpt3.WireID().Wire};
549  float const tempXFeature3 = detProp.ConvertTicksToX(endpt3.DriftTime(), endpt3_wireid);
550 
551  if (std::abs(tempXFeature3 - tempXFeature2) >= 1.0 ||
552  std::abs(tempXFeature3 - tempXFeature1) >= 1.0) {
553  continue;
554  }
555 
556  // Make sure our third feature has an intersecting channel with our other two channels.
557  if (!geom->ChannelsIntersect(geom->PlaneWireToChannel(endpt3_wireid),
558  geom->PlaneWireToChannel(endpt1_wireid),
559  yy3,
560  zz3) ||
561  !geom->ChannelsIntersect(geom->PlaneWireToChannel(endpt3_wireid),
562  geom->PlaneWireToChannel(endpt2_wireid),
563  yy2,
564  zz2)) {
565  continue;
566  }
567  candidate_x.push_back(
568  detProp.ConvertTicksToX(endpt1.DriftTime(), endpt1_wireid.asPlaneID()));
569 
570  // Finding intersection points
571  geom->IntersectionPoint(endpt1_wireid, endpt2_wireid, y, z);
572 
573  candidate_y.push_back(y);
574  candidate_z.push_back(z);
575  candidate_strength.push_back(endpt1.Strength() + endpt2.Strength() + endpt3.Strength());
576 
577  // Note: If I've made it here I have a matched triplet...since I don't want to use any
578  // of these features again I am going to iterate each of the counters so we move to the
579  // next one.
580  if (iendpt1 < numEndPoints) { ++iendpt1; }
581  if (iendpt2 < numEndPoints) { ++iendpt2; }
582  if (iendpt3 < numEndPoints) { ++iendpt3; }
583  } //<---End iendpt3
584  } //<---End iendpt2 loop
585  } //<---End iendpt1 loop
586  } //<---End TPC loop
587  } //<---End Get3dVertexCandidates
588 
589  // -----------------------------------------------------------------------------
590  // Get 2d Vertex Candidates from clusters
591  // -----------------------------------------------------------------------------
593  detinfo::DetectorPropertiesData const& detProp,
594  art::PtrVector<recob::Cluster> RawClusters,
596  {
598 
599  int nClustersFound = 0;
600 
601  // Initialize Cls
602  for (auto& c : Cls)
603  c.clear();
604 
605  for (size_t iclu = 0; iclu < RawClusters.size(); ++iclu) {
606  // Gathering the hits associated with the current cluster
607  std::vector<art::Ptr<recob::Hit>> hit = fmhit.at(iclu);
608 
609  // I only want to consider this cluster if it has a sufficient number
610  // of hits
611  if (hit.size() < 15) { continue; }
612 
613  // Determine which view the current cluster is in
614 
615  const geo::View_t view = RawClusters.at(iclu)->View();
616  if (view >= Cls.size()) {
617  std::cerr << __PRETTY_FUNCTION__ << "\033[93m Logic error in my code ... view " << view
618  << " not supported ! \033[00m" << std::endl;
619  throw std::exception();
620  }
621 
622  Cls.at(RawClusters.at(iclu)->View()).push_back(iclu);
623 
624  // Filling wires and times into a TGraph for the cluster
625 
626  std::vector<double> wires;
627  std::vector<double> times;
628 
629  // Counting the number of hits in the current cluster (n)
630  int n = 0;
631 
632  // Looping over the hits in the cluster
633 
634  for (size_t i = 0; i < hit.size(); ++i) {
635  wires.push_back(hit[i]->WireID().Wire);
636  times.push_back(hit[i]->PeakTime());
637  ++n;
638  }
639 
640  // If there are 2 or more hits in the cluster fill a TGraph and
641  // fit a from a polynomial or order 1
642 
643  if (n >= 15) {
644 
645  // Push the hits into a TGraph
646 
647  TGraph* the2Dtrack = new TGraph(n, &wires[0], &times[0]);
648  // === Try to fit the TGraph with a 1st order polynomial ===
649  try {
650  the2Dtrack->Fit("pol1", "Q");
651  TF1* pol1 = (TF1*)the2Dtrack->GetFunction("pol1");
652  double par[2];
653  double parerror[2];
654  pol1->GetParameters(par);
655  parerror[1] = pol1->GetParError(1);
656  double FitChi2 = pol1->GetChisquare();
657  double FitNDF = pol1->GetNDF();
658 
659  double fitGoodness = FitChi2 / FitNDF;
660 
661  // Skipping the fitted slope if the Chi2/NDF < 5
662  if (fitGoodness > 10) {
663  dtdwstart.push_back(std::tan(RawClusters[iclu]->StartAngle()));
664  continue;
665  }
666 
667  // Take change in time tick vs change in wire (dT/dW) from the fit
668 
669  dtdwstart.push_back(par[1]);
670  dtdwstartError.push_back(parerror[1]);
671  } //<---End try to fit with a polynomial order 1
672 
673  // If the fitter fails just take dT/dW from the cluster
674  catch (...) {
675  mf::LogWarning("FeatureVertexFinder")
676  << "Fitter failed, using the clusters default dTdW()";
677  delete the2Dtrack;
678  dtdwstart.push_back(std::tan(RawClusters[iclu]->StartAngle()));
679  continue;
680  }
681 
682  delete the2Dtrack;
683  } //<---End if the cluster has 2 or more hits
684 
685  // If the cluster has fewer than 2 hits just take the dT/dW from
686  // the cluster
687 
688  else {
689  dtdwstart.push_back(std::tan(RawClusters[iclu]->StartAngle()));
690  }
691  } //<---End loop over clusters iclu
692 
693  // Now that I have slopes for all the clusters move on to finding
694  // intersection points
695 
696  // Looping over cryostats
697 
698  for (auto const& tpc : geom->Iterate<geo::TPCGeo>()) {
699  for (unsigned int i = 0; i < tpc.Nplanes(); ++i) {
700 
701  // If there is at least one cluster found
702 
703  if (Cls[i].size() >= 1) {
704 
705  // Loop over each cluster
706 
707  for (unsigned int j = 0; j < Cls[i].size(); ++j) {
708  // === Current Clusters Plane ===
709  Clu_Plane.push_back(RawClusters.at(Cls.at(i).at(j))->View());
710  // === Current Clusters StartPos ===
711  Clu_StartPos_Wire.push_back(RawClusters.at(Cls.at(i).at(j))->StartWire());
712  Clu_StartPos_TimeTick.push_back(RawClusters.at(Cls.at(i).at(j))->StartTick());
713  // === Current Clusters EndPos ===
714  Clu_EndPos_Wire.push_back(RawClusters.at(Cls.at(i).at(j))->EndWire());
715  Clu_EndPos_TimeTick.push_back(RawClusters.at(Cls.at(i).at(j))->EndTick());
716  // Current Clusters Slope (In Wire and Time Tick)
717  Clu_Slope.push_back(dtdwstart[Cls[i][j]]);
718  Clu_Length.push_back(std::sqrt(pow((RawClusters.at(Cls.at(i).at(j))->StartWire() -
719  RawClusters.at(Cls.at(i).at(j))->EndWire()) *
720  13.5,
721  2) +
722  pow(RawClusters.at(Cls.at(i).at(j))->StartTick() -
723  RawClusters.at(Cls.at(i).at(j))->EndTick(),
724  2)));
725 
726  // Given a slope and a point find the y-intercept
727  // c = y-mx
728 
729  Clu_Yintercept.push_back(
730  RawClusters.at(Cls.at(i).at(j))->StartTick() -
731  (dtdwstart[Cls[i][j]] * RawClusters.at(Cls.at(i).at(j))->StartWire()));
732 
733  // Also calculating the y-intercept but using the end
734  // time of the cluster correct for the possibility that
735  // the clustering didn't get start and end points right
736 
737  Clu_Yintercept2.push_back(
738  RawClusters.at(Cls.at(i).at(j))->EndTick() -
739  (dtdwstart[Cls[i][j]] * RawClusters.at(Cls.at(i).at(j))->EndWire()));
740 
741  // Iterating on the total number of clusters found
742 
743  nClustersFound++;
744  } //<---End loop over all clusters
745 
746  } //<---End check if we have at least one cluster
747 
748  // If no clusters were found then put in dummy vertex values
749  else {
750  TwoDvtx_wire.push_back(-1);
751  TwoDvtx_time.push_back(-1);
752  TwoDvtx_plane.push_back(-1);
753  }
754 
755  } //<---End loop over planes (i)
756  } //<---End loop over tpc's
757 
758  // Now loop over all the clusters found and establish a
759  // preliminary set of 2d-verticies based on the slope/intercept of
760  // those clusters
761 
762  for (unsigned int n = nClustersFound; n > 0; n--) {
763 
764  // Looping over the clusters starting from the first cluster and
765  // checking against the nCluster
766 
767  for (unsigned int m = 0; m < n; m++) {
768 
769  // Checking to make sure clusters are in the same view
770 
771  if (Clu_Plane[n] == Clu_Plane[m]) {
772  // --- Skip the vertex if the lines slope don't intercept ---
773  if (Clu_Slope[m] - Clu_Slope[n] == 0) { break; }
774 
775  // === X intersection = (yInt2 - yInt1) / (slope1 - slope2) ===
776  float intersection_X =
777  (Clu_Yintercept[n] - Clu_Yintercept[m]) / (Clu_Slope[m] - Clu_Slope[n]);
778 
779  // === Y intersection = (slope1 * XInt) + yInt1 ===
780  float intersection_Y = (Clu_Slope[m] * intersection_X) + Clu_Yintercept[m];
781 
782  // === X intersection = (yInt2 - yInt1) / (slope1 - slope2) ===
783  float intersection_X2 =
785 
786  // === Y intersection = (slope1 * XInt) + yInt1 ===
787  float intersection_Y2 = (Clu_Slope[m] * intersection_X2) + Clu_Yintercept2[m];
788 
789  // Skipping crap intersection points
790  geo::PlaneID const planeid(0, 0, Clu_Plane[m]);
791  if (intersection_X2 < 1) { intersection_X2 = -999; }
792  if (intersection_X2 > geom->Nwires(planeid)) { intersection_X2 = -999; }
793  if (intersection_Y2 < 0) { intersection_Y2 = -999; }
794  if (intersection_Y2 > detProp.NumberTimeSamples()) { intersection_Y2 = -999; }
795  if (intersection_X < 1) { intersection_X = -999; }
796  if (intersection_X > geom->Nwires(planeid)) { intersection_X = -999; }
797  if (intersection_Y < 0) { intersection_Y = -999; }
798  if (intersection_Y > detProp.NumberTimeSamples()) { intersection_Y = -999; }
799 
800  // Putting in a protection for the findManyHit function
801 
802  try {
803  // Gathering the hits associated with the current cluster
804  std::vector<art::Ptr<recob::Hit>> hitClu1 = fmhit.at(m);
805  std::vector<art::Ptr<recob::Hit>> hitClu2 = fmhit.at(n);
806 
807  // If the intersection point is 80 or more wires away from
808  // either cluster and one of the clusters has fewer than 8 hits the
809  // intersection is likely a crap one and we won't save this point
810  if ((abs(Clu_EndPos_Wire[m] - intersection_X2) > 80 && hitClu1.size() < 8) ||
811  (abs(Clu_EndPos_Wire[n] - intersection_X2) > 80 && hitClu2.size() < 8)) {
812  intersection_X2 = -999;
813  intersection_Y2 = -999;
814  }
815 
816  if ((abs(Clu_StartPos_Wire[m] - intersection_X) > 80 && hitClu1.size() < 8) ||
817  (abs(Clu_StartPos_Wire[n] - intersection_X) > 80 && hitClu2.size() < 8)) {
818  intersection_X = -999;
819  intersection_Y = -999;
820  }
821 
822  // If the intersection point is 50 or more wires away from
823  // either cluster and the one of the clusters has fewer than 3 hits
824  // the intersection is likely a crap one and we won't save this
825  // point
826  if ((abs(Clu_EndPos_Wire[m] - intersection_X2) > 50 && hitClu1.size() < 4) ||
827  (abs(Clu_EndPos_Wire[n] - intersection_X2) > 50 && hitClu2.size() < 4)) {
828  intersection_X2 = -999;
829  intersection_Y2 = -999;
830  }
831 
832  if ((abs(Clu_StartPos_Wire[m] - intersection_X) > 50 && hitClu1.size() < 4) ||
833  (abs(Clu_StartPos_Wire[n] - intersection_X) > 50 && hitClu2.size() < 4)) {
834  intersection_X = -999;
835  intersection_Y = -999;
836  }
837  }
838  catch (...) {
839  mf::LogWarning("FeatureVertexFinder") << "FindManyHit Function faild";
840  intersection_X = -999;
841  intersection_Y = -999;
842  intersection_X2 = -999;
843  intersection_Y2 = -999;
844  continue;
845  }
846 
847  // Push back a candidate 2dClusterVertex if it is inside the
848  // detector
849 
850  if (intersection_X2 > 1 && intersection_Y2 > 0 &&
851  (intersection_X2 < geom->Nwires(planeid)) &&
852  (intersection_Y2 < detProp.NumberTimeSamples())) {
853 
854  TwoDvtx_wire.push_back(intersection_X2);
855  TwoDvtx_time.push_back(intersection_Y2);
856  TwoDvtx_plane.push_back(Clu_Plane[m]);
857  } //<---End saving a "good 2d vertex" candidate
858 
859  // Push back a candidate 2dClusterVertex if it is inside the
860  // detector
861 
862  if (intersection_X > 1 && intersection_Y > 0 &&
863  (intersection_X < geom->Nwires(planeid)) &&
864  (intersection_Y < detProp.NumberTimeSamples())) {
865  TwoDvtx_wire.push_back(intersection_X);
866  TwoDvtx_time.push_back(intersection_Y);
867  TwoDvtx_plane.push_back(Clu_Plane[m]);
868  } //<---End saving a "good 2d vertex" candidate
869 
870  } //<---End check that they are in differnt planes
871  } //<---End m loop
872  } //<---End n loop
873 
874  } //<---End 2dClusterVertexCandidates
875 
876  // -----------------------------------------------------------------------------
877  // Get 3d Vertex Candidates from clusters 2d Vertex candidates
878  // -----------------------------------------------------------------------------
880  detinfo::DetectorPropertiesData const& detProp,
881  std::vector<double> const& Wire_2dvtx,
882  std::vector<double> const& Time_2dvtx,
883  std::vector<double> const& Plane_2dvtx)
884  {
885 
886  std::vector<double> vtx_wire_merged;
887  std::vector<double> vtx_time_merged;
888  std::vector<double> vtx_plane_merged;
889 
890  double y_coord = 0, z_coord = 0;
891 
892  bool merged = false;
893 
895 
896  // MERGING THE LONG LIST OF 2D CANDIDATES
897 
898  // Looping over 2d-verticies (loop1)
899 
900  for (size_t vtxloop1 = 0; vtxloop1 < Wire_2dvtx.size(); vtxloop1++) {
901  if (Wire_2dvtx[vtxloop1] < 0) { continue; }
902 
903  merged = false;
904 
905  // Looping over 2d-verticies (loop2)
906 
907  for (size_t vtxloop2 = vtxloop1 + 1; vtxloop2 < Wire_2dvtx.size(); vtxloop2++) {
908  if (Wire_2dvtx[vtxloop2] < 0) { continue; }
909 
910  // Make sure the 2d-Verticies are in the same plane
911 
912  if (Plane_2dvtx[vtxloop1] == Plane_2dvtx[vtxloop2]) {
913 
914  // Considering merging 2d vertices if they
915  // are within 3 wires of each other
916 
917  if (fabs(Wire_2dvtx[vtxloop1] - Wire_2dvtx[vtxloop2]) < 4) {
918 
919  // Merge the verticies if they are within 10 time ticks
920 
921  if (fabs(Time_2dvtx[vtxloop1] - Time_2dvtx[vtxloop2]) < 10) {
922  vtx_wire_merged.push_back(((Wire_2dvtx[vtxloop2] + Wire_2dvtx[vtxloop1]) / 2));
923  vtx_time_merged.push_back(((Time_2dvtx[vtxloop2] + Time_2dvtx[vtxloop1]) / 2));
924  vtx_plane_merged.push_back(Plane_2dvtx[vtxloop1]);
925 
926  merged = true;
927  if (vtxloop2 < Wire_2dvtx.size()) { vtxloop2++; }
928  if (vtxloop1 < Wire_2dvtx.size()) { vtxloop1++; }
929  } //<---End the check within 10 time ticks for merging
930  } //<---Looking at vertices that are within 3 wires of each other
931  } //<---Only looking at vertices that are in the same plane
932  } //<---End vtxloop2
933  if (!merged) {
934  vtx_wire_merged.push_back(Wire_2dvtx[vtxloop1]);
935  vtx_time_merged.push_back(Time_2dvtx[vtxloop1]);
936  vtx_plane_merged.push_back(Plane_2dvtx[vtxloop1]);
937  } //<---end saving unmerged verticies
938  } //<---End vtxloop1
939 
940  // Variables for channel numbers
941 
942  uint32_t vtx1_channel = 0;
943  uint32_t vtx2_channel = 0;
944 
945  // Having now found a very long list of potential 2-d end points
946  // we need to check if any of them match between planes and only
947  // keep those that have matches
948 
949  for (auto const& tpcid : geom->Iterate<geo::TPCID>()) {
950  for (unsigned int vtx = vtx_wire_merged.size(); vtx > 0; --vtx) {
951  for (unsigned int vtx1 = 0; vtx1 < vtx; ++vtx1) {
952 
953  // Make sure we are comparing verticies from different planes
954  if (vtx_plane_merged[vtx1] == vtx_plane_merged[vtx]) continue;
955 
956  // To figure out if these two verticies are from a common
957  // point we need to check if the channels intersect and if
958  // they are close in time ticks as well...to do this we have
959  // to do some converting to use geom->PlaneWireToChannel(...)
960  geo::PlaneID const vtx1_planeid(tpcid, vtx_plane_merged[vtx1]);
961  geo::WireID const vtx1_wireid(vtx1_planeid, vtx_wire_merged[vtx1]);
962  try {
963  vtx1_channel = geom->PlaneWireToChannel(vtx1_wireid);
964  }
965  catch (...) {
966  mf::LogWarning("FeatureVertexFinder") << "PlaneWireToChannel Failed";
967  continue;
968  }
969 
970  geo::PlaneID const vtx2_planeid(tpcid, vtx_plane_merged[vtx]);
971  geo::WireID const vtx2_wireid(vtx2_planeid, vtx_wire_merged[vtx]);
972  try {
973  vtx2_channel = geom->PlaneWireToChannel(vtx2_wireid);
974  }
975  catch (...) {
976  mf::LogWarning("FeatureVertexFinder") << "PlaneWireToChannel Failed";
977  continue;
978  }
979 
980  // Check to see if the channels intersect and save the y and z coordinate
981  if (!geom->ChannelsIntersect(vtx1_channel, vtx2_channel, y_coord, z_coord)) {
982  mf::LogWarning("FeatureVertexFinder") << "match failed for some reason";
983  continue;
984  }
985 
986  // If the channels intersect establish if they are close in "X"
987  float tempXCluster1 = detProp.ConvertTicksToX(vtx_time_merged[vtx1], vtx1_planeid);
988  float tempXCluster2 = detProp.ConvertTicksToX(vtx_time_merged[vtx], vtx2_planeid);
989 
990  // Now check if the matched channels are within 0.5 cm when
991  // projected in X and that we have less than 100 of these
992  // candidates...because more than that seems silly
993 
994  if (std::abs(tempXCluster1 - tempXCluster2) < 0.5 && candidate_x.size() < 101) {
995  candidate_x.push_back(detProp.ConvertTicksToX(vtx_time_merged[vtx1], vtx1_planeid));
996  candidate_y.push_back(y_coord);
997  candidate_z.push_back(z_coord);
998  candidate_strength.push_back(
999  10); //<--For cluster verticies I give it a strength of "10"
1000  // arbitrarily for now
1001 
1002  } //<---End Checking if the vertices agree "well enough" in time tick
1003  } //<---end vtx1 for loop
1004  } //<---End vtx for loop
1005  } //<---End loop over TPC's
1006 
1007  } //<---End Find3dVtxFrom2dClusterVtxCand
1008 
1009  // -----------------------------------------------------------------------------
1010  // Get 3d Vertex Candidates from clusters 2d Vertex candidates
1011  // -----------------------------------------------------------------------------
1013  std::vector<double> const& merge_vtxX,
1014  std::vector<double> const& merge_vtxY,
1015  std::vector<double> const& merge_vtxZ,
1016  std::vector<double> const& merge_vtxStgth)
1017  {
1018 
1019  std::vector<double> x_3dVertex_dupRemoved = {0.};
1020  std::vector<double> y_3dVertex_dupRemoved = {0.};
1021  std::vector<double> z_3dVertex_dupRemoved = {0.};
1022  std::vector<double> strength_dupRemoved = {0.};
1023 
1024  // Looping over the 3d candidates found thus far
1025 
1026  for (size_t dup = 0; dup < merge_vtxX.size(); dup++) {
1027  // Temperary storing the current vertex
1028  float tempX_dup = merge_vtxX[dup];
1029  float tempY_dup = merge_vtxY[dup];
1030  float tempZ_dup = merge_vtxZ[dup];
1031  float tempStgth = merge_vtxStgth[dup];
1032 
1033  // Setting a boolian to see if this is a duplicate
1034 
1035  bool duplicate_found = false;
1036 
1037  // Looping over the rest of the list for duplicate check
1038 
1039  for (size_t check = dup + 1; check < merge_vtxX.size(); check++) {
1040 
1041  // I am going to call a duplicate vertex one that matches in x, y,
1042  // and z within 0.1 cm for all 3 coordinates simultaneously
1043 
1044  if (std::abs(merge_vtxX[check] - tempX_dup) < 0.1 &&
1045  std::abs(merge_vtxY[check] - tempY_dup) < 0.1 &&
1046  std::abs(merge_vtxZ[check] - tempZ_dup) < 0.1) {
1047  duplicate_found = true;
1048  } //<---End checking to see if this is a duplicate vertex
1049 
1050  } //<---End check for loop
1051 
1052  // If we didn't find a duplicate then lets save this 3d vertex
1053  // as a real candidate for consideration
1054 
1055  if (!duplicate_found && tempX_dup > 0) {
1056  x_3dVertex_dupRemoved.push_back(tempX_dup);
1057  y_3dVertex_dupRemoved.push_back(tempY_dup);
1058  z_3dVertex_dupRemoved.push_back(tempZ_dup);
1059  strength_dupRemoved.push_back(tempStgth);
1060  } //<---End storing only non-duplicates
1061 
1062  } //<---End dup for loop
1063 
1064  // Sorting the verticies I have found such that the first in the
1065  // list is the vertex with the highest vertex strength and the
1066  // lowest z location
1067 
1068  int flag = 1;
1069  double tempX, tempY, tempZ, tempS;
1070 
1071  // Looping over all duplicate removed candidates
1072 
1073  for (size_t npri = 0; (npri < x_3dVertex_dupRemoved.size()) && flag; npri++) {
1074  flag = 0;
1075  for (size_t mpri = 0; mpri < x_3dVertex_dupRemoved.size() - 1; mpri++) {
1076  // Swap the order of the two elements
1077  if (strength_dupRemoved[mpri + 1] > strength_dupRemoved[mpri] ||
1078  (strength_dupRemoved[mpri + 1] == strength_dupRemoved[mpri] &&
1079  z_3dVertex_dupRemoved[mpri] > z_3dVertex_dupRemoved[mpri + 1])) {
1080  tempX = x_3dVertex_dupRemoved[mpri];
1081  x_3dVertex_dupRemoved[mpri] = x_3dVertex_dupRemoved[mpri + 1];
1082  x_3dVertex_dupRemoved[mpri + 1] = tempX;
1083 
1084  tempY = y_3dVertex_dupRemoved[mpri];
1085  y_3dVertex_dupRemoved[mpri] = y_3dVertex_dupRemoved[mpri + 1];
1086  y_3dVertex_dupRemoved[mpri + 1] = tempY;
1087 
1088  tempZ = z_3dVertex_dupRemoved[mpri];
1089  z_3dVertex_dupRemoved[mpri] = z_3dVertex_dupRemoved[mpri + 1];
1090  z_3dVertex_dupRemoved[mpri + 1] = tempZ;
1091 
1092  tempS = strength_dupRemoved[mpri];
1093  strength_dupRemoved[mpri] = strength_dupRemoved[mpri + 1];
1094  strength_dupRemoved[mpri + 1] = tempS;
1095 
1096  flag = 1;
1097 
1098  } //<---Inside swap loop
1099  } //<---End mpri
1100  } //<---End npri loop
1101 
1102  // Pushing into a vector of merged and sorted verticies
1103 
1104  for (size_t count = 0; count < x_3dVertex_dupRemoved.size(); count++) {
1105  MergeSort3dVtx_xpos.push_back(x_3dVertex_dupRemoved[count]);
1106  MergeSort3dVtx_ypos.push_back(y_3dVertex_dupRemoved[count]);
1107  MergeSort3dVtx_zpos.push_back(z_3dVertex_dupRemoved[count]);
1108  MergeSort3dVtx_strength.push_back(strength_dupRemoved[count]);
1109  }
1110 
1111  } // End MergeAndSort3dVtxCandidate
1112 
1114 
1115  // -----------------------------------------------------------------------------
1116 
1117 } //<---End namespace vertex
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Encapsulate the construction of a single cyostat.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Double_t z
Definition: plot.C:276
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
Geometry information for a single TPC.
Definition: TPCGeo.h:36
void produce(art::Event &evt) override
constexpr auto abs(T v)
Returns the absolute value of the argument.
bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double &y, double &z) const
Returns an intersection point of two channels.
Double_t zz
Definition: plot.C:276
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
void MergeAndSort3dVtxCandidate(std::vector< double > const &merge_vtxX, std::vector< double > const &merge_vtxY, std::vector< double > const &merge_vtxZ, std::vector< double > const &merge_vtxStgth)
std::vector< double > MergeSort3dVtx_strength
Definition of vertex object for LArSoft.
Definition: Vertex.h:35
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::vector< std::vector< int > > Cls
void Find3dVtxFrom2dClusterVtxCand(detinfo::DetectorPropertiesData const &detProp, std::vector< double > const &Wire_2dvtx, std::vector< double > const &Time_2dvtx, std::vector< double > const &Plane_2dvtx)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
void Find2dClusterVertexCandidates(detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Cluster > RawClusters, art::FindManyP< recob::Hit > fmhit)
Provides recob::Track data product.
reference at(size_type n)
Definition: PtrVector.h:359
FeatureVertexFinder(fhicl::ParameterSet const &pset)
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Declaration of cluster object.
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
double ConvertTicksToX(double ticks, int p, int t, int c) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
WireID NearestWireID(Point_t const &point, PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
raw::ChannelID_t NearestChannel(Point_t const &worldLoc, PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
constexpr WireID()=default
Default constructor: an invalid TPC ID.
Char_t n[5]
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Definition: fwd.h:26
void Get3dVertexCandidates(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::EndPoint2D >> const &EndPoints, bool PlaneDet)
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane.
vertex reconstruction