LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CCTrackMaker_module.cc
Go to the documentation of this file.
1 //
3 // CCTrackMaker
4 //
5 // Make tracks using ClusterCrawler clusters and vertex info
6 //
7 // baller@fnal.gov, October 2014
8 // May 2015: Major re-write
9 //
11 
12 // C++ includes
13 #include <algorithm>
14 #include <cmath>
15 #include <iomanip>
16 #include <iostream>
17 #include <string>
18 
19 // Framework includes
26 #include "fhiclcpp/ParameterSet.h"
28 
29 // LArSoft includes
42 
43 struct CluLen {
44  int index;
45  int length;
46 };
47 
48 bool greaterThan(CluLen c1, CluLen c2)
49 {
50  return (c1.length > c2.length);
51 }
52 bool lessThan(CluLen c1, CluLen c2)
53 {
54  return (c1.length < c2.length);
55 }
56 
57 namespace trkf {
58 
59  class CCTrackMaker : public art::EDProducer {
60  public:
61  explicit CCTrackMaker(fhicl::ParameterSet const& pset);
62 
63  private:
64  void produce(art::Event& evt) override;
65 
66  std::string fHitModuleLabel;
67  std::string fClusterModuleLabel;
68  std::string fVertexModuleLabel;
69 
70  // services
72 
75 
76  // Track matching parameters
77  unsigned short algIndex;
78  std::vector<short> fMatchAlgs;
79  std::vector<float> fXMatchErr;
80  std::vector<float> fAngleMatchErr;
81  std::vector<float> fChgAsymFactor;
82  std::vector<float> fMatchMinLen;
83  std::vector<bool> fMakeAlgTracks;
84 
85  // Cluster merging parameters
86  float fMaxDAng;
87  float fChainMaxdX;
88  float fChainVtxAng;
92 
93  float fChgWindow;
94  float fWirePitch;
95  // cosmic ray tagging
96  float fFiducialCut;
97  float fDeltaRayCut;
98 
99  bool fMakePFPs;
100 
101  // vertex fitting
102  unsigned short fNVtxTrkHitsFit;
104 
105  // temp
106  bool fuBCode;
107 
108  // debugging inputs
109  short fDebugAlg;
110  short fDebugPlane;
113  bool prt;
114 
115  // hit indexing info
116  std::array<unsigned int, 3> firstWire;
117  std::array<unsigned int, 3> lastWire;
118  std::array<unsigned int, 3> firstHit;
119  std::array<unsigned int, 3> lastHit;
120  std::array<std::vector<std::pair<int, int>>, 3> WireHitRange;
121 
122  std::vector<art::Ptr<recob::Hit>> allhits;
123 
124  // Cluster parameters
125  struct clPar {
126  std::array<float, 2> Wire; // Begin/End Wire
127  std::array<float, 2> X; // Begin/End X
128  std::array<short, 2> Time; // Begin/End Time
129  std::array<float, 2> Slope; // Begin/End slope
130  std::array<float, 2> Charge; // Begin/End charge
131  std::array<float, 2> ChgNear; // Charge near the cluster at each end
132  std::array<float, 2> Angle; // Begin/End angle (radians)
133  std::array<short, 2> Dir; // Product of end * slope
134  std::array<short, 2> VtxIndex; // Vertex index
135  std::array<short, 2> mVtxIndex; // "Maybe" Vertex index
136  std::array<short, 2> BrkIndex; // Broken cluster index
137  std::array<float, 2> MergeError; // Broken cluster merge error (See MakeClusterChains)
138  unsigned short EvtIndex; // index of the cluster in clusterlist
139  short InTrack; // cluster -> chain index
140  unsigned short Length; // cluster length (wires)
141  float TotChg; // total charge of cluster (or series of clusters)
142  };
143  // vector of cluster parameters in each plane
144  std::array<std::vector<clPar>, 3> cls;
145 
146  // cluster chain parameters
147  struct ClsChainPar {
148  std::array<float, 2> Wire; // Begin/End Wire
149  std::array<float, 2> X; // Begin/End X
150  std::array<float, 2> Time; // Begin/End Time
151  std::array<float, 2> Slope; // Begin/End slope
152  std::array<float, 2> Angle; // Begin/End angle (radians)
153  std::array<short, 2> VtxIndex; // Vertex index
154  std::array<short, 2> Dir;
155  std::array<float, 2> ChgNear; // Charge near the cluster at each end
156  std::array<short, 2> mBrkIndex; // a "Maybe" cluster index
157  unsigned short Length;
158  float TotChg;
159  std::vector<unsigned short> ClsIndex;
160  std::vector<unsigned short> Order;
161  short InTrack; // cluster -> track ID (-1 if none, 0 if under construction)
162  };
163  // vector of cluster parameters in each plane
164  std::array<std::vector<ClsChainPar>, 3> clsChain;
165 
166  // 3D Vertex info
167  struct vtxPar {
168  short ID;
169  unsigned short EvtIndex;
170  float X;
171  float Y;
172  float Z;
173  std::array<unsigned short, 3> nClusInPln;
174  bool Neutrino;
175  };
176 
177  std::vector<vtxPar> vtx;
178 
179  // cluster indices assigned to one vertex. Filled in VtxMatch
180  std::array<std::vector<unsigned short>, 3> vxCls;
181 
182  struct TrkPar {
183  short ID;
184  unsigned short Proc; // 1 = VtxMatch, 2 = ...
185  std::array<std::vector<art::Ptr<recob::Hit>>, 3> TrkHits;
186  std::vector<TVector3> TrjPos;
187  std::vector<TVector3> TrjDir;
188  std::array<short, 2> VtxIndex;
189  std::vector<unsigned short> ClsEvtIndices;
190  float Length;
191  short ChgOrder;
192  short MomID;
193  std::array<bool, 2> EndInTPC;
194  std::array<bool, 2> GoodEnd; // set true if there are hits in all planes at the end point
195  std::vector<short> DtrID;
196  short PDGCode;
197  };
198  std::vector<TrkPar> trk;
199 
200  // Array of pointers to hits in each plane for one track
201  std::array<std::vector<art::Ptr<recob::Hit>>, 3> trkHits;
202  // and for one seed
203  std::array<std::vector<art::Ptr<recob::Hit>>, 3> seedHits;
204  // relative charge normalization between planes
205  std::array<float, 3> ChgNorm;
206 
207  // Vector of PFParticle -> track IDs. The first element will be
208  // track ID = 0 indicating that this is a neutrino PFParticle and
209  // there is no associated track
210  std::vector<unsigned short> pfpToTrkID;
211 
212  // characterize the match between clusters in 2 or 3 planes
213  struct MatchPars {
214  std::array<short, 3> Cls;
215  std::array<unsigned short, 3> End;
216  std::array<float, 3> Chg;
217  short Vtx;
218  float dWir; // wire difference at the matching end
219  float dAng; // angle difference at the matching end
220  float dX; // X difference
221  float Err; // Wire,Angle,Time match error
222  short oVtx;
223  float odWir; // wire difference at the other end
224  float odAng; // angle difference at the other end
225  float odX; // time difference at the other end
226  float dSP; // space point difference
227  float oErr; // dAngle dX match error
228  };
229  // vector of many match combinations
230  std::vector<MatchPars> matcomb;
231 
233  geo::TPCID const& tpcid) const;
234 
235  void PrintTracks() const;
236 
237  void MakeClusterChains(detinfo::DetectorPropertiesData const& detProp,
238  art::FindManyP<recob::Hit> const& fmCluHits,
239  geo::TPCID const& tpcid);
240  float dXClTraj(art::FindManyP<recob::Hit> const& fmCluHits,
241  unsigned short ipl,
242  unsigned short icl1,
243  unsigned short end1);
244  void FillChgNear(detinfo::DetectorPropertiesData const& detProp, geo::TPCID const& tpcid);
245  void FillWireHitRange(geo::TPCID const& tpcid);
246 
247  // Find clusters that point to vertices but do not have a
248  // cluster-vertex association made by ClusterCrawler
249  void FindMaybeVertices(geo::TPCID const& tpcid);
250  // match clusters associated with vertices
251  void VtxMatch(detinfo::DetectorPropertiesData const& detProp,
252  art::FindManyP<recob::Hit> const& fmCluHits,
253  geo::TPCID const& tpcid);
254  // match clusters in all planes
255  void PlnMatch(detinfo::DetectorPropertiesData const& detProp, geo::TPCID const& tpcid);
256  // match clusters in all planes
257  void AngMatch(art::FindManyP<recob::Hit> const& fmCluHits);
258 
259  // Make the track/vertex and mother/daughter relationships
260  void MakeFamily(geo::TPCID const& tpcid);
261  void TagCosmics(geo::TPCID const& tpcid);
262 
263  void FitVertices(detinfo::DetectorPropertiesData const& detProp, geo::TPCID const& tpcid);
264 
265  // fill the end matching parameters in the MatchPars struct
266  void FillEndMatch(detinfo::DetectorPropertiesData const& detProp,
267  geo::TPCID const& tpcid,
268  MatchPars& match);
269  // 2D version
270  void FillEndMatch2(MatchPars& match);
271 
272  float ChargeAsym(std::array<float, 3>& mChg);
273 
274  bool FindMissingCluster(unsigned short kpl,
275  short& kcl,
276  unsigned short& kend,
277  float kWir,
278  float kX,
279  float okWir,
280  float okX);
281 
282  bool DupMatch(MatchPars& match, unsigned short nplanes);
283 
284  void SortMatches(detinfo::DetectorPropertiesData const& detProp,
285  art::FindManyP<recob::Hit> const& fmCluHits,
286  unsigned short procCode,
287  geo::TPCID const& tpcid);
288 
289  // fill the trkHits array using information.
290  void FillTrkHits(art::FindManyP<recob::Hit> const& fmCluHits,
291  unsigned short imat,
292  unsigned short nplanes);
293 
294  // store the track in the trk vector
295  void StoreTrack(detinfo::DetectorPropertiesData const& detProp,
296  art::FindManyP<recob::Hit> const& fmCluHits,
297  unsigned short imat,
298  unsigned short procCode,
299  geo::TPCID const& tpcid);
300 
301  // returns the charge along the line between (wire1, time1) and (wire2, time2)
302  float ChargeNear(geo::PlaneID const& planeid,
303  unsigned short wire1,
304  float time1,
305  unsigned short wire2,
306  float time2);
307 
308  // inflate cuts at large angle
309  float AngleFactor(float slope);
310 
311  }; // class CCTrackMaker
312 
313  //-------------------------------------------------
314  CCTrackMaker::CCTrackMaker(fhicl::ParameterSet const& pset) : EDProducer{pset}
315  {
316  fHitModuleLabel = pset.get<std::string>("HitModuleLabel");
317  fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
318  fVertexModuleLabel = pset.get<std::string>("VertexModuleLabel");
319  // track matching
320  fMatchAlgs = pset.get<std::vector<short>>("MatchAlgs");
321  fXMatchErr = pset.get<std::vector<float>>("XMatchErr");
322  fAngleMatchErr = pset.get<std::vector<float>>("AngleMatchErr");
323  fChgAsymFactor = pset.get<std::vector<float>>("ChgAsymFactor");
324  fMatchMinLen = pset.get<std::vector<float>>("MatchMinLen");
325  fMakeAlgTracks = pset.get<std::vector<bool>>("MakeAlgTracks");
326  // Cluster merging
327  fMaxDAng = pset.get<float>("MaxDAng");
328  fChainMaxdX = pset.get<float>("ChainMaxdX");
329  fChainVtxAng = pset.get<float>("ChainVtxAng");
330  fMergeChgAsym = pset.get<float>("MergeChgAsym");
331  // Cosmic ray tagging
332  fFiducialCut = pset.get<float>("FiducialCut");
333  fDeltaRayCut = pset.get<float>("DeltaRayCut");
334  // make PFParticles
335  fMakePFPs = pset.get<bool>("MakePFPs");
336  // vertex fitting
337  fNVtxTrkHitsFit = pset.get<unsigned short>("NVtxTrkHitsFit");
338  fHitFitErrFac = pset.get<float>("HitFitErrFac");
339  // uB code
340  fuBCode = pset.get<bool>("uBCode");
341  // debugging inputs
342  fDebugAlg = pset.get<short>("DebugAlg");
343  fDebugPlane = pset.get<short>("DebugPlane");
344  fDebugCluster = pset.get<short>("DebugCluster");
345  fPrintAllClusters = pset.get<bool>("PrintAllClusters");
346 
347  // Consistency check
348  if (fMatchAlgs.size() > fXMatchErr.size() || fMatchAlgs.size() > fAngleMatchErr.size() ||
349  fMatchAlgs.size() > fChgAsymFactor.size() || fMatchAlgs.size() > fMatchMinLen.size() ||
350  fMatchAlgs.size() > fMakeAlgTracks.size()) {
351  mf::LogError("CCTM") << "Incompatible fcl input vector sizes";
352  return;
353  }
354  // Reality check
355  for (unsigned short ii = 0; ii < fMatchAlgs.size(); ++ii) {
356  if (fAngleMatchErr[ii] <= 0 || fXMatchErr[ii] <= 0) {
357  mf::LogError("CCTM") << "Invalid matching parameters " << fAngleMatchErr[ii] << " "
358  << fXMatchErr[ii];
359  return;
360  }
361  } // ii
362 
363  produces<std::vector<recob::PFParticle>>();
364  produces<art::Assns<recob::PFParticle, recob::Track>>();
365  produces<art::Assns<recob::PFParticle, recob::Cluster>>();
366  produces<art::Assns<recob::PFParticle, recob::Seed>>();
367  produces<art::Assns<recob::PFParticle, recob::Vertex>>();
368  produces<std::vector<recob::Vertex>>();
369  produces<std::vector<recob::Track>>();
370  produces<art::Assns<recob::Track, recob::Hit>>();
371  produces<std::vector<recob::Seed>>();
372  }
373 
374  //------------------------------------------------------------------------------------//
376  {
378 
379  fChgWindow = 40; // window (ticks) for identifying shower-like clusters
380 
381  auto tcol = std::make_unique<std::vector<recob::Track>>();
382  auto thassn = std::make_unique<art::Assns<recob::Track, recob::Hit>>();
383  auto vcol = std::make_unique<std::vector<recob::Vertex>>();
384  auto pcol = std::make_unique<std::vector<recob::PFParticle>>();
385  auto ptassn = std::make_unique<art::Assns<recob::PFParticle, recob::Track>>();
386  auto pcassn = std::make_unique<art::Assns<recob::PFParticle, recob::Cluster>>();
387  auto psassn = std::make_unique<art::Assns<recob::PFParticle, recob::Seed>>();
388  auto pvassn = std::make_unique<art::Assns<recob::PFParticle, recob::Vertex>>();
389 
390  // seed collection
391  auto scol = std::make_unique<std::vector<recob::Seed>>();
392  auto shassn = std::make_unique<art::Assns<recob::Seed, recob::Hit>>();
393 
394  // all hits
395  art::Handle<std::vector<recob::Hit>> allhitsListHandle;
396  // cluster list
397  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
398  std::vector<art::Ptr<recob::Cluster>> clusterlist;
399  // ClusterCrawler Vertices
401  std::vector<art::Ptr<recob::Vertex>> vtxlist;
402 
403  // get Hits
404  allhits.clear();
405  if (evt.getByLabel(fHitModuleLabel, allhitsListHandle))
406  art::fill_ptr_vector(allhits, allhitsListHandle);
407 
408  // get Clusters
409  if (evt.getByLabel(fClusterModuleLabel, clusterListHandle))
410  art::fill_ptr_vector(clusterlist, clusterListHandle);
411  if (clusterlist.size() == 0) return;
412  // get cluster - hit associations
413  art::FindManyP<recob::Hit> fmCluHits(clusterListHandle, evt, fClusterModuleLabel);
414 
415  // get Vertices
416  if (evt.getByLabel(fVertexModuleLabel, VtxListHandle))
417  art::fill_ptr_vector(vtxlist, VtxListHandle);
419 
420  std::vector<CluLen> clulens;
421 
422  unsigned short ipl, icl, end, itr, tID, tIndex;
423 
424  // maximum error (see MakeClusterChains) for considering clusters broken
425  fMaxMergeError = 30;
426  // Cut on the error for a definitely broken cluster. Clusters with fMergeErrorCut < MergeError < fMaxMergeError
427  // are possibly broken clusters but we will consider these when matching between planes
428  fMergeErrorCut = 10;
429 
430  std::vector<art::Ptr<recob::Hit>> tmpHits;
431  std::vector<art::Ptr<recob::Cluster>> tmpCls;
432  std::vector<art::Ptr<recob::Vertex>> tmpVtx;
433 
434  // vector for PFParticle constructor
435  std::vector<size_t> dtrIndices;
436 
437  // some vectors for recob::Seed
438  double sPos[3], sDir[3];
439  double sErr[3] = {0, 0, 0};
440 
441  auto const detProp =
443 
444  // check consistency between clusters and associated hits
445  std::vector<art::Ptr<recob::Hit>> clusterhits;
446  for (icl = 0; icl < clusterlist.size(); ++icl) {
447  ipl = clusterlist[icl]->Plane().Plane;
448  clusterhits = fmCluHits.at(icl);
449  if (clusterhits[0]->WireID().Wire != std::nearbyint(clusterlist[icl]->EndWire())) {
450  std::cout << "CCTM Cluster-Hit End wire mis-match " << clusterhits[0]->WireID().Wire
451  << " vs " << std::nearbyint(clusterlist[icl]->EndWire()) << " Bail out! \n";
452  return;
453  }
454  for (unsigned short iht = 0; iht < clusterhits.size(); ++iht) {
455  if (clusterhits[iht]->WireID().Plane != ipl) {
456  std::cout << "CCTM Cluster-Hit plane mis-match " << ipl << " vs "
457  << clusterhits[iht]->WireID().Plane << " on hit " << iht << " Bail out! \n";
458  return;
459  } // hit-cluster plane mis-match
460  } // iht
461  } // icl
462  // end check consistency
463 
464  vtx.clear();
465  trk.clear();
466  for (auto const& tpcgeom : geom->Iterate<geo::TPCGeo>()) {
467  unsigned int const nplanes = tpcgeom.Nplanes();
468  if (nplanes > 3) continue;
469  for (ipl = 0; ipl < 3; ++ipl) {
470  cls[ipl].clear();
471  clsChain[ipl].clear();
472  trkHits[ipl].clear();
473  } // ipl
474  // FillWireHitRange also calculates the charge in each plane
475  FillWireHitRange(tpcgeom.ID());
476  for (auto const& planeid : geom->Iterate<geo::PlaneID>(tpcgeom.ID())) {
477  clulens.clear();
478  // sort clusters by increasing End wire number
479  for (icl = 0; icl < clusterlist.size(); ++icl) {
480  if (clusterlist[icl]->Plane() != planeid) continue;
481  CluLen clulen;
482  clulen.index = icl;
483  clulen.length = clusterlist[icl]->EndWire();
484  clulens.push_back(clulen);
485  }
486  if (clulens.empty()) continue;
487  // sort clusters
488  std::sort(clulens.begin(), clulens.end(), lessThan);
489  if (clulens.empty()) continue;
490  for (unsigned short ii = 0; ii < clulens.size(); ++ii) {
491  const unsigned short icl = clulens[ii].index;
492  clPar clstr;
493  clstr.EvtIndex = icl;
494  recob::Cluster const& cluster = *(clusterlist[icl]);
495  // Begin info -> end index 1 (DS)
496  clstr.Wire[1] = cluster.StartWire();
497  clstr.Time[1] = cluster.StartTick();
498  clstr.X[1] = (float)detProp.ConvertTicksToX(cluster.StartTick(), planeid);
499  clstr.Angle[1] = cluster.StartAngle();
500  clstr.Slope[1] = std::tan(cluster.StartAngle());
501  clstr.Dir[1] = 0;
502  clstr.Charge[1] = ChgNorm[ipl] * cluster.StartCharge();
503  // this will be filled later
504  clstr.ChgNear[1] = 0;
505  clstr.VtxIndex[1] = -1;
506  clstr.mVtxIndex[1] = -1;
507  clstr.BrkIndex[1] = -1;
508  clstr.MergeError[1] = fMaxMergeError;
509  // End info -> end index 0 (US)
510  clstr.Wire[0] = cluster.EndWire();
511  clstr.Time[0] = cluster.EndTick();
512  clstr.X[0] = (float)detProp.ConvertTicksToX(cluster.EndTick(), planeid);
513  clstr.Angle[0] = cluster.EndAngle();
514  clstr.Slope[0] = std::tan(cluster.EndAngle());
515  clstr.Dir[0] = 0;
516  if (clstr.Time[1] > clstr.Time[0]) {
517  clstr.Dir[0] = 1;
518  clstr.Dir[1] = -1;
519  }
520  else {
521  clstr.Dir[0] = -1;
522  clstr.Dir[1] = 1;
523  }
524  clstr.Charge[0] = ChgNorm[ipl] * cluster.EndCharge();
525  // this will be filled later
526  clstr.ChgNear[1] = 0;
527  clstr.VtxIndex[0] = -1;
528  clstr.mVtxIndex[0] = -1;
529  clstr.BrkIndex[0] = -1;
530  clstr.MergeError[0] = fMaxMergeError;
531  // other info
532  clstr.InTrack = -1;
533  clstr.Length = (unsigned short)(0.5 + clstr.Wire[1] - clstr.Wire[0]);
534  clstr.TotChg = ChgNorm[ipl] * cluster.Integral();
535  if (clstr.TotChg <= 0) clstr.TotChg = 1;
536  clusterhits = fmCluHits.at(icl);
537  if (clusterhits.size() == 0) {
538  mf::LogError("CCTM") << "No associated cluster hits " << icl;
539  continue;
540  }
541  // correct charge for missing cluster hits
542  clstr.TotChg *= clstr.Length / (float)clusterhits.size();
543  cls[ipl].push_back(clstr);
544  } // ii (icl)
545  } // planeid
546 
547  FillChgNear(detProp, tpcgeom.ID());
548 
549  // and finally the vertices
550  double xyz[3];
551  for (unsigned short ivx = 0; ivx < vtxlist.size(); ++ivx) {
552  vtxPar aVtx;
553  aVtx.ID = ivx + 1;
554  aVtx.EvtIndex = ivx;
555  vtxlist[ivx]->XYZ(xyz);
556  aVtx.X = xyz[0];
557  aVtx.Y = xyz[1];
558  aVtx.Z = xyz[2];
559  aVtx.nClusInPln[0] = 0;
560  aVtx.nClusInPln[1] = 0;
561  aVtx.nClusInPln[2] = 0;
562  std::vector<art::Ptr<recob::Cluster>> const& vtxCls = fmVtxCls.at(ivx);
563  std::vector<const unsigned short*> const& vtxClsEnd = fmVtxCls.data(ivx);
564  for (unsigned short vcass = 0; vcass < vtxCls.size(); ++vcass) {
565  icl = vtxCls[vcass].key();
566  // the cluster plane
567  ipl = vtxCls[vcass]->Plane().Plane;
568  end = *vtxClsEnd[vcass];
569  if (end > 1)
570  throw cet::exception("CCTM")
571  << "Invalid end data from vertex - cluster association" << end;
572  bool gotit = false;
573  // CCTM end is opposite of CC end
574  end = 1 - end;
575  for (unsigned short jcl = 0; jcl < cls[ipl].size(); ++jcl) {
576  if (cls[ipl][jcl].EvtIndex == icl) {
577  cls[ipl][jcl].VtxIndex[end] = ivx;
578  ++aVtx.nClusInPln[ipl];
579  gotit = true;
580  break;
581  } // index check
582  } // jcl
583  if (!gotit)
584  throw cet::exception("CCTM") << "Bad index from vertex - cluster association" << icl;
585  } // icl
586  vtx.push_back(aVtx);
587  } // ivx
588  // Find broken clusters
589  MakeClusterChains(detProp, fmCluHits, tpcgeom.ID());
590  FindMaybeVertices(tpcgeom.ID());
591 
592  // call algorithms in the specified order
593  matcomb.clear();
594  for (algIndex = 0; algIndex < fMatchAlgs.size(); ++algIndex) {
595  if (fMatchAlgs[algIndex] == 1) {
596  prt = (fDebugAlg == 1);
597  VtxMatch(detProp, fmCluHits, tpcgeom.ID());
598  if (fMakeAlgTracks[algIndex]) SortMatches(detProp, fmCluHits, 1, tpcgeom.ID());
599  }
600  else if (fMatchAlgs[algIndex] == 2) {
601  prt = (fDebugAlg == 2);
602  PlnMatch(detProp, tpcgeom.ID());
603  if (fMakeAlgTracks[algIndex]) SortMatches(detProp, fmCluHits, 2, tpcgeom.ID());
604  }
605  if (prt) PrintClusters(detProp, tpcgeom.ID());
606  } // algIndex
607  prt = false;
608  pfpToTrkID.clear();
609  // Determine the vertex/track hierarchy
610  if (fMakePFPs) {
611  TagCosmics(tpcgeom.ID());
612  MakeFamily(tpcgeom.ID());
613  }
614  FitVertices(detProp, tpcgeom.ID());
615 
616  // Make PFParticles -> tracks
617  for (unsigned short ipf = 0; ipf < pfpToTrkID.size(); ++ipf) {
618  // trackID of this PFP
619  tID = pfpToTrkID[ipf];
620  if (tID > trk.size() + 1) {
621  mf::LogWarning("CCTM") << "Bad track ID";
622  continue;
623  }
624  dtrIndices.clear();
625  // load the daughter PFP indices
626  mf::LogVerbatim("CCTM") << "PFParticle " << ipf << " tID " << tID;
627  for (unsigned short jpf = 0; jpf < pfpToTrkID.size(); ++jpf) {
628  itr = pfpToTrkID[jpf] - 1; // convert from track ID to track index
629  if (trk[itr].MomID == tID) dtrIndices.push_back(jpf);
630  if (trk[itr].MomID == tID)
631  mf::LogVerbatim("CCTM") << " dtr jpf " << jpf << " itr " << itr;
632  } // jpf
633  unsigned short parentIndex = USHRT_MAX;
634  if (tID == 0) {
635  // make neutrino PFP USHRT_MAX == primary PFP
636  recob::PFParticle pfp(14, ipf, parentIndex, dtrIndices);
637  pcol->emplace_back(std::move(pfp));
638  for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
639  if (!vtx[ivx].Neutrino) continue;
640  // make the vertex
641  xyz[0] = vtx[ivx].X;
642  xyz[1] = vtx[ivx].Y;
643  xyz[2] = vtx[ivx].Z;
644  size_t vStart = vcol->size();
645  recob::Vertex vertex(xyz, vtx[ivx].ID);
646  vcol->emplace_back(std::move(vertex));
647  size_t vEnd = vcol->size();
648  // PFParticle - vertex association
649  util::CreateAssn(evt, *pcol, *vcol, *pvassn, vStart, vEnd);
650  vtx[ivx].ID = -vtx[ivx].ID;
651  break;
652  } // ivx
653  }
654  else if (tID > 0) {
655  // make non-neutrino PFP. Need to find the parent PFP index
656  // Track index of this PFP
657  tIndex = tID - 1;
658  // trackID of this PFP parent
659  for (unsigned short ii = 0; ii < pfpToTrkID.size(); ++ii) {
660  if (pfpToTrkID[ii] == trk[tIndex].MomID) {
661  parentIndex = ii;
662  break;
663  }
664  } // ii
665  // PFParticle
666  mf::LogVerbatim("CCTM") << "daughters tID " << tID << " pdg " << trk[tIndex].PDGCode
667  << " ipf " << ipf << " parentIndex " << parentIndex
668  << " dtr size " << dtrIndices.size();
669  recob::PFParticle pfp(trk[tIndex].PDGCode, ipf, parentIndex, dtrIndices);
670  pcol->emplace_back(std::move(pfp));
671  // track
672  // make the track
673  size_t tStart = tcol->size();
677  recob::Track::Flags_t(trk[itr].TrjPos.size()),
678  false),
679  0,
680  -1.,
681  0,
684  tID);
685  tcol->emplace_back(std::move(track));
686  size_t tEnd = tcol->size();
687  // PFParticle - track association
688  util::CreateAssn(evt, *pcol, *tcol, *ptassn, tStart, tEnd);
689  // flag this track as already put in the event
690  trk[tIndex].ID = -trk[tIndex].ID;
691  // track -> hits association
692  tmpHits.clear();
693  for (ipl = 0; ipl < nplanes; ++ipl)
694  tmpHits.insert(
695  tmpHits.end(), trk[tIndex].TrkHits[ipl].begin(), trk[tIndex].TrkHits[ipl].end());
696  util::CreateAssn(evt, *tcol, tmpHits, *thassn);
697  // Find seed hits and the end of the track that is best
698  end = 0;
699  unsigned short itj = 0;
700  if (end > 0) itj = trk[tIndex].TrjPos.size() - 1;
701  for (unsigned short ii = 0; ii < 3; ++ii) {
702  sPos[ii] = trk[tIndex].TrjPos[itj](ii);
703  sDir[ii] = trk[tIndex].TrjDir[itj](ii);
704  } // ii
705  size_t sStart = scol->size();
706  recob::Seed seed(sPos, sDir, sErr, sErr);
707  scol->emplace_back(std::move(seed));
708  size_t sEnd = scol->size();
709  // PFP-seed association
710  util::CreateAssn(evt, *pcol, *scol, *psassn, sStart, sEnd);
711  // Seed-hit association
712  tmpHits.clear();
713  for (ipl = 0; ipl < nplanes; ++ipl)
714  tmpHits.insert(tmpHits.end(), seedHits[ipl].begin(), seedHits[ipl].end());
715  util::CreateAssn(evt, *scol, tmpHits, *shassn);
716  // cluster association
717  // PFP-cluster association
718  tmpCls.clear();
719  for (unsigned short ii = 0; ii < trk[tIndex].ClsEvtIndices.size(); ++ii) {
720  icl = trk[tIndex].ClsEvtIndices[ii];
721  tmpCls.push_back(clusterlist[icl]);
722  } // ii
723  util::CreateAssn(evt, *pcol, tmpCls, *pcassn);
724  } // non-neutrino PFP
725  } // ipf
726  // make non-PFP tracks
727  for (unsigned short itr = 0; itr < trk.size(); ++itr) {
728  // ignore already saved tracks
729  if (trk[itr].ID < 0) continue;
733  recob::Track::Flags_t(trk[itr].TrjPos.size()),
734  false),
735  0,
736  -1.,
737  0,
740  trk[itr].ID);
741  tcol->emplace_back(std::move(track));
742  tmpHits.clear();
743  for (ipl = 0; ipl < nplanes; ++ipl)
744  tmpHits.insert(tmpHits.end(), trk[itr].TrkHits[ipl].begin(), trk[itr].TrkHits[ipl].end());
745  util::CreateAssn(evt, *tcol, tmpHits, *thassn);
746  } // itr
747  // make remnant vertices
748  for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
749  if (vtx[ivx].ID < 0) continue;
750  xyz[0] = vtx[ivx].X;
751  xyz[1] = vtx[ivx].Y;
752  xyz[2] = vtx[ivx].Z;
753  recob::Vertex vertex(xyz, vtx[ivx].ID);
754  vcol->emplace_back(std::move(vertex));
755  }
756  if (fDebugAlg > 0) PrintTracks();
757 
758  double orphanLen = 0;
759  for (ipl = 0; ipl < nplanes; ++ipl) {
760  for (icl = 0; icl < cls[ipl].size(); ++icl) {
761  if (cls[ipl][icl].Length > 40 && cls[ipl][icl].InTrack < 0) {
762  orphanLen += cls[ipl][icl].Length;
763  // unused cluster
764  mf::LogVerbatim("CCTM")
765  << "Orphan long cluster " << ipl << ":" << icl << ":" << cls[ipl][icl].Wire[0] << ":"
766  << (int)cls[ipl][icl].Time[0] << " length " << cls[ipl][icl].Length;
767  }
768  } // icl
769 
770  cls[ipl].clear();
771  clsChain[ipl].clear();
772  } // ipl
773  std::cout << "Total orphan length " << orphanLen << "\n";
774  trkHits[ipl].clear();
775  seedHits[ipl].clear();
776  vxCls[ipl].clear();
777  } // tpcs
778 
779  evt.put(std::move(pcol));
780  evt.put(std::move(ptassn));
781  evt.put(std::move(pcassn));
782  evt.put(std::move(pvassn));
783  evt.put(std::move(psassn));
784  evt.put(std::move(tcol));
785  evt.put(std::move(thassn));
786  evt.put(std::move(scol));
787  evt.put(std::move(vcol));
788 
789  // final cleanup
790  vtx.clear();
791  trk.clear();
792  allhits.clear();
793  matcomb.clear();
794  pfpToTrkID.clear();
795 
796  } // produce
797 
800  geo::TPCID const& tpcid)
801  {
802  std::vector<std::vector<geo::WireID>> hitWID;
803  std::vector<std::vector<double>> hitX;
804  std::vector<std::vector<double>> hitXErr;
805  TVector3 pos, posErr;
806  std::vector<TVector3> trkDir;
807  std::vector<TVector3> trkDirErr;
808  float ChiDOF;
809 
810  if (fNVtxTrkHitsFit == 0) return;
811 
812  unsigned short indx, indx2, iht, nHitsFit;
813 
814  for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
815  if (!vtx[ivx].Neutrino) continue;
816  hitWID.clear();
817  hitX.clear();
818  hitXErr.clear();
819  trkDir.clear();
820  // find the tracks associated with this vertex
821  unsigned int thePln, theTPC, theCst;
822  for (unsigned short itk = 0; itk < trk.size(); ++itk) {
823  for (unsigned short end = 0; end < 2; ++end) {
824  if (trk[itk].VtxIndex[end] != ivx) continue;
825  unsigned short itj = 0;
826  if (end == 1) itj = trk[itk].TrjPos.size() - 1;
827  // increase the size of the hit vectors by 1 for this track
828  indx = hitX.size();
829  hitWID.resize(indx + 1);
830  hitX.resize(indx + 1);
831  hitXErr.resize(indx + 1);
832  trkDir.resize(indx + 1);
833  trkDir[indx] = trk[itk].TrjDir[itj];
834  auto const nplanes = geom->Nplanes(tpcid);
835  for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
836  if (trk[itk].TrkHits[ipl].size() < 2) continue;
837  // make slots for the hits on this track in this plane
838  nHitsFit = trk[itk].TrkHits[ipl].size();
839  if (nHitsFit > fNVtxTrkHitsFit) nHitsFit = fNVtxTrkHitsFit;
840  indx2 = hitWID[indx].size();
841  hitWID[indx].resize(indx2 + nHitsFit);
842  hitX[indx].resize(indx2 + nHitsFit);
843  hitXErr[indx].resize(indx2 + nHitsFit);
844  for (unsigned short ii = 0; ii < nHitsFit; ++ii) {
845  if (end == 0) { iht = ii; }
846  else {
847  iht = trk[itk].TrkHits[ipl].size() - ii - 1;
848  }
849  hitWID[indx][indx2 + ii] = trk[itk].TrkHits[ipl][iht]->WireID();
850  thePln = trk[itk].TrkHits[ipl][iht]->WireID().Plane;
851  theTPC = trk[itk].TrkHits[ipl][iht]->WireID().TPC;
852  theCst = trk[itk].TrkHits[ipl][iht]->WireID().Cryostat;
853  hitX[indx][indx2 + ii] = detProp.ConvertTicksToX(
854  trk[itk].TrkHits[ipl][iht]->PeakTime(), thePln, theTPC, theCst);
855  hitXErr[indx][indx2 + ii] = fHitFitErrFac * trk[itk].TrkHits[ipl][iht]->RMS();
856  } // ii
857  } // ipl
858  } // end
859  } // itk
860  if (hitX.size() < 2) {
861  mf::LogVerbatim("CCTM") << "Not enough hits to fit vtx " << ivx;
862  continue;
863  } // hitX.size() < 2
864  pos(0) = vtx[ivx].X;
865  pos(1) = vtx[ivx].Y;
866  pos(2) = vtx[ivx].Z;
867  fVertexFitAlg.VertexFit(hitWID, hitX, hitXErr, pos, posErr, trkDir, trkDirErr, ChiDOF);
868  if (ChiDOF > 3000) continue;
869  // update the vertex position
870  vtx[ivx].X = pos(0);
871  vtx[ivx].Y = pos(1);
872  vtx[ivx].Z = pos(2);
873  // and the track trajectory
874  unsigned short fitTrk = 0;
875  for (unsigned short itk = 0; itk < trk.size(); ++itk) {
876  for (unsigned short end = 0; end < 2; ++end) {
877  if (trk[itk].VtxIndex[end] != ivx) continue;
878  unsigned short itj = 0;
879  if (end == 1) itj = trk[itk].TrjPos.size() - 1;
880  trk[itk].TrjDir[itj] = trkDir[fitTrk];
881  ++fitTrk;
882  } // end
883  } // itk
884  } // ivx
885  } // FitVertices
886 
889  geo::TPCID const& tpcid)
890  {
891  // fills the CHgNear array
892 
893  unsigned short wire, nwires, indx;
894  float dir, ctime, cx, chgWinLo, chgWinHi;
895  float cnear;
896 
897  for (auto const& planeid : geom->Iterate<geo::PlaneID>(tpcid)) {
898  auto const ipl = planeid.Plane;
899  for (unsigned short icl = 0; icl < cls[ipl].size(); ++icl) {
900  // find the nearby charge at the US and DS ends
901  nwires = cls[ipl][icl].Length / 2;
902  if (nwires < 2) continue;
903  // maximum of 30 wires for long clusters
904  if (nwires > 30) nwires = 30;
905  for (unsigned short end = 0; end < 2; ++end) {
906  cnear = 0;
907  // direction for adding/subtracting wire numbers
908  dir = 1 - 2 * end;
909  for (unsigned short w = 0; w < nwires; ++w) {
910  wire = cls[ipl][icl].Wire[end] + dir * w;
911  cx = cls[ipl][icl].X[end] + dir * w * cls[ipl][icl].Slope[end] * fWirePitch;
912  ctime = detProp.ConvertXToTicks(cx, planeid);
913  chgWinLo = ctime - fChgWindow;
914  chgWinHi = ctime + fChgWindow;
915  indx = wire - firstWire[ipl];
916  if (WireHitRange[ipl][indx].first < 0) continue;
917  unsigned int firhit = WireHitRange[ipl][indx].first;
918  unsigned int lashit = WireHitRange[ipl][indx].second;
919  for (unsigned int hit = firhit; hit < lashit; ++hit) {
920  if (hit > allhits.size() - 1) {
921  mf::LogError("CCTM")
922  << "FillChgNear bad lashit " << lashit << " size " << allhits.size() << "\n";
923  continue;
924  }
925  if (allhits[hit]->PeakTime() < chgWinLo) continue;
926  if (allhits[hit]->PeakTime() > chgWinHi) continue;
927  cnear += allhits[hit]->Integral();
928  } // hit
929  } // w
930  cnear /= (float)(nwires - 1);
931  if (cnear > cls[ipl][icl].Charge[end]) {
932  cls[ipl][icl].ChgNear[end] = ChgNorm[ipl] * cnear / cls[ipl][icl].Charge[end];
933  }
934  else {
935  cls[ipl][icl].ChgNear[end] = 0;
936  }
937  } // end
938  } // icl
939  } // ipl
940 
941  } // FillChgNear
942 
945  {
946  // define the track/vertex and mother/daughter relationships
947 
948  unsigned short ivx, itr, ipl, ii, jtr;
949  unsigned short nus, nds, nuhs, ndhs;
950  float longUSTrk, longDSTrk, qual;
951 
952  // min distance^2 between a neutrino vertex candidate and a through
953  // going muon
954  float tgMuonCut2 = 9;
955 
956  // struct for neutrino vertex candidates
957  struct NuVtx {
958  unsigned short VtxIndex;
959  unsigned short nUSTk;
960  unsigned short nDSTk;
961  unsigned short nUSHit;
962  unsigned short nDSHit;
963  float longUSTrk;
964  float longDSTrk;
965  float Qual;
966  };
967  std::vector<NuVtx> nuVtxCand;
968 
969  NuVtx aNuVtx;
970 
971  // analyze all of the vertices
972  float best = 999, dx, dy, dz, dr;
973  short imbest = -1;
974  bool skipVtx;
975  unsigned short itj;
976  auto const nplanes = geom->Nplanes(tpcid);
977  for (ivx = 0; ivx < vtx.size(); ++ivx) {
978  vtx[ivx].Neutrino = false;
979  nus = 0;
980  nds = 0;
981  nuhs = 0;
982  ndhs = 0;
983  longUSTrk = 0;
984  longDSTrk = 0;
985  skipVtx = false;
986  // skip vertices that are close to through-going muons
987  for (itr = 0; itr < trk.size(); ++itr) {
988  if (trk[itr].ID < 0) continue;
989  if (trk[itr].PDGCode != 13) continue;
990  for (itj = 0; itj < trk[itr].TrjPos.size(); ++itj) {
991  dx = trk[itr].TrjPos[itj](0) - vtx[ivx].X;
992  dy = trk[itr].TrjPos[itj](1) - vtx[ivx].Y;
993  dz = trk[itr].TrjPos[itj](2) - vtx[ivx].Z;
994  dr = dx * dx + dy * dy + dz * dz;
995  if (dr < tgMuonCut2) {
996  skipVtx = true;
997  break;
998  }
999  if (skipVtx) break;
1000  } // itj
1001  if (skipVtx) break;
1002  } // itr
1003  if (skipVtx) continue;
1004  for (itr = 0; itr < trk.size(); ++itr) {
1005  if (trk[itr].ID < 0) continue;
1006  if (trk[itr].VtxIndex[0] == ivx) {
1007  // DS-going track
1008  ++nds;
1009  if (trk[itr].Length > longDSTrk) longDSTrk = trk[itr].Length;
1010  for (ipl = 0; ipl < nplanes; ++ipl)
1011  ndhs += trk[itr].TrkHits[ipl].size();
1012  } // DS-going track
1013  // Reject this vertex as a neutrino candidate if the track being
1014  // considered has a starting vertex
1015  if (trk[itr].VtxIndex[1] == ivx && trk[itr].VtxIndex[0] >= 0) {
1016  skipVtx = true;
1017  break;
1018  } // trk[itr].VtxIndex[0] > 0
1019  if (trk[itr].VtxIndex[1] == ivx && trk[itr].VtxIndex[0] < 0) {
1020  // US-going track w no US vertex
1021  ++nus;
1022  if (trk[itr].Length > longUSTrk) longUSTrk = trk[itr].Length;
1023  for (ipl = 0; ipl < nplanes; ++ipl)
1024  nuhs += trk[itr].TrkHits[ipl].size();
1025  } // US-going track
1026  } // itr
1027  if (skipVtx) continue;
1028  if (nds == 0) continue;
1029  qual = 1 / (float)nds;
1030  qual /= (float)ndhs;
1031  if (nus > 0) qual *= (float)nuhs / (float)ndhs;
1032  if (qual < best) {
1033  best = qual;
1034  imbest = ivx;
1035  }
1036  if (nds > 0 && longDSTrk > 5) {
1037  // at least one longish track going DS
1038  aNuVtx.VtxIndex = ivx;
1039  aNuVtx.nUSTk = nus;
1040  aNuVtx.nDSTk = nds;
1041  aNuVtx.nUSHit = nuhs;
1042  aNuVtx.nDSHit = ndhs;
1043  aNuVtx.longUSTrk = longUSTrk;
1044  aNuVtx.longDSTrk = longDSTrk;
1045  aNuVtx.Qual = qual;
1046  nuVtxCand.push_back(aNuVtx);
1047  }
1048  } // ivx
1049  if (imbest < 0) return;
1050 
1051  // Found the neutrino interaction vertex
1052  ivx = imbest;
1053  vtx[ivx].Neutrino = true;
1054  // Put a 0 into the PFParticle -> track vector to identify a neutrino
1055  // track with no trajectory or hits
1056  pfpToTrkID.push_back(0);
1057 
1058  // list of DS-going current generation daughters so we can fix next
1059  // generation daughter assignments. This code section assumes that there
1060  // are no decays or 2ndry interactions with US-going primary tracks.
1061  std::vector<unsigned short> dtrGen;
1062  std::vector<unsigned short> dtrNextGen;
1063  for (itr = 0; itr < trk.size(); ++itr) {
1064  if (trk[itr].ID < 0) continue;
1065  if (trk[itr].VtxIndex[0] == ivx) {
1066  // DS-going primary track
1067  trk[itr].MomID = 0;
1068  // call every track coming from a neutrino vertex a proton
1069  trk[itr].PDGCode = 2212;
1070  pfpToTrkID.push_back(trk[itr].ID);
1071  dtrGen.push_back(itr);
1072  } // DS-going primary track
1073  if (trk[itr].VtxIndex[1] == ivx) {
1074  // US-going primary track
1075  trk[itr].MomID = 0;
1076  // call every track coming from a neutrino vertex a proton
1077  trk[itr].PDGCode = 2212;
1078  pfpToTrkID.push_back(trk[itr].ID);
1079  // reverse the track trajectory
1080  std::reverse(trk[itr].TrjPos.begin(), trk[itr].TrjPos.end());
1081  for (ii = 0; ii < trk[itr].TrjDir.size(); ++ii)
1082  trk[itr].TrjDir[ii] = -trk[itr].TrjDir[ii];
1083  } // DS-going primary track
1084  } // itr
1085 
1086  if (dtrGen.size() == 0) return;
1087 
1088  unsigned short tmp, indx;
1089  unsigned short nit = 0;
1090 
1091  // follow daughters through all generations (< 10). Daughter tracks
1092  // may go US or DS
1093  while (nit < 10) {
1094  ++nit;
1095  dtrNextGen.clear();
1096  // Look for next generation daughters
1097  for (ii = 0; ii < dtrGen.size(); ++ii) {
1098  itr = dtrGen[ii];
1099  if (trk[itr].VtxIndex[1] >= 0) {
1100  // found a DS vertex
1101  ivx = trk[itr].VtxIndex[1];
1102  // look for a track associated with this vertex
1103  for (jtr = 0; jtr < trk.size(); ++jtr) {
1104  if (jtr == itr) continue;
1105  if (trk[jtr].VtxIndex[0] == ivx) {
1106  // DS-going track
1107  indx = trk[itr].DtrID.size();
1108  trk[itr].DtrID.resize(indx + 1);
1109  trk[itr].DtrID[indx] = jtr;
1110  trk[jtr].MomID = trk[itr].ID;
1111  // call all daughters pions
1112  trk[jtr].PDGCode = 211;
1113  dtrNextGen.push_back(jtr);
1114  pfpToTrkID.push_back(trk[jtr].ID);
1115  } // DS-going track
1116  if (trk[jtr].VtxIndex[1] == ivx) {
1117  // US-going track
1118  indx = trk[itr].DtrID.size();
1119  trk[itr].DtrID.resize(indx + 1);
1120  trk[itr].DtrID[indx] = jtr;
1121  trk[jtr].MomID = trk[itr].ID;
1122  // call all daughters pions
1123  trk[jtr].PDGCode = 211;
1124  dtrNextGen.push_back(jtr);
1125  pfpToTrkID.push_back(trk[jtr].ID);
1126  // reverse the track trajectory
1127  std::reverse(trk[jtr].TrjPos.begin(), trk[jtr].TrjPos.end());
1128  for (unsigned short jj = 0; jj < trk[jtr].TrjDir.size(); ++jj)
1129  trk[jtr].TrjDir[jj] = -trk[jtr].TrjDir[jj];
1130  // interchange the trk - vtx assignments
1131  tmp = trk[jtr].VtxIndex[0];
1132  trk[jtr].VtxIndex[0] = trk[jtr].VtxIndex[1];
1133  trk[jtr].VtxIndex[1] = tmp;
1134  } // DS-going track
1135  } // jtr
1136  } // trk[itr].VtxIndex[0] >= 0
1137  } // ii (itr)
1138  // break out if no next gen daughters found
1139  if (dtrNextGen.size() == 0) break;
1140  dtrGen = dtrNextGen;
1141  } // nit
1142 
1143  } // MakeFamily
1144 
1147  {
1148  // Make cosmic ray PFParticles
1149  unsigned short ipf, itj;
1150  bool skipit = true;
1151 
1152  // Y,Z limits of the detector
1153  const geo::TPCGeo& thetpc = geom->TPC(tpcid);
1154  auto const world = thetpc.GetCenter();
1155 
1156  float XLo = world.X() - thetpc.HalfWidth() + fFiducialCut;
1157  float XHi = world.X() + thetpc.HalfWidth() - fFiducialCut;
1158  float YLo = world.Y() - thetpc.HalfHeight() + fFiducialCut;
1159  float YHi = world.Y() + thetpc.HalfHeight() - fFiducialCut;
1160  float ZLo = world.Z() - thetpc.Length() / 2 + fFiducialCut;
1161  float ZHi = world.Z() + thetpc.Length() / 2 - fFiducialCut;
1162 
1163  bool startsIn, endsIn;
1164 
1165  for (unsigned short itk = 0; itk < trk.size(); ++itk) {
1166  // ignore already used tracks
1167  if (trk[itk].ID < 0) continue;
1168  // ignore short tracks (< 10 cm)
1169  if (trk[itk].Length < 10) continue;
1170  // check for already identified PFPs
1171  skipit = false;
1172  for (ipf = 0; ipf < pfpToTrkID.size(); ++ipf) {
1173  if (pfpToTrkID[ipf] == trk[itk].ID) {
1174  skipit = true;
1175  break;
1176  }
1177  } // ipf
1178  if (skipit) continue;
1179  startsIn = true;
1180  if (trk[itk].TrjPos[0](0) < XLo || trk[itk].TrjPos[0](0) > XHi) startsIn = false;
1181  if (trk[itk].TrjPos[0](1) < YLo || trk[itk].TrjPos[0](1) > YHi) startsIn = false;
1182  if (trk[itk].TrjPos[0](2) < ZLo || trk[itk].TrjPos[0](2) > ZHi) startsIn = false;
1183  if (startsIn) continue;
1184  endsIn = true;
1185  itj = trk[itk].TrjPos.size() - 1;
1186  if (trk[itk].TrjPos[itj](0) < XLo || trk[itk].TrjPos[itj](0) > XHi) endsIn = false;
1187  if (trk[itk].TrjPos[itj](1) < YLo || trk[itk].TrjPos[itj](1) > YHi) endsIn = false;
1188  if (trk[itk].TrjPos[itj](2) < ZLo || trk[itk].TrjPos[itj](2) > ZHi) endsIn = false;
1189  if (endsIn) continue;
1190  // call it a cosmic muon
1191  trk[itk].PDGCode = 13;
1192  pfpToTrkID.push_back(trk[itk].ID);
1193  } // itk
1194 
1195  if (fDeltaRayCut <= 0) return;
1196 
1197  for (unsigned short itk = 0; itk < trk.size(); ++itk) {
1198  // find a tagged cosmic ray
1199  if (trk[itk].PDGCode != 13) continue;
1200 
1201  } // itk
1202 
1203  } // TagCosmics
1204 
1207  art::FindManyP<recob::Hit> const& fmCluHits,
1208  geo::TPCID const& tpcid)
1209  {
1210  // Use vertex assignments to match clusters
1211  unsigned short ivx, ii, ipl, icl, jj, jpl, jcl, kk, kpl, kcl;
1212  short idir, iend, jdir, jend, kdir, kend, ioend;
1213 
1214  auto const nplanes = geom->Nplanes(tpcid);
1215  for (ivx = 0; ivx < vtx.size(); ++ivx) {
1216 
1217  // list of cluster chains associated with this vertex in each plane
1218  for (ipl = 0; ipl < nplanes; ++ipl) {
1219  vxCls[ipl].clear();
1220  for (icl = 0; icl < clsChain[ipl].size(); ++icl) {
1221  if (clsChain[ipl][icl].InTrack >= 0) continue;
1222  for (iend = 0; iend < 2; ++iend) {
1223  if (clsChain[ipl][icl].VtxIndex[iend] == vtx[ivx].EvtIndex) vxCls[ipl].push_back(icl);
1224  } // end
1225  } // icl
1226  } // ipl
1227 
1228  if (prt) {
1229  mf::LogVerbatim myprt("CCTM");
1230  myprt << "VtxMatch: Vertex ID " << vtx[ivx].EvtIndex << "\n";
1231  for (ipl = 0; ipl < nplanes; ++ipl) {
1232  myprt << "ipl " << ipl << " cls";
1233  for (unsigned short ii = 0; ii < vxCls[ipl].size(); ++ii)
1234  myprt << " " << vxCls[ipl][ii];
1235  myprt << "\n";
1236  } // ipl
1237  } // prt
1238  // match between planes, requiring clusters matched to this vertex
1239  iend = 0;
1240  jend = 0;
1241  bool gotkcl;
1242  float totErr;
1243  for (ipl = 0; ipl < nplanes; ++ipl) {
1244  if (nplanes == 2 && ipl > 0) continue;
1245  for (ii = 0; ii < vxCls[ipl].size(); ++ii) {
1246  icl = vxCls[ipl][ii];
1247  // ignore used clusters
1248  if (clsChain[ipl][icl].InTrack >= 0) continue;
1249  jpl = (ipl + 1) % nplanes;
1250  kpl = (jpl + 1) % nplanes;
1251  for (jj = 0; jj < vxCls[jpl].size(); ++jj) {
1252  jcl = vxCls[jpl][jj];
1253  if (clsChain[jpl][jcl].InTrack >= 0) continue;
1254  for (iend = 0; iend < 2; ++iend) {
1255  if (clsChain[ipl][icl].VtxIndex[iend] != vtx[ivx].EvtIndex) continue;
1256  ioend = 1 - iend;
1257  idir = clsChain[ipl][icl].Dir[iend];
1258  for (jend = 0; jend < 2; ++jend) {
1259  if (clsChain[jpl][jcl].VtxIndex[jend] != vtx[ivx].EvtIndex) continue;
1260  jdir = clsChain[jpl][jcl].Dir[jend];
1261  if (idir != 0 && jdir != 0 && idir != jdir) continue;
1262  // ignore outrageously bad other end X matches
1263  if (fabs(clsChain[jpl][jcl].X[1 - jend] - clsChain[ipl][icl].X[ioend]) > 50)
1264  continue;
1265  MatchPars match;
1266  match.Cls[ipl] = icl;
1267  match.End[ipl] = iend;
1268  match.Cls[jpl] = jcl;
1269  match.End[jpl] = jend;
1270  match.Vtx = ivx;
1271  match.oVtx = -1;
1272  // set large so that DupMatch doesn't get confused when called before FillEndMatch
1273  match.Err = 1E6;
1274  match.oErr = 1E6;
1275  if (nplanes == 2) {
1276  FillEndMatch2(match);
1277  if (prt)
1278  mf::LogVerbatim("CCTM")
1279  << "FillEndMatch2: Err " << match.Err << " oErr " << match.oErr;
1280  if (match.Err + match.oErr > 100) continue;
1281  if (DupMatch(match, nplanes)) continue;
1282  matcomb.push_back(match);
1283  continue;
1284  }
1285  match.Cls[kpl] = -1;
1286  match.End[kpl] = 0;
1287  if (prt)
1288  mf::LogVerbatim("CCTM")
1289  << "VtxMatch: check " << ipl << ":" << icl << ":" << iend << " and " << jpl
1290  << ":" << jcl << ":" << jend << " for cluster in kpl " << kpl;
1291  gotkcl = false;
1292  for (kk = 0; kk < vxCls[kpl].size(); ++kk) {
1293  kcl = vxCls[kpl][kk];
1294  if (clsChain[kpl][kcl].InTrack >= 0) continue;
1295  for (kend = 0; kend < 2; ++kend) {
1296  kdir = clsChain[kpl][kcl].Dir[kend];
1297  if (idir != 0 && kdir != 0 && idir != kdir) continue;
1298  if (clsChain[kpl][kcl].VtxIndex[kend] != vtx[ivx].EvtIndex) continue;
1299  // rough check of other end match
1300  // TODO: SHOWER-LIKE CLUSTER CHECK
1301  match.Cls[kpl] = kcl;
1302  match.End[kpl] = kend;
1303  // first call to ignore redundant matches
1304  if (DupMatch(match, nplanes)) continue;
1305  FillEndMatch(detProp, tpcid, match);
1306  // ignore if no signal at the other end
1307  if (match.Chg[kpl] <= 0) continue;
1308  if (match.Err + match.oErr > 100) continue;
1309  // second call to keep matches with better error
1310  if (DupMatch(match, nplanes)) continue;
1311  matcomb.push_back(match);
1312  gotkcl = true;
1313  // break;
1314  } // kend
1315  } // kk -> kcl
1316  if (gotkcl) continue;
1317  // look for a cluster that missed the vertex assignment
1318  float best = 10;
1319  short kbst = -1;
1320  unsigned short kbend = 0;
1321  if (prt)
1322  mf::LogVerbatim("CCTM") << "VtxMatch: look for missed cluster chain in kpl";
1323  for (kcl = 0; kcl < clsChain[kpl].size(); ++kcl) {
1324  if (clsChain[kpl][kcl].InTrack >= 0) continue;
1325  for (kend = 0; kend < 2; ++kend) {
1326  kdir = clsChain[kpl][kcl].Dir[kend];
1327  if (idir != 0 && kdir != 0 && idir != kdir) continue;
1328  if (clsChain[kpl][kcl].VtxIndex[kend] >= 0) continue;
1329  // make a rough dX cut at the match end
1330  if (fabs(clsChain[kpl][kcl].X[kend] - vtx[ivx].X) > 5) continue;
1331  // and at the other end
1332  if (fabs(clsChain[kpl][kcl].X[1 - kend] - clsChain[ipl][icl].X[ioend]) > 50)
1333  continue;
1334  // check the error
1335  match.Cls[kpl] = kcl;
1336  match.End[kpl] = kend;
1337  if (DupMatch(match, nplanes)) continue;
1338  FillEndMatch(detProp, tpcid, match);
1339  totErr = match.Err + match.oErr;
1340  if (prt) {
1341  mf::LogVerbatim myprt("CCTM");
1342  myprt << "VtxMatch: Chk missing cluster match ";
1343  for (unsigned short ii = 0; ii < nplanes; ++ii)
1344  myprt << " " << ii << ":" << match.Cls[ii] << ":" << match.End[ii];
1345  myprt << " Err " << match.Err << "\n";
1346  }
1347  if (totErr > 100) continue;
1348  if (totErr < best) {
1349  best = totErr;
1350  kbst = kcl;
1351  kbend = kend;
1352  }
1353  } // kend
1354  } // kcl
1355  if (kbst >= 0) {
1356  // found a decent match
1357  match.Cls[kpl] = kbst;
1358  match.End[kpl] = kbend;
1359  FillEndMatch(detProp, tpcid, match);
1360  matcomb.push_back(match);
1361  // assign the vertex to this cluster
1362  clsChain[kpl][kbst].VtxIndex[kbend] = ivx;
1363  // and update vxCls
1364  vxCls[kpl].push_back(kbst);
1365  }
1366  else {
1367  // Try a 2 plane match if a 3 plane match didn't work
1368  match.Cls[kpl] = -1;
1369  match.End[kpl] = 0;
1370  if (DupMatch(match, nplanes)) continue;
1371  FillEndMatch(detProp, tpcid, match);
1372  if (match.Err + match.oErr < 100) matcomb.push_back(match);
1373  }
1374  } // jend
1375  } // iend
1376  } // jj
1377  } // ii -> icl
1378  } // ipl
1379 
1380  if (matcomb.size() == 0) continue;
1381  SortMatches(detProp, fmCluHits, 1, tpcid);
1382 
1383  } // ivx
1384 
1385  for (ipl = 0; ipl < 3; ++ipl)
1386  vxCls[ipl].clear();
1387 
1388  } // VtxMatch
1389 
1392  {
1393  // Project clusters to vertices and fill mVtxIndex. No requirement is
1394  // made that charge exists on the line between the Begin (End) of the
1395  // cluster and the vertex
1396 
1397  if (vtx.empty()) return;
1398 
1399  for (auto const& planeID : geom->Iterate<geo::PlaneID>(tpcid)) {
1400  unsigned int const ipl = planeID.Cryostat;
1401  for (unsigned int icl = 0; icl < cls[ipl].size(); ++icl) {
1402  for (unsigned int end = 0; end < 2u; ++end) {
1403  // ignore already attached clusters
1404  if (cls[ipl][icl].VtxIndex[end] >= 0) continue;
1405  short ibstvx = -1;
1406  float best = 1.;
1407  // index of the other end
1408  unsigned int oend = 1 - end;
1409  for (std::size_t ivx = 0; ivx < vtx.size(); ++ivx) {
1410  // ignore if the other end is attached to this vertex (which can happen with short clusters)
1411  if (cls[ipl][icl].VtxIndex[oend] == static_cast<short>(ivx)) continue;
1412  float const dWire =
1413  geom->WireCoordinate(geo::Point_t{0, vtx[ivx].Y, vtx[ivx].Z}, planeID) -
1414  cls[ipl][icl].Wire[end];
1415  if (end == 0) {
1416  if (dWire > 30 || dWire < -2) continue;
1417  }
1418  else {
1419  if (dWire < -30 || dWire > 2) continue;
1420  }
1421  // project the cluster to the vertex wire
1422  float const dX = fabs(cls[ipl][icl].X[end] +
1423  cls[ipl][icl].Slope[end] * fWirePitch * dWire - vtx[ivx].X);
1424  if (dX < best) {
1425  best = dX;
1426  ibstvx = ivx;
1427  }
1428  } // ivx
1429  if (ibstvx >= 0) {
1430  // attach
1431  cls[ipl][icl].VtxIndex[end] = ibstvx;
1432  cls[ipl][icl].mVtxIndex[end] = ibstvx;
1433  }
1434  } // end
1435  } // icl
1436  } // ipl
1437 
1438  } // FindMaybeVertices
1439 
1442  art::FindManyP<recob::Hit> const& fmCluHits,
1443  geo::TPCID const& tpcid)
1444  {
1445  unsigned short icl, icl1, icl2;
1446  float dw, dx, dWCut, dw1Max, dw2Max;
1447  float dA, dA2, dACut = fMaxDAng, chgAsymCut;
1448  float dXCut, chgasym, mrgErr;
1449  // long straight clusters
1450  bool ls1, ls2;
1451  bool gotprt = false;
1452  auto const nplanes = geom->Nplanes(tpcid);
1453  for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
1454  if (cls[ipl].size() > 1) {
1455  for (icl1 = 0; icl1 < cls[ipl].size() - 1; ++icl1) {
1456  prt = (fDebugAlg == 666 && ipl == fDebugPlane && icl1 == fDebugCluster);
1457  if (prt) gotprt = true;
1458  // maximum delta Wire overlap is 10% of the total length
1459  dw1Max = 0.6 * cls[ipl][icl1].Length;
1460  ls1 = (cls[ipl][icl1].Length > 100 &&
1461  fabs(cls[ipl][icl1].Angle[0] - cls[ipl][icl1].Angle[1]) < 0.04);
1462  for (icl2 = icl1 + 1; icl2 < cls[ipl].size(); ++icl2) {
1463  ls2 = (cls[ipl][icl2].Length > 100 &&
1464  fabs(cls[ipl][icl2].Angle[0] - cls[ipl][icl2].Angle[1]) < 0.04);
1465  dw2Max = 0.6 * cls[ipl][icl2].Length;
1466  // set overlap cut to be the shorter of the two
1467  dWCut = dw1Max;
1468  if (dw2Max < dWCut) dWCut = dw2Max;
1469  // but not exceeding 20 for very long clusters
1470  if (dWCut > 100) dWCut = 100;
1471  if (dWCut < 2) dWCut = 2;
1472  chgAsymCut = fMergeChgAsym;
1473  // Compare end 1 of icl1 with end 0 of icl2
1474 
1475  if (prt)
1476  mf::LogVerbatim("CCTM")
1477  << "MCC P:C:W icl1 " << ipl << ":" << icl1 << ":" << cls[ipl][icl1].Wire[1]
1478  << " vtx " << cls[ipl][icl1].VtxIndex[1] << " ls1 " << ls1 << " icl2 " << ipl << ":"
1479  << icl2 << ":" << cls[ipl][icl2].Wire[0] << " vtx " << cls[ipl][icl2].VtxIndex[0]
1480  << " ls2 " << ls2 << " dWCut " << dWCut;
1481  if (std::abs(cls[ipl][icl1].Wire[1] - cls[ipl][icl2].Wire[0]) > dWCut) continue;
1482  // ignore if the clusters begin/end on the same wire
1483  // if(cls[ipl][icl1].Wire[1] == cls[ipl][icl2].Wire[0]) continue;
1484  // or if the angle exceeds the cut
1485  float af = AngleFactor(cls[ipl][icl1].Slope[1]);
1486  dACut = fMaxDAng * af;
1487  dXCut = fChainMaxdX * 5 * af;
1488  dA = fabs(cls[ipl][icl1].Angle[1] - cls[ipl][icl2].Angle[0]);
1489  // compare the match angle at the opposite ends of the clusters.
1490  // May have a bad end/begin angle if there is a delta-ray in the middle
1491  dA2 = fabs(cls[ipl][icl1].Angle[0] - cls[ipl][icl2].Angle[1]);
1492 
1493  if (prt)
1494  mf::LogVerbatim("CCTM")
1495  << " dA " << dA << " dA2 " << dA2 << " DACut " << dACut << " dXCut " << dXCut;
1496 
1497  if (dA2 < dA) dA = dA2;
1498  // ignore vertices that have only two associated clusters and
1499  // the angle is small
1500  if (dA < fChainVtxAng && cls[ipl][icl1].VtxIndex[1] >= 0) {
1501  dw = fWirePitch * (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1]);
1502  dx = cls[ipl][icl1].X[1] + cls[ipl][icl1].Slope[1] * dw * fWirePitch -
1503  cls[ipl][icl2].X[0];
1504  unsigned short ivx = cls[ipl][icl1].VtxIndex[1];
1505  if (vtx[ivx].nClusInPln[ipl] == 2 && fabs(dx) < 1) {
1506  cls[ipl][icl1].VtxIndex[1] = -2;
1507  cls[ipl][icl2].VtxIndex[0] = -2;
1508  vtx[ivx].nClusInPln[ipl] = 0;
1509  if (prt) mf::LogVerbatim("CCTM") << " clobbered vertex " << ivx;
1510  } // vertices match
1511  } // dA < 0.1 && ...
1512 
1513  // don't merge if a vertex exists at these ends
1514  if (cls[ipl][icl1].VtxIndex[1] >= 0) continue;
1515  if (cls[ipl][icl2].VtxIndex[0] >= 0) continue;
1516 
1517  // expand the angle match cut for clusters that appear to be stopping
1518  if (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1] < 3 &&
1519  (cls[ipl][icl1].Length < 3 || cls[ipl][icl2].Length < 3)) {
1520  if (prt) mf::LogVerbatim("CCTM") << "Stopping cluster";
1521  dACut *= 1.5;
1522  chgAsymCut *= 1.5;
1523  dXCut *= 3;
1524  } // stopping US cluster
1525 
1526  // find the angle made by the endpoints of the clusters
1527  dw = fWirePitch * (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1]);
1528  if (dw != 0) {
1529  dx = cls[ipl][icl2].X[0] - cls[ipl][icl1].X[1];
1530  float dA3 = std::abs(atan(dx / dw) - cls[ipl][icl1].Angle[1]);
1531  if (prt) mf::LogVerbatim("CCTM") << " dA3 " << dA3;
1532  if (dA3 > dA) dA = dA3;
1533  }
1534 
1535  // angle matching
1536  if (dA > dACut) continue;
1537 
1538  if (prt)
1539  mf::LogVerbatim("CCTM")
1540  << " rough dX " << fabs(cls[ipl][icl1].X[1] - cls[ipl][icl2].X[0]) << " cut = 20";
1541 
1542  // make a rough dX cut
1543  if (fabs(cls[ipl][icl1].X[1] - cls[ipl][icl2].X[0]) > 20) continue;
1544 
1545  // handle cosmic ray clusters that are broken at delta rays
1546  if (ls1 || ls2) {
1547  // tighter angle cuts but no charge cut
1548  if (dA > fChainVtxAng) continue;
1549  }
1550  else {
1551  chgasym = fabs(cls[ipl][icl1].Charge[1] - cls[ipl][icl2].Charge[0]);
1552  chgasym /= cls[ipl][icl1].Charge[1] + cls[ipl][icl2].Charge[0];
1553  if (prt) mf::LogVerbatim("CCTM") << " chgasym " << chgasym << " cut " << chgAsymCut;
1554  if (chgasym > chgAsymCut) continue;
1555  } // ls1 || ls2
1556  // project the longer cluster to the end of the shorter one
1557  if (cls[ipl][icl1].Length > cls[ipl][icl2].Length) {
1558  dw = fWirePitch * (cls[ipl][icl2].Wire[0] - cls[ipl][icl1].Wire[1]);
1559  dx = cls[ipl][icl1].X[1] + cls[ipl][icl1].Slope[1] * dw * fWirePitch -
1560  cls[ipl][icl2].X[0];
1561  }
1562  else {
1563  dw = fWirePitch * (cls[ipl][icl1].Wire[1] - cls[ipl][icl2].Wire[0]);
1564  dx = cls[ipl][icl2].X[0] + cls[ipl][icl2].Slope[0] * dw * fWirePitch -
1565  cls[ipl][icl1].X[1];
1566  }
1567 
1568  // handle overlapping clusters
1569  if (dA2 < 0.01 && abs(dx) > dXCut && dx < -1) {
1570  dx = dXClTraj(fmCluHits, ipl, icl1, 1);
1571  if (prt) mf::LogVerbatim("CCTM") << " new dx from dXClTraj " << dx;
1572  }
1573 
1574  if (prt)
1575  mf::LogVerbatim("CCTM")
1576  << " X0 " << cls[ipl][icl1].X[1] << " slp " << cls[ipl][icl1].Slope[1] << " dw "
1577  << dw << " oX " << cls[ipl][icl2].X[0] << " dx " << dx << " cut " << dXCut;
1578 
1579  if (fabs(dx) > dXCut) continue;
1580 
1581  // calculate a merge error that will be used to adjudicate between multiple merge attempts AngleFactor
1582  float xerr = dx / dXCut;
1583  float aerr = dA / dACut;
1584  mrgErr = xerr * xerr + aerr * aerr;
1585 
1586  if (prt)
1587  mf::LogVerbatim("CCTM")
1588  << "icl1 mrgErr " << mrgErr << " MergeError " << cls[ipl][icl1].MergeError[1]
1589  << " icl2 MergeError " << cls[ipl][icl2].MergeError[0];
1590 
1591  // this merge better than a previous one?
1592  if (mrgErr > cls[ipl][icl1].MergeError[1]) continue;
1593  if (mrgErr > cls[ipl][icl2].MergeError[0]) continue;
1594 
1595  // un-merge icl1 - this should always be true but check anyway
1596  if (cls[ipl][icl1].BrkIndex[1] >= 0) {
1597  unsigned short ocl = cls[ipl][icl1].BrkIndex[1];
1598  if (prt) mf::LogVerbatim("CCTM") << "clobber old icl1 BrkIndex " << ocl;
1599  if (cls[ipl][ocl].BrkIndex[0] == icl1) {
1600  cls[ipl][ocl].BrkIndex[0] = -1;
1601  cls[ipl][ocl].MergeError[0] = fMaxMergeError;
1602  }
1603  if (cls[ipl][ocl].BrkIndex[1] == icl1) {
1604  cls[ipl][ocl].BrkIndex[1] = -1;
1605  cls[ipl][ocl].MergeError[1] = fMaxMergeError;
1606  }
1607  } // cls[ipl][icl1].BrkIndex[1] >= 0
1608  cls[ipl][icl1].BrkIndex[1] = icl2;
1609  cls[ipl][icl1].MergeError[1] = mrgErr;
1610 
1611  // un-merge icl2
1612  if (cls[ipl][icl2].BrkIndex[0] >= 0) {
1613  unsigned short ocl = cls[ipl][icl2].BrkIndex[0];
1614  if (prt) mf::LogVerbatim("CCTM") << "clobber old icl2 BrkIndex " << ocl;
1615  if (cls[ipl][ocl].BrkIndex[0] == icl1) {
1616  cls[ipl][ocl].BrkIndex[0] = -1;
1617  cls[ipl][ocl].MergeError[0] = fMaxMergeError;
1618  }
1619  if (cls[ipl][ocl].BrkIndex[1] == icl1) {
1620  cls[ipl][ocl].BrkIndex[1] = -1;
1621  cls[ipl][ocl].MergeError[1] = fMaxMergeError;
1622  }
1623  } // cls[ipl][icl2].BrkIndex[0] >= 0
1624  cls[ipl][icl2].BrkIndex[0] = icl1;
1625  cls[ipl][icl2].MergeError[0] = mrgErr;
1626  if (prt) mf::LogVerbatim("CCTM") << " merge " << icl1 << " and " << icl2;
1627 
1628  } // icl2
1629  } // icl1
1630 
1631  // look for broken clusters in which have a C shape similar to a Begin-Begin vertex. The clusters
1632  // will have large and opposite sign angles
1633  bool gotone;
1634  for (icl1 = 0; icl1 < cls[ipl].size() - 1; ++icl1) {
1635  gotone = false;
1636  for (icl2 = icl1 + 1; icl2 < cls[ipl].size(); ++icl2) {
1637  // check both ends
1638  for (unsigned short end = 0; end < 2; ++end) {
1639  // Ignore already identified broken clusters
1640  if (cls[ipl][icl1].BrkIndex[end] >= 0) continue;
1641  if (cls[ipl][icl2].BrkIndex[end] >= 0) continue;
1642  // require a large angle cluster
1643  if (fabs(cls[ipl][icl1].Angle[end]) < 1) continue;
1644  // and a second large angle cluster
1645  if (fabs(cls[ipl][icl2].Angle[end]) < 1) continue;
1646  if (prt)
1647  mf::LogVerbatim("CCTM")
1648  << "BrokenC: clusters " << cls[ipl][icl1].Wire[end] << ":"
1649  << (int)cls[ipl][icl1].Time[end] << " " << cls[ipl][icl2].Wire[end] << ":"
1650  << (int)cls[ipl][icl2].Time[end] << " angles " << cls[ipl][icl1].Angle[end] << " "
1651  << cls[ipl][icl2].Angle[end] << " dWire "
1652  << fabs(cls[ipl][icl1].Wire[end] - cls[ipl][icl2].Wire[end]);
1653  if (fabs(cls[ipl][icl1].Wire[end] - cls[ipl][icl2].Wire[end]) > 5) continue;
1654  // This is really crude but maybe OK
1655  // project 1 -> 2
1656  float dsl = cls[ipl][icl2].Slope[end] - cls[ipl][icl1].Slope[end];
1657  float fvw =
1658  (cls[ipl][icl1].X[end] - cls[ipl][icl1].Wire[end] * cls[ipl][icl1].Slope[end] -
1659  cls[ipl][icl2].X[end] + cls[ipl][icl2].Wire[end] * cls[ipl][icl2].Slope[end]) /
1660  dsl;
1661  if (prt) mf::LogVerbatim("CCTM") << " fvw " << fvw;
1662  if (fabs(cls[ipl][icl1].Wire[end] - fvw) > 4) continue;
1663  if (fabs(cls[ipl][icl2].Wire[end] - fvw) > 4) continue;
1664  cls[ipl][icl1].BrkIndex[end] = icl2;
1665  // TODO This could use some improvement if necessary
1666  cls[ipl][icl1].MergeError[end] = 1;
1667  cls[ipl][icl2].BrkIndex[end] = icl1;
1668  cls[ipl][icl2].MergeError[end] = 1;
1669  gotone = true;
1670  dx = fabs(cls[ipl][icl1].X[end] - cls[ipl][icl2].X[end]);
1671  if (prt)
1672  mf::LogVerbatim("CCTM")
1673  << "BrokenC: icl1:W " << icl1 << ":" << cls[ipl][icl1].Wire[end] << " icl2:W "
1674  << icl2 << ":" << cls[ipl][icl2].Wire[end] << " end " << end << " dx " << dx;
1675  } // end
1676  if (gotone) break;
1677  } // icl2
1678  } // icl1
1679 
1680  } // cls[ipl].size() > 1
1681 
1682  // follow mother-daughter broken clusters and put them in the cluster chain array
1683  unsigned short end, mom, momBrkEnd, dtrBrkEnd, nit;
1684  short dtr;
1685 
1686  std::vector<bool> gotcl(cls[ipl].size());
1687  for (icl = 0; icl < cls[ipl].size(); ++icl)
1688  gotcl[icl] = false;
1689  if (prt)
1690  mf::LogVerbatim("CCTM") << "ipl " << ipl << " cls.size() " << cls[ipl].size() << "\n";
1691 
1692  std::vector<unsigned short> sCluster;
1693  std::vector<unsigned short> sOrder;
1694  for (icl = 0; icl < cls[ipl].size(); ++icl) {
1695  sCluster.clear();
1696  sOrder.clear();
1697  if (gotcl[icl]) continue;
1698  // don't start with a cluster broken at both ends
1699  if (cls[ipl][icl].BrkIndex[0] >= 0 && cls[ipl][icl].BrkIndex[1] >= 0) continue;
1700  for (end = 0; end < 2; ++end) {
1701  if (cls[ipl][icl].BrkIndex[end] < 0) continue;
1702  if (cls[ipl][icl].MergeError[end] > fMergeErrorCut) continue;
1703  gotcl[icl] = true;
1704  mom = icl;
1705  // end where the mom is broken
1706  momBrkEnd = end;
1707  sCluster.push_back(mom);
1708  if (momBrkEnd == 1) {
1709  // typical case - broken at the DS end
1710  sOrder.push_back(0);
1711  }
1712  else {
1713  // broken at the US end
1714  sOrder.push_back(1);
1715  }
1716  dtr = cls[ipl][icl].BrkIndex[end];
1717  nit = 0;
1718  while (dtr >= 0 && dtr < (short)cls[ipl].size() && nit < cls[ipl].size()) {
1719  // determine which end of the dtr should be attached to mom
1720  for (dtrBrkEnd = 0; dtrBrkEnd < 2; ++dtrBrkEnd)
1721  if (cls[ipl][dtr].BrkIndex[dtrBrkEnd] == mom) break;
1722  if (dtrBrkEnd == 2) {
1723  gotcl[icl] = false;
1724  break;
1725  }
1726  // check for reasonable merge error
1727  if (cls[ipl][dtr].MergeError[dtrBrkEnd] < fMergeErrorCut) {
1728  sCluster.push_back(dtr);
1729  sOrder.push_back(dtrBrkEnd);
1730  gotcl[dtr] = true;
1731  }
1732  ++nit;
1733  // set up to check the new mom
1734  mom = dtr;
1735  // at the other end
1736  momBrkEnd = 1 - dtrBrkEnd;
1737  // with the new dtr
1738  dtr = cls[ipl][mom].BrkIndex[momBrkEnd];
1739  } // dtr >= 0 ...
1740  if (dtrBrkEnd == 2) continue;
1741  } // end
1742 
1743  if (!gotcl[icl]) {
1744  // a single unbroken cluster
1745  sCluster.push_back(icl);
1746  sOrder.push_back(0);
1747  }
1748 
1749  if (sCluster.size() == 0) {
1750  mf::LogError("CCTM") << "MakeClusterChains error in plane " << ipl << " cluster " << icl;
1751  return;
1752  }
1753 
1754  ClsChainPar ccp;
1755  // fill the struct parameters assuming this is a cluster chain containing only one cluster
1756  unsigned short jcl = sCluster[0];
1757  if (jcl > cls[ipl].size()) std::cout << "oops MCC\n";
1758  unsigned short oend;
1759  for (end = 0; end < 2; ++end) {
1760  oend = end;
1761  if (sOrder[0] > 0) oend = 1 - end;
1762  ccp.Wire[end] = cls[ipl][jcl].Wire[oend];
1763  ccp.Time[end] = cls[ipl][jcl].Time[oend];
1764  ccp.X[end] = cls[ipl][jcl].X[oend];
1765  ccp.Slope[end] = cls[ipl][jcl].Slope[oend];
1766  ccp.Angle[end] = cls[ipl][jcl].Angle[oend];
1767  ccp.Dir[end] = cls[ipl][icl].Dir[oend];
1768  ccp.VtxIndex[end] = cls[ipl][jcl].VtxIndex[oend];
1769  ccp.ChgNear[end] = cls[ipl][jcl].ChgNear[oend];
1770  ccp.mBrkIndex[end] = cls[ipl][jcl].BrkIndex[oend];
1771  } // end
1772  ccp.Length = cls[ipl][icl].Length;
1773  ccp.TotChg = cls[ipl][icl].TotChg;
1774  ccp.InTrack = -1;
1775 
1776  for (unsigned short ii = 1; ii < sCluster.size(); ++ii) {
1777  jcl = sCluster[ii];
1778  if (jcl > cls[ipl].size()) std::cout << "oops MCC\n";
1779  // end is the end where the break is being mended
1780  end = sOrder[ii];
1781  if (end > 1) std::cout << "oops2 MCC\n";
1782  oend = 1 - end;
1783  // update the parameters at the other end of the chain
1784  ccp.Wire[1] = cls[ipl][jcl].Wire[oend];
1785  ccp.Time[1] = cls[ipl][jcl].Time[oend];
1786  ccp.X[1] = cls[ipl][jcl].X[oend];
1787  ccp.Slope[1] = cls[ipl][jcl].Slope[oend];
1788  ccp.Angle[1] = cls[ipl][jcl].Angle[oend];
1789  ccp.Dir[1] = cls[ipl][jcl].Dir[oend];
1790  ccp.VtxIndex[1] = cls[ipl][jcl].VtxIndex[oend];
1791  ccp.ChgNear[1] = cls[ipl][jcl].ChgNear[oend];
1792  ccp.mBrkIndex[1] = cls[ipl][jcl].BrkIndex[oend];
1793  ccp.Length += cls[ipl][jcl].Length;
1794  ccp.TotChg += cls[ipl][jcl].TotChg;
1795  } // ii
1796  ccp.ClsIndex = sCluster;
1797  ccp.Order = sOrder;
1798  // redo the direction
1799  if (ccp.Time[1] > ccp.Time[0]) {
1800  ccp.Dir[0] = 1;
1801  ccp.Dir[1] = -1;
1802  }
1803  else {
1804  ccp.Dir[0] = -1;
1805  ccp.Dir[1] = 1;
1806  }
1807  clsChain[ipl].push_back(ccp);
1808 
1809  } // icl
1810 
1811  // re-index mBrkIndex to point to a cluster chain, not a cluster
1812  unsigned short brkCls;
1813  bool gotit;
1814  for (unsigned short ccl = 0; ccl < clsChain[ipl].size(); ++ccl) {
1815  for (unsigned short end = 0; end < 2; ++end) {
1816  if (clsChain[ipl][ccl].mBrkIndex[end] < 0) continue;
1817  brkCls = clsChain[ipl][ccl].mBrkIndex[end];
1818  gotit = false;
1819  // find this cluster index in a cluster chain
1820  for (unsigned short ccl2 = 0; ccl2 < clsChain[ipl].size(); ++ccl2) {
1821  if (ccl2 == ccl) continue;
1822  if (std::find(clsChain[ipl][ccl2].ClsIndex.begin(),
1823  clsChain[ipl][ccl2].ClsIndex.end(),
1824  brkCls) == clsChain[ipl][ccl2].ClsIndex.end())
1825  continue;
1826  // found it
1827  clsChain[ipl][ccl].mBrkIndex[end] = ccl2;
1828  gotit = true;
1829  break;
1830  } // ccl2
1831  if (!gotit)
1832  mf::LogError("CCTM") << "MCC: Cluster chain " << ccl << " end " << end
1833  << " Failed to find brkCls " << brkCls << " in plane " << ipl;
1834  } // end
1835  } // ccl
1836 
1837  } // ipl
1838 
1839  if (gotprt) PrintClusters(detProp, tpcid);
1840  prt = false;
1841 
1842  } // MakeClusterChains
1843 
1846  unsigned short ipl,
1847  unsigned short icl1,
1848  unsigned short end1)
1849  {
1850  // project cluster icl1 at end1 to find the best intersection with icl2
1851  float dw, dx, best = 999;
1852  std::vector<art::Ptr<recob::Hit>> clusterhits = fmCluHits.at(cls[ipl][icl1].EvtIndex);
1853  for (unsigned short hit = 0; hit < clusterhits.size(); ++hit) {
1854  dw = clusterhits[hit]->WireID().Wire - cls[ipl][icl1].Wire[end1];
1855  dx = fabs(cls[ipl][icl1].Time[end1] + dw * fWirePitch * cls[ipl][icl1].Slope[end1] -
1856  clusterhits[hit]->PeakTime());
1857  if (dx < best) best = dx;
1858  if (dx < 0.01) break;
1859  } // hit
1860  return best;
1861  } // dXClTraj
1862 
1865  art::FindManyP<recob::Hit> const& fmCluHits,
1866  unsigned short imat,
1867  unsigned short procCode,
1868  geo::TPCID const& tpcid)
1869  {
1870  // store the current "under construction" track in the trk vector
1871 
1872  TrkPar newtrk;
1873 
1874  if (imat > matcomb.size() - 1) {
1875  mf::LogError("CCTM") << "Bad imat in StoreTrack";
1876  return;
1877  }
1878 
1879  // ensure there are at least 2 hits in at least 2 planes
1880  unsigned short nhitinpl = 0;
1881  auto const nplanes = geom->Nplanes(tpcid);
1882  for (unsigned short ipl = 0; ipl < nplanes; ++ipl)
1883  if (trkHits[ipl].size() > 1) ++nhitinpl;
1884  if (nhitinpl < 2) {
1885  mf::LogError("CCTM") << "StoreTrack: Not enough hits in each plane\n";
1886  return;
1887  }
1888  if (prt)
1889  mf::LogVerbatim("CCTM") << "In StoreTrack: matcomb " << imat << " cluster chains "
1890  << matcomb[imat].Cls[0] << " " << matcomb[imat].Cls[1] << " "
1891  << matcomb[imat].Cls[2];
1892 
1893  // Track hit vectors for fitting the trajectory
1894  std::array<std::vector<geo::WireID>, 3> trkWID;
1895  std::array<std::vector<double>, 3> trkX;
1896  std::array<std::vector<double>, 3> trkXErr;
1897 
1898  // track trajectory for a track
1899  std::vector<TVector3> trkPos;
1900  std::vector<TVector3> trkDir;
1901 
1902  newtrk.ID = trk.size() + 1;
1903  newtrk.Proc = procCode;
1904  newtrk.TrkHits = trkHits;
1905  // BUG the double brace syntax is required to work around clang bug 21629
1906  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
1907  newtrk.VtxIndex = {{-1, -1}};
1908  newtrk.ChgOrder = 0;
1909  newtrk.MomID = -1;
1910  // BUG the double brace syntax is required to work around clang bug 21629
1911  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
1912  newtrk.EndInTPC = {{false, false}};
1913  newtrk.GoodEnd = {{false, false}};
1914  newtrk.DtrID = {0};
1915  newtrk.PDGCode = -1;
1916 
1917  unsigned short icl, iht;
1918 
1919  if (prt)
1920  mf::LogVerbatim("CCTM") << "CCTM: Make traj for track " << newtrk.ID << " procCode "
1921  << procCode << " nhits in planes " << trkHits[0].size() << " "
1922  << trkHits[1].size() << " " << trkHits[2].size();
1923  // make the track trajectory
1924  if (nplanes == 2) {
1925  trkWID[2].resize(0);
1926  trkX[2].resize(0);
1927  trkXErr[2].resize(0);
1928  }
1929  for (auto const& planeid : geom->Iterate<geo::PlaneID>(tpcid)) {
1930  auto const ipl = planeid.Plane;
1931  trkWID[ipl].resize(trkHits[ipl].size());
1932  trkX[ipl].resize(trkHits[ipl].size());
1933  trkXErr[ipl].resize(trkHits[ipl].size());
1934  for (iht = 0; iht < trkHits[ipl].size(); ++iht) {
1935  trkWID[ipl][iht] = trkHits[ipl][iht]->WireID();
1936  trkX[ipl][iht] = detProp.ConvertTicksToX(trkHits[ipl][iht]->PeakTime(), planeid);
1937  trkXErr[ipl][iht] =
1938  fHitFitErrFac * trkHits[ipl][iht]->RMS() * trkHits[ipl][iht]->Multiplicity();
1939  } // iht
1940  } // ipl
1941  fTrackTrajectoryAlg.TrackTrajectory(trkWID, trkX, trkXErr, trkPos, trkDir);
1942  if (trkPos.size() < 2) {
1943  mf::LogError("CCTM") << "StoreTrack: No trajectory points on failed track " << newtrk.ID
1944  << " in StoreTrack: matcomb " << imat << " cluster chains "
1945  << matcomb[imat].Cls[0] << " " << matcomb[imat].Cls[1] << " "
1946  << matcomb[imat].Cls[2];
1947  // make a garbage trajectory
1948  trkPos.resize(2);
1949  trkPos[1](2) = 1;
1950  trkDir.resize(2);
1951  trkDir[1](2) = 1;
1952  }
1953  newtrk.TrjPos = trkPos;
1954  newtrk.TrjDir = trkDir;
1955 
1956  if (prt) mf::LogVerbatim("CCTM") << " number of traj points " << trkPos.size();
1957 
1958  // determine if each end is good in the sense that there are hits in each plane
1959  // that are consistent in time and are presumed to form a good 3D space point
1960  unsigned short end, nClose, indx, jndx;
1961  float xErr;
1962  for (end = 0; end < 2; ++end) {
1963  nClose = 0;
1964  for (unsigned short ipl = 0; ipl < nplanes - 1; ++ipl) {
1965  if (trkX[ipl].size() == 0) continue;
1966  for (unsigned short jpl = ipl + 1; jpl < nplanes; ++jpl) {
1967  if (trkX[jpl].size() == 0) continue;
1968  if (end == 0) {
1969  indx = 0;
1970  jndx = 0;
1971  }
1972  else {
1973  indx = trkXErr[ipl].size() - 1;
1974  jndx = trkXErr[jpl].size() - 1;
1975  }
1976  xErr = 3 * (trkXErr[ipl][indx] + trkXErr[jpl][jndx]);
1977  if (std::abs(trkX[ipl][indx] - trkX[jpl][jndx]) < xErr) ++nClose;
1978  } // jpl
1979  } // ipl
1980  if (nClose == nplanes) newtrk.GoodEnd[end] = true;
1981  } // end
1982 
1983  // set trajectory end points to a vertex if one exists
1984  unsigned short ivx, itj, ccl;
1985  float dx, dy, dz, dr0, dr1;
1986  unsigned short attachEnd;
1987  for (end = 0; end < 2; ++end) {
1988  ivx = USHRT_MAX;
1989  if (end == 0 && matcomb[imat].Vtx >= 0) ivx = matcomb[imat].Vtx;
1990  if (end == 1 && matcomb[imat].oVtx >= 0) ivx = matcomb[imat].oVtx;
1991  if (ivx == USHRT_MAX) continue;
1992  // determine the proper end using the TrjPos order and brute force
1993  itj = 0;
1994  dx = vtx[ivx].X - newtrk.TrjPos[itj](0);
1995  dy = vtx[ivx].Y - newtrk.TrjPos[itj](1);
1996  dz = vtx[ivx].Z - newtrk.TrjPos[itj](2);
1997  dr0 = dx * dx + dy * dy + dz * dz;
1998  itj = newtrk.TrjPos.size() - 1;
1999  dx = vtx[ivx].X - newtrk.TrjPos[itj](0);
2000  dy = vtx[ivx].Y - newtrk.TrjPos[itj](1);
2001  dz = vtx[ivx].Z - newtrk.TrjPos[itj](2);
2002  dr1 = dx * dx + dy * dy + dz * dz;
2003  attachEnd = 1;
2004  if (dr0 < dr1) {
2005  itj = 0;
2006  attachEnd = 0;
2007  // a really bad match to the vertex
2008  if (dr0 > 5) return;
2009  }
2010  else {
2011  // a really bad match to the vertex
2012  if (dr1 > 5) return;
2013  }
2014  newtrk.TrjPos[itj](0) = vtx[ivx].X;
2015  newtrk.TrjPos[itj](1) = vtx[ivx].Y;
2016  newtrk.TrjPos[itj](2) = vtx[ivx].Z;
2017  newtrk.VtxIndex[attachEnd] = ivx;
2018  // correct the trajectory direction
2019  TVector3 dir;
2020  if (itj == 0) {
2021  dir = newtrk.TrjPos[1] - newtrk.TrjPos[0];
2022  newtrk.TrjDir[0] = dir.Unit();
2023  }
2024  else {
2025  dir = newtrk.TrjPos[itj - 1] - newtrk.TrjPos[itj];
2026  newtrk.TrjDir[itj] = dir.Unit();
2027  }
2028  } // end
2029 
2030  if (newtrk.VtxIndex[0] >= 0 && newtrk.VtxIndex[0] == newtrk.VtxIndex[1]) {
2031  mf::LogError("CCTM")
2032  << "StoreTrack: Trying to attach a vertex to both ends of a track. imat = " << imat;
2033  return;
2034  }
2035 
2036  // calculate the length
2037  newtrk.Length = 0;
2038  float norm;
2039  double X, Y, Z;
2040  for (unsigned short itj = 1; itj < newtrk.TrjPos.size(); ++itj) {
2041  X = newtrk.TrjPos[itj](0) - newtrk.TrjPos[itj - 1](0);
2042  Y = newtrk.TrjPos[itj](1) - newtrk.TrjPos[itj - 1](1);
2043  Z = newtrk.TrjPos[itj](2) - newtrk.TrjPos[itj - 1](2);
2044  norm = sqrt(X * X + Y * Y + Z * Z);
2045  newtrk.Length += norm;
2046  }
2047 
2048  // store the cluster -> track assignment
2049  newtrk.ClsEvtIndices.clear();
2050  for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
2051  if (matcomb[imat].Cls[ipl] < 0) continue;
2052  ccl = matcomb[imat].Cls[ipl];
2053  if (ccl > clsChain[ipl].size()) std::cout << "oops StoreTrack\n";
2054  clsChain[ipl][ccl].InTrack = newtrk.ID;
2055  for (unsigned short icc = 0; icc < clsChain[ipl][ccl].ClsIndex.size(); ++icc) {
2056  icl = clsChain[ipl][ccl].ClsIndex[icc];
2057  if (icl > cls[ipl].size()) std::cout << "oops StoreTrack\n";
2058  cls[ipl][icl].InTrack = newtrk.ID;
2059  if (cls[ipl][icl].EvtIndex > fmCluHits.size() - 1) {
2060  std::cout << "ooops2 store track EvtIndex " << cls[ipl][icl].EvtIndex << " size "
2061  << fmCluHits.size() << " icl " << icl << "\n";
2062  continue;
2063  }
2064  newtrk.ClsEvtIndices.push_back(cls[ipl][icl].EvtIndex);
2065  } // icc
2066  } // ipl
2067 
2068  if (prt) mf::LogVerbatim("CCTM") << " track ID " << newtrk.ID << " stored in StoreTrack";
2069 
2070  trk.push_back(newtrk);
2071  } // StoreTrack
2072 
2075  geo::TPCID const& tpcid)
2076  {
2077  // Match clusters in all planes
2078  bool ignoreSign;
2079  float kSlp, kAng, kX, kWir, okWir;
2080  short idir, ioend, jdir, joend, kdir;
2081 
2082  double yp, zp;
2083  float tpcSizeY = geom->DetHalfWidth();
2084  float tpcSizeZ = geom->DetLength();
2085 
2086  float dxcut = 2;
2087  float dxkcut;
2088  float dwcut = 6;
2089  if (fuBCode) {
2090  dxcut = 20;
2091  dwcut = 60;
2092  }
2093 
2094  // temp array for making a rough charge asymmetry cut
2095  std::array<float, 3> mchg;
2096  auto const nplanes = geom->Nplanes(tpcid);
2097  for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
2098  geo::PlaneID const iplane_id{tpcid, ipl};
2099  for (unsigned short icl = 0; icl < clsChain[ipl].size(); ++icl) {
2100  if (clsChain[ipl][icl].InTrack >= 0) continue;
2101  // skip short clusters
2102  if (clsChain[ipl][icl].Length < fMatchMinLen[algIndex]) continue;
2103  unsigned short jpl = (ipl + 1) % nplanes;
2104  unsigned short kpl = (jpl + 1) % nplanes;
2105  geo::PlaneID const jplane_id{tpcid, jpl};
2106  for (unsigned short jcl = 0; jcl < clsChain[jpl].size(); ++jcl) {
2107  if (clsChain[jpl][jcl].InTrack >= 0) continue;
2108  // skip short clusters
2109  if (clsChain[jpl][jcl].Length < fMatchMinLen[algIndex]) continue;
2110  // make first charge asymmetry cut
2111  mchg[0] = clsChain[ipl][icl].TotChg;
2112  mchg[1] = clsChain[jpl][jcl].TotChg;
2113  mchg[2] = mchg[1];
2114  if (fChgAsymFactor[algIndex] > 0 && ChargeAsym(mchg) > 0.5) continue;
2115  for (unsigned short iend = 0; iend < 2; ++iend) {
2116  idir = clsChain[ipl][icl].Dir[iend];
2117  for (unsigned short jend = 0; jend < 2; ++jend) {
2118  jdir = clsChain[jpl][jcl].Dir[jend];
2119  if (idir != 0 && jdir != 0 && idir != jdir) continue;
2120  // make an X cut
2121  if (fabs(clsChain[ipl][icl].X[iend] - clsChain[jpl][jcl].X[jend]) > dxcut) continue;
2122  ioend = 1 - iend;
2123  joend = 1 - jend;
2124  // Find the expected third (k) plane parameters
2125  kSlp = geom->ThirdPlaneSlope(
2126  ipl, clsChain[ipl][icl].Slope[iend], jpl, clsChain[jpl][jcl].Slope[jend], tpcid);
2127  kAng = atan(kSlp);
2128  // Ensure the match end is within the TPC
2130  geo::WireID{iplane_id, (unsigned int)(0.5 + clsChain[ipl][icl].Wire[iend])},
2131  geo::WireID{jplane_id, (unsigned int)(0.5 + clsChain[jpl][jcl].Wire[jend])},
2132  yp,
2133  zp);
2134  if (yp > tpcSizeY || yp < -tpcSizeY) continue;
2135  if (zp < 0 || zp > tpcSizeZ) continue;
2136  kX = 0.5 * (clsChain[ipl][icl].X[iend] + clsChain[jpl][jcl].X[jend]);
2137  kWir = geom->WireCoordinate(geo::Point_t{0, yp, zp}, geo::PlaneID{tpcid, kpl});
2138  // now look at the other end
2140  geo::WireID{iplane_id, (unsigned int)(0.5 + clsChain[ipl][icl].Wire[ioend])},
2141  geo::WireID{jplane_id, (unsigned int)(0.5 + clsChain[jpl][jcl].Wire[joend])},
2142  yp,
2143  zp);
2144  if (yp > tpcSizeY || yp < -tpcSizeY) continue;
2145  if (zp < 0 || zp > tpcSizeZ) continue;
2146  okWir = geom->WireCoordinate(geo::Point_t{0, yp, zp}, geo::PlaneID{tpcid, kpl});
2147  if (prt)
2148  mf::LogVerbatim("CCTM")
2149  << "PlnMatch: chk i " << ipl << ":" << icl << ":" << iend << " idir " << idir
2150  << " X " << clsChain[ipl][icl].X[iend] << " j " << jpl << ":" << jcl << ":"
2151  << jend << " jdir " << jdir << " X " << clsChain[jpl][jcl].X[jend];
2152 
2153  if (prt)
2154  mf::LogVerbatim("CCTM")
2155  << "PlnMatch: chk j " << ipl << ":" << icl << ":" << iend << " " << jpl << ":"
2156  << jcl << ":" << jend << " iSlp " << std::setprecision(2)
2157  << clsChain[ipl][icl].Slope[iend] << " jSlp " << std::setprecision(2)
2158  << clsChain[jpl][jcl].Slope[jend] << " kWir " << (int)kWir << " okWir "
2159  << (int)okWir << " kSlp " << std::setprecision(2) << kSlp << " kAng "
2160  << std::setprecision(2) << kAng << " kX " << std::setprecision(1) << kX;
2161 
2162  // handle the case near pi/2, where the errors on large slopes
2163  // could result in a wrong-sign kAng
2164  ignoreSign = (fabs(kSlp) > 1.5);
2165  if (ignoreSign) kAng = fabs(kAng);
2166  dxkcut = dxcut * AngleFactor(kSlp);
2167  bool gotkcl = false;
2168  for (unsigned short kcl = 0; kcl < clsChain[kpl].size(); ++kcl) {
2169  if (clsChain[kpl][kcl].InTrack >= 0) continue;
2170  // make second charge asymmetry cut
2171  mchg[0] = clsChain[ipl][icl].TotChg;
2172  mchg[1] = clsChain[jpl][jcl].TotChg;
2173  mchg[2] = clsChain[kpl][kcl].TotChg;
2174  if (fChgAsymFactor[algIndex] > 0 && ChargeAsym(mchg) > 0.5) continue;
2175  for (unsigned short kend = 0; kend < 2; ++kend) {
2176  kdir = clsChain[kpl][kcl].Dir[kend];
2177  if (idir != 0 && kdir != 0 && idir != kdir) continue;
2178  if (prt)
2179  mf::LogVerbatim("CCTM")
2180  << " kcl " << kcl << " kend " << kend << " dx "
2181  << std::abs(clsChain[kpl][kcl].X[kend] - kX) << " dxkcut " << dxkcut;
2182  if (std::abs(clsChain[kpl][kcl].X[kend] - kX) > dxkcut) continue;
2183  // rough dWire cut
2184  if (prt)
2185  mf::LogVerbatim("CCTM")
2186  << " kcl " << kcl << " kend " << kend << " dw "
2187  << (clsChain[kpl][kcl].Wire[kend] - kWir) << " ignoreSign " << ignoreSign;
2188  if (fabs(clsChain[kpl][kcl].Wire[kend] - kWir) > dwcut) continue;
2189  if (prt) mf::LogVerbatim("CCTM") << " chk k " << kpl << ":" << kcl << ":" << kend;
2190  MatchPars match;
2191  match.Cls[ipl] = icl;
2192  match.End[ipl] = iend;
2193  match.Cls[jpl] = jcl;
2194  match.End[jpl] = jend;
2195  match.Cls[kpl] = kcl;
2196  match.End[kpl] = kend;
2197  match.Err = 100;
2198  if (DupMatch(match, nplanes)) continue;
2199  match.Chg[ipl] = 0;
2200  match.Chg[jpl] = 0;
2201  match.Chg[kpl] = 0;
2202  match.Vtx = clsChain[ipl][icl].VtxIndex[iend];
2203  match.oVtx = -1;
2204  FillEndMatch(detProp, tpcid, match);
2205  if (prt)
2206  mf::LogVerbatim("CCTM") << " PlnMatch: match k " << kpl << ":" << match.Cls[kpl]
2207  << ":" << match.End[kpl] << " oChg " << match.Chg[kpl]
2208  << " mErr " << match.Err << " oErr " << match.oErr;
2209  if (match.Chg[kpl] == 0) continue;
2210  if (match.Err > 10 || match.oErr > 10) continue;
2211  if (prt) mf::LogVerbatim("CCTM") << " dup? ";
2212  if (DupMatch(match, nplanes)) continue;
2213  matcomb.push_back(match);
2214  gotkcl = true;
2215  } // kend
2216  } // kcl
2217  if (prt) mf::LogVerbatim("CCTM") << " PlnMatch: gotkcl " << gotkcl;
2218  if (!gotkcl) {
2219  // make a 2-plane match and try again
2220  MatchPars match;
2221  match.Cls[ipl] = icl;
2222  match.End[ipl] = iend;
2223  match.Cls[jpl] = jcl;
2224  match.End[jpl] = jend;
2225  match.Cls[kpl] = -1;
2226  match.End[kpl] = 0;
2227  match.Err = 100;
2228  if (DupMatch(match, nplanes)) continue;
2229  match.Chg[ipl] = 0;
2230  match.Chg[jpl] = 0;
2231  match.Chg[kpl] = 0;
2232  match.Vtx = clsChain[ipl][icl].VtxIndex[iend];
2233  match.oVtx = -1;
2234  FillEndMatch(detProp, tpcid, match);
2235  if (prt)
2236  mf::LogVerbatim("CCTM")
2237  << " Tried 2-plane match"
2238  << " k " << kpl << ":" << match.Cls[kpl] << ":" << match.End[kpl] << " Chg "
2239  << match.Chg[kpl] << " Err " << match.Err << " match.oErr " << match.oErr;
2240  if (match.Chg[kpl] <= 0) continue;
2241  if (match.Err > 10 || match.oErr > 10) continue;
2242  matcomb.push_back(match);
2243  } // !gotkcl
2244  } // jend
2245  } // iend
2246  } // jcl
2247  } // icl
2248  } // ipl
2249  } // PlnMatch
2250 
2252  bool CCTrackMaker::DupMatch(MatchPars& match, unsigned short const nplanes)
2253  {
2254  unsigned short nMatCl, nMiss;
2255  float toterr = match.Err + match.oErr;
2256  for (unsigned int imat = 0; imat < matcomb.size(); ++imat) {
2257  // check for exact matches
2258  if (match.Cls[0] == matcomb[imat].Cls[0] && match.Cls[1] == matcomb[imat].Cls[1] &&
2259  match.Cls[2] == matcomb[imat].Cls[2]) {
2260 
2261  // compare the error
2262  if (toterr < matcomb[imat].Err + matcomb[imat].oErr) {
2263  // keep the better one
2264  matcomb[imat].End[0] = match.End[0];
2265  matcomb[imat].End[1] = match.End[1];
2266  matcomb[imat].End[2] = match.End[2];
2267  matcomb[imat].Vtx = match.Vtx;
2268  matcomb[imat].dWir = match.dWir;
2269  matcomb[imat].dAng = match.dAng;
2270  matcomb[imat].dX = match.dX;
2271  matcomb[imat].Err = match.Err;
2272  matcomb[imat].oVtx = match.oVtx;
2273  matcomb[imat].odWir = match.odWir;
2274  matcomb[imat].odAng = match.odAng;
2275  matcomb[imat].odX = match.odX;
2276  matcomb[imat].oErr = match.oErr;
2277  }
2278  return true;
2279  } // test
2280  // check for a 3-plane match vs 2-plane match
2281  nMatCl = 0;
2282  nMiss = 0;
2283  for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
2284  if (match.Cls[ipl] >= 0) {
2285  if (match.Cls[ipl] == matcomb[imat].Cls[ipl] &&
2286  (match.End[0] == matcomb[imat].End[0] || match.End[1] == matcomb[imat].End[1]))
2287  ++nMatCl;
2288  }
2289  else {
2290  ++nMiss;
2291  }
2292  } // ipl
2293  if (nMatCl == 2 && nMiss == 1) return true;
2294  } // imat
2295  return false;
2296  } // DupMatch
2297 
2300  art::FindManyP<recob::Hit> const& fmCluHits,
2301  unsigned short const procCode,
2302  geo::TPCID const& tpcid)
2303  {
2304  // sort cluster matches by increasing total match error. Find the minimum total error of all
2305  // cluster match combinations and make tracks from them
2306  CluLen merr;
2307  std::vector<CluLen> materr;
2308  unsigned int ii, im;
2309 
2310  if (matcomb.size() == 0) return;
2311 
2312  // sort by decreasing error
2313  for (ii = 0; ii < matcomb.size(); ++ii) {
2314  merr.index = ii;
2315  merr.length = matcomb[ii].Err + matcomb[ii].oErr;
2316  materr.push_back(merr);
2317  } // ii
2318  std::sort(materr.begin(), materr.end(), lessThan);
2319 
2320  auto const nplanes = geom->Nplanes(tpcid);
2321  if (prt) {
2322  mf::LogVerbatim myprt("CCTM");
2323  myprt << "SortMatches\n";
2324  myprt << " ii im Vx Err dW dA dX oVx oErr odW odA odX Asym "
2325  " icl jcl kcl \n";
2326  for (ii = 0; ii < materr.size(); ++ii) {
2327  im = materr[ii].index;
2328  float asym =
2329  fabs(matcomb[im].Chg[0] - matcomb[im].Chg[1]) / (matcomb[im].Chg[0] + matcomb[im].Chg[1]);
2330  asym *=
2331  fabs(matcomb[im].Chg[1] - matcomb[im].Chg[2]) / (matcomb[im].Chg[1] + matcomb[im].Chg[2]);
2332  myprt << std::fixed << std::right << std::setw(5) << ii << std::setw(5) << im
2333  << std::setw(4) << matcomb[im].Vtx << std::setw(7) << std::setprecision(2)
2334  << matcomb[im].Err << std::setw(7) << std::setprecision(1) << matcomb[im].dWir
2335  << std::setw(7) << std::setprecision(2) << matcomb[im].dAng << std::setw(7)
2336  << std::setprecision(2) << matcomb[im].dX << std::setw(4) << matcomb[im].oVtx
2337  << std::setw(7) << std::setprecision(2) << matcomb[im].oErr << std::setw(7)
2338  << std::setprecision(1) << matcomb[im].odWir << std::setw(7) << std::setprecision(2)
2339  << matcomb[im].odAng << std::setw(7) << std::setprecision(2) << matcomb[im].odX
2340  << std::setw(7) << std::setprecision(3) << asym << " 0:" << matcomb[im].Cls[0] << ":"
2341  << matcomb[im].End[0] << " 1:" << matcomb[im].Cls[1] << ":" << matcomb[im].End[1];
2342  if (nplanes > 2) myprt << " 2:" << matcomb[im].Cls[2] << ":" << matcomb[im].End[2];
2343  myprt << "\n";
2344  } // ii
2345  } // prt
2346 
2347  // define an array to ensure clusters are only used once
2348  std::array<std::vector<bool>, 3> pclUsed;
2349  unsigned short ipl;
2350  for (ipl = 0; ipl < nplanes; ++ipl) {
2351  pclUsed[ipl].resize(clsChain[ipl].size());
2352  }
2353 
2354  // count the total number of clusters and length used in matcomb
2355  unsigned short matcombTotCl = 0;
2356  float matcombTotLen = 0;
2357  unsigned short icl;
2358  for (ii = 0; ii < matcomb.size(); ++ii) {
2359  for (ipl = 0; ipl < nplanes; ++ipl) {
2360  if (matcomb[ii].Cls[ipl] < 0) continue;
2361  icl = matcomb[ii].Cls[ipl];
2362  ++matcombTotCl;
2363  matcombTotLen += clsChain[ipl][icl].Length;
2364  }
2365  }
2366 
2367  if (prt)
2368  mf::LogVerbatim("CCTM") << "Number of clusters to match " << matcombTotCl << " total length "
2369  << matcombTotLen;
2370 
2371  if (matcombTotLen <= 0) {
2372  mf::LogError("CCTM") << "SortMatches: bad matcomb total length " << matcombTotLen;
2373  return;
2374  }
2375 
2376  // vector of matcomb indices of unique cluster matches
2377  std::vector<unsigned short> matIndex;
2378  // vector of matcomb indices of unique cluster matches that have the best total error
2379  std::vector<unsigned short> bestMatIndex;
2380  float totLen, totErr, bestTotErr = 9999;
2381  // start with the best match
2382  unsigned short jj, jm, nused, jcl;
2383  // fraction of the length of all clustters in matcomb that are used in a match
2384  float fracLen;
2385 
2386  for (ii = 0; ii < materr.size(); ++ii) {
2387  im = materr[ii].index;
2388  matIndex.clear();
2389  // skip really bad matches
2390  if (matcomb[im].Err > bestTotErr) continue;
2391  totLen = 0;
2392  // initialize pclUsed and flag the clusters in this match
2393  for (ipl = 0; ipl < nplanes; ++ipl) {
2394  // initialize to no clusters used
2395  std::fill(pclUsed[ipl].begin(), pclUsed[ipl].end(), false);
2396  // check for 2 plane match
2397  if (matcomb[im].Cls[ipl] < 0) continue;
2398  icl = matcomb[im].Cls[ipl];
2399  pclUsed[ipl][icl] = true;
2400  totLen += clsChain[ipl][icl].Length;
2401  } // ipl
2402  // Initialize the error sum
2403  totErr = matcomb[im].Err;
2404  // Save the index
2405  matIndex.push_back(im);
2406  // look for matches in the rest of the list that are not already matched.
2407  for (jj = 0; jj < materr.size(); ++jj) {
2408  if (jj == ii) continue;
2409  jm = materr[jj].index;
2410  // skip really bad matches
2411  if (matcomb[jm].Err > bestTotErr) continue;
2412  // check for non-unique cluster indices
2413  nused = 0;
2414  for (ipl = 0; ipl < nplanes; ++ipl) {
2415  if (matcomb[jm].Cls[ipl] < 0) continue;
2416  jcl = matcomb[jm].Cls[ipl];
2417  if (pclUsed[ipl][jcl]) ++nused;
2418  // This cluster chain was used in a previous match
2419  if (nused > 0) break;
2420  totLen += clsChain[ipl][jcl].Length;
2421  } // ipl
2422  // at least one of the clusters in this match have been used
2423  if (nused != 0) continue;
2424  // found a match with an unmatched set of clusters. Update the total error and flag them
2425  totErr += matcomb[jm].Err;
2426  matIndex.push_back(jm);
2427  // Flag the clusters used and see if all of them are used
2428  for (ipl = 0; ipl < nplanes; ++ipl) {
2429  if (matcomb[jm].Cls[ipl] < 0) continue;
2430  jcl = matcomb[jm].Cls[ipl];
2431  pclUsed[ipl][jcl] = true;
2432  } // ipl
2433  } // jm
2434  if (totLen == 0) continue;
2435  nused = 0;
2436  for (ipl = 0; ipl < nplanes; ++ipl) {
2437  for (unsigned short indx = 0; indx < pclUsed[ipl].size(); ++indx)
2438  if (pclUsed[ipl][indx]) ++nused;
2439  } // ipl
2440  if (totLen > matcombTotLen) std::cout << "Oops " << totLen << " " << matcombTotLen << "\n";
2441  // weight the total error by the total length of all clusters
2442  fracLen = totLen / matcombTotLen;
2443  totErr /= fracLen;
2444  if (prt) {
2445  mf::LogVerbatim myprt("CCTM");
2446  myprt << "match " << im << " totErr " << totErr << " nused " << nused << " fracLen "
2447  << fracLen << " totLen " << totLen << " mat: ";
2448  for (unsigned short indx = 0; indx < matIndex.size(); ++indx)
2449  myprt << " " << matIndex[indx];
2450  } // prt
2451  // check for more used clusters and a better total error
2452  if (totErr < bestTotErr) {
2453  bestTotErr = totErr;
2454  bestMatIndex = matIndex;
2455  if (nused == matcombTotCl) break;
2456  if (prt) {
2457  mf::LogVerbatim myprt("CCTM");
2458  myprt << "bestTotErr " << bestTotErr << " nused " << nused << " matcombTotCl "
2459  << matcombTotCl << " mat: ";
2460  for (unsigned short indx = 0; indx < bestMatIndex.size(); ++indx)
2461  myprt << " " << bestMatIndex[indx];
2462  } // prt
2463  // stop looking if we have found everything
2464  if (fracLen > 0.999) break;
2465  } // totErr < bestTotErr
2466  } // im
2467 
2468  if (bestTotErr > 9000) return;
2469 
2470  for (ii = 0; ii < bestMatIndex.size(); ++ii) {
2471  im = bestMatIndex[ii];
2472  FillTrkHits(fmCluHits, im, nplanes);
2473  // look for missing clusters?
2474  // store this track with processor code 1
2475  StoreTrack(detProp, fmCluHits, im, procCode, tpcid);
2476  } // ii
2477 
2478  } // SortMatches
2479 
2482  {
2483  // 2D version of FillEndMatch
2484 
2485  match.Err = 100;
2486  match.oErr = 100;
2487  match.Chg[2] = 0;
2488  match.dWir = 0;
2489  match.dAng = 0;
2490  match.odWir = 0;
2491  match.odAng = 0;
2492 
2493  unsigned short ipl = 0;
2494  unsigned short jpl = 1;
2495 
2496  if (match.Cls[0] < 0 || match.Cls[1] < 0) return;
2497 
2498  unsigned short icl = match.Cls[ipl];
2499  unsigned short iend = match.End[ipl];
2500  match.Chg[ipl] = clsChain[ipl][icl].TotChg;
2501  // cluster i match end
2502  float miX = clsChain[ipl][icl].X[iend];
2503  // cluster i other end
2504  unsigned short oiend = 1 - iend;
2505  float oiX = clsChain[ipl][icl].X[oiend];
2506 
2507  unsigned short jcl = match.Cls[jpl];
2508  unsigned short jend = match.End[jpl];
2509  match.Chg[jpl] = clsChain[jpl][jcl].TotChg;
2510  // cluster j match end
2511  float mjX = clsChain[jpl][jcl].X[jend];
2512  // cluster j other end
2513  unsigned short ojend = 1 - jend;
2514  float ojX = clsChain[jpl][jcl].X[ojend];
2515 
2516  // look for a match end vertex match
2517  match.Vtx = -1;
2518  if (clsChain[ipl][icl].VtxIndex[iend] >= 0 &&
2519  clsChain[ipl][icl].VtxIndex[iend] == clsChain[jpl][jcl].VtxIndex[jend]) {
2520  match.Vtx = clsChain[ipl][icl].VtxIndex[iend];
2521  miX = vtx[match.Vtx].X;
2522  mjX = vtx[match.Vtx].X;
2523  }
2524 
2525  // look for an other end vertex match
2526  match.oVtx = -1;
2527  if (clsChain[ipl][icl].VtxIndex[oiend] >= 0 &&
2528  clsChain[ipl][icl].VtxIndex[oiend] == clsChain[jpl][jcl].VtxIndex[ojend]) {
2529  match.oVtx = clsChain[ipl][icl].VtxIndex[oiend];
2530  oiX = vtx[match.oVtx].X;
2531  ojX = vtx[match.oVtx].X;
2532  }
2533 
2534  // find the charge asymmetry
2535  float chgAsym = 1;
2536  if (fChgAsymFactor[algIndex] > 0) {
2537  chgAsym = fabs(match.Chg[ipl] - match.Chg[jpl]) / (match.Chg[ipl] + match.Chg[jpl]);
2538  if (chgAsym > 0.5) return;
2539  chgAsym = 1 + fChgAsymFactor[algIndex] * chgAsym;
2540  }
2541 
2542  // find the error at the match end
2543  float maxSlp = fabs(clsChain[ipl][icl].Slope[iend]);
2544  if (fabs(clsChain[jpl][jcl].Slope[jend]) > maxSlp)
2545  maxSlp = fabs(clsChain[jpl][jcl].Slope[jend]);
2546  float sigmaX = fXMatchErr[algIndex] + std::max(maxSlp, (float)20);
2547  match.dX = fabs(miX - mjX);
2548  match.Err = chgAsym * match.dX / sigmaX;
2549 
2550  // find the error at the other end
2551  maxSlp = fabs(clsChain[ipl][icl].Slope[oiend]);
2552  if (fabs(clsChain[jpl][jcl].Slope[ojend]) > maxSlp)
2553  maxSlp = fabs(clsChain[jpl][jcl].Slope[ojend]);
2554  sigmaX = fXMatchErr[algIndex] + std::max(maxSlp, (float)20);
2555  match.odX = fabs(oiX - ojX);
2556  match.oErr = chgAsym * match.odX / sigmaX;
2557 
2558  if (prt)
2559  mf::LogVerbatim("CCTM") << "FEM2: m " << ipl << ":" << icl << ":" << iend << " miX " << miX
2560  << " - " << jpl << ":" << jcl << ":" << jend << " mjX " << mjX
2561  << " match.dX " << match.dX << " match.Err " << match.Err
2562  << " chgAsym " << chgAsym << " o "
2563  << " oiX " << oiX << " ojX " << ojX << " match.odX " << match.odX
2564  << " match.oErr " << match.oErr << "\n";
2565 
2566  } // FillEndMatch2
2567 
2570  geo::TPCID const& tpcid,
2571  MatchPars& match)
2572  {
2573  // fill the matching parameters for this cluster match. The calling routine
2574  // should set the match end vertex ID (if applicable) as well as the
2575  // cluster IDs and matching ends in each plane. Note that the matching variables
2576  // Note that dWir, dAng and dTim are not filled if there is a vertex (match.Vtx >= 0).
2577  // Likewise, odWir, odAng and odX are not filled if there is a vertex match
2578  // at the other end
2579  auto const nplanes = geom->Nplanes(tpcid);
2580 
2581  if (nplanes == 2) {
2582  FillEndMatch2(match);
2583  return;
2584  }
2585 
2586  std::array<short, 3> mVtx;
2587  std::array<short, 3> oVtx;
2588  std::array<float, 3> oWir;
2589  std::array<float, 3> oSlp;
2590  std::array<float, 3> oAng;
2591  std::array<float, 3> oX;
2592 
2593  std::array<float, 3> mChg;
2594 
2595  unsigned short ii, ipl, iend, jpl, jend, kpl, kend, oend;
2596  short icl, jcl, kcl;
2597 
2598  for (ipl = 0; ipl < 3; ++ipl) {
2599  mVtx[ipl] = -1;
2600  oVtx[ipl] = -1;
2601  oWir[ipl] = -66;
2602  oSlp[ipl] = -66;
2603  oAng[ipl] = -66;
2604  oX[ipl] = -66;
2605  mChg[ipl] = -1;
2606  } // ipl
2607 
2608  // initialize parameters that shouldn't have been set by the calling routine
2609  match.dWir = 0;
2610  match.dAng = 0;
2611  match.dX = 0;
2612  match.Err = 100;
2613  match.odWir = 0;
2614  match.odAng = 0;
2615  match.odX = 0;
2616  match.oErr = 100;
2617  match.oVtx = -1;
2618 
2619  if (prt) {
2620  mf::LogVerbatim myprt("CCTM");
2621  myprt << "FEM ";
2622  for (ipl = 0; ipl < nplanes; ++ipl) {
2623  myprt << " " << ipl << ":" << match.Cls[ipl] << ":" << match.End[ipl];
2624  }
2625  }
2626 
2627  short missingPlane = -1;
2628  unsigned short nClInPln = 0;
2629  // number of vertex matches at each end
2630  short aVtx = -1;
2631  unsigned short novxmat = 0;
2632  short aoVtx = -1;
2633  unsigned short nvxmat = 0;
2634  unsigned short nShortCl = 0;
2635  // fill the other end parameters in each plane
2636  for (ipl = 0; ipl < nplanes; ++ipl) {
2637  if (match.Cls[ipl] < 0) {
2638  missingPlane = ipl;
2639  continue;
2640  }
2641  ++nClInPln;
2642  icl = match.Cls[ipl];
2643  match.Chg[ipl] = clsChain[ipl][icl].TotChg;
2644  mChg[ipl] = clsChain[ipl][icl].TotChg;
2645  iend = match.End[ipl];
2646  mVtx[ipl] = clsChain[ipl][icl].VtxIndex[iend];
2647  if (clsChain[ipl][icl].Length < 6) ++nShortCl;
2648  if (mVtx[ipl] >= 0) {
2649  if (aVtx < 0) aVtx = mVtx[ipl];
2650  if (mVtx[ipl] == aVtx) ++nvxmat;
2651  }
2652  if (prt)
2653  mf::LogVerbatim("CCTM") << "FEM: m " << ipl << ":" << icl << ":" << iend << " Vtx "
2654  << mVtx[ipl] << " Wir " << clsChain[ipl][icl].Wire[iend]
2655  << std::fixed << std::setprecision(3) << " Slp "
2656  << clsChain[ipl][icl].Slope[iend] << std::fixed
2657  << std::setprecision(1) << " X " << clsChain[ipl][icl].X[iend];
2658 
2659  oend = 1 - iend;
2660  oWir[ipl] = clsChain[ipl][icl].Wire[oend];
2661  oAng[ipl] = clsChain[ipl][icl].Angle[oend];
2662  oSlp[ipl] = clsChain[ipl][icl].Slope[oend];
2663  oX[ipl] = clsChain[ipl][icl].X[oend];
2664  oVtx[ipl] = clsChain[ipl][icl].VtxIndex[oend];
2665  if (oVtx[ipl] >= 0) {
2666  if (aoVtx < 0) aoVtx = oVtx[ipl];
2667  if (oVtx[ipl] == aoVtx) ++novxmat;
2668  }
2669 
2670  if (prt)
2671  mf::LogVerbatim("CCTM") << " o " << ipl << ":" << icl << ":" << oend << " oVtx "
2672  << oVtx[ipl] << " oWir " << oWir[ipl] << std::fixed
2673  << std::setprecision(3) << " oSlp " << oSlp[ipl] << std::fixed
2674  << std::setprecision(1) << " oX " << oX[ipl] << " Chg "
2675  << (int)mChg[ipl];
2676 
2677  } // ipl
2678 
2679  bool isShort = (nShortCl > 1);
2680 
2681  if (nClInPln < 2) {
2682  mf::LogWarning("CCTM") << "Not enough matched planes supplied";
2683  return;
2684  }
2685 
2686  if (prt)
2687  mf::LogVerbatim("CCTM") << "FEM: Vtx m " << aVtx << " count " << nvxmat << " o " << aoVtx
2688  << " count " << novxmat << " missingPlane " << missingPlane
2689  << " nClInPln " << nClInPln;
2690 
2691  // perfect match
2692  if (nvxmat == 3 && novxmat == 3) {
2693  match.Vtx = aVtx;
2694  match.Err = 0;
2695  match.oVtx = aoVtx;
2696  match.oErr = 0;
2697  return;
2698  }
2699 
2700  // 2-plane vertex match?
2701  // factors applied to error = 1 (no vtx), 0.5 (2 pln vtx), 0.33 (3 pln vtx)
2702  float vxFactor = 1;
2703  float ovxFactor = 1;
2704  if (nClInPln == 3) {
2705  // a cluster in all 3 planes
2706  if (nvxmat == 3) {
2707  // and all vertex assignments agree at the match end
2708  match.Vtx = aVtx;
2709  vxFactor = 0.33;
2710  }
2711  if (novxmat == 3) {
2712  // and all vertex assignments agree at the other end
2713  match.oVtx = aoVtx;
2714  ovxFactor = 0.33;
2715  }
2716  }
2717  else {
2718  // a cluster in 2 planes
2719  if (nvxmat == 2) {
2720  match.Vtx = aVtx;
2721  vxFactor = 0.5;
2722  }
2723  if (novxmat == 2) {
2724  match.oVtx = aoVtx;
2725  ovxFactor = 0.5;
2726  }
2727  } // nClInPln
2728 
2729  // find wire, X and Time at both ends
2730 
2731  // Find the "other end" matching parameters with
2732  // two cases: a 3-plane match or a 2-plane match
2733  // and with/without an other end vertex
2734 
2735  double ypos, zpos;
2736  float kWir, okWir;
2737  float kSlp, okSlp, kAng, okAng, okX, kX, kTim, okTim;
2738 
2739  if (nClInPln == 3) {
2740  ipl = 0;
2741  jpl = 1;
2742  kpl = 2;
2743  }
2744  else {
2745  // 2-plane match
2746  kpl = missingPlane;
2747  if (kpl == 0) {
2748  ipl = 1;
2749  jpl = 2;
2750  }
2751  else if (kpl == 1) {
2752  ipl = 2;
2753  jpl = 0;
2754  }
2755  else {
2756  ipl = 0;
2757  jpl = 1;
2758  } // kpl test
2759  } // missing plane
2760  iend = match.End[ipl];
2761  jend = match.End[jpl];
2762  icl = match.Cls[ipl];
2763  jcl = match.Cls[jpl];
2764  if (nplanes > 2) {
2765  kcl = match.Cls[kpl];
2766  kend = match.End[kpl];
2767  }
2768 
2769  geo::PlaneID const iplaneid{tpcid, ipl};
2770  geo::PlaneID const jplaneid{tpcid, jpl};
2771  geo::PlaneID const kplaneid{tpcid, kpl};
2773  okSlp = geom->ThirdPlaneSlope(ipl, oSlp[ipl], jpl, oSlp[jpl], tpcid);
2774  okAng = atan(okSlp);
2775  // handle the case near pi/2, where the errors on large slopes could result in
2776  // a wrong-sign kAng
2777  bool ignoreSign = (fabs(okSlp) > 10);
2778  if (ignoreSign) okAng = fabs(okAng);
2779  if (match.oVtx >= 0) {
2780  // a vertex exists at the other end
2781  okWir = geom->WireCoordinate(geo::Point_t{0, vtx[match.oVtx].Y, vtx[match.oVtx].Z}, kplaneid);
2782  okX = vtx[match.oVtx].X;
2783  }
2784  else {
2785  // no vertex at the other end
2787  geo::WireID(iplaneid, oWir[ipl]), geo::WireID(jplaneid, oWir[jpl]), ypos, zpos);
2788  okWir = (0.5 + geom->WireCoordinate(geo::Point_t{0, ypos, zpos}, kplaneid));
2789  okX = 0.5 * (oX[ipl] + oX[jpl]);
2790  }
2791  okTim = detProp.ConvertXToTicks(okX, kplaneid);
2792  if (prt)
2793  mf::LogVerbatim("CCTM") << "FEM: oEnd"
2794  << " kpl " << kpl << " okSlp " << okSlp << " okAng " << okAng
2795  << " okWir " << (int)okWir << " okX " << okX;
2796 
2798 
2799  kSlp = geom->ThirdPlaneSlope(
2800  ipl, clsChain[ipl][icl].Slope[iend], jpl, clsChain[jpl][jcl].Slope[jend], kplaneid);
2801  kAng = atan(kSlp);
2802  if (ignoreSign) kAng = fabs(kAng);
2803  if (match.Vtx >= 0) {
2804  if (vtx.size() == 0 || (unsigned int)match.Vtx > vtx.size() - 1) {
2805  mf::LogError("CCTM") << "FEM: Bad match.Vtx " << match.Vtx << " vtx size " << vtx.size();
2806  return;
2807  }
2808  // a vertex exists at the match end
2809  kWir = geom->WireCoordinate(geo::Point_t{0, vtx[match.Vtx].Y, vtx[match.Vtx].Z}, kplaneid);
2810  kX = vtx[match.Vtx].X;
2811  }
2812  else {
2813  // no vertex at the match end
2814  geom->IntersectionPoint(geo::WireID(iplaneid, clsChain[ipl][icl].Wire[iend]),
2815  geo::WireID(jplaneid, clsChain[jpl][jcl].Wire[jend]),
2816  ypos,
2817  zpos);
2818  kWir = (0.5 + geom->WireCoordinate(geo::Point_t{0, ypos, zpos}, kplaneid));
2819  kX = 0.5 * (clsChain[ipl][icl].X[iend] + clsChain[jpl][jcl].X[jend]);
2820  }
2821  kTim = detProp.ConvertXToTicks(kX, kplaneid);
2822  if (prt)
2823  mf::LogVerbatim("CCTM") << "FEM: mEnd"
2824  << " kpl " << kpl << " kSlp " << kSlp << " kAng " << kAng << " kX "
2825  << kX;
2826 
2827  // try to find a 3-plane match using this information
2828  if (nClInPln < 3 && FindMissingCluster(kpl, kcl, kend, kWir, kX, okWir, okX)) {
2829  nClInPln = 3;
2830  // update local variables
2831  match.Cls[kpl] = kcl;
2832  match.End[kpl] = kend;
2833  match.Chg[kpl] = clsChain[kpl][kcl].TotChg;
2834  mChg[kpl] = clsChain[kpl][kcl].TotChg;
2835  oend = 1 - kend;
2836  oWir[kpl] = clsChain[kpl][kcl].Wire[oend];
2837  oX[kpl] = clsChain[kpl][kcl].X[oend];
2838  oAng[kpl] = clsChain[kpl][kcl].Angle[oend];
2839  oSlp[kpl] = clsChain[kpl][kcl].Slope[oend];
2840  } // FindMissingCluster
2841 
2842  // decide whether to continue with a 2-plane match. The distance between match and other end should
2843  // be large enough to create a cluster
2844  if (nClInPln == 2 && fabs(okWir - kWir) > 3) return;
2845 
2846  // Calculate the cluster charge asymmetry. This factor will be applied
2847  // to the error of the end matches
2848  float chgAsym = 1;
2849  // Get the charge in the plane without a matching cluster
2850  if (nClInPln < 3 && mChg[missingPlane] <= 0) {
2851  if (missingPlane != kpl)
2852  mf::LogError("CCTM") << "FEM bad missingPlane " << missingPlane << " " << kpl << "\n";
2853  mChg[kpl] = ChargeNear(kplaneid, (unsigned short)kWir, kTim, (unsigned short)okWir, okTim);
2854  match.Chg[kpl] = mChg[kpl];
2855  if (prt)
2856  mf::LogVerbatim("CCTM") << "FEM: Missing cluster in " << kpl << " ChargeNear " << (int)kWir
2857  << ":" << (int)kTim << " " << (int)okWir << ":" << (int)okTim
2858  << " chg " << mChg[kpl];
2859  if (mChg[kpl] <= 0) return;
2860  }
2861 
2862  if (fChgAsymFactor[algIndex] > 0) {
2863  chgAsym = ChargeAsym(mChg);
2864  if (chgAsym > 0.5) return;
2865  chgAsym = 1 + fChgAsymFactor[algIndex] * chgAsym;
2866  }
2867 
2868  if (prt) mf::LogVerbatim("CCTM") << "FEM: charge asymmetry factor " << chgAsym;
2869  float sigmaX, sigmaA;
2870  float da, dx, dw;
2871 
2873  // check for vertex consistency at the match end
2874  aVtx = -1;
2875  bool allPlnVtxMatch = false;
2876  if (nClInPln == 3) {
2877  unsigned short nmvtx = 0;
2878  for (ii = 0; ii < nplanes; ++ii) {
2879  if (mVtx[ii] >= 0) {
2880  if (aVtx < 0) aVtx = mVtx[ii];
2881  ++nmvtx;
2882  }
2883  } // ii
2884  // same vertex in all planes
2885  if (nmvtx) allPlnVtxMatch = true;
2886  } // nClInPln
2887 
2888  // inflate the X match error to allow for missing one wire on the end of a cluster
2889  sigmaX = fXMatchErr[algIndex] + std::max(kSlp, (float)20);
2890  sigmaA = fAngleMatchErr[algIndex] * AngleFactor(kSlp);
2891  if (prt)
2892  mf::LogVerbatim("CCTM") << "bb " << algIndex << " " << fXMatchErr[algIndex] << " "
2893  << fAngleMatchErr[algIndex] << " kslp " << kSlp;
2894 
2895  if (nClInPln == 3) {
2896  kcl = match.Cls[kpl];
2897  kend = match.End[kpl];
2898  dw = kWir - clsChain[kpl][kcl].Wire[kend];
2899  match.dWir = dw;
2900  if (fabs(match.dWir) > 100) return;
2901  if (match.Vtx >= 0) { match.dX = kX - clsChain[kpl][kcl].X[kend]; }
2902  else {
2903  match.dX = std::abs(clsChain[ipl][icl].X[iend] - clsChain[jpl][jcl].X[jend]) +
2904  std::abs(clsChain[ipl][icl].X[iend] - clsChain[kpl][kcl].X[kend]);
2905  }
2906  if (prt) mf::LogVerbatim("CCTM") << " dw " << dw << " dx " << match.dX;
2907  // TODO: Angle matching has problems with short clusters
2908  if (!isShort) {
2909  if (ignoreSign) { match.dAng = kAng - fabs(clsChain[kpl][kcl].Angle[kend]); }
2910  else {
2911  match.dAng = kAng - clsChain[kpl][kcl].Angle[kend];
2912  }
2913  } // !isShort
2914  da = fabs(match.dAng) / sigmaA;
2915  dx = fabs(match.dX) / sigmaX;
2916  if (allPlnVtxMatch) {
2917  // matched vertex. Use angle for match error
2918  match.Err = vxFactor * chgAsym * da / 3;
2919  if (prt) mf::LogVerbatim("CCTM") << " 3-pln w Vtx match.Err " << match.Err;
2920  }
2921  else {
2922  dw /= 2;
2923  // divide by 9
2924  match.Err = vxFactor * chgAsym * sqrt(dx * dx + da * da + dw * dw) / 9;
2925  if (prt) mf::LogVerbatim("CCTM") << " 3-pln match.Err " << match.Err;
2926  }
2927  }
2928  else {
2929  // 2-plane match
2930  match.dWir = -1;
2931  match.dAng = -1;
2932  match.dX = clsChain[ipl][icl].X[iend] - clsChain[jpl][jcl].X[jend];
2933  // degrade error by 3 for 2-plane matches
2934  match.Err = 3 + vxFactor * chgAsym * fabs(match.dX) / sigmaX;
2935  if (prt) mf::LogVerbatim("CCTM") << " 2-pln Err " << match.Err;
2936  } // !(nClInPln == 3)
2937 
2939  if (nClInPln == 3) {
2940  // A cluster in all 3 planes
2941  dw = okWir - oWir[kpl];
2942  match.odWir = dw;
2943  if (match.oVtx >= 0) { match.odX = okX - oX[kpl]; }
2944  else {
2945  match.odX = std::abs(oX[ipl] - oX[jpl]) + std::abs(oX[ipl] - oX[kpl]);
2946  }
2947  if (prt)
2948  mf::LogVerbatim("CCTM") << " odw " << match.odWir << " odx " << match.odX << " sigmaX "
2949  << sigmaX;
2950  // TODO: CHECK FOR SHOWER-LIKE CLUSTER OTHER END MATCH
2951  if (!isShort) {
2952  if (ignoreSign) { match.odAng = okAng - fabs(oAng[kpl]); }
2953  else {
2954  match.odAng = okAng - oAng[kpl];
2955  }
2956  } // !isShort
2957  da = match.odAng / sigmaA;
2958  dx = fabs(match.odX) / sigmaX;
2959  // error for wire number match
2960  dw /= 2;
2961  // divide by number of planes with clusters * 3 for dx, da and dw
2962  match.oErr = ovxFactor * chgAsym * sqrt(dx * dx + da * da + dw * dw) / 9;
2963  if (prt) mf::LogVerbatim("CCTM") << " 3-pln match.oErr " << match.oErr;
2964  }
2965  else {
2966  // Only 2 clusters in 3 planes
2967  match.odX = (oX[ipl] - oX[jpl]) / sigmaX;
2968  match.oErr = 3 + ovxFactor * chgAsym * fabs(match.odX);
2969  if (prt) mf::LogVerbatim("CCTM") << " 2-pln match.oErr " << match.oErr;
2970  }
2971  } // FillEndMatch
2972 
2974  bool CCTrackMaker::FindMissingCluster(unsigned short kpl,
2975  short& kcl,
2976  unsigned short& kend,
2977  float kWir,
2978  float kX,
2979  float okWir,
2980  float okX)
2981  {
2982  // try to attach a missing cluster to the cluster chain kcl. kend is the "match end"
2983 
2984  unsigned short okend;
2985  float dxcut;
2986 
2987  if (kcl >= 0) return false;
2988 
2989  // Look for a missing cluster with loose cuts
2990  float kslp = fabs((okX - kX) / (okWir - kWir));
2991  if (kslp > 20) kslp = 20;
2992  // expand dX cut assuming there is a missing hit on the end of a cluster => 1 wire
2993  dxcut = 3 * fXMatchErr[algIndex] + kslp;
2994  unsigned short nfound = 0;
2995  unsigned short foundCl = 0, foundEnd = 0;
2996  for (unsigned short ccl = 0; ccl < clsChain[kpl].size(); ++ccl) {
2997  if (clsChain[kpl][ccl].InTrack >= 0) continue;
2998  // require a match at both ends
2999  for (unsigned short end = 0; end < 2; ++end) {
3000  okend = 1 - end;
3001  if (fabs(clsChain[kpl][ccl].Wire[end] - kWir) > 4) continue;
3002  if (fabs(clsChain[kpl][ccl].Wire[okend] - okWir) > 4) continue;
3003  // require at least one end to match
3004  if (fabs(clsChain[kpl][ccl].X[end] - kX) > dxcut &&
3005  fabs(clsChain[kpl][ccl].X[okend] - okX) > dxcut)
3006  continue;
3007  ++nfound;
3008  foundCl = ccl;
3009  foundEnd = end;
3010  } // end
3011  } // ccl
3012  if (nfound == 0) return false;
3013  if (nfound > 1) {
3014  mf::LogVerbatim("CCTM") << "FindMissingCluster: Found too many matches. Write some code "
3015  << nfound;
3016  return false;
3017  }
3018  kcl = foundCl;
3019  kend = foundEnd;
3020  return true;
3021 
3022  } // FindMissingCluster
3023 
3025  float CCTrackMaker::ChargeAsym(std::array<float, 3>& mChg)
3026  {
3027  // find charge asymmetry between the cluster with the highest (lowest)
3028  // charge
3029  float big = 0, small = 1.E9;
3030  for (unsigned short ii = 0; ii < 3; ++ii) {
3031  if (mChg[ii] < small) small = mChg[ii];
3032  if (mChg[ii] > big) big = mChg[ii];
3033  }
3034  // chgAsym varies between 0 (perfect charge match) and 1 (dreadfull)
3035  return (big - small) / (big + small);
3036  } // CalculateChargeAsym
3037 
3040  unsigned short imat,
3041  unsigned short nplanes)
3042  {
3043  // Fills the trkHits vector using cluster hits associated with the match combo imat
3044 
3045  unsigned short ipl;
3046 
3047  for (ipl = 0; ipl < 3; ++ipl)
3048  trkHits[ipl].clear();
3049 
3050  if (imat > matcomb.size()) return;
3051 
3052  unsigned short indx;
3053  std::vector<art::Ptr<recob::Hit>> clusterhits;
3054  unsigned short icc, ccl, icl, ecl, iht, ii;
3055  short endOrder, fillOrder;
3056 
3057  if (prt)
3058  mf::LogVerbatim("CCTM") << "In FillTrkHits: matcomb " << imat << " cluster chains "
3059  << matcomb[imat].Cls[0] << " " << matcomb[imat].Cls[1] << " "
3060  << matcomb[imat].Cls[2];
3061 
3062  for (ipl = 0; ipl < nplanes; ++ipl) {
3063  if (matcomb[imat].Cls[ipl] < 0) continue;
3064  // ccl is the cluster chain index
3065  ccl = matcomb[imat].Cls[ipl];
3066  // endOrder = 1 for normal order (hits added from US to DS), and -1 for reverse order
3067  endOrder = 1 - 2 * matcomb[imat].End[ipl];
3068  // re-order the sequence of cluster indices for reverse order
3069  if (endOrder < 0) {
3070  std::reverse(clsChain[ipl][ccl].ClsIndex.begin(), clsChain[ipl][ccl].ClsIndex.end());
3071  std::reverse(clsChain[ipl][ccl].Order.begin(), clsChain[ipl][ccl].Order.end());
3072  for (ii = 0; ii < clsChain[ipl][ccl].Order.size(); ++ii)
3073  clsChain[ipl][ccl].Order[ii] = 1 - clsChain[ipl][ccl].Order[ii];
3074  }
3075  if (ccl > clsChain[ipl].size() - 1) {
3076  mf::LogError("CCTM") << "Bad cluster chain index " << ccl << " in plane " << ipl;
3077  continue;
3078  }
3079  // loop over all the clusters in the chain
3080  for (icc = 0; icc < clsChain[ipl][ccl].ClsIndex.size(); ++icc) {
3081  icl = clsChain[ipl][ccl].ClsIndex[icc];
3082  if (icl > fmCluHits.size() - 1) {
3083  std::cout << "oops in FTH " << icl << " clsChain size " << clsChain[ipl].size() << "\n";
3084  exit(1);
3085  }
3086  ecl = cls[ipl][icl].EvtIndex;
3087  if (ecl > fmCluHits.size() - 1) {
3088  std::cout << "FTH bad EvtIndex " << ecl << " fmCluHits size " << fmCluHits.size() << "\n";
3089  continue;
3090  }
3091  clusterhits = fmCluHits.at(ecl);
3092  if (clusterhits.size() == 0) {
3093  std::cout << "FTH no cluster hits for EvtIndex " << cls[ipl][icl].EvtIndex << "\n";
3094  continue;
3095  }
3096  indx = trkHits[ipl].size();
3097  trkHits[ipl].resize(indx + clusterhits.size());
3098  // ensure the hit fill ordering is consistent
3099  fillOrder = 1 - 2 * clsChain[ipl][ccl].Order[icc];
3100  if (fillOrder == 1) {
3101  for (iht = 0; iht < clusterhits.size(); ++iht) {
3102  if (indx + iht > trkHits[ipl].size() - 1) std::cout << "FTH oops3\n";
3103  trkHits[ipl][indx + iht] = clusterhits[iht];
3104  }
3105  }
3106  else {
3107  for (ii = 0; ii < clusterhits.size(); ++ii) {
3108  iht = clusterhits.size() - 1 - ii;
3109  if (indx + ii > trkHits[ipl].size() - 1) std::cout << "FTH oops4\n";
3110  trkHits[ipl][indx + ii] = clusterhits[iht];
3111  } // ii
3112  }
3113  } // icc
3114  ii = trkHits[ipl].size() - 1;
3115  if (prt)
3116  mf::LogVerbatim("CCTM") << "plane " << ipl << " first p " << trkHits[ipl][0]->WireID().Plane
3117  << " w " << trkHits[ipl][0]->WireID().Wire << ":"
3118  << (int)trkHits[ipl][0]->PeakTime() << " last p "
3119  << trkHits[ipl][ii]->WireID().Plane << " w "
3120  << trkHits[ipl][ii]->WireID().Wire << ":"
3121  << (int)trkHits[ipl][ii]->PeakTime();
3122  } // ipl
3123 
3124  // TODO Check the ends of trkHits to see if there are missing hits that should have been included
3125  // in a cluster
3126 
3127  } // FillTrkHits
3128 
3131  {
3132  mf::LogVerbatim myprt("CCTM");
3133  myprt << "********* PrintTracks \n";
3134  myprt << "vtx Index X Y Z\n";
3135  for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
3136  myprt << std::right << std::setw(4) << ivx << std::setw(4) << vtx[ivx].EvtIndex;
3137  myprt << std::fixed;
3138  myprt << std::right << std::setw(10) << std::setprecision(1) << vtx[ivx].X;
3139  myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Y;
3140  myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Z;
3141  if (vtx[ivx].Neutrino) myprt << " Neutrino vertex";
3142  myprt << "\n";
3143  } // ivx
3144 
3145  myprt << ">>>>>>>>>> Tracks \n";
3146  myprt << "trk ID Proc nht nTrj sX sY sZ eX eY eZ sVx eVx sGd eGd "
3147  "ChgOrd dirZ Mom PDG ClsIndices\n";
3148  for (unsigned short itr = 0; itr < trk.size(); ++itr) {
3149  myprt << std::right << std::setw(3) << itr << std::setw(4) << trk[itr].ID;
3150  myprt << std::right << std::setw(5) << std::setw(4) << trk[itr].Proc;
3151  unsigned short nht = 0;
3152  for (unsigned short ii = 0; ii < 3; ++ii)
3153  nht += trk[itr].TrkHits[ii].size();
3154  myprt << std::right << std::setw(5) << nht;
3155  myprt << std::setw(4) << trk[itr].TrjPos.size();
3156  myprt << std::fixed;
3157  myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[0](0);
3158  myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[0](1);
3159  myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[0](2);
3160  unsigned short itj = trk[itr].TrjPos.size() - 1;
3161  myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[itj](0);
3162  myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[itj](1);
3163  myprt << std::right << std::setw(7) << std::setprecision(1) << trk[itr].TrjPos[itj](2);
3164  myprt << std::setw(4) << trk[itr].VtxIndex[0] << std::setw(4) << trk[itr].VtxIndex[1];
3165  myprt << std::setw(4) << trk[itr].GoodEnd[0];
3166  myprt << std::setw(4) << trk[itr].GoodEnd[1];
3167  myprt << std::setw(4) << trk[itr].ChgOrder;
3168  myprt << std::right << std::setw(10) << std::setprecision(3) << trk[itr].TrjDir[itj](2);
3169  myprt << std::right << std::setw(4) << trk[itr].MomID;
3170  myprt << std::right << std::setw(5) << trk[itr].PDGCode << " ";
3171  for (unsigned short ii = 0; ii < trk[itr].ClsEvtIndices.size(); ++ii)
3172  myprt << " " << trk[itr].ClsEvtIndices[ii];
3173  myprt << "\n";
3174  } // itr
3175 
3176  } // PrintTracks
3177 
3180  geo::TPCID const& tpcid) const
3181  {
3182 
3183  unsigned short iTime;
3184  mf::LogVerbatim myprt("CCTM");
3185  myprt << "******* PrintClusters ********* Num_Clusters_in Wire:Time\n";
3186  myprt << "vtx Index X Y Z Pln0 Pln1 Pln2 Pln0 Pln1 Pln2\n";
3187  for (unsigned short ivx = 0; ivx < vtx.size(); ++ivx) {
3188  myprt << std::right << std::setw(3) << ivx << std::setw(7) << ivx;
3189  myprt << std::fixed;
3190  myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].X;
3191  myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Y;
3192  myprt << std::right << std::setw(7) << std::setprecision(1) << vtx[ivx].Z;
3193  myprt << std::right << std::setw(5) << vtx[ivx].nClusInPln[0];
3194  myprt << std::right << std::setw(5) << vtx[ivx].nClusInPln[1];
3195  myprt << std::right << std::setw(5) << vtx[ivx].nClusInPln[2];
3196  myprt << " ";
3197  for (auto const& planeID : geom->Iterate<geo::PlaneID>(tpcid)) {
3198  int time = (0.5 + detProp.ConvertXToTicks(vtx[ivx].X, planeID));
3199  int wire = geom->WireCoordinate(geo::Point_t{0, vtx[ivx].Y, vtx[ivx].Z}, planeID);
3200  myprt << std::right << std::setw(7) << wire << ":" << time;
3201  }
3202 
3203  myprt << "\n";
3204  } // ivx
3205 
3206  auto const nplanes = geom->Nplanes(tpcid);
3207  for (unsigned short ipl = 0; ipl < nplanes; ++ipl) {
3208  myprt << ">>>>>>>>>> Cluster chains in Plane " << ipl << "\n";
3209  myprt << "ipl ccl Len Chg W0:T0 Ang0 Dir0 Vx0 mBk0 W1:T1 Ang1 Dir1 Vx1 "
3210  " mBk1 InTk cls:Order \n";
3211  for (unsigned short ccl = 0; ccl < clsChain[ipl].size(); ++ccl) {
3212  myprt << std::right << std::setw(3) << ipl;
3213  myprt << std::right << std::setw(5) << ccl;
3214  myprt << std::right << std::setw(6) << clsChain[ipl][ccl].Length;
3215  myprt << std::right << std::setw(8) << (int)clsChain[ipl][ccl].TotChg;
3216  for (unsigned short end = 0; end < 2; ++end) {
3217  iTime = clsChain[ipl][ccl].Time[end];
3218  myprt << std::right << std::setw(5) << (int)clsChain[ipl][ccl].Wire[end] << ":"
3219  << std::setprecision(1) << iTime;
3220  if (iTime < 10) { myprt << " "; }
3221  else if (iTime < 100) {
3222  myprt << " ";
3223  }
3224  else if (iTime < 1000)
3225  myprt << " ";
3226  myprt << std::right << std::setw(7) << std::setprecision(2)
3227  << clsChain[ipl][ccl].Angle[end];
3228  myprt << std::right << std::setw(5) << clsChain[ipl][ccl].Dir[end];
3229  myprt << std::right << std::setw(5) << clsChain[ipl][ccl].VtxIndex[end];
3230  myprt << std::fixed << std::right << std::setw(6) << std::setprecision(1)
3231  << clsChain[ipl][ccl].mBrkIndex[end];
3232  }
3233  myprt << std::right << std::setw(7) << clsChain[ipl][ccl].InTrack;
3234  myprt << " ";
3235  for (unsigned short ii = 0; ii < clsChain[ipl][ccl].ClsIndex.size(); ++ii)
3236  myprt << " " << clsChain[ipl][ccl].ClsIndex[ii] << ":" << clsChain[ipl][ccl].Order[ii];
3237  myprt << "\n";
3238  } // ccl
3239  if (fPrintAllClusters) {
3240  myprt << ">>>>>>>>>> Clusters in Plane " << ipl << "\n";
3241  myprt << "ipl icl Evt Len Chg W0:T0 Ang0 Dir0 Vx0 CN0 W1:T1 Ang1 "
3242  "Dir1 Vx1 CN1 InTk Brk0 MrgEr0 Brk1 MrgEr1\n";
3243  for (unsigned short icl = 0; icl < cls[ipl].size(); ++icl) {
3244  myprt << std::right << std::setw(3) << ipl;
3245  myprt << std::right << std::setw(5) << icl;
3246  myprt << std::right << std::setw(5) << cls[ipl][icl].EvtIndex;
3247  myprt << std::right << std::setw(6) << cls[ipl][icl].Length;
3248  myprt << std::right << std::setw(8) << (int)cls[ipl][icl].TotChg;
3249  for (unsigned short end = 0; end < 2; ++end) {
3250  iTime = cls[ipl][icl].Time[end];
3251  myprt << std::right << std::setw(5) << (int)cls[ipl][icl].Wire[end] << ":" << iTime;
3252  if (iTime < 10) { myprt << " "; }
3253  else if (iTime < 100) {
3254  myprt << " ";
3255  }
3256  else if (iTime < 1000)
3257  myprt << " ";
3258  myprt << std::right << std::setw(7) << std::setprecision(2) << cls[ipl][icl].Angle[end];
3259  myprt << std::right << std::setw(5) << cls[ipl][icl].Dir[end];
3260  myprt << std::right << std::setw(5) << cls[ipl][icl].VtxIndex[end];
3261  myprt << std::fixed << std::right << std::setw(5) << std::setprecision(1)
3262  << cls[ipl][icl].ChgNear[end];
3263  }
3264  myprt << std::fixed;
3265  myprt << std::right << std::setw(5) << cls[ipl][icl].InTrack;
3266  myprt << std::right << std::setw(5) << (int)cls[ipl][icl].BrkIndex[0];
3267  myprt << std::right << std::setw(7) << std::setprecision(1)
3268  << cls[ipl][icl].MergeError[0];
3269  myprt << std::right << std::setw(5) << (int)cls[ipl][icl].BrkIndex[1];
3270  myprt << std::right << std::setw(7) << std::setprecision(1)
3271  << cls[ipl][icl].MergeError[1];
3272  myprt << "\n";
3273  } // icl
3274  } // fPrintAllClusters
3275  } // ipl
3276  } // PrintClusters
3277 
3279  float CCTrackMaker::AngleFactor(float slope)
3280  {
3281  float slp = fabs(slope);
3282  if (slp > 10.) slp = 30.;
3283  // return a value between 1 and 46
3284  return 1 + 0.05 * slp * slp;
3285  } // AngleFactor
3286 
3289  unsigned short wire1,
3290  float time1,
3291  unsigned short wire2,
3292  float time2)
3293  {
3294  // returns the hit charge along a line between (wire1, time1) and
3295  // (wire2, time2)
3296 
3297  // put in increasing wire order (wire2 > wire1)
3298  unsigned short w1 = wire1;
3299  unsigned short w2 = wire2;
3300  double t1 = time1;
3301  double t2 = time2;
3302  double slp, prtime;
3303  if (w1 == w2) { slp = 0; }
3304  else {
3305  if (w1 > w2) {
3306  w1 = wire2;
3307  w2 = wire1;
3308  t1 = time2;
3309  t2 = time1;
3310  }
3311  slp = (t2 - t1) / (w2 - w1);
3312  }
3313 
3314  unsigned short wire;
3315 
3316  float chg = 0;
3317  for (unsigned short hit = 0; hit < allhits.size(); ++hit) {
3318  if (allhits[hit]->WireID().asPlaneID() != planeid) continue;
3319  wire = allhits[hit]->WireID().Wire;
3320  if (wire < w1) continue;
3321  if (wire > w2) continue;
3322  prtime = t1 + (wire - w1) * slp;
3323  if (prtime > allhits[hit]->PeakTimePlusRMS(3)) continue;
3324  if (prtime < allhits[hit]->PeakTimeMinusRMS(3)) continue;
3325  chg += ChgNorm[planeid.Plane] * allhits[hit]->Integral();
3326  } // hit
3327  return chg;
3328  } // ChargeNear
3329 
3332  {
3333  // fills the WireHitRange vector. Slightly modified version of the one in ClusterCrawlerAlg
3334 
3335  unsigned short ipl;
3336 
3337  // initialize everything
3338  for (ipl = 0; ipl < 3; ++ipl) {
3339  firstWire[ipl] = INT_MAX;
3340  lastWire[ipl] = 0;
3341  firstHit[ipl] = INT_MAX;
3342  lastHit[ipl] = 0;
3343  WireHitRange[ipl].clear();
3344  ChgNorm[ipl] = 0;
3345  }
3346 
3347  // find the first and last wire with a hit
3348  unsigned short oldipl = 0;
3349  for (unsigned int hit = 0; hit < allhits.size(); ++hit) {
3350  if (static_cast<geo::TPCID const&>(allhits[hit]->WireID()) != tpcid) continue;
3351  ipl = allhits[hit]->WireID().Plane;
3352  if (allhits[hit]->WireID().Wire > geom->Nwires(allhits[hit]->WireID().asPlaneID())) {
3353  if (lastWire[ipl] < firstWire[ipl]) {
3354  mf::LogError("CCTM") << "Invalid WireID().Wire " << allhits[hit]->WireID().Wire;
3355  return;
3356  }
3357  }
3358  // ChgNorm[ipl] += allhits[hit]->Integral();
3359  if (ipl < oldipl) {
3360  mf::LogError("CCTM") << "Hits are not in increasing-plane order\n";
3361  return;
3362  }
3363  oldipl = ipl;
3364  if (firstHit[ipl] == INT_MAX) {
3365  firstHit[ipl] = hit;
3366  firstWire[ipl] = allhits[hit]->WireID().Wire;
3367  }
3368  lastHit[ipl] = hit;
3369  lastWire[ipl] = allhits[hit]->WireID().Wire;
3370  } // hit
3371 
3372  // xxx
3373  auto const nplanes = geom->Nplanes(tpcid);
3374  for (ipl = 0; ipl < nplanes; ++ipl) {
3375  if (lastWire[ipl] < firstWire[ipl]) {
3376  mf::LogError("CCTM") << "Invalid first/last wire in plane " << ipl;
3377  return;
3378  }
3379  } // ipl
3380 
3381  // normalize charge in induction planes to the collection plane
3382  // for(ipl = 0; ipl < nplanes; ++ipl) ChgNorm[ipl] = ChgNorm[nplanes - 1] / ChgNorm[ipl];
3383  for (ipl = 0; ipl < nplanes; ++ipl)
3384  ChgNorm[ipl] = 1;
3385 
3386  // get the service to learn about channel status
3387  //lariov::ChannelStatusProvider const& channelStatus
3388  // = art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider();
3389 
3390  // now we can define the WireHitRange vector.
3391  int sflag, nwires, wire;
3392  unsigned int indx, thisWire, thisHit, lastFirstHit;
3393  //uint32_t chan;
3394  for (ipl = 0; ipl < nplanes; ++ipl) {
3395  nwires = lastWire[ipl] - firstWire[ipl] + 1;
3396  WireHitRange[ipl].resize(nwires);
3397  // start by defining the "no hits on wire" condition
3398  sflag = -2;
3399  for (wire = 0; wire < nwires; ++wire)
3400  WireHitRange[ipl][wire] = std::make_pair(sflag, sflag);
3401  // overwrite with the "dead wires" condition
3402  sflag = -1;
3403  // next overwrite with the index of the first/last hit on each wire
3404  lastWire[ipl] = firstWire[ipl];
3405  thisHit = firstHit[ipl];
3406  lastFirstHit = firstHit[ipl];
3407  for (unsigned int hit = firstHit[ipl]; hit <= lastHit[ipl]; ++hit) {
3408  thisWire = allhits[hit]->WireID().Wire;
3409  if (thisWire > lastWire[ipl]) {
3410  indx = lastWire[ipl] - firstWire[ipl];
3411  int tmp1 = lastFirstHit;
3412  int tmp2 = thisHit;
3413  WireHitRange[ipl][indx] = std::make_pair(tmp1, tmp2);
3414  lastWire[ipl] = thisWire;
3415  lastFirstHit = thisHit;
3416  }
3417  else if (thisWire < lastWire[ipl]) {
3418  mf::LogError("CCTM") << "Hit not in proper order in plane " << ipl;
3419  exit(1);
3420  }
3421  ++thisHit;
3422  } // hit
3423  // define the last wire
3424  indx = lastWire[ipl] - firstWire[ipl];
3425  int tmp1 = lastFirstHit;
3426  ++lastHit[ipl];
3427  int tmp2 = lastHit[ipl];
3428  WireHitRange[ipl][indx] = std::make_pair(tmp1, tmp2);
3429  // add one to lastWire and lastHit for more standard indexing
3430  ++lastWire[ipl];
3431  } // ipl
3432 
3433  // error checking
3434  for (ipl = 0; ipl < nplanes; ++ipl) {
3435  if (firstWire[ipl] < INT_MAX) continue;
3436  if (lastWire[ipl] > 0) continue;
3437  if (firstHit[ipl] < INT_MAX) continue;
3438  if (lastHit[ipl] > 0) continue;
3439  std::cout << "FWHR problem\n";
3440  exit(1);
3441  } // ipl
3442 
3443  unsigned int nht = 0;
3444  std::vector<bool> hchk(allhits.size());
3445  for (unsigned int ii = 0; ii < hchk.size(); ++ii)
3446  hchk[ii] = false;
3447  for (ipl = 0; ipl < nplanes; ++ipl) {
3448  for (unsigned int w = firstWire[ipl]; w < lastWire[ipl]; ++w) {
3449  indx = w - firstWire[ipl];
3450  if (indx > lastWire[ipl]) {
3451  std::cout << "FWHR bad index " << indx << "\n";
3452  exit(1);
3453  }
3454  // no hit on this wire
3455  if (WireHitRange[ipl][indx].first < 0) continue;
3456  unsigned int firhit = WireHitRange[ipl][indx].first;
3457  unsigned int lashit = WireHitRange[ipl][indx].second;
3458  for (unsigned int hit = firhit; hit < lashit; ++hit) {
3459  ++nht;
3460  if (hit > allhits.size() - 1) {
3461  std::cout << "FWHR: Bad hchk " << hit << " size " << allhits.size() << "\n";
3462  continue;
3463  }
3464  hchk[hit] = true;
3465  if (allhits[hit]->WireID().Plane != ipl || allhits[hit]->WireID().Wire != w) {
3466  std::cout << "FWHR bad plane " << allhits[hit]->WireID().Plane << " " << ipl
3467  << " or wire " << allhits[hit]->WireID().Wire << " " << w << "\n";
3468  exit(1);
3469  }
3470  } // hit
3471  } // w
3472  } // ipl
3473 
3474  if (nht != allhits.size()) {
3475  std::cout << "FWHR hit count problem " << nht << " " << allhits.size() << "\n";
3476  for (unsigned int ii = 0; ii < hchk.size(); ++ii)
3477  if (!hchk[ii])
3478  std::cout << " " << ii << " " << allhits[ii]->WireID().Plane << " "
3479  << allhits[ii]->WireID().Wire << " " << (int)allhits[ii]->PeakTime() << "\n";
3480  exit(1);
3481  }
3482 
3483  } // FillWireHitRange
3485 
3486 } // namespace
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
float Length(const PFPStruct &pfp)
Definition: PFPUtils.cxx:3269
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
std::vector< MatchPars > matcomb
float ChargeNear(geo::PlaneID const &planeid, unsigned short wire1, float time1, unsigned short wire2, float time2)
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
art::ServiceHandle< geo::Geometry const > geom
TTree * t1
Definition: plottest35.C:26
std::string fVertexModuleLabel
void StoreTrack(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, unsigned short imat, unsigned short procCode, geo::TPCID const &tpcid)
std::array< float, 3 > ChgNorm
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:78
std::vector< float > fChgAsymFactor
std::array< unsigned short, 3 > End
Declaration of signal hit object.
void PrintClusters()
std::array< unsigned int, 3 > firstWire
void MakeClusterChains(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, geo::TPCID const &tpcid)
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
Float_t Y
Definition: plot.C:37
Planes which measure X direction.
Definition: geo_types.h:140
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
Geometry information for a single TPC.
Definition: TPCGeo.h:36
std::vector< art::Ptr< recob::Hit > > allhits
void FindMaybeVertices(geo::TPCID const &tpcid)
std::array< short, 2 > Dir
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:67
constexpr auto abs(T v)
Returns the absolute value of the argument.
float StartWire() const
Returns the wire coordinate of the start of the cluster.
Definition: Cluster.h:276
float dXClTraj(art::FindManyP< recob::Hit > const &fmCluHits, unsigned short ipl, unsigned short icl1, unsigned short end1)
std::array< std::vector< art::Ptr< recob::Hit > >, 3 > seedHits
std::string fClusterModuleLabel
void FitVertices(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)
TrackTrajectoryAlg fTrackTrajectoryAlg
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
void FillWireHitRange(geo::TPCID inTPCID)
Definition: Utils.cxx:4323
Set of hits with a 2D structure.
Definition: Cluster.h:69
std::array< short, 2 > BrkIndex
float EndTick() const
Returns the tick coordinate of the end of the cluster.
Definition: Cluster.h:331
Float_t tmp
Definition: plot.C:35
std::vector< unsigned short > pfpToTrkID
void FillTrkHits(art::FindManyP< recob::Hit > const &fmCluHits, unsigned short imat, unsigned short nplanes)
data_const_reference data(size_type i) const
Definition: FindManyP.h:446
std::vector< TVector3 > TrjPos
std::vector< TrkPar > trk
bool FindMissingCluster(unsigned short kpl, short &kcl, unsigned short &kend, float kWir, float kX, float okWir, float okX)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Cluster finding and building.
std::vector< float > fAngleMatchErr
std::array< float, 2 > ChgNear
std::array< std::vector< clPar >, 3 > cls
float StartAngle() const
Returns the starting angle of the cluster.
Definition: Cluster.h:461
Definition of vertex object for LArSoft.
Definition: Vertex.h:35
std::array< float, 2 > Slope
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
std::array< bool, 2 > EndInTPC
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:104
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
Float_t Z
Definition: plot.C:37
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::array< std::vector< art::Ptr< recob::Hit > >, 3 > trkHits
std::array< std::vector< ClsChainPar >, 3 > clsChain
float EndCharge() const
Returns the charge on the last wire of the cluster.
Definition: Cluster.h:483
Point_t GetCenter() const
Returns the center of the TPC volume in world coordinates [cm].
Definition: TPCGeo.h:247
TCanvas * c1
Definition: plotHisto.C:7
TCanvas * c2
Definition: plot_hist.C:75
void TagCosmics(geo::TPCID const &tpcid)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
long seed
Definition: chem4.cc:67
A trajectory in space reconstructed from hits.
void FillEndMatch(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid, MatchPars &match)
std::vector< float > fMatchMinLen
std::vector< unsigned short > Order
float AngleFactor(float slope)
std::array< short, 2 > VtxIndex
std::array< short, 2 > mVtxIndex
double ConvertXToTicks(double X, int p, int t, int c) const
std::vector< vtxPar > vtx
std::array< float, 2 > MergeError
Provides recob::Track data product.
void produce(art::Event &evt) override
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
std::array< bool, 2 > GoodEnd
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
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.
void TrackTrajectory(std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
bool DupMatch(MatchPars &match, unsigned short nplanes)
std::array< unsigned int, 3 > firstHit
std::array< float, 2 > Charge
void VtxMatch(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, geo::TPCID const &tpcid)
std::vector< unsigned short > ClsIndex
std::vector< unsigned short > ClsEvtIndices
Detector simulation of raw signals on wires.
std::array< unsigned int, 3 > lastWire
TTree * t2
Definition: plottest35.C:36
Class defining a plane for tracking. It provides various functionalities to convert track parameters ...
Definition: TrackingPlane.h:37
void MakeFamily(geo::TPCID const &tpcid)
double ConvertTicksToX(double ticks, int p, int t, int c) const
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:100
bool greaterThan(CluLen c1, CluLen c2)
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
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
void VertexFit(std::vector< std::vector< geo::WireID >> const &hitWID, std::vector< std::vector< double >> const &hitX, std::vector< std::vector< double >> const &hitXErr, TVector3 &VtxPos, TVector3 &VtxPosErr, std::vector< TVector3 > &TrkDir, std::vector< TVector3 > &TrkDirErr, float &ChiDOF) const
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
void FillEndMatch2(MatchPars &match)
Utility object to perform functions of association.
Float_t norm
float StartCharge() const
Returns the charge on the first wire of the cluster.
Definition: Cluster.h:440
TDirectory * dir
Definition: macro.C:5
std::array< short, 2 > VtxIndex
void PrintClusters(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid) const
std::array< unsigned int, 3 > lastHit
std::array< float, 2 > X
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double ThirdPlaneSlope(PlaneID const &pid1, double slope1, PlaneID const &pid2, double slope2, PlaneID const &output_plane) const
Returns the slope on the third plane, given it in the other two.
void PlnMatch(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
bool lessThan(CluLen c1, CluLen c2)
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
std::vector< TVector3 > TrjDir
void FillWireHitRange(geo::TPCID const &tpcid)
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
std::vector< bool > fMakeAlgTracks
std::vector< short > fMatchAlgs
std::vector< float > fXMatchErr
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
float EndAngle() const
Returns the ending angle of the cluster.
Definition: Cluster.h:504
recob::tracking::Plane Plane
Definition: TrackState.h:17
float ChargeAsym(std::array< float, 3 > &mChg)
Float_t X
Definition: plot.C:37
Float_t track
Definition: plot.C:35
std::array< float, 2 > Wire
Float_t w
Definition: plot.C:20
float StartTick() const
Returns the tick coordinate of the start of the cluster.
Definition: Cluster.h:287
std::array< std::vector< art::Ptr< recob::Hit > >, 3 > TrkHits
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
void FillChgNear(detinfo::DetectorPropertiesData const &detProp, geo::TPCID const &tpcid)
std::array< short, 2 > Time
art framework interface to geometry description
double HalfWidth() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:96
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
std::array< float, 2 > Angle
vec_iX clear()
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane.
std::array< std::vector< std::pair< int, int > >, 3 > WireHitRange
std::array< std::vector< unsigned short >, 3 > vxCls
std::array< unsigned short, 3 > nClusInPln
unsigned short fNVtxTrkHitsFit
float Integral() const
Returns the total charge of the cluster from hit shape.
Definition: Cluster.h:581
float EndWire() const
Returns the wire coordinate of the end of the cluster.
Definition: Cluster.h:318
void SortMatches(detinfo::DetectorPropertiesData const &detProp, art::FindManyP< recob::Hit > const &fmCluHits, unsigned short procCode, geo::TPCID const &tpcid)
vertex reconstruction