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