LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
TCVertex.cxx
Go to the documentation of this file.
1 #include "larreco/RecoAlg/TCAlg/TCVertex.h"
2 
3 namespace tca {
4 
5  struct SortEntry{
6  unsigned int index;
7  float val;
8  };
9 
10  bool valDecreasing (SortEntry c1, SortEntry c2) { return (c1.val > c2.val);}
11  bool valIncreasing (SortEntry c1, SortEntry c2) { return (c1.val < c2.val);}
12 
14  void MakeJunkVertices(TCSlice& slc, const CTP_t& inCTP)
15  {
16  // Vertices between poorly reconstructed tjs (especially junk slc) and normal
17  // tjs can fail because the junk tj trajectory parameters are inaccurate. This function
18  // uses proximity and not pointing to make junk vertices
19  // Don't use this if standard vertex reconstruction is disabled
20  if(tcc.vtx2DCuts[0] <= 0) return;
21  if(!tcc.useAlg[kJunkVx]) return;
22  if(slc.tjs.size() < 2) return;
23 
24  // Look for tjs that are within maxSep of the end of a Tj
25  constexpr float maxSep = 4;
26 
27  geo::PlaneID planeID = DecodeCTP(inCTP);
28  bool prt = (tcc.dbgVxJunk && tcc.dbgSlc);
29  if(prt) {
30  mf::LogVerbatim("TC")<<"MakeJunkVertices: prt set for plane "<<planeID.Plane<<" maxSep btw tjs "<<maxSep;
31 // PrintAllTraj("MJTi", slc, debug, USHRT_MAX, slc.tjs.size());
32  }
33 
34  // make a template vertex
35  VtxStore junkVx;
36  junkVx.CTP = inCTP;
37  junkVx.Topo = 9;
38  junkVx.Stat[kJunkVx] = true;
39  junkVx.Stat[kFixed] = true;
40  // set an invalid ID
41  junkVx.ID = USHRT_MAX;
42  // put in generous errors
43  // BUG the double brace syntax is required to work around clang bug 21629
44  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
45  junkVx.PosErr = {{2.0, 2.0}};
46  // define a minimal score so it won't get clobbered
47  junkVx.Score = tcc.vtx2DCuts[7] + 0.1;
48 
49  // look at both ends of long tjs
50  for(unsigned short it1 = 0; it1 < slc.tjs.size() - 1; ++it1) {
51  auto& tj1 = slc.tjs[it1];
52  if(tj1.AlgMod[kKilled] || tj1.AlgMod[kHaloTj]) continue;
53  if(tj1.SSID > 0 || tj1.AlgMod[kShowerLike]) continue;
54  if(tj1.CTP != inCTP) continue;
55  if(tj1.AlgMod[kJunkTj]) continue;
56  if(TrajLength(tj1) < 10) continue;
57  if(tj1.MCSMom < 100) continue;
58  for(unsigned short end1 = 0; end1 < 2; ++end1) {
59  // existing vertex?
60  if(tj1.VtxID[end1] > 0) continue;
61  auto& tp1 = tj1.Pts[tj1.EndPt[end1]];
62  // get a list of tjs in this vicinity
63  auto tjlist = FindCloseTjs(slc, tp1, tp1, maxSep);
64  if(tjlist.empty()) continue;
65  // set to an invalid ID
66  junkVx.ID = USHRT_MAX;
67  for(auto tj2id : tjlist) {
68  auto& tj2 = slc.tjs[tj2id - 1];
69  if(tj2.CTP != inCTP) continue;
70  if(tj2id == tj1.ID) continue;
71  if(tj2.SSID > 0 || tj2.AlgMod[kShowerLike]) continue;
72  float close = maxSep;
73  unsigned short closeEnd = USHRT_MAX;
74  for(unsigned short end2 = 0; end2 < 2; ++end2) {
75  auto& tp2 = tj2.Pts[tj2.EndPt[end2]];
76  float sep = PosSep(tp1.Pos, tp2.Pos);
77  if(sep < close) {
78  close = sep;
79  closeEnd = end2;
80  } // sep
81  } // end2
82  if(closeEnd > 1) continue;
83  auto& tp2 = tj2.Pts[tj2.EndPt[closeEnd]];
84  bool signalBetween = SignalBetween(slc, tp1, tp2, 0.8);
85  if(!signalBetween) continue;
86  if(junkVx.ID == USHRT_MAX) {
87  // define the new vertex
88  junkVx.ID = slc.vtxs.size() + 1;
89  junkVx.Pos = tp1.Pos;
90  } // new vertex
91  tj2.VtxID[closeEnd] = junkVx.ID;
92  tj1.VtxID[end1] = junkVx.ID;
93  } // tjid
94  if(junkVx.ID == USHRT_MAX) continue;
95  if(!StoreVertex(slc, junkVx)) {
96  mf::LogVerbatim("TC")<<"MJV: StoreVertex failed";
97  for(auto& tj : slc.tjs) {
98  if(tj.AlgMod[kKilled]) continue;
99  if(tj.VtxID[0] == junkVx.ID) tj.VtxID[0] = 0;
100  if(tj.VtxID[1] == junkVx.ID) tj.VtxID[1] = 0;
101  } // tj
102  continue;
103  } // StoreVertex failed
104  if(prt) {
105  mf::LogVerbatim("TC")<<" New junk 2V"<<junkVx.ID<<" at "<<std::fixed<<std::setprecision(1)<<junkVx.Pos[0]<<":"<<junkVx.Pos[1]/tcc.unitsPerTick;
106  } // prt
107  junkVx.ID = USHRT_MAX;
108  } // end1
109  } // it1
110 
111  } // MakeJunkVertices
112 
114  void Find2DVertices(TCSlice& slc, const CTP_t& inCTP, unsigned short pass)
115  {
116  // Find 2D vertices between pairs of tjs that have a same-end topology. Using an example
117  // where StepDir = 1 (end 0 is at small wire number) vertices will be found with Topo = 0
118  // with a vertex US of the ends (<) or Topo = 2 with a vertex DS of the ends (>). This is reversed
119  // if StepDir = -1. Vertices with Topo = 1 (/\) and (\/) are found in EndMerge.
120 
121  // tcc.vtx2DCuts fcl input usage
122  // 0 = maximum length of a short trajectory
123  // 1 = max vertex - trajectory separation for short trajectories
124  // 2 = max vertex - trajectory separation for long trajectories
125  // 3 = max position pull for adding TJs to a vertex
126  // 4 = max allowed vertex position error
127  // 5 = min MCSMom
128  // 6 = min Pts/Wire fraction
129  // 7 min Score
130  // 8 Min charge fraction near a merge point (not a vertex)
131  // 9 max MCSmom asymmetry for a merge
132  // 10 Require charge on wires between a vtx and the start of the tjs in induction planes? (1 = yes)
133 
134  if(tcc.vtx2DCuts[0] <= 0) return;
135  if(slc.tjs.size() < 2) return;
136 
137  bool firstPassCuts = (tcc.useAlg[kNewVtxCuts] && pass == 0);
138 
139  geo::PlaneID planeID = DecodeCTP(inCTP);
140 
141  // require charge between the vertex and the tj start points?
142  bool requireVtxTjChg = true;
143  if(tcc.vtx2DCuts[10] == 0 && int(planeID.Plane) < slc.nPlanes - 1) requireVtxTjChg = false;
144 
145  bool prt = (tcc.dbg2V && tcc.dbgSlc && debug.Plane == (int)planeID.Plane);
146  if(prt) {
147  mf::LogVerbatim("TC")<<"prt set for plane "<<planeID.Plane<<" in Find2DVertices. firstPassCuts? "<<firstPassCuts<<" requireVtxTjChg "<<requireVtxTjChg;
148  PrintAllTraj("F2DVi", slc, USHRT_MAX, slc.tjs.size());
149  }
150 
151  unsigned short maxShortTjLen = tcc.vtx2DCuts[0];
152  for(unsigned short it1 = 0; it1 < slc.tjs.size() - 1; ++it1) {
153  auto& tj1 = slc.tjs[it1];
154  if(tj1.AlgMod[kKilled] || tj1.AlgMod[kHaloTj]) continue;
155  if(tj1.SSID > 0 || tj1.AlgMod[kShowerLike]) continue;
156  if(tj1.CTP != inCTP) continue;
157  bool tj1Short = (TrajLength(tj1) < maxShortTjLen);
158  for(unsigned short end1 = 0; end1 < 2; ++end1) {
159  // vertex assignment exists?
160  if(tj1.VtxID[end1] > 0) continue;
161  // wrong end of a high energy electron?
162  if(tj1.PDGCode == 111 && end1 != tj1.StartEnd) continue;
163  // default condition is to use the end point to define the trajectory and direction
164  // at the end
165  short endPt1 = tj1.EndPt[end1];
166  float wire1 = tj1.Pts[endPt1].Pos[0];
167  // unless there are few points fitted, indicating that the trajectory fit
168  // may have been biased by the presence of another trajectory at the vertex or by
169  // other close unresolved tracks
170  if(tj1.Pts.size() > 6 && tj1.Pts[endPt1].NTPsFit < 4) {
171  if(end1 == 0 && endPt1 < int(tj1.Pts.size()) - 3) {
172  endPt1 += 3;
173  } else if (end1 == 1 && endPt1 >=3 ) {
174  endPt1 -= 3;
175  }
176  if(tj1.Pts[endPt1].Chg == 0) endPt1 = NearestPtWithChg(slc, tj1, endPt1);
177  } // few points fit at end1
178  TrajPoint tp1 = tj1.Pts[endPt1];
179  MoveTPToWire(tp1, wire1);
180  // re-purpose endPt1 to reference the end point. This will be used the find the point on
181  // tj1 that is closest to the vertex position
182  endPt1 = tj1.EndPt[end1];
183  short oendPt1 = tj1.EndPt[1-end1];
184  // reference to the other end of tj1
185  auto& otp1 = tj1.Pts[oendPt1];
186  for(unsigned short it2 = it1 + 1; it2 < slc.tjs.size(); ++it2) {
187  auto& tj2 = slc.tjs[it2];
188  if(tj2.AlgMod[kKilled] || tj2.AlgMod[kHaloTj]) continue;
189  if(tj2.SSID > 0 || tj2.AlgMod[kShowerLike]) continue;
190  if(tj2.CTP != inCTP) continue;
191  if(tj1.VtxID[end1] > 0) continue;
192  if(tj1.MCSMom < tcc.vtx2DCuts[5] && tj2.MCSMom < tcc.vtx2DCuts[5]) continue;
193  bool tj2Short = (TrajLength(tj2) < maxShortTjLen);
194  // find the end that is closer to tp1
195  unsigned short end2 = 0;
196  if(PosSep2(tj2.Pts[tj2.EndPt[1]].Pos, tp1.Pos) < PosSep2(tj2.Pts[tj2.EndPt[0]].Pos, tp1.Pos)) end2 = 1;
197  if(tj2.VtxID[end2] > 0) continue;
198  // wrong end of a high energy electron?
199  if(tj2.PDGCode == 111 && end2 != tj2.StartEnd) continue;
200  // check for a vertex between these tjs at the other ends
201  if(tj1.VtxID[1 - end1] > 0 && tj1.VtxID[1 - end1] == tj2.VtxID[1 - end2]) continue;
202  // see if the other ends are closer
203  unsigned short oendPt2 = tj2.EndPt[1-end2];
204  auto& otp2 = tj2.Pts[oendPt2];
205  if(PosSep2(otp1.Pos, otp2.Pos) < PosSep2(tp1.Pos, tj2.Pts[tj2.EndPt[end2]].Pos)) continue;
206  short endPt2 = tj2.EndPt[end2];
207  float wire2 = tj2.Pts[endPt2].Pos[0];
208  if(tj2.Pts.size() > 6 && tj2.Pts[endPt2].NTPsFit < 4) {
209  if(end2 == 0 && endPt2 < int(tj2.Pts.size()) - 3) {
210  endPt2 += 3;
211  } else if (end2 == 1 && endPt2 >= 3){
212  endPt2 -= 3;
213  }
214  if(tj2.Pts[endPt2].Chg == 0) endPt2 = NearestPtWithChg(slc, tj2, endPt2);
215  } // few points fit at end1
216  TrajPoint tp2 = tj2.Pts[endPt2];
217  MoveTPToWire(tp2, wire2);
218  // re-purpose endPt2
219  endPt2 = tj2.EndPt[end2];
220  // Rough first cut on the separation between the end points of the
221  // two trajectories
222  float sepCut = 100;
223  if(std::abs(tp1.Pos[0] - tp2.Pos[0]) > sepCut) continue;
224  if(std::abs(tp1.Pos[1] - tp2.Pos[1]) > sepCut) continue;
225  float wint, tint;
226  TrajIntersection(tp1, tp2, wint, tint);
227  // make sure this is inside the TPC.
228  if(wint < 0 || wint > tcc.maxPos0[planeID.Plane] - 3) continue;
229  if(tint < 0 || tint > tcc.maxPos1[planeID.Plane]) continue;
230  // Next cut on separation between the TPs and the intersection point
231  if(tj1Short || tj2Short) { sepCut = tcc.vtx2DCuts[1]; } else { sepCut = tcc.vtx2DCuts[2]; }
232  // NewVtxCuts: require close separation on the first pass
233  if(firstPassCuts) sepCut = tcc.vtx2DCuts[1];
234  Point2_t vPos {{wint, tint}};
235  float vt1Sep = PosSep(vPos, tp1.Pos);
236  float vt2Sep = PosSep(vPos, tp2.Pos);
237  float dwc1 = DeadWireCount(slc, wint, tp1.Pos[0], tp1.CTP);
238  float dwc2 = DeadWireCount(slc, wint, tp2.Pos[0], tp1.CTP);
239  vt1Sep -= dwc1;
240  vt2Sep -= dwc2;
241  bool vtxOnDeadWire = (DeadWireCount(slc, wint, wint, tp1.CTP) == 1);
242  if(prt && vt1Sep < 200 && vt2Sep < 200) {
243  mf::LogVerbatim myprt("TC");
244  myprt<<"F2DV candidate T"<<tj1.ID<<"_"<<end1<<"-T"<<tj2.ID<<"_"<<end2;
245  myprt<<" vtx pos "<<(int)wint<<":"<<(int)(tint/tcc.unitsPerTick)<<" tp1 "<<PrintPos(slc, tp1)<<" tp2 "<<PrintPos(slc, tp2);
246  myprt<<" dwc1 "<<dwc1<<" dwc2 "<<dwc2<<" on dead wire? "<<vtxOnDeadWire;
247  myprt<<" vt1Sep "<<vt1Sep<<" vt2Sep "<<vt2Sep<<" sepCut "<<sepCut;
248  }
249  if(vt1Sep > sepCut || vt2Sep > sepCut) continue;
250  // make sure that the other end isn't closer
251  if(PosSep(vPos, slc.tjs[it1].Pts[oendPt1].Pos) < vt1Sep) {
252  if(prt) mf::LogVerbatim("TC")<<" tj1 other end "<<PrintPos(slc, tj1.Pts[oendPt1])<<" is closer to the vertex";
253  continue;
254  }
255  if(PosSep(vPos, slc.tjs[it2].Pts[oendPt2].Pos) < vt2Sep) {
256  if(prt) mf::LogVerbatim("TC")<<" tj2 other end "<<PrintPos(slc, tj2.Pts[oendPt2])<<" is closer to the vertex";
257  continue;
258  }
259  // Ensure that the vertex position is close to the end of each Tj
260  unsigned short closePt1;
261  float doca1 = sepCut;
262  if(!TrajClosestApproach(tj1, wint, tint, closePt1, doca1)) continue;
263  // dpt1 (and dpt2) will be 0 if the vertex is at the end
264  short stepDir = -1;
265  if(tcc.modes[kStepDir]) stepDir = 1;
266  short dpt1 = stepDir * (closePt1 - endPt1);
267  if(prt) mf::LogVerbatim("TC")<<" endPt1 "<<endPt1<<" closePt1 "<<closePt1<<" dpt1 "<<dpt1<<" doca1 "<<doca1;
268  // BB April 19, 2018: require vertex to be near the end
269  if(dpt1 < -1) continue;
270  if(slc.tjs[it1].EndPt[1] > 4) {
271  if(dpt1 > 3) continue;
272  } else {
273  // tighter cut for short trajectories
274  if(dpt1 > 2) continue;
275  }
276  unsigned short closePt2;
277  float doca2 = sepCut;
278  if(!TrajClosestApproach(tj2, wint, tint, closePt2, doca2)) continue;
279  short dpt2 = stepDir * (closePt2 - endPt2);
280  if(prt) mf::LogVerbatim("TC")<<" endPt2 "<<endPt2<<" closePt2 "<<closePt2<<" dpt2 "<<dpt2<<" doca2 "<<doca2;
281  // BB April 19, 2018: require vertex to be near the end
282  if(dpt2 < -1) continue;
283  if(slc.tjs[it2].EndPt[1] > 4) {
284  if(dpt2 > 3) continue;
285  } else {
286  // tighter cut for short trajectories
287  if(dpt2 > 2) continue;
288  }
289  if(prt) mf::LogVerbatim("TC")<<" wint:tint "<<(int)wint<<":"<<(int)(tint/tcc.unitsPerTick);
290  bool fixVxPos = false;
291  if(requireVtxTjChg) {
292  // ensure that there is a signal between these TPs and the vertex on most of the wires
293  bool signalBetween = true;
294  short dpt = abs(wint - tp1.Pos[0]);
295  if(dpt > 2 && !SignalBetween(slc, tp1, wint, tcc.vtx2DCuts[6])) {
296  if(prt) mf::LogVerbatim("TC")<<" Fails SignalBetween for tp1 "<<dpt;
297  signalBetween = false;
298  }
299  dpt = abs(wint - tp2.Pos[0]);
300  if(dpt > 2 && !SignalBetween(slc, tp2, wint, tcc.vtx2DCuts[6])) {
301  if(prt) mf::LogVerbatim("TC")<<" Fails SignalBetween for tp2 "<<dpt;
302  signalBetween = false;
303  }
304  // consider the case where the intersection point is wrong because the
305  // end TP angles are screwed up but the Tjs are close to each other near the end
306  if(!signalBetween) {
307  unsigned short ipt1, ipt2;
308  float maxSep = 3;
309  bool isClose = TrajTrajDOCA(slc, tj1, tj2, ipt1, ipt2, maxSep, false);
310  // require that they are close at the correct end
311  if(isClose) isClose = (abs(ipt1 - endPt1) < 4 && abs(ipt2 - endPt2) < 4);
312  if(isClose) {
313  if(prt) mf::LogVerbatim("TC")<<" TrajTrajDOCA are close with minSep "<<maxSep<<" near "<<PrintPos(slc, tj1.Pts[ipt1].Pos)<<" "<<PrintPos(slc, tj2.Pts[ipt2].Pos);
314  // put the vertex at the TP that is closest to the intersection point
315  Point2_t vpos = {{wint, tint}};
316  if(PosSep2(tp1.Pos, vpos) < PosSep2(tp2.Pos, vpos)) {
317  wint = tp1.Pos[0];
318  tint = tp1.Pos[1];
319  } else {
320  wint = tp2.Pos[0];
321  tint = tp2.Pos[1];
322  }
323  fixVxPos = true;
324  if(prt) mf::LogVerbatim("TC")<<" new wint:tint "<<(int)wint<<":"<<(int)(tint/tcc.unitsPerTick);
325  } else {
326  // closest approach > 3
327  continue;
328  }
329  } // no signal between
330  } // requireVtxTjChg
331  // make a new temporary vertex
332  VtxStore aVtx;
333  aVtx.Pos[0] = wint;
334  aVtx.Pos[1] = tint;
335  aVtx.NTraj = 0;
336  aVtx.Pass = tj1.Pass;
337  // Topo 0 has this topology (<) and Topo 2 has this (>)
338  aVtx.Topo = 2 * end1;
339  aVtx.ChiDOF = 0;
340  aVtx.CTP = inCTP;
341  aVtx.Stat[kOnDeadWire] = vtxOnDeadWire;
342  // fix the vertex position if we needed to move it significantly, or if it is on a dead wire
343  aVtx.Stat[kFixed] = fixVxPos;
344  aVtx.Stat[kVtxIndPlnNoChg] = !requireVtxTjChg;
345  // try to fit it. We need to give it an ID to do that. Take the next
346  // available ID
347  unsigned short newVtxID = slc.vtxs.size() + 1;
348  aVtx.ID = newVtxID;
349  tj1.VtxID[end1] = newVtxID;
350  tj2.VtxID[end2] = newVtxID;
351  if(!FitVertex(slc, aVtx, prt)) {
352  tj1.VtxID[end1] = 0;
353  tj2.VtxID[end2] = 0;
354  continue;
355  }
356  // check proximity to nearby vertices
357  unsigned short mergeMeWithVx = IsCloseToVertex(slc, aVtx);
358  if(mergeMeWithVx > 0 && MergeWithVertex(slc, aVtx, mergeMeWithVx)) {
359  if(prt) mf::LogVerbatim("TC")<<" Merged with close vertex "<<mergeMeWithVx;
360  continue;
361  }
362  // Save it
363  if(!StoreVertex(slc, aVtx)) continue;
364  if(prt) {
365  mf::LogVerbatim myprt("TC");
366  myprt<<" New vtx 2V"<<aVtx.ID;
367  myprt<<" Tjs "<<tj1.ID<<"_"<<end1<<"-"<<tj2.ID<<"_"<<end2;
368  myprt<<" at "<<std::fixed<<std::setprecision(1)<<aVtx.Pos[0]<<":"<<aVtx.Pos[1]/tcc.unitsPerTick;
369  }
370  AttachAnyTrajToVertex(slc, slc.vtxs.size() - 1, prt);
371 /*
372  if(tcc.useAlg[kNewVtxCuts]) {
373  // try stepping away from the vertex
374  auto& newVx2 = slc.vtxs[slc.vtxs.size() - 1];
375  auto tjlist = GetAssns(slc, "2V", newVx2.ID, "T");
376  // count the number of long straight Tjs
377  unsigned short nlong = 0;
378  for(auto tjid : tjlist) {
379  auto& vtj = slc.tjs[tjid - 1];
380  unsigned short npts = vtj.EndPt[1] - vtj.EndPt[0] + 1;
381  if(npts > 20 && vtj.MCSMom > 100) ++nlong;
382  } // tjid
383  if(nlong > 1) {
384  std::cout<<"Try 2V"<<newVx2.ID<<" nlong "<<nlong<<"\n";
385  FindVtxTjs(slc, newVx2);
386  }
387  } // new cuts
388 */
389  SetVx2Score(slc);
390  } // it2
391  } // end1
392  } // it1
393 
394  // check the consistency of the Tjs for the newly added vertices
395  ChkVxTjs(slc, inCTP, prt);
396 
397  // only call these on the last pass
398  if(pass == USHRT_MAX) {
399  FindHammerVertices(slc, inCTP);
400  FindHammerVertices2(slc, inCTP);
401  }
402 
403  if(prt) PrintAllTraj("F2DVo", slc, USHRT_MAX, USHRT_MAX);
404 
405  } // Find2DVertices
406 
407 
409  void FindVtxTjs(TCSlice& slc, VtxStore& vx2)
410  {
411  // Look for available hits in the vicinity of this vertex and try to make
412  // a vertex trajectory from them
413 
414  if(!tcc.useAlg[kVtxTj]) return;
415 
416  if(vx2.Stat[kVtxTrjTried]) return;
417  // ignore low score
418  if(vx2.Score < tcc.vtx2DCuts[7]) return;
419 
420  bool prt = (tcc.dbg2V || tcc.dbgAlg[kVtxTj]);
421 
422  std::array<int, 2> wireWindow;
423  std::array<float, 2> timeWindow;
424 
425  // on the first try we look for small angle trajectories which will have hits
426  // with a large wire window and a small time window
427  // Vertex2DCuts fcl input usage
428  // 0 User definition of a short Tj => max number of Tj points
429  // 1 max separation between a vertex and the start of a trajectory for a user-defined short Tj
430  // 2 max separation for a user-defined long Tj
431  // 3 max position pull when attaching a Tj to a vertex
432  // 4 max position error for creating a Tj or attaching Tjs to an existing vertex
433  // 5 Min MCSMom of Tjs that can be used to create a vertex
434  // 6 min frac of Points/Wire between a vtx and a Tj. Ideally one if the efficiency is good
435  // 7 min Score
436  // 8 ID of a vertex for printing special debugging information
437  wireWindow[0] = std::nearbyint(vx2.Pos[0] - tcc.vtx2DCuts[1]);
438  wireWindow[1] = std::nearbyint(vx2.Pos[0] + tcc.vtx2DCuts[1]);
439  timeWindow[0] = vx2.Pos[1] - tcc.vtx2DCuts[1];
440  timeWindow[1] = vx2.Pos[1] + tcc.vtx2DCuts[1];
441 
442  geo::PlaneID planeID = DecodeCTP(vx2.CTP);
443  unsigned short ipl = planeID.Plane;
444 
445  if(prt) mf::LogVerbatim("TC")<<"inside FindVtxTjs 2v"<<vx2.ID<<" Window "<<wireWindow[0]<<" "<<wireWindow[1]<<" "<<timeWindow[0]/tcc.unitsPerTick<<" "<<timeWindow[1]/tcc.unitsPerTick<<" in plane "<<ipl;
446 
447  // find nearby available hits
448  bool hitsNear;
449  std::vector<unsigned int> closeHits = FindCloseHits(slc, wireWindow, timeWindow, ipl, kUnusedHits, true, hitsNear);
450  if(prt) {
451  mf::LogVerbatim myprt("TC");
452  myprt<<"closeHits";
453  for(auto& iht : closeHits) myprt<<" "<<PrintHit(slc.slHits[iht]);
454  }
455  if(closeHits.empty()) return;
456  // sort by distance from the vertex
457  std::vector<SortEntry> sortVec(closeHits.size());
458  SortEntry sortEntry;
459  for(unsigned short ii = 0; ii < closeHits.size(); ++ii) {
460  unsigned int iht = closeHits[ii];
461  auto& hit = (*evt.allHits)[slc.slHits[iht].allHitsIndex];
462  float dw = hit.WireID().Wire - vx2.Pos[0];
463  float dt = tcc.unitsPerTick * hit.PeakTime() - vx2.Pos[1];
464  float d2 = dw * dw + dt * dt;
465  sortEntry.index = ii;
466  sortEntry.val = d2;
467  sortVec[ii] = sortEntry;
468  } // ii
469  std::sort(sortVec.begin(), sortVec.end(), valDecreasing);
470  unsigned int vWire = std::nearbyint(vx2.Pos[0]);
471  int vTick = vx2.Pos[1]/tcc.unitsPerTick;
472  if(prt) PrintHeader("FVT");
473  for(unsigned short ii = 0; ii < closeHits.size(); ++ii) {
474  unsigned int iht = closeHits[sortVec[ii].index];
475  if(slc.slHits[iht].InTraj > 0) continue;
476  auto& hit = (*evt.allHits)[slc.slHits[iht].allHitsIndex];
477  // the direction will be poorly defined if a hit is very close to the vertex and it is in this list.
478  // Ignore these hits
479  if(hit.WireID().Wire == vWire) {
480  // on the vertex wire. Check for a close time
481  if(abs(hit.PeakTime() - vTick) < 10) continue;
482  } // hit on vtx wire
483  float toWire = hit.WireID().Wire;
484  float toTick = hit.PeakTime();
485  // assume the last pass and fix it later after the angle is calculated
486  unsigned short pass = tcc.minPts.size() - 1;
487  Trajectory tj;
488  if(!StartTraj(slc, tj, vx2.Pos[0], vx2.Pos[1]/tcc.unitsPerTick, toWire, toTick, vx2.CTP, pass)) continue;
489  // ensure that the first TP is good
490  if(tj.Pts[0].Pos[0] < 0) continue;
491  tj.VtxID[0] = vx2.ID;
492  TrajPoint tp = tj.Pts[0];
493  // Move the Pt to the hit
494  MoveTPToWire(tp, toWire);
495  // attach the hit
496  tp.Hits.push_back(iht);
497  tp.UseHit[tp.Hits.size()-1] = true;
498  slc.slHits[iht].InTraj = tj.ID;
499  tp.UseHit[tp.Hits.size()-1] = false;
500  if(prt) PrintTrajPoint("FVT", slc, 0, tj.StepDir, tj.Pass, tp);
501  // Step away and see what happens
502  StepAway(slc, tj);
503  // check for a major failure
504  if(!slc.isValid) return;
505  // Check the quality of the trajectory
506  CheckTraj(slc, tj);
507  if(!tj.IsGood || NumPtsWithCharge(slc, tj, true) < 2) {
508  if(prt) mf::LogVerbatim("TC")<<" xxxxxxx Not enough points "<<NumPtsWithCharge(slc, tj, true)<<" minimum "<<tcc.minPts[tj.Pass]<<" or !IsGood";
509  ReleaseHits(slc, tj);
510  continue;
511  }
512  tj.AlgMod[kVtxTj] = true;
513  slc.isValid = StoreTraj(slc, tj);
514  if(tcc.useAlg[kChkInTraj]) {
515  slc.isValid = InTrajOK(slc, "FVT");
516  if(!slc.isValid) {
517  mf::LogVerbatim("TC")<<"InTrajOK failed in FindVtxTjs";
518  return;
519  }
520  }
521  if(prt) mf::LogVerbatim("TC")<<"FindVtxTjs: calling StoreTraj with npts "<<tj.EndPt[1];
522  } // ii
523 
524  // Flag this as tried so we don't try again
525  vx2.Stat[kVtxTrjTried] = true;
526  } // FindVtxTjs
527 
529  void FindNeutralVertices(TCSlice& slc)
530  {
531  // Look for 2D neutral vertices between Tjs
532  if(!tcc.useAlg[kVxNeutral]) return;
533  if(tcc.neutralVxCuts.size() < 4) return;
534  if(slc.nPlanes < 3) return;
535  if(slc.pfps.size() < 2) return;
536 
537  /* TODO This algorithm needs work.
538 
539  bool prt = tcc.dbgVxNeutral && tcc.dbgSlc;
540 
541  struct CandVx {
542  unsigned short ip1;
543  unsigned short end1;
544  unsigned short ip2;
545  unsigned short end2;
546  Point3_t intersect;
547  float sepSum;
548  bool isValid;
549  };
550  std::vector<CandVx> candVxs;
551 
552  for(unsigned short ip1 = 0; ip1 < slc.pfps.size() - 1; ++ip1) {
553  auto& p1 = slc.pfps[ip1];
554  if(p1.ID == 0) continue;
555  if(p1.Tp3s.empty()) continue;
556  float len1 = PosSep(p1.XYZ[0], p1.XYZ[1]);
557  if(len1 < tcc.neutralVxCuts[3]) continue;
558  for(unsigned short end1 = 0; end1 < 2; ++end1) {
559  float cfne1 = ChgFracNearEnd(slc, p1, end1);
560  if(cfne1 < tcc.neutralVxCuts[0]) continue;
561  for(unsigned short ip2 = ip1 + 1; ip2 < slc.pfps.size(); ++ip2) {
562  auto& p2 = slc.pfps[ip2];
563  if(p2.ID == 0) continue;
564  if(p2.Tp3s.empty()) continue;
565  float len2 = PosSep(p2.XYZ[0], p2.XYZ[1]);
566  if(len2 < tcc.neutralVxCuts[3]) continue;
567  for(unsigned short end2 = 0; end2 < 2; ++end2) {
568  float cfne2 = ChgFracNearEnd(slc, p2, end2);
569  if(cfne2 < tcc.neutralVxCuts[0]) continue;
570  float vxDOCA = 1E6;
571  Point3_t intersect;
572  if(!PointDirIntersect(p1.XYZ[end1], p1.Dir[end1], p2.XYZ[end2], p2.Dir[end2], intersect, vxDOCA)) continue;
573  if(intersect[0] < slc.xLo || intersect[0] > slc.xHi) continue;
574  if(intersect[1] < slc.yLo || intersect[1] > slc.yHi) continue;
575  if(intersect[2] < slc.zLo || intersect[2] > slc.zHi) continue;
576  // ensure that the pfp end and the vertex are consistent
577  float sep1 = PosSep(intersect, p1.XYZ[end1]);
578  if(PosSep(intersect, p1.XYZ[1-end1]) < sep1) continue;
579  float sep2 = PosSep(intersect, p2.XYZ[end2]);
580  if(PosSep(intersect, p2.XYZ[1-end2]) < sep2) continue;
581  if(vxDOCA > tcc.neutralVxCuts[1]) continue;
582  // ensure that there isn't a lot of charge between the end of each pfp and the intersection point
583  float cfb1 = ChgFracBetween(slc, p1.XYZ[end1], intersect);
584  float cfb2 = ChgFracBetween(slc, p2.XYZ[end2], intersect);
585  if(prt) {
586  mf::LogVerbatim myprt("TC");
587  myprt<<"FNV: P"<<p1.ID<<"_"<<end1<<" sep1 "<<std::fixed<<std::setprecision(2)<<sep1;
588  myprt<<" cfne1 "<<cfne1<<" cfb1 "<<cfb1;
589  myprt<<" P"<<p2.ID<<"_"<<end2<<" sep2 "<<sep2<<" cfne2 "<<cfne2<<" cfb2 "<<cfb2;
590  myprt<<" intersect "<<intersect[0]<<" "<<intersect[1]<<" "<<intersect[2];
591  myprt<<" vxDOCA "<<vxDOCA;
592  } // prt
593  if(cfb1 > tcc.neutralVxCuts[2]) continue;
594  if(cfb2 > tcc.neutralVxCuts[2]) continue;
595  // check existing candidates
596  float sepSum = sep1 + sep2;
597  bool skipit = false;
598  for(auto& candVx : candVxs) {
599  if(!candVx.isValid) continue;
600  if(candVx.ip1 != ip1 && candVx.ip2 != ip2) continue;
601  // see if the separation sum is smaller
602  if(sepSum < candVx.sepSum) {
603  // flag the saved one as not valid
604  candVx.isValid = false;
605  } else {
606  skipit = true;
607  break;
608  }
609  } // candVx
610  if(skipit) continue;
611  CandVx candVx;
612  candVx.ip1 = ip1;
613  candVx.end1 = end1;
614  candVx.ip2 = ip2;
615  candVx.end2 = end2;
616  candVx.intersect = intersect;
617  candVx.sepSum = sep1 + sep2;
618  candVx.isValid = true;
619  candVxs.push_back(candVx);
620  if(prt) mf::LogVerbatim("TC")<<" candidate";
621  } // end2
622  } // ip2
623  } // end1
624  } // ip1
625 
626  if(candVxs.empty()) return;
627 
628  // Make vertices with the valid candidates
629  for(auto& candVx : candVxs) {
630  if(!candVx.isValid) continue;
631  Vtx3Store vx3;
632  auto& p1 = slc.pfps[candVx.ip1];
633  auto& p2 = slc.pfps[candVx.ip2];
634  vx3.TPCID = p1.TPCID;
635  vx3.Vx2ID.resize(slc.nPlanes);
636  vx3.ID = slc.vtx3s.size() + 1;
637  // try to improve the vertex position by fitting the Tjs in 2D
638  std::vector<TrajPoint> vxTps;
639  Vector3_t dir = {{0, 0, 0}};
640  for(unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
641  vxTps.clear();
642  CTP_t inCTP = EncodeCTP(p1.TPCID.Cryostat, p1.TPCID.TPC, plane);
643  // Need the position in this plane to find the closest end
644  auto vtp = MakeBareTP(slc, candVx.intersect, dir, inCTP);
645  if(vtp.Pos[0] < 0) continue;
646  for(auto tid : p1.TjIDs) {
647  auto& tj = slc.tjs[tid - 1];
648  if(tj.CTP != vtp.CTP) continue;
649  unsigned short end = CloseEnd(slc, tj, vtp.Pos);
650  auto vxTp = tj.Pts[tj.EndPt[end]];
651  vxTp.Step = tid;
652  vxTp.AngleCode = end;
653  vxTps.push_back(vxTp);
654  } // tid
655  for(auto tid : p2.TjIDs) {
656  auto& tj = slc.tjs[tid - 1];
657  if(tj.CTP != vtp.CTP) continue;
658  unsigned short end = CloseEnd(slc, tj, vtp.Pos);
659  auto vxTp = tj.Pts[tj.EndPt[end]];
660  vxTp.Step = tid;
661  vxTp.AngleCode = end;
662  vxTps.push_back(vxTp);
663  } // tid
664  if(vxTps.size() < 2) continue;
665  VtxStore vx2;
666  vx2.CTP = inCTP;
667  vx2.Topo = 11;
668  if(!FitVertex(slc, vx2, vxTps, prt)) {
669  if(prt) mf::LogVerbatim("TC")<<"FNV: vertex fit failed";
670  continue;
671  } // fit failed
672  vx2.ID = slc.vtxs.size() + 1;
673  for(auto& vxTp : vxTps) {
674  int tid = vxTp.Step;
675  auto& tj = slc.tjs[tid - 1];
676  tj.VtxID[vxTp.AngleCode] = vx2.ID;
677  }
678  if(prt) mf::LogVerbatim("TC")<<"FNV: 2V fit inCTP "<<inCTP<<" pos "<<PrintPos(slc, vx2.Pos);
679  vx2.Vx3ID = vx3.ID;
680  vx3.Vx2ID[plane] = vx2.ID;
681  SetVx2Score(slc, vx2);
682  if(!StoreVertex(slc, vx2)) {
683  if(prt) mf::LogVerbatim("TC")<<"FNV: store vertex failed";
684  continue;
685  } // StoreVertex failed
686  } // plane
687  vx3.X = candVx.intersect[0];
688  vx3.Y = candVx.intersect[1];
689  vx3.Z = candVx.intersect[2];
690  vx3.Primary = true;
691  SetVx3Score(slc, vx3);
692  ++evt.globalS3ID;
693  vx3.UID = evt.globalS3ID;
694  slc.vtx3s.push_back(vx3);
695  if(prt) mf::LogVerbatim("TC")<<"FNV: P"<<p1.ID<<"_"<<candVx.end1<<" P"<<p2.ID<<"_"<<candVx.end2<<" -> 3V"<<vx3.ID;
696  } // candVx
697 */
698  } // FindNeutralVertices
699 
701  bool MergeWithVertex(TCSlice& slc, VtxStore& vx, unsigned short oVxID)
702  {
703  // Attempts to merge the trajectories attached to vx with an existing 2D vertex
704  // referenced by existingVxID. This function doesn't use the existing end0/end1 vertex association.
705  // It returns true if the merging was successful in which case the calling function should
706  // not store vx. The calling function needs to have set VtxID to vx.ID for tjs that are currently attached
707  // to vx. It assumed that vx hasn't yet been pushed onto slc.vtxs
708 
709  if(!tcc.useAlg[kVxMerge]) return false;
710 
711  bool prt = tcc.dbgVxMerge && tcc.dbgSlc;
712 
713  if(oVxID > slc.vtxs.size()) return false;
714  auto& oVx = slc.vtxs[oVxID - 1];
715  if(vx.CTP != oVx.CTP) return false;
716 
717  // get a list of tjs attached to both vertices
718  std::vector<int> tjlist = GetVtxTjIDs(slc, vx);
719 // auto tjlist = GetAssns(slc, "2V", vx.ID, "T");
720  if(tjlist.empty()) return false;
721  std::vector<int> tmp = GetVtxTjIDs(slc, oVx);
722 // auto tmp = GetAssns(slc, "2V", oVx.ID, "T");
723  if(tmp.empty()) return false;
724  for(auto tjid : tmp) {
725  if(std::find(tjlist.begin(), tjlist.end(), tjid) == tjlist.end()) tjlist.push_back(tjid);
726  } // tjid
727  if(tjlist.size() < 2) return false;
728  // handle the simple case
729  if(tjlist.size() == 2) {
730  // Unset the fixed bit
731  vx.Stat[kFixed] = false;
732  oVx.Stat[kFixed] = false;
733  // assign the vx tjs to oVx
734  for(auto tjid : tjlist) {
735  auto& tj = slc.tjs[tjid - 1];
736  for(unsigned short end = 0; end < 2; ++end) {
737  if(tj.VtxID[end] == vx.ID) tj.VtxID[end] = oVx.ID;
738  } // end
739  } // tjid
740  if(!FitVertex(slc, oVx, prt)) {
741  if(prt) mf::LogVerbatim("TC")<<"MWV: merge failed "<<vx.ID<<" and existing "<<oVx.ID;
742  return false;
743  }
744  return true;
745  } // size = 2
746 
747  // sort by decreasing length
748  std::vector<SortEntry> sortVec(tjlist.size());
749  for(unsigned int indx = 0; indx < sortVec.size(); ++indx) {
750  sortVec[indx].index = indx;
751  auto& tj = slc.tjs[tjlist[indx] - 1];
752  sortVec[indx].val = tj.Pts.size();
753  } // indx
754  std::sort(sortVec.begin(), sortVec.end(), valDecreasing);
755  // re-order the list of Tjs
756  auto ttl = tjlist;
757  for(unsigned short ii = 0; ii < sortVec.size(); ++ii) tjlist[ii] = ttl[sortVec[ii].index];
758  // Create a local vertex using the two longest slc, then add the shorter ones
759  // until the pull reaches the cut
760  VtxStore aVx;
761  aVx.CTP = vx.CTP;
762  std::vector<TrajPoint> tjpts(tjlist.size());
763  // determine which point on each Tj that will be used in the vertex fit and stash it in
764  // the traj point Step variable. This requires knowing the real position of the merged vertex
765  // which we estimate by averaging
766  std::array<float, 2> vpos;
767  vpos[0] = 0.5 * (vx.Pos[0] + oVx.Pos[0]);
768  vpos[1] = 0.5 * (vx.Pos[1] + oVx.Pos[1]);
769  for(unsigned short ii = 0; ii < tjpts.size(); ++ii) {
770  auto& tj = slc.tjs[tjlist[ii] - 1];
771  unsigned short npwc = NumPtsWithCharge(slc, tj, false);
772  unsigned short end = CloseEnd(slc, tj, vpos);
773  // assume that we will use the end point of the tj
774  unsigned short endPt = tj.EndPt[end];
775  if(npwc > 6 && tj.Pts[endPt].NTPsFit < 4) {
776  if(end == 0) {
777  endPt += 3;
778  } else {
779  endPt -= 3;
780  }
781  endPt = NearestPtWithChg(slc, tj, endPt);
782  } // few points fit at the end
783  if(endPt < tj.EndPt[0]) endPt = tj.EndPt[0];
784  if(endPt > tj.EndPt[1]) endPt = tj.EndPt[1];
785  // define tjpts
786  tjpts[ii].CTP = tj.CTP;
787  tjpts[ii].Pos = tj.Pts[endPt].Pos;
788  tjpts[ii].Dir = tj.Pts[endPt].Dir;
789  tjpts[ii].Ang = tj.Pts[endPt].Ang;
790  tjpts[ii].AngErr = tj.Pts[endPt].AngErr;
791  // stash the point in Step
792  tjpts[ii].Step = endPt;
793  // and the end in AngleCode
794  tjpts[ii].AngleCode = end;
795  // stash the ID in Hits
796  tjpts[ii].Hits.resize(1, tj.ID);
797  } // tjid
798  if(prt) {
799  mf::LogVerbatim myprt("TC");
800  myprt<<"MWV: "<<oVxID;
801  myprt<<" Fit TPs";
802  for(unsigned short ii = 0; ii < tjpts.size(); ++ii) {
803  auto& tjpt = tjpts[ii];
804  myprt<<" "<<tjlist[ii]<<"_"<<tjpt.Step<<"_"<<PrintPos(slc, tjpt.Pos);
805  }
806  } // prt
807  // create a subset of the first two for the first fit
808  auto fitpts = tjpts;
809  fitpts.resize(2);
810  if(!FitVertex(slc, aVx, fitpts, prt)) {
811  if(prt) mf::LogVerbatim("TC")<<"MWV: first fit failed ";
812  return false;
813  }
814  // Fit and add tjs to the vertex
815  bool needsUpdate = false;
816  for(unsigned short ii = 2; ii < tjlist.size(); ++ii) {
817  fitpts.push_back(tjpts[ii]);
818  if(FitVertex(slc, aVx, fitpts, prt)) {
819  needsUpdate = false;
820  } else {
821  // remove the last Tj point and keep going
822  fitpts.pop_back();
823  needsUpdate = true;
824  }
825  } // ii
826 
827  if(needsUpdate) FitVertex(slc, aVx, fitpts, prt);
828  if(prt) mf::LogVerbatim("TC")<<"MWV: done "<<vx.ID<<" and existing "<<oVx.ID;
829 
830  // update. Remove old associations
831  for(auto& tj : slc.tjs) {
832  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
833  if(tj.CTP != vx.CTP) continue;
834  for(unsigned short end = 0; end < 2; ++end) {
835  if(tj.VtxID[end] == vx.ID) tj.VtxID[end] = 0;
836  if(tj.VtxID[end] == oVxID) tj.VtxID[end] = 0;
837  }
838  } // tj
839  // set the new associations
840  for(unsigned short ii = 0; ii < fitpts.size(); ++ii) {
841  auto& tjpt = fitpts[ii];
842  unsigned short end = tjpt.AngleCode;
843  auto& tj = slc.tjs[tjpt.Hits[0] - 1];
844  if(tj.VtxID[end] != 0) {
845  std::cout<<"MWV: coding error. tj "<<tj.ID<<" end "<<end<<" VtxID "<<tj.VtxID[end]<<" != 0\n";
846  return false;
847  }
848  tj.VtxID[end] = oVxID;
849  } // ii
850 
851  // Update oVx
852  oVx.Pos = aVx.Pos;
853  oVx.PosErr = aVx.PosErr;
854  oVx.ChiDOF = aVx.ChiDOF;
855  oVx.NTraj = fitpts.size();
856  // Update the score and the charge fraction
857  SetVx2Score(slc, oVx);
858  oVx.Stat[kVtxMerged] = true;
859  oVx.Stat[kFixed] = false;
860  if(prt) {
861  mf::LogVerbatim myprt("TC");
862  myprt<<"MWV: "<<oVxID;
863  myprt<<" Done TPs";
864  for(unsigned short ii = 0; ii < fitpts.size(); ++ii) {
865  auto& tjpt = fitpts[ii];
866  myprt<<" "<<tjpt.Hits[0]<<"_"<<tjpt.AngleCode<<"_"<<PrintPos(slc, tjpt.Pos);
867  }
868  } // prt
869 
870  return true;
871  } // MergeWithVertex
872 
874  void ChkVxTjs(TCSlice& slc, const CTP_t& inCTP, bool prt)
875  {
876  //
877 
878  if(!tcc.useAlg[kChkVxTj]) return;
879 
880  for(unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
881  auto& vx2 = slc.vtxs[ivx];
882  if(vx2.ID == 0) continue;
883  if(vx2.CTP != inCTP) continue;
884  auto vxtjs = GetAssns(slc, "2V", vx2.ID, "T");
885  if(vxtjs.size() < 2) continue;
886  if(tcc.useAlg[kNewVtxCuts] && vx2.Stat[kOnDeadWire]) continue;
887  // BB added June 4, 2018
888  // find the closest separation between the vertex and the ends of the tjs
889  float close = 200;
890  int closeID = 0;
891  unsigned short closeEnd = 0;
892  for(auto tid : vxtjs) {
893  auto& tj = slc.tjs[tid - 1];
894  unsigned short nearEnd = 1 - FarEnd(slc, tj, vx2.Pos);
895  if(tj.VtxID[nearEnd] != vx2.ID) continue;
896  float sep = PosSep(tj.Pts[tj.EndPt[nearEnd]].Pos, vx2.Pos);
897  if(sep > close) continue;
898  close = sep;
899  closeID = tid;
900  closeEnd = nearEnd;
901  } // tid
902  if(close > 1.5 && closeID > 0) {
903  auto& closeTj = slc.tjs[closeID - 1];
904  auto& closeTP = closeTj.Pts[closeTj.EndPt[closeEnd]];
905  if(prt) mf::LogVerbatim("TC")<<"CVTjs: moved 2V"<<vx2.ID<<" from "<<PrintPos(slc, vx2.Pos)<<" to "<<PrintPos(slc, closeTP);
906  vx2.Pos = closeTP.Pos;
907  vx2.Stat[kFixed] = true;
908  }
909  for(unsigned short it1 = 0; it1 < vxtjs.size() - 1; ++it1) {
910  auto& tj1 = slc.tjs[vxtjs[it1] - 1];
911  if(tj1.AlgMod[kKilled] || tj1.AlgMod[kHaloTj]) continue;
912  unsigned short end1 = 0;
913  if(tj1.VtxID[1] == vx2.ID) end1 = 1;
914  auto& vtp1 = tj1.Pts[tj1.EndPt[end1]];
915  auto& otp1 = tj1.Pts[tj1.EndPt[1 - end1]];
916  float tj1sep = PosSep(vtp1.Pos, vx2.Pos);
917  for(unsigned short it2 = it1 + 1; it2 < vxtjs.size(); ++it2) {
918  auto& tj2 = slc.tjs[vxtjs[it2] - 1];
919  if(tj2.AlgMod[kKilled] || tj2.AlgMod[kHaloTj]) continue;
920  unsigned short end2 = 0;
921  if(tj2.VtxID[2] == vx2.ID) end2 = 1;
922  auto& vtp2 = tj2.Pts[tj2.EndPt[end2]];
923  auto& otp2 = tj2.Pts[tj2.EndPt[1 - end2]];
924  float tj2sep = PosSep(vtp2.Pos, vx2.Pos);
925  float otj1tj2 = PosSep(otp1.Pos, vtp2.Pos);
926  float delta12 = PointTrajDOCA(slc, otp1.Pos[0], otp1.Pos[1], vtp2);
927  float dang12 = DeltaAngle(otp1.Ang, vtp2.Ang);
928  if(otj1tj2 < tj2sep && delta12 < 1 && otj1tj2 < 4) {
929  if(prt) {
930  mf::LogVerbatim myprt("TC");
931  myprt<<"CVTjs: "<<vx2.ID<<" tj1 "<<tj1.ID<<" tj2 "<<tj2.ID;
932  myprt<<" otj1tj2 "<<otj1tj2;
933  myprt<<" delta12 "<<delta12;
934  myprt<<" dang12 "<<dang12;
935  myprt<<" Try to merge";
936  }
937  // End 1 of tj1 is closer to end0 of tj2 than tj2 is to the vertex
938  tj2.VtxID[end2] = 0;
939  if(CompatibleMerge(slc, tj1, tj2, prt) && MergeAndStore(slc, vxtjs[it1], vxtjs[it2], prt)) {
940  auto& newTj = slc.tjs[slc.tjs.size()-1];
941  newTj.AlgMod[kChkVxTj] = true;
942  if(prt) mf::LogVerbatim("TC")<<"CVTjs: Merged tjs "<<tj1.ID<<" and "<<tj2.ID<<" -> "<<newTj.ID;
943  } else {
944  if(prt) mf::LogVerbatim("TC")<<"CVTjs: Merge failed";
945  }
946  continue;
947  } // other end is closer
948  // now check the other end of tj2
949  float tj1otj2 = PosSep(vtp1.Pos, otp2.Pos);
950  if(tj1otj2 < tj1sep && delta12 < 1 && tj1otj2 < 4) {
951  // End 1 of tj1 is closer to end0 of tj2 than tj2 is to the vertex
952  tj1.VtxID[end1] = 0;
953  if(CompatibleMerge(slc, tj2, tj1, prt) && MergeAndStore(slc, vxtjs[it2], vxtjs[it1], prt)) {
954  auto& newTj = slc.tjs[slc.tjs.size()-1];
955  newTj.AlgMod[kChkVxTj] = true;
956  if(prt) mf::LogVerbatim("TC")<<"CVTjs: Merged tjs "<<tj1.ID<<" and "<<tj2.ID<<" -> "<<newTj.ID;
957  } else {
958  if(prt) mf::LogVerbatim("TC")<<"CVTjs: Merge failed";
959  }
960  } // the other end is closer
961  } // it2
962  } // it1
963  // Check for delta-rays that have a vertex when they should have been merged
964  if(vx2.Topo == 1 && vxtjs.size() == 2) {
965  auto& tj1 = slc.tjs[vxtjs[0] - 1];
966  auto& tj2 = slc.tjs[vxtjs[1] - 1];
967  // ensure that these weren't killed above
968  if(tj1.AlgMod[kKilled] || tj2.AlgMod[kKilled]) continue;
969  if(tj1.AlgMod[kDeltaRay] || tj2.AlgMod[kDeltaRay]) {
970  if(prt) mf::LogVerbatim("TC")<<"CVTjs: Merge delta rays "<<tj1.ID<<" and "<<tj2.ID<<" CompatibleMerge? "<<CompatibleMerge(slc, tj1, tj2, prt);
971  MakeVertexObsolete("CVTjs", slc, vx2, true);
972  MergeAndStore(slc, vxtjs[0] - 1, vxtjs[1] - 1, prt);
973  } // one is a tagged delta-ray
974  } // delta-ray check
975  } // ivx
976  } // ChkVxTjs
977 
979  void FindHammerVertices2(TCSlice& slc, const CTP_t& inCTP)
980  {
981  // Variant of FindHammerVertices with slightly different requirements:
982  // 1) tj1 is a straight trajectory with most of the points fit
983  // 2) No angle requirement between tj1 and tj2
984  // 3) Large charge near the intersection point X on tj2
985  // tj2 ---X---
986  // tj1 /
987  // tj1 /
988  // tj1 /
989  // minimum^2 DOCA of tj1 endpoint with tj2
990 
991  if(!tcc.useAlg[kHamVx2]) return;
992 
993  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kHamVx2]);
994  if(prt) mf::LogVerbatim("TC")<<"Inside HamVx2";
995 
996  for(unsigned short it1 = 0; it1 < slc.tjs.size(); ++it1) {
997  if(slc.tjs[it1].CTP != inCTP) continue;
998  if(slc.tjs[it1].AlgMod[kKilled] || slc.tjs[it1].AlgMod[kHaloTj]) continue;
999  if(slc.tjs[it1].AlgMod[kHamVx]) continue;
1000  if(slc.tjs[it1].AlgMod[kHamVx2]) continue;
1001  if(slc.tjs[it1].AlgMod[kJunkTj]) continue;
1002  if(slc.tjs[it1].PDGCode == 111) continue;
1003  unsigned short numPtsWithCharge1 = NumPtsWithCharge(slc, slc.tjs[it1], false);
1004  if(numPtsWithCharge1 < 6) continue;
1005  // Check each end of tj1
1006  bool didaSplit = false;
1007  for(unsigned short end1 = 0; end1 < 2; ++end1) {
1008  // vertex assignment exists?
1009  if(slc.tjs[it1].VtxID[end1] > 0) continue;
1010  unsigned short endPt1 = slc.tjs[it1].EndPt[end1];
1011  for(unsigned short it2 = 0; it2 < slc.tjs.size(); ++it2) {
1012  if(it1 == it2) continue;
1013  if(slc.tjs[it2].AlgMod[kKilled] || slc.tjs[it2].AlgMod[kHaloTj]) continue;
1014  if(slc.tjs[it2].AlgMod[kHamVx]) continue;
1015  if(slc.tjs[it2].AlgMod[kHamVx2]) continue;
1016  // require that both be in the same CTP
1017  if(slc.tjs[it2].CTP != inCTP) continue;
1018  if(slc.tjs[it2].AlgMod[kShowerLike]) continue;
1019  if(slc.tjs[it2].AlgMod[kJunkTj]) continue;
1020  if(slc.tjs[it2].PDGCode == 111) continue;
1021  unsigned short numPtsWithCharge2 = NumPtsWithCharge(slc, slc.tjs[it2], true);
1022  if(numPtsWithCharge2 < 6) continue;
1023  // ignore muon-like tjs
1024  if(tcc.useAlg[kNewVtxCuts] && numPtsWithCharge2 > 100 && slc.tjs[it2].MCSMom > 500) continue;
1025  // Find the minimum separation between tj1 and tj2
1026  float minDOCA = 5;
1027  float doca = minDOCA;
1028  unsigned short closePt2 = 0;
1029  TrajPointTrajDOCA(slc, slc.tjs[it1].Pts[endPt1], slc.tjs[it2], closePt2, doca);
1030  if(doca == minDOCA) continue;
1031  if(prt) {
1032  mf::LogVerbatim myprt("TC");
1033  auto& tj1 = slc.tjs[it1];
1034  auto& tj2 = slc.tjs[it2];
1035  myprt<<" FHV2 CTP"<<tj1.CTP;
1036  myprt<<" t"<<tj1.ID<<"_"<<end1<<" MCSMom "<<tj1.MCSMom<<" ChgRMS "<<tj1.ChgRMS;
1037  myprt<<" split t"<<tj2.ID<<"? MCSMom "<<tj2.MCSMom<<" ChgRMS "<<tj2.ChgRMS;
1038  myprt<<" doca "<<doca<<" tj2.EndPt[0] "<<tj2.EndPt[0]<<" closePt2 "<<closePt2;
1039  myprt<<" tj2.EndPt[1] "<<tj2.EndPt[1];
1040  } // prt
1041  // ensure that the closest point is not near an end
1042  if(closePt2 < slc.tjs[it2].EndPt[0] + 3) continue;
1043  if(closePt2 > slc.tjs[it2].EndPt[1] - 3) continue;
1044  // Find the intersection point between the tj1 end and tj2 closest Pt
1045  float wint, tint;
1046  TrajIntersection(slc.tjs[it1].Pts[endPt1], slc.tjs[it2].Pts[closePt2], wint, tint);
1047  // make an angle cut
1048  float dang = DeltaAngle(slc.tjs[it1].Pts[endPt1].Ang, slc.tjs[it2].Pts[closePt2].Ang);
1049  if(dang < 0.2) continue;
1050  // ensure that tj1 doesn't cross tj2 but ensuring that the closest point on tj1 is at closePt1
1051  doca = 5;
1052  unsigned short closePt1 = 0;
1053  TrajPointTrajDOCA(slc, slc.tjs[it2].Pts[closePt2], slc.tjs[it1], closePt1, doca);
1054  if(closePt1 != endPt1) continue;
1055  if(prt) mf::LogVerbatim("TC")<<" intersection W:T "<<(int)wint<<":"<<(int)(tint/tcc.unitsPerTick)<<" dang "<<dang;
1056  // Find the point on tj2 that is closest to this point, overwriting closePt
1057  doca = minDOCA;
1058  // the point on tj2 that is closest to the intersection point
1059  unsigned short intPt2;
1060  TrajClosestApproach(slc.tjs[it2], wint, tint, intPt2, doca);
1061  if(prt) mf::LogVerbatim("TC")<<" intPt2 on tj2 "<<intPt2<<" at Pos "<<PrintPos(slc, slc.tjs[it2].Pts[intPt2])<<" doca "<<doca;
1062  if(doca == minDOCA) continue;
1063  // find the MCSMom for both sections of tj2
1064  float mcsmom = slc.tjs[it2].MCSMom;
1065  float mcsmom1 = MCSMom(slc, slc.tjs[it2], slc.tjs[it2].EndPt[0], intPt2);
1066  float mcsmom2 = MCSMom(slc, slc.tjs[it2], intPt2, slc.tjs[it2].EndPt[1]);
1067  // require that the both MCSMoms be greater than
1068  if(prt) mf::LogVerbatim("TC")<<" Check MCSMom after split: mcsmom1 "<<mcsmom1<<" mcsmom2 "<<mcsmom2;
1069  if(mcsmom1 < mcsmom || mcsmom2 < mcsmom) continue;
1070  // start scanning for the point on tj2 that has the best IP with the end of tj1 in the direction
1071  // from closePt2 -> endPt2
1072  short dir = 1;
1073  if(intPt2 < closePt2) dir = -1;
1074  unsigned short nit = 0;
1075  unsigned short ipt = intPt2;
1076  float mostChg = slc.tjs[it2].Pts[ipt].Chg;
1077  if(prt) mf::LogVerbatim("TC")<<" ipt "<<ipt<<" at Pos "<<PrintPos(slc, slc.tjs[it2].Pts[ipt])<<" chg "<<mostChg;
1078  while(nit < 20) {
1079  ipt += dir;
1080  if(ipt < 3 || ipt > slc.tjs[it2].Pts.size() - 4) break;
1081  float delta = PointTrajDOCA(slc, slc.tjs[it2].Pts[ipt].Pos[0], slc.tjs[it2].Pts[ipt].Pos[1], slc.tjs[it1].Pts[endPt1]);
1082  float sep = PosSep(slc.tjs[it2].Pts[ipt].Pos, slc.tjs[it1].Pts[endPt1].Pos);
1083  float dang = delta / sep;
1084  float pull = dang / slc.tjs[it1].Pts[endPt1].DeltaRMS;
1085  // if(prt) mf::LogVerbatim("TC")<<" intPt2 "<<intPt2<<" at Pos "<<PrintPos(slc, slc.tjs[it2].Pts[ipt])<<" delta "<<delta<<" chg "<<slc.tjs[it2].Pts[ipt].Chg<<" pull "<<pull;
1086  if(pull < 2 && slc.tjs[it2].Pts[ipt].Chg > mostChg) {
1087  mostChg = slc.tjs[it2].Pts[ipt].Chg;
1088  intPt2 = ipt;
1089  }
1090  }
1091  // require a lot of charge in tj2 in this vicinity compared with the average.
1092  float chgPull = (mostChg - slc.tjs[it2].AveChg) / slc.tjs[it2].ChgRMS;
1093  if(prt) mf::LogVerbatim("TC")<<" chgPull at intersection point "<<chgPull;
1094  if(chgPull < 10) continue;
1095  // require a signal on every wire between it1 and intPt2
1096  TrajPoint ltp = slc.tjs[it1].Pts[endPt1];
1097  if(slc.tjs[it1].Pts[endPt1].Pos[0] < -0.4) continue;
1098  unsigned int fromWire = std::nearbyint(slc.tjs[it1].Pts[endPt1].Pos[0]);
1099  if(slc.tjs[it2].Pts[intPt2].Pos[0] < -0.4) continue;
1100  unsigned int toWire = std::nearbyint(slc.tjs[it2].Pts[intPt2].Pos[0]);
1101  if(fromWire > toWire) {
1102  unsigned int tmp = fromWire;
1103  fromWire = toWire;
1104  toWire = tmp;
1105  }
1106  bool skipIt = false;
1107  for(unsigned int wire = fromWire + 1; wire < toWire; ++wire) {
1108  MoveTPToWire(ltp, (float)wire);
1109  if(!SignalAtTp(slc, ltp)) {
1110  skipIt = true;
1111  break;
1112  }
1113  } // wire
1114  if(skipIt) continue;
1115  // we have a winner
1116  // create a new vertex
1117  VtxStore aVtx;
1118  aVtx.Pos = slc.tjs[it2].Pts[intPt2].Pos;
1119  aVtx.NTraj = 3;
1120  aVtx.Pass = slc.tjs[it2].Pass;
1121  aVtx.Topo = 6;
1122  aVtx.ChiDOF = 0;
1123  aVtx.CTP = inCTP;
1124  aVtx.ID = slc.vtxs.size() + 1;
1125  unsigned short ivx = slc.vtxs.size();
1126  if(!StoreVertex(slc, aVtx)) continue;
1127  if(!SplitTraj(slc, it2, intPt2, ivx, prt)) {
1128  if(prt) mf::LogVerbatim("TC")<<"FHV2: Failed to split trajectory";
1129  // we can just remove the vertex since no Tj VtxID association has been made yet
1130  slc.vtxs.pop_back();
1131  continue;
1132  }
1133  slc.tjs[it1].VtxID[end1] = slc.vtxs[ivx].ID;
1134  slc.tjs[it1].AlgMod[kHamVx2] = true;
1135  slc.tjs[it2].AlgMod[kHamVx2] = true;
1136  unsigned short newTjIndex = slc.tjs.size() - 1;
1137  slc.tjs[newTjIndex].AlgMod[kHamVx2] = true;
1138  AttachAnyTrajToVertex(slc, ivx, prt);
1139  SetVx2Score(slc);
1140  // Update the PDGCode for the chopped trajectory
1141  SetPDGCode(slc, it2, true);
1142  // and for the new trajectory
1143  SetPDGCode(slc, newTjIndex, true);
1144  if(prt) mf::LogVerbatim("TC")<<" FHV2: New vtx 2V"<<slc.vtxs[ivx].ID<<" Score "<<slc.vtxs[ivx].Score;
1145  didaSplit = true;
1146  break;
1147  } // it2
1148  if(didaSplit) break;
1149  } // end1
1150  } // it1
1151  } // FindHammerVertices2
1152 
1154  void FindHammerVertices(TCSlice& slc, const CTP_t& inCTP)
1155  {
1156  // Look for a trajectory that intersects another. Split
1157  // the trajectory and make a vertex. The convention used
1158  // is shown pictorially here. Trajectory tj1 must be longer
1159  // than tj2
1160  // tj2 ------
1161  // tj1 /
1162  // tj1 /
1163  // tj1 /
1164 
1165  if(!tcc.useAlg[kHamVx]) return;
1166 
1167  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kHamVx]);
1168 
1169  for(unsigned short it1 = 0; it1 < slc.tjs.size(); ++it1) {
1170  if(slc.tjs[it1].CTP != inCTP) continue;
1171  if(slc.tjs[it1].AlgMod[kKilled] || slc.tjs[it1].AlgMod[kHaloTj]) continue;
1172  if(slc.tjs[it1].AlgMod[kShowerLike]) continue;
1173  if(slc.tjs[it1].AlgMod[kJunkTj]) continue;
1174  if(slc.tjs[it1].PDGCode == 111) continue;
1175  // minimum length requirements
1176  unsigned short tj1len = slc.tjs[it1].EndPt[1] - slc.tjs[it1].EndPt[0] + 1;
1177  if(tj1len < 5) continue;
1178  // Check each end of tj1
1179  bool didaSplit = false;
1180  for(unsigned short end1 = 0; end1 < 2; ++end1) {
1181  // vertex assignment exists?
1182  if(slc.tjs[it1].VtxID[end1] > 0) continue;
1183  unsigned short endPt1 = slc.tjs[it1].EndPt[end1];
1184  for(unsigned short it2 = 0; it2 < slc.tjs.size(); ++it2) {
1185  if(slc.tjs[it2].CTP != inCTP) continue;
1186  if(it1 == it2) continue;
1187  if(slc.tjs[it2].AlgMod[kKilled] || slc.tjs[it2].AlgMod[kHaloTj]) continue;
1188  if(slc.tjs[it2].AlgMod[kJunkTj]) continue;
1189  if(slc.tjs[it2].PDGCode == 111) continue;
1190  // length of tj2 cut
1191  unsigned short tj2len = slc.tjs[it2].EndPt[1] - slc.tjs[it2].EndPt[0] + 1;
1192  if(tj2len < 6) continue;
1193  // ignore if tj1 is a lot shorter than tj2
1194 // if(tj1len < 0.5 * tj2len) continue;
1195  // ignore very long straight trajectories (probably a cosmic muon)
1196  unsigned short end20 = slc.tjs[it2].EndPt[0];
1197  unsigned short end21 = slc.tjs[it2].EndPt[1];
1198  if(tcc.useAlg[kNewVtxCuts]) {
1199  if(tj2len > 200 && slc.tjs[it2].PDGCode == 13) continue;
1200  } else {
1201  if(tj2len > 100 && DeltaAngle(slc.tjs[it2].Pts[end20].Ang, slc.tjs[it2].Pts[end21].Ang) < 0.2) continue;
1202  }
1203  // Require no vertex associated with itj2
1204 // if(slc.tjs[it2].VtxID[0] > 0 || slc.tjs[it2].VtxID[1] > 0) continue;
1205  float minDOCA = 5;
1206  float doca = minDOCA;
1207  unsigned short closePt2 = 0;
1208  TrajPointTrajDOCA(slc, slc.tjs[it1].Pts[endPt1], slc.tjs[it2], closePt2, doca);
1209  if(doca == minDOCA) continue;
1210  // ensure that the closest point is not near an end
1211  if(prt) mf::LogVerbatim("TC")<<"FHV: Candidate t"<<slc.tjs[it1].ID<<" t"<<slc.tjs[it2].ID<<" doca "<<doca<<" tj2.EndPt[0] "<<slc.tjs[it2].EndPt[0]<<" closePt2 "<<closePt2<<" tj2.EndPt[1] "<<slc.tjs[it2].EndPt[1];
1212  if(closePt2 < slc.tjs[it2].EndPt[0] + 3) continue;
1213  if(closePt2 > slc.tjs[it2].EndPt[1] - 3) continue;
1214  // make an angle cut
1215  float dang = DeltaAngle(slc.tjs[it1].Pts[endPt1].Ang, slc.tjs[it2].Pts[closePt2].Ang);
1216  if(prt) mf::LogVerbatim("TC")<<" dang "<<dang<<" imposing a hard cut of 0.4 for now ";
1217  if(dang < 0.4) continue;
1218  // check the cleanliness in this area
1219  std::vector<int> tjids(2);
1220  tjids[0] = slc.tjs[it1].ID;
1221  tjids[1] = slc.tjs[it2].ID;
1222  float chgFrac = ChgFracNearPos(slc, slc.tjs[it2].Pts[closePt2].Pos, tjids);
1223  if(prt) mf::LogVerbatim("TC")<<" chgFrac "<<chgFrac;
1224  if(chgFrac < 0.9) continue;
1225  // create a new vertex
1226  VtxStore aVtx;
1227  aVtx.Pos = slc.tjs[it2].Pts[closePt2].Pos;
1228  aVtx.NTraj = 3;
1229  aVtx.Pass = slc.tjs[it2].Pass;
1230  aVtx.Topo = 5;
1231  aVtx.ChiDOF = 0;
1232  aVtx.CTP = inCTP;
1233  aVtx.ID = slc.vtxs.size() + 1;
1234  if(!StoreVertex(slc, aVtx)) continue;
1235  unsigned short ivx = slc.vtxs.size() - 1;
1236  if(!SplitTraj(slc, it2, closePt2, ivx, prt)) {
1237  if(prt) mf::LogVerbatim("TC")<<"FHV: Failed to split trajectory";
1238  slc.vtxs.pop_back();
1239  continue;
1240  }
1241  slc.tjs[it1].VtxID[end1] = slc.vtxs[ivx].ID;
1242  slc.tjs[it1].AlgMod[kHamVx] = true;
1243  slc.tjs[it2].AlgMod[kHamVx] = true;
1244  unsigned short newTjIndex = slc.tjs.size() - 1;
1245  slc.tjs[newTjIndex].AlgMod[kHamVx] = true;
1246  SetVx2Score(slc);
1247  // Update the PDGCode for the chopped trajectory
1248  SetPDGCode(slc, it2, true);
1249  // and for the new trajectory
1250  SetPDGCode(slc, newTjIndex, true);
1251  if(prt) {
1252  auto& vx2 = slc.vtxs[ivx];
1253  mf::LogVerbatim myprt("TC");
1254  myprt<<" new 2V"<<vx2.ID<<" Score "<<vx2.Score<<" Tjs";
1255  auto tjlist = GetAssns(slc, "2V", vx2.ID, "T");
1256  for(auto tid : tjlist) myprt<<" t"<<tid;
1257  }
1258  didaSplit = true;
1259  break;
1260  } // tj2
1261  if(didaSplit) break;
1262  } // end1
1263  } // tj1
1264 
1265  } // FindHammerVertices
1266 
1268  void SplitTrajCrossingVertices(TCSlice& slc, CTP_t inCTP)
1269  {
1270  // This is kind of self-explanatory...
1271 
1272  if(!tcc.useAlg[kSplitTjCVx]) return;
1273 
1274  if(slc.vtxs.empty()) return;
1275  if(slc.tjs.empty()) return;
1276 
1277  constexpr float docaCut = 4;
1278 
1279  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kSplitTjCVx]);
1280  if(prt) mf::LogVerbatim("TC")<<"Inside SplitTrajCrossingVertices inCTP "<<inCTP;
1281 
1282  geo::PlaneID planeID = DecodeCTP(inCTP);
1283 
1284  unsigned short nTraj = slc.tjs.size();
1285  for(unsigned short itj = 0; itj < nTraj; ++itj) {
1286  // NOTE: Don't use a reference variable because it may get lost if the tj is split
1287 // auto& tj = slc.tjs[itj];
1288  if(slc.tjs[itj].CTP != inCTP) continue;
1289  // obsolete trajectory
1290  if(slc.tjs[itj].AlgMod[kKilled] || slc.tjs[itj].AlgMod[kHaloTj]) continue;
1291  if(slc.tjs[itj].AlgMod[kSplitTjCVx]) continue;
1292  // too short
1293  if(slc.tjs[itj].EndPt[1] < 6) continue;
1294  for(unsigned short iv = 0; iv < slc.vtxs.size(); ++iv) {
1295  auto& vx2 = slc.vtxs[iv];
1296  // obsolete vertex
1297  if(vx2.NTraj == 0) continue;
1298  // trajectory already associated with vertex
1299  if(slc.tjs[itj].VtxID[0] == vx2.ID ||
1300  slc.tjs[itj].VtxID[1] == vx2.ID) continue;
1301  // not in the cryostat/tpc/plane
1302  if(slc.tjs[itj].CTP != vx2.CTP) continue;
1303  // poor quality
1304  if(vx2.Score < tcc.vtx2DCuts[7]) continue;
1305  float doca = docaCut;
1306  // make the cut significantly larger if the vertex is in a dead
1307  // wire gap to get the first TP that is just outside the gap.
1308  if(vx2.Stat[kOnDeadWire]) doca = 100;
1309  unsigned short closePt = 0;
1310  if(!TrajClosestApproach(slc.tjs[itj], vx2.Pos[0], vx2.Pos[1], closePt, doca)) continue;
1311  if(vx2.Stat[kOnDeadWire]) {
1312  // special handling for vertices in dead wire regions. Find the IP between
1313  // the closest point on the Tj and the vertex
1314  doca = PointTrajDOCA(slc, vx2.Pos[0], vx2.Pos[1], slc.tjs[itj].Pts[closePt]);
1315  }
1316  if(doca > docaCut) continue;
1317  if(prt) mf::LogVerbatim("TC")<<" doca "<<doca<<" btw T"<<slc.tjs[itj].ID<<" and 2V"<<slc.vtxs[iv].ID<<" closePt "<<closePt<<" in plane "<<planeID.Plane;
1318  // compare the length of the Tjs used to make the vertex with the length of the
1319  // Tj that we want to split. Don't allow a vertex using very short Tjs to split a long
1320  // Tj in the 3rd plane
1321  auto vxtjs = GetVtxTjIDs(slc, vx2);
1322  if(vxtjs.empty()) continue;
1323  unsigned short maxPts = 0;
1324  // ensure that there is a large angle between a Tj already attached to the vertex and the
1325  // tj that we want to split. We might be considering a delta-ray here
1326  float maxdang = 0;
1327  float tjAng = slc.tjs[itj].Pts[closePt].Ang;
1328  for(auto tjid : vxtjs) {
1329  auto& vtj = slc.tjs[tjid - 1];
1330  if(vtj.AlgMod[kDeltaRay]) continue;
1331  unsigned short npwc = NumPtsWithCharge(slc, vtj, false);
1332  if(npwc > maxPts) maxPts = npwc;
1333  unsigned short end = 0;
1334  if(vtj.VtxID[1] == slc.vtxs[iv].ID) end = 1;
1335  auto& vtp = vtj.Pts[vtj.EndPt[end]];
1336  float dang = DeltaAngle(vtp.Ang, tjAng);
1337  if(dang > maxdang) maxdang = dang;
1338  } // tjid
1339  // skip this operation if any of the Tjs in the split list are > 3 * maxPts
1340  maxPts *= 3;
1341  bool skipit = false;
1342  if(NumPtsWithCharge(slc, slc.tjs[itj], false) > maxPts && maxPts < 100) skipit = true;
1343  if(!skipit && maxdang < tcc.kinkCuts[0]) skipit = true;
1344  if(prt) mf::LogVerbatim("TC")<<" maxPts "<<maxPts<<" vxtjs[0] "<<vxtjs[0]<<" maxdang "<<maxdang<<" skipit? "<<skipit;
1345  if(skipit) {
1346  // kill the vertex?
1347  if(doca < 1) MakeVertexObsolete("STCV", slc, vx2, true);
1348  continue;
1349  }
1350 
1351  // make some adjustments to closePt
1352  if(vx2.Stat[kOnDeadWire]) {
1353  // ensure that the tj will be split at the gap. The closePt point may be
1354  // on the wrong side of it
1355  auto& closeTP = slc.tjs[itj].Pts[closePt];
1356  if(slc.tjs[itj].StepDir > 0 && closePt > slc.tjs[itj].EndPt[0]) {
1357  if(closeTP.Pos[0] > vx2.Pos[0]) --closePt;
1358  } else if(slc.tjs[itj].StepDir < 0 && closePt < slc.tjs[itj].EndPt[1]) {
1359  if(closeTP.Pos[0] < vx2.Pos[0]) ++closePt;
1360  }
1361  } else {
1362  // improve closePt based on vertex position
1363  // check if closePt and EndPt[1] are the two sides of vertex
1364  // take dot product of closePt-vtx and EndPt[1]-vtx
1365  if ((slc.tjs[itj].Pts[closePt].Pos[0]-vx2.Pos[0])*(slc.tjs[itj].Pts[slc.tjs[itj].EndPt[1]].Pos[0]-vx2.Pos[0]) + (slc.tjs[itj].Pts[closePt].Pos[1]-vx2.Pos[1])*(slc.tjs[itj].Pts[slc.tjs[itj].EndPt[1]].Pos[1]-vx2.Pos[1]) <0 && closePt < slc.tjs[itj].EndPt[1] - 1) ++closePt;
1366  else if ((slc.tjs[itj].Pts[closePt].Pos[0]-vx2.Pos[0])*(slc.tjs[itj].Pts[slc.tjs[itj].EndPt[0]].Pos[0]-vx2.Pos[0]) + (slc.tjs[itj].Pts[closePt].Pos[1]-vx2.Pos[1])*(slc.tjs[itj].Pts[slc.tjs[itj].EndPt[0]].Pos[1]-slc.vtxs[iv].Pos[1]) <0 && closePt > slc.tjs[itj].EndPt[0] + 1) --closePt;
1367  }
1368 
1369 
1370  if(prt) {
1371  mf::LogVerbatim("TC")<<"Good doca "<<doca<<" btw T"<<slc.tjs[itj].ID<<" and 2V"<<vx2.ID<<" closePt "<<closePt<<" in plane "<<planeID.Plane<<" CTP "<<slc.vtxs[iv].CTP;
1372  PrintTrajPoint("STCV", slc, closePt, 1, slc.tjs[itj].Pass, slc.tjs[itj].Pts[closePt]);
1373  }
1374  // ensure that the closest point is not near an end
1375  if(closePt < slc.tjs[itj].EndPt[0] + 3) continue;
1376  if(closePt > slc.tjs[itj].EndPt[1] - 3) continue;
1377  if(!SplitTraj(slc, itj, closePt, iv, prt)) {
1378  if(prt) mf::LogVerbatim("TC")<<"SplitTrajCrossingVertices: Failed to split trajectory";
1379  continue;
1380  }
1381  slc.tjs[itj].AlgMod[kSplitTjCVx] = true;
1382  unsigned short newTjIndex = slc.tjs.size() - 1;
1383  slc.tjs[newTjIndex].AlgMod[kSplitTjCVx] = true;
1384  // re-fit the vertex position
1385  FitVertex(slc, vx2, prt);
1386  } // iv
1387  } // itj
1388 
1389  } // SplitTrajCrossingVertices
1390 
1392  void Find3DVertices(TCSlice& slc)
1393  {
1394  // Create 3D vertices from 2D vertices. 3D vertices that are matched
1395  // in all three planes have Vtx2ID > 0 for all planes. This function re-scores all
1396  // 2D and 3D vertices and flags Tjs that have high-score 3D vertices
1397 
1398  if(tcc.vtx3DCuts[0] < 0) return;
1399  if(slc.vtxs.size() < 2) return;
1400  bool newCuts = (tcc.vtx3DCuts.size() > 2);
1401 
1402  // create a array/vector of 2D vertex indices in each plane
1403  std::vector<std::vector<unsigned short>> vIndex(3);
1404  for(unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
1405  // obsolete vertex
1406  if(slc.vtxs[ivx].ID == 0) continue;
1407  geo::PlaneID planeID = DecodeCTP(slc.vtxs[ivx].CTP);
1408  unsigned short plane = planeID.Plane;
1409  if(plane > 2) continue;
1410  vIndex[plane].push_back(ivx);
1411  }
1412 
1413  unsigned short vtxInPln = 0;
1414  for(unsigned short plane = 0; plane < slc.nPlanes; ++plane) if(vIndex[plane].size() > 0) ++vtxInPln;
1415  if(vtxInPln < 2) return;
1416 
1417  float thirdPlanedXCut = 2 * tcc.vtx3DCuts[0];
1418  bool prt = (tcc.dbg3V && tcc.dbgSlc);
1419  if(prt) {
1420  mf::LogVerbatim("TC")<<"Inside Find3DVertices. dX cut "<<tcc.vtx3DCuts[0]<<" thirdPlanedXCut "<<thirdPlanedXCut;
1421 // PrintAllTraj("F3DV", slc, USHRT_MAX, slc.tjs.size());
1422  }
1423 
1424  size_t vsize = slc.vtxs.size();
1425  // vector of 2D vertices -> 3D vertices.
1426  std::vector<short> vPtr(vsize, -1);
1427  // fill temp vectors of 2D vertex X and X errors
1428  std::vector<float> vX(vsize, -100);
1429 
1430  for(unsigned short ivx = 0; ivx < vsize; ++ivx) {
1431  if(slc.vtxs[ivx].ID <= 0) continue;
1432  if(tcc.useAlg[kNewVtxCuts] && slc.vtxs[ivx].Score < tcc.vtx2DCuts[7]) continue;
1433  geo::PlaneID planeID = DecodeCTP(slc.vtxs[ivx].CTP);
1434  if(slc.vtxs[ivx].Pos[0] < -0.4) continue;
1435  unsigned int wire = std::nearbyint(slc.vtxs[ivx].Pos[0]);
1436  if(!tcc.geom->HasWire(geo::WireID(planeID, wire))) continue;
1437  // Convert 2D vertex time error to X error
1438  double ticks = slc.vtxs[ivx].Pos[1] / tcc.unitsPerTick;
1439  vX[ivx] = tcc.detprop->ConvertTicksToX(ticks, planeID);
1440  } // ivx
1441 
1442  // temp vector of all 2D vertex matches
1443  std::vector<Vtx3Store> v3temp;
1444 
1445  unsigned int cstat = slc.TPCID.Cryostat;
1446  unsigned int tpc = slc.TPCID.TPC;
1447 
1448  TrajPoint tp;
1449  float maxScore = 0;
1450  constexpr float maxSep = 4;
1451  // i, j, k indicates 3 different wire planes
1452  // compare vertices in each view
1453  for(unsigned short ipl = 0; ipl < 2; ++ipl) {
1454  for(unsigned short ii = 0; ii < vIndex[ipl].size(); ++ii) {
1455  unsigned short ivx = vIndex[ipl][ii];
1456  if(vX[ivx] < 0) continue;
1457  auto& ivx2 = slc.vtxs[ivx];
1458  if(ivx2.Pos[0] < -0.4) continue;
1459  unsigned int iWire = std::nearbyint(ivx2.Pos[0]);
1460  for(unsigned short jpl = ipl + 1; jpl < 3; ++jpl) {
1461  for(unsigned short jj = 0; jj < vIndex[jpl].size(); ++jj) {
1462  unsigned short jvx = vIndex[jpl][jj];
1463  if(vX[jvx] < 0) continue;
1464  auto& jvx2 = slc.vtxs[jvx];
1465  if(jvx2.Pos[0] < -0.4) continue;
1466  unsigned int jWire = std::nearbyint(jvx2.Pos[0]);
1467  float dX = std::abs(vX[ivx] - vX[jvx]);
1468  if(dX > tcc.vtx3DCuts[0]) continue;
1469  // see if this pair is already in the list
1470  //bool gotit = false;
1471  for(auto& v3t : v3temp) {
1472  unsigned short cnt = 0;
1473  if(std::find(v3t.Vx2ID.begin(), v3t.Vx2ID.end(), ivx2.ID) != v3t.Vx2ID.end()) ++cnt;
1474  if(std::find(v3t.Vx2ID.begin(), v3t.Vx2ID.end(), jvx2.ID) != v3t.Vx2ID.end()) ++cnt;
1475  if(cnt == slc.nPlanes) {
1476  //gotit = true;
1477  break;
1478  } // cnt == slc.nPlanes
1479  } // v3t
1480  if(prt) {
1481  mf::LogVerbatim("TC")<<"F3DV: ipl "<<ipl<<" i2V"<<ivx2.ID<<" iX "<<vX[ivx]
1482  <<" jpl "<<jpl<<" j2V"<<jvx2.ID<<" jvX "<<vX[jvx]<<" W:T "<<(int)jvx2.Pos[0]<<":"<<(int)jvx2.Pos[1]<<" dX "<<dX;
1483  }
1484  double y = -1000, z = -1000;
1485  tcc.geom->IntersectionPoint(iWire, jWire, ipl, jpl, cstat, tpc, y, z);
1486  if(y < slc.yLo || y > slc.yHi || z < slc.zLo || z > slc.zHi) continue;
1487  unsigned short kpl = 3 - ipl - jpl;
1488  float kX = 0.5 * (vX[ivx] + vX[jvx]);
1489  int kWire = -1;
1490  if(slc.nPlanes > 2) {
1491  kWire = (int)(tcc.geom->WireCoordinate(y, z, kpl, tpc, cstat) + 0.5);
1492  if(kWire < 0 || (unsigned int)kWire > slc.nWires[kpl]) continue;
1493  if(!tcc.geom->HasWire(geo::WireID(cstat, tpc, kpl, kWire))) continue;
1494  std::array<int, 2> wireWindow;
1495  std::array<float, 2> timeWindow;
1496  wireWindow[0] = kWire - maxSep;
1497  wireWindow[1] = kWire + maxSep;
1498  float time = tcc.detprop->ConvertXToTicks(kX, kpl, (int)tpc, (int)cstat) * tcc.unitsPerTick;
1499  timeWindow[0] = time - maxSep;
1500  timeWindow[1] = time + maxSep;
1501  bool hitsNear;
1502  std::vector<unsigned int> closeHits = FindCloseHits(slc, wireWindow, timeWindow, kpl, kAllHits, true, hitsNear);
1503  if(prt) {
1504  mf::LogVerbatim myprt("TC");
1505  myprt<<" Hits near "<<kpl<<":"<<kWire<<":"<<(int)(time/tcc.unitsPerTick)<<" = ";
1506  for(auto iht : closeHits) myprt<<" "<<PrintHit(slc.slHits[iht]);
1507  }
1508  if(!hitsNear) continue;
1509  } // 3-plane TPC
1510  // save this incomplete 3D vertex
1511  Vtx3Store v3d;
1512  // give it a non-zero ID so that SetVx3Score returns a valid score
1513  v3d.ID = 666;
1514  v3d.Vx2ID.resize(slc.nPlanes);
1515  v3d.Vx2ID[ipl] = ivx2.ID;
1516  v3d.Vx2ID[jpl] = jvx2.ID;
1517  if(slc.nPlanes == 2) v3d.Vx2ID[2] = -1;
1518  v3d.X = kX;
1519  // Use XErr to store dX
1520  v3d.XErr = dX;
1521  v3d.Y = y;
1522  v3d.Z = z;
1523  v3d.Wire = kWire;
1524  float posError = dX / tcc.vtx3DCuts[0];
1525  float vxScoreWght = 0;
1526  if(newCuts) {
1527  SetVx3Score(slc, v3d);
1528  vxScoreWght = tcc.vtx3DCuts[2] / v3d.Score;
1529  if(posError < 0.5) posError = 0;
1530  }
1531  v3d.Score = posError + vxScoreWght;
1532  v3d.TPCID = slc.TPCID;
1533  // push the incomplete vertex onto the list
1534  v3temp.push_back(v3d);
1535 
1536  if(prt) {
1537  mf::LogVerbatim myprt("TC");
1538  myprt<<"F3DV: 2 Plane match i2V";
1539  myprt<<slc.vtxs[ivx].ID<<" P:W:T "<<ipl<<":"<<(int)slc.vtxs[ivx].Pos[0]<<":"<<(int)slc.vtxs[ivx].Pos[1];
1540  myprt<<" j2V"<<slc.vtxs[jvx].ID<<" P:W:T "<<jpl<<":"<<(int)slc.vtxs[jvx].Pos[0]<<":"<<(int)slc.vtxs[jvx].Pos[1];
1541  myprt<<std::fixed<<std::setprecision(3);
1542  myprt<<" dX "<<dX<<" posError "<<posError<<" vxScoreWght "<<vxScoreWght<<" Score "<<v3d.Score;
1543  }
1544 
1545  if(slc.nPlanes == 2) continue;
1546 
1547  // look for a 3 plane match
1548  for(unsigned short kk = 0; kk < vIndex[kpl].size(); ++kk) {
1549  unsigned short kvx = vIndex[kpl][kk];
1550  if(vX[kvx] < 0) continue;
1551  float dX = std::abs(vX[kvx] - v3d.X);
1552  // Wire difference error
1553  float dW = tcc.wirePitch * std::abs(slc.vtxs[kvx].Pos[0] - kWire);
1554  if(prt) mf::LogVerbatim("TC")<<" k2V"<<kvx+1<<" dX "<<dX<<" dW "<<dW;
1555  if(dX > thirdPlanedXCut) continue;
1556  if(dW > tcc.vtx3DCuts[1]) continue;
1557  // put the Y,Z difference in YErr and ZErr
1558  double y = -1000, z = -1000;
1559  tcc.geom->IntersectionPoint(iWire, kWire, ipl, kpl, cstat, tpc, y, z);
1560  v3d.YErr = y - v3d.Y;
1561  v3d.ZErr = z - v3d.Z;
1562  v3d.Vx2ID[kpl] = slc.vtxs[kvx].ID;
1563  v3d.Wire = -1;
1564  // hijack the Score variable to hold the separation^2, weighted by the
1565  // vertex3DCuts
1566  dX = (vX[kvx] - v3d.X) / tcc.vtx3DCuts[0];
1567  float dY = v3d.YErr / tcc.vtx3DCuts[1];
1568  float dZ = v3d.ZErr / tcc.vtx3DCuts[1];
1569  posError = dX * dX + dY * dY + dZ * dZ;
1570  vxScoreWght = 0;
1571  if(newCuts) {
1572  SetVx3Score(slc, v3d);
1573  vxScoreWght = tcc.vtx3DCuts[2] / v3d.Score;
1574  if(posError < 0.5) posError = 0;
1575  } // newCuts
1576  v3d.Score = posError + vxScoreWght;
1577  if(v3d.Score > maxScore) maxScore = v3d.Score;
1578  v3temp.push_back(v3d);
1579  } // kk
1580  } // jj
1581  } // jpl
1582  } // ii
1583  } // ipl
1584 
1585  if(v3temp.empty()) return;
1586 
1587  // We will sort this list by increasing score. First add the maxScore for 2-plane matches so that
1588  // they are considered after the 3-plane matches
1589  maxScore += 1;
1590  for(auto& v3 : v3temp) if(v3.Wire >= 0) v3.Score += maxScore;
1591 
1592  if(prt) {
1593  mf::LogVerbatim("TC")<<"v3temp list";
1594  for(auto& v3 : v3temp) {
1595  if(slc.nPlanes == 2) {
1596  mf::LogVerbatim("TC")<<"2V"<<v3.Vx2ID[0]<<" 2V"<<v3.Vx2ID[1]<<" wire "<<v3.Wire<<" Score "<<v3.Score;
1597  } else {
1598  mf::LogVerbatim("TC")<<"2V"<<v3.Vx2ID[0]<<" 2V"<<v3.Vx2ID[1]<<" 2V"<<v3.Vx2ID[2]<<" wire "<<v3.Wire<<" Score "<<v3.Score;
1599  }
1600  } // v3
1601  }
1602  SortEntry sEntry;
1603  std::vector<SortEntry> sortVec(v3temp.size());
1604  for(unsigned short ivx = 0; ivx < v3temp.size(); ++ivx) {
1605  sEntry.index = ivx;
1606  sEntry.val = v3temp[ivx].Score;
1607  sortVec[ivx] = sEntry;
1608  } // ivx
1609  if(sortVec.size() > 1) std::sort(sortVec.begin(), sortVec.end(), valIncreasing);
1610  // create a new vector of selected 3D vertices
1611  std::vector<Vtx3Store> v3sel;
1612  for(unsigned short ii = 0; ii < sortVec.size(); ++ii) {
1613  unsigned short ivx = sortVec[ii].index;
1614  // ensure that all 2D vertices are unique
1615  bool skipit = false;
1616  for(auto& v3 : v3sel) {
1617  for(unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
1618  if(v3temp[ivx].Vx2ID[ipl] == 0) continue;
1619  if(v3temp[ivx].Vx2ID[ipl] == v3.Vx2ID[ipl]) {
1620  skipit = true;
1621  break;
1622  }
1623  } // ipl
1624  if(skipit) break;
1625  } // v3
1626  if(skipit) continue;
1627  v3sel.push_back(v3temp[ivx]);
1628  } // ii
1629  v3temp.clear();
1630 
1631  if(prt) {
1632  mf::LogVerbatim myprt("TC");
1633  myprt<<"v3sel list\n";
1634  for(auto& v3d : v3sel) {
1635  for(auto vx2id : v3d.Vx2ID) if(vx2id > 0) myprt<<" 2V"<<vx2id;
1636  myprt<<" wire "<<v3d.Wire<<" Score "<<v3d.Score;
1637  myprt<<"\n";
1638  } // v3d
1639  } // prt
1640 
1641  // Count the number of incomplete vertices and store
1642  unsigned short ninc = 0;
1643  for(auto& vx3 : v3sel) {
1644  if(slc.nPlanes == 3 && vx3.Wire >= 0) ++ninc;
1645  vx3.ID = slc.vtx3s.size() + 1;
1646  ++evt.globalS3ID;
1647  vx3.UID = evt.globalS3ID;
1648  if(prt) {
1649  mf::LogVerbatim myprt("TC");
1650  myprt<<" 3V"<<vx3.ID;
1651  for(auto v2id : vx3.Vx2ID) myprt<<" 2V"<<v2id;
1652  myprt<<" wire "<<vx3.Wire;
1653  } // prt
1654  slc.vtx3s.push_back(vx3);
1655  // make the 2D -> 3D associations
1656  for(unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
1657  if(vx3.Vx2ID[ipl] == 0) continue;
1658  VtxStore& vx2 = slc.vtxs[vx3.Vx2ID[ipl]-1];
1659  vx2.Vx3ID = vx3.ID;
1660  } // ipl
1661  } // ivx
1662 
1663  // Try to complete incomplete vertices
1664  if(ninc > 0) {
1665  CompleteIncomplete3DVerticesInGaps(slc);
1666  CompleteIncomplete3DVertices(slc);
1667  }
1668 
1669  // Score and flag Tjs that are attached to high-score vertices
1670  // First remove Tj vertex flags
1671  for(auto& tj : slc.tjs) {
1672  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1673  geo::PlaneID planeID = DecodeCTP(tj.CTP);
1674  if(planeID.TPC != tpc || planeID.Cryostat != cstat) continue;
1675  tj.AlgMod[kTjHiVx3Score] = false;
1676  } // tj
1677  // Score the 2D vertices
1678  for(auto& vx2 : slc.vtxs) {
1679  if(vx2.ID == 0) continue;
1680  geo::PlaneID planeID = DecodeCTP(vx2.CTP);
1681  if(planeID.TPC != tpc || planeID.Cryostat != cstat) continue;
1682  SetVx2Score(slc, vx2);
1683  } // vx2
1684  // and the 3D vertices
1685  for(auto& vx3 : slc.vtx3s) {
1686  if(vx3.ID == 0) continue;
1687  SetVx3Score(slc, vx3);
1688  } // vx3
1689 
1690  } // Find3DVertices
1691 
1693  void Match3DVtxTjs(TCSlice& slc, bool prt)
1694  {
1695  // Matches Tjs that are attached to 2D vertices that are matched in 3D. This function does not attempt
1696  // to determine the appropriate ends of matched Tjs when there is a 3D vertex at both ends.
1697 
1698  if(!tcc.useAlg[kM3DVxTj]) return;
1699  if(slc.vtx3s.empty()) return;
1700  if(slc.matchVec.empty()) return;
1701 
1702  // sort the vertices by decreasing score
1703  std::vector<SortEntry> sortVec;
1704  for(unsigned short iv3 = 0; iv3 < slc.vtx3s.size(); ++iv3) {
1705  auto& vx3 = slc.vtx3s[iv3];
1706  if(vx3.ID == 0) continue;
1707  // put the TjIDs into a vector for matching
1708  float score = 0;
1709  auto v3TjIDs = GetVtxTjIDs(slc, vx3, score);
1710  if(v3TjIDs.empty()) continue;
1711  if(score < tcc.vtx2DCuts[7]) continue;
1712  SortEntry se;
1713  se.index = iv3;
1714  se.val = score;
1715  sortVec.push_back(se);
1716  } // vx3
1717  if(sortVec.empty()) return;
1718  if(sortVec.size() > 1) std::sort(sortVec.begin(), sortVec.end(), valDecreasing);
1719 
1720  for(unsigned short ii = 0; ii < sortVec.size(); ++ii) {
1721  auto& vx3 = slc.vtx3s[sortVec[ii].index];
1722  float score = 0;
1723  std::vector<int> v3TjIDs = GetVtxTjIDs(slc, vx3, score);
1724  if(prt) {
1725  mf::LogVerbatim myprt("TC");
1726  myprt<<"M3DVTj 3V"<<vx3.ID<<" score "<<score<<" Tjs";
1727  for(auto& tjID : v3TjIDs) myprt<<" T"<<tjID;
1728  }
1729  for(unsigned int ims = 0; ims < slc.matchVec.size(); ++ims) {
1730  auto& ms = slc.matchVec[ims];
1731  if(ms.Count == 0) continue;
1732  bool skipit = false;
1733  for(auto tjid : ms.TjIDs) {
1734  auto& tj = slc.tjs[tjid - 1];
1735  if(tj.AlgMod[kMat3D]) skipit = true;
1736  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) skipit = true;
1737  }
1738  if(skipit) continue;
1739  std::vector<int> shared = SetIntersection(ms.TjIDs, v3TjIDs);
1740  if(shared.size() < 2) continue;
1741  if(prt) mf::LogVerbatim("TC")<<" ims "<<ims<<" shared size "<<shared.size();
1742  PFPStruct pfp = CreatePFP(slc);
1743  pfp.TjIDs = ms.TjIDs;
1744  pfp.Vx3ID[0] = vx3.ID;
1745  // note that the ms position is the average position of all 3D matched Tp3s at this point.
1746  // It is not the start position. This will be determined in DefinePFP.
1747  pfp.XYZ[0] = ms.Pos;
1748  pfp.Dir[0] = ms.Dir;
1749  pfp.MatchVecIndex = ims;
1750  if(prt) mf::LogVerbatim("TC")<<"M3DVTj: pfp P"<<pfp.ID<<" 3V"<<vx3.ID<<" ims "<<ims;
1751  // Set the PDGCode so DefinePFP can ignore incompatible matches
1752  pfp.PDGCode = PDGCodeVote(slc, pfp.TjIDs, prt);
1753  // make a copy of the TjIDs to test for changes
1754  auto saveTjIDs = pfp.TjIDs;
1755  if(!DefinePFP("M3DVTj1", slc, pfp, prt)) continue;
1756  // separation distance (cm) for kink detection.
1757  double sep = 1;
1758  bool didSplit = Split3DKink(slc, pfp, sep, prt);
1759  if(prt) PrintPFP("M3D", slc, pfp, true);
1760  if(!didSplit && shared.size() != ms.TjIDs.size()) {
1761  // Try to repair the PFParticle by merging the Tj that was in the match list but
1762  // wasn't attached to the vertex. Hopefully there aren't more than one...
1763  auto tjNotInVx = SetDifference(ms.TjIDs, shared);
1764  }
1765  AnalyzePFP(slc, pfp, prt);
1766  if(!StorePFP(slc, pfp)) continue;
1767  if(pfp.TjIDs != saveTjIDs) {
1768  // v3TjIDs is wrong if Tjs were merged so re-make it.
1769  auto tmp = GetVtxTjIDs(slc, vx3, score);
1770  v3TjIDs.clear();
1771  // then re-build it
1772  for(auto tjid : tmp) {
1773  auto& tj = slc.tjs[tjid - 1];
1774  if(!tj.AlgMod[kMat3D]) v3TjIDs.push_back(tjid);
1775  } // tjid
1776  if(v3TjIDs.size() < 2) break;
1777  } // pfp.TjIDs != saveTjIDs
1778  ms.Count = 0;
1779  // clobber MatchStructs that use the Tjs in this pfp
1780  for(auto& allms : slc.matchVec) {
1781  auto mfpfp = SetIntersection(allms.TjIDs, pfp.TjIDs);
1782  if(!mfpfp.empty()) allms.Count = 0;
1783  } // allms
1784  auto leftover = SetDifference(v3TjIDs, pfp.TjIDs);
1785  if(prt) {
1786  mf::LogVerbatim myprt("TC");
1787  myprt<<" leftover";
1788  for(auto tjid : leftover) myprt<<" T"<<tjid;
1789  }
1790  if(leftover.size() < 2) break;
1791  // keep looking using the leftovers
1792  v3TjIDs = leftover;
1793  } // ims
1794  if(v3TjIDs.size() > 1 && v3TjIDs.size() <= slc.nPlanes) {
1795  // a last-ditch attempt
1796  PFPStruct pfp = CreatePFP(slc);
1797  // put the available leftovers in
1798  for(auto& tjid : v3TjIDs) {
1799  auto& tj = slc.tjs[tjid - 1];
1800  if(tj.AlgMod[kMat3D]) continue;
1801  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1802  pfp.TjIDs.push_back(tjid);
1803  } // tjid
1804  if(pfp.TjIDs.size() < 2) continue;
1805  pfp.Vx3ID[0] = vx3.ID;
1806  if(!DefinePFP("M3DVTj2", slc, pfp, prt)) continue;
1807  Split3DKink(slc, pfp, 1, prt);
1808  AnalyzePFP(slc, pfp, prt);
1809  if(prt) PrintPFP("left", slc, pfp, true);
1810  if(!StorePFP(slc, pfp)) continue;
1811  }
1812  } // ims
1813  } // Match3DVtxTjs
1814 
1816  unsigned short TPNearVertex(TCSlice& slc, const TrajPoint& tp)
1817  {
1818  // Returns the index of a vertex if tp is nearby
1819  for(unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
1820  if(slc.vtxs[ivx].ID == 0) continue;
1821  if(slc.vtxs[ivx].CTP != tp.CTP) continue;
1822  if(std::abs(slc.vtxs[ivx].Pos[0] - tp.Pos[0]) > 1.2) continue;
1823  if(std::abs(slc.vtxs[ivx].Pos[1] - tp.Pos[1]) > 1.2) continue;
1824  return ivx;
1825  } // ivx
1826  return USHRT_MAX;
1827  } // TPNearVertex
1828 
1830  bool AttachPFPToVertex(TCSlice& slc, PFPStruct& pfp, unsigned short end, unsigned short vx3ID, bool prt)
1831  {
1832  if(vx3ID > int(slc.vtx3s.size())) return false;
1833  if(pfp.ID > int(slc.pfps.size())) return false;
1834  if(pfp.PDGCode == 22) return false;
1835  if(end > 1) return false;
1836 
1837  auto& vx3 = slc.vtx3s[vx3ID - 1];
1838 
1839  pfp.Vx3ID[end] = vx3.ID;
1840 
1841  // We are done if this a PFP-only vertex
1842  if(vx3.Wire == -2) return true;
1843 
1844  // Update the 2D and 3D vertex and tj associations
1845  for(auto tjid : pfp.TjIDs) {
1846  auto& tj = slc.tjs[tjid - 1];
1847  unsigned short plane = DecodeCTP(tj.CTP).Plane;
1848  // TODO: Check to see if the Tjs have been ordered correctly?
1849  if(tj.VtxID[end] == 0) {
1850  // tj is available to be attached to a 2D vertex. See if the 3D vertex is matched to
1851  // an existing 2D vertex in this plane
1852  if(vx3.Vx2ID[plane] == 0) {
1853  // not matched. Look for one
1854  std::array<float, 2> pos;
1855  PosInPlane(slc, vx3, plane, pos);
1856 // if(prt) std::cout<<" tj "<<tj.ID<<" has no 2D vertex. Look for one vertex near "<<tj.CTP<<":"<<PrintPos(slc, pos)<<" Events processed "<<slc.EventsProcessed<<"\n";
1857  } else {
1858  // Existing 2D vertex matched to the 3D vertex
1859 // if(prt) std::cout<<" tj "<<tj.ID<<" has no 2D vertex in CTP "<<tj.CTP<<" but vx3 is matched to 2D vertex"<<vx3.Vx2ID[plane]<<". Attach it? Events processed "<<slc.EventsProcessed<<"\n";
1860  }
1861  }
1862  } // tjid
1863 
1864  return true;
1865  } // AttachPFPToVertex
1866 
1868  bool AttachAnyTrajToVertex(TCSlice& slc, unsigned short ivx, bool prt)
1869  {
1870 
1871  if(ivx > slc.vtxs.size() - 1) return false;
1872  if(slc.vtxs[ivx].ID == 0) return false;
1873  if(tcc.vtx2DCuts[0] < 0) return false;
1874 
1875  VtxStore& vx = slc.vtxs[ivx];
1876  // Hammer vertices should be isolated and clean
1877  if(vx.Topo == 5 || vx.Topo == 6) return false;
1878 
1879  unsigned short nadd = 0;
1880  for(auto& tj : slc.tjs) {
1881  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1882  if(tj.CTP != vx.CTP) continue;
1883  if(tj.VtxID[0] == vx.ID || tj.VtxID[1] == vx.ID) continue;
1884  if(AttachTrajToVertex(slc, tj, vx, prt)) ++nadd;
1885  } // tj
1886  if(prt) mf::LogVerbatim("TC")<<" AttachAnyTrajToVertex: nadd "<<nadd;
1887  if(nadd == 0) return false;
1888  return true;
1889 
1890  } // AttachAnyTrajToVertex
1891 
1893  bool AttachTrajToVertex(TCSlice& slc, Trajectory& tj, VtxStore& vx, bool prt)
1894  {
1895  // Note that this function does not require a signal between the end of the Tj and the vertex
1896 
1897  // tcc.vtx2DCuts fcl input usage
1898  // 0 = maximum length of a short trajectory
1899  // 1 = max vertex - trajectory separation for short trajectories
1900  // 2 = max vertex - trajectory separation for long trajectories
1901  // 3 = max position pull for adding TJs to a vertex
1902  // 4 = max allowed vertex position error
1903  // 5 = min MCSMom
1904  // 6 = min Pts/Wire fraction
1905 
1906  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) return false;
1907  if(tj.CTP != vx.CTP) return false;
1908  // already attached?
1909  if(tj.VtxID[0] == vx.ID || tj.VtxID[1] == vx.ID) return false;
1910 
1911  unsigned short maxShortTjLen = tcc.vtx2DCuts[0];
1912  // square the separation cut to simplify testing in the loop
1913  float maxSepCutShort2 = tcc.vtx2DCuts[1] * tcc.vtx2DCuts[1];
1914  float maxSepCutLong2 = tcc.vtx2DCuts[2] * tcc.vtx2DCuts[2];
1915 
1916  // assume that end 0 is closest to the vertex
1917  unsigned short end = 0;
1918  float vtxTjSep2 = PosSep2(vx.Pos, tj.Pts[tj.EndPt[0]].Pos);
1919  float sep1 = PosSep2(vx.Pos, tj.Pts[tj.EndPt[1]].Pos);
1920  if(sep1 < vtxTjSep2) {
1921  // End 1 is closer
1922  end = 1;
1923  vtxTjSep2 = sep1;
1924  }
1925  // is there a vertex already assigned to this end?
1926  if(tj.VtxID[end] > 0) return false;
1927 
1928  // is the trajectory short?
1929  bool tjShort = (tj.EndPt[1] - tj.EndPt[0] < maxShortTjLen);
1930  // use the short Tj cut if the trajectory looks like an electron
1931  if(tcc.useAlg[kNewVtxCuts] && !tjShort && tj.ChgRMS > 0.5) tjShort = true;
1932  float closestApproach;
1933  // ignore bad separation between the closest tj end and the vertex
1934  if(tjShort) {
1935  if(vtxTjSep2 > maxSepCutShort2) return false;
1936  closestApproach = tcc.vtx2DCuts[1];
1937  } else {
1938  closestApproach = tcc.vtx2DCuts[2];
1939  if(vtxTjSep2 > maxSepCutLong2) return false;
1940  }
1941 
1942  // Calculate the pull on the vertex
1943  TrajPoint& tp = tj.Pts[tj.EndPt[end]];
1944  float tpVxPull = TrajPointVertexPull(slc, tp, vx);
1945  bool signalBetween = SignalBetween(slc, tp, vx.Pos[0], 0.8);
1946 
1947  // See if the vertex position is close to an end
1948  unsigned short closePt;
1949  TrajClosestApproach(tj, vx.Pos[0], vx.Pos[1], closePt, closestApproach);
1950  // count the number of points between the end of the trajectory and the vertex.
1951  // tj ------------- tj ------------
1952  // vx * >> dpt = 0 vx * >> dpt = 2
1953  short dpt;
1954  if(end == 0) {
1955  dpt = closePt - tj.EndPt[end];
1956  } else {
1957  dpt = tj.EndPt[end] - closePt;
1958  }
1959 
1960  float length = TrajLength(tj);
1961  // don't attach it if the tj length is shorter than the separation distance
1962  if(length > 4 && length < closestApproach) return false;
1963 
1964  float pullCut = tcc.vtx2DCuts[3];
1965  // Dec 21, 2017 Loosen up the pull cut for short close slc. These are likely to
1966  // be poorly reconstructed. It is better to have them associated with the vertex
1967  // than not. Oct 5, 2018 This was too loose.
1968  if(tcc.useAlg[kNewVtxCuts]) {
1969  if(tjShort) pullCut *= 2;
1970  } else {
1971  if(tjShort) pullCut = 10;
1972  }
1973 
1974  if(prt) {
1975  mf::LogVerbatim myprt("TC");
1976  myprt<<"ATTV: 2V"<<vx.ID;
1977  myprt<<" NTraj "<<vx.NTraj;
1978  myprt<<" oldTJs";
1979  for(unsigned short itj = 0; itj < slc.tjs.size(); ++itj) {
1980  Trajectory& tj = slc.tjs[itj];
1981  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1982  if(tj.CTP != vx.CTP) continue;
1983  if(tj.VtxID[0] == vx.ID) myprt<<" T"<<tj.ID<<"_0";
1984  if(tj.VtxID[1] == vx.ID) myprt<<" T"<<tj.ID<<"_1";
1985  }
1986  myprt<<" + T"<<tj.ID<<"_"<<end<<" vtxTjSep "<<sqrt(vtxTjSep2)<<" tpVxPull "<<tpVxPull<<" pullCut "<<pullCut<<" dpt "<<dpt;
1987  }
1988 // if(tpVxPull > tcc.vtx2DCuts[3]) return false;
1989  if(tpVxPull > pullCut) return false;
1990  if(dpt > 2) return true;
1991 
1992  // remove the fixed position flag if there are more than 2 tjs
1993  bool fixedBit = vx.Stat[kFixed];
1994  if(fixedBit && vx.NTraj < 2) vx.Stat[kFixed] = false;
1995 
1996  // don't allow a short Tj with a large pull to bias the fit
1997  if(tjShort && tpVxPull > tcc.vtx2DCuts[3]) tj.AlgMod[kNoFitToVx] = true;
1998 
1999  // Passed all the cuts. Attach it to the vertex and try a fit
2000  tj.VtxID[end] = vx.ID;
2001  // flag as a photon Tj so it isn't included in the fit
2002  tj.AlgMod[kPhoton] = !signalBetween;
2003  // make a copy of the vertex and fit it
2004  auto vxTmp = vx;
2005  if(FitVertex(slc, vxTmp, prt)) {
2006  SetVx2Score(slc, vxTmp);
2007  if(prt) mf::LogVerbatim("TC")<<" Success";
2008  vx = vxTmp;
2009  return true;
2010  }
2011 
2012  // fit failed so remove the tj -> vx assignment if it is long and
2013  // set noFitToVtx if it is short
2014  if(tjShort) {
2015  tj.AlgMod[kNoFitToVx] = true;
2016  if(prt) mf::LogVerbatim("TC")<<" Poor fit. Keep short Tj "<<tj.ID<<" with kNoFitToVx";
2017  return true;
2018  } else {
2019  tj.VtxID[end] = 0;
2020  // restore the fixed flag
2021  vx.Stat[kFixed] = fixedBit;
2022  if(prt) mf::LogVerbatim("TC")<<" Poor fit. Removed Tj "<<tj.ID;
2023  return false;
2024  }
2025 
2026  } // AttachTrajToVertex
2027 
2029  float TrajPointVertexPull(TCSlice& slc, const TrajPoint& tp, const VtxStore& vx)
2030  {
2031  // Calculates the position pull between a trajectory point and a vertex
2032 
2033  // impact parameter between them
2034  double ip = PointTrajDOCA(slc, vx.Pos[0], vx.Pos[1], tp);
2035  // separation^2
2036  double sep2 = PosSep2(vx.Pos, tp.Pos);
2037 
2038  // Find the projection of the vertex error ellipse in a coordinate system
2039  // along the TP direction
2040  double vxErrW = vx.PosErr[0] * tp.Dir[1];
2041  double vxErrT = vx.PosErr[1] * tp.Dir[0];
2042  double vxErr2 = vxErrW * vxErrW + vxErrT * vxErrT;
2043  // add the TP position error^2
2044  vxErr2 += tp.HitPosErr2;
2045 
2046  // close together so ignore the TP projection error and return
2047  // the pull using the vertex error and TP position error
2048  if(sep2 < 1) return (float)(ip/sqrt(vxErr2));
2049 
2050  double dang = ip / sqrt(sep2);
2051 
2052  // calculate the angle error.
2053  // Start with the vertex error^2
2054  double angErr = vxErr2 / sep2;
2055  // Add the TP angle error^2
2056  angErr += tp.AngErr * tp.AngErr;
2057  if(angErr == 0) return 999;
2058  angErr = sqrt(angErr);
2059  return (float)(dang / angErr);
2060 
2061  } // TrajPointVertexPull
2062 
2064  float VertexVertexPull(TCSlice& slc, const Vtx3Store& vx1, const Vtx3Store& vx2)
2065  {
2066  // Calculates the position pull between two vertices
2067  double dx = vx1.X - vx2.X;
2068  double dy = vx1.Y - vx2.Y;
2069  double dz = vx1.Z - vx2.Z;
2070  double dxErr2 = (vx1.XErr * vx1.XErr + vx2.XErr * vx2.XErr) / 2;
2071  double dyErr2 = (vx1.YErr * vx1.YErr + vx2.YErr * vx2.YErr) / 2;
2072  double dzErr2 = (vx1.ZErr * vx1.ZErr + vx2.ZErr * vx2.ZErr) / 2;
2073  dx = dx * dx / dxErr2;
2074  dy = dy * dy / dyErr2;
2075  dz = dz * dz / dzErr2;
2076  return (float)(sqrt(dx + dy + dz)/3);
2077  }
2078 
2080  float VertexVertexPull(TCSlice& slc, const VtxStore& vx1, const VtxStore& vx2)
2081  {
2082  // Calculates the position pull between two vertices
2083  double dw = vx1.Pos[0] - vx2.Pos[0];
2084  double dt = vx1.Pos[1] - vx2.Pos[1];
2085  double dwErr2 = (vx1.PosErr[0] * vx1.PosErr[0] + vx2.PosErr[0] * vx2.PosErr[0]) / 2;
2086  double dtErr2 = (vx1.PosErr[1] * vx1.PosErr[1] + vx2.PosErr[1] * vx2.PosErr[1]) / 2;
2087  dw = dw * dw / dwErr2;
2088  dt = dt * dt / dtErr2;
2089  return (float)sqrt(dw + dt);
2090  }
2091 
2093  bool StoreVertex(TCSlice& slc, VtxStore& vx)
2094  {
2095  // jacket around the push to ensure that the Tj and vtx CTP is consistent.
2096  // The calling function should score the vertex after the trajectories are attached
2097 
2098  if(vx.ID != int(slc.vtxs.size() + 1)) {
2099  mf::LogVerbatim("TC")<<"StoreVertex: Invalid ID "<<vx.ID<<" It should be "<<slc.vtxs.size() + 1;
2100  return false;
2101  }
2102 
2103  ++evt.globalS2ID;
2104  vx.UID = evt.globalS2ID;
2105 
2106  unsigned short nvxtj = 0;
2107  unsigned short nok = 0;
2108  for(auto& tj : slc.tjs) {
2109  if(tj.AlgMod[kKilled]) continue;
2110  if(vx.ID == tj.VtxID[0] || vx.ID == tj.VtxID[1]) ++nvxtj;
2111  if(vx.CTP != tj.CTP) continue;
2112  if(vx.ID == tj.VtxID[0] || vx.ID == tj.VtxID[1]) ++nok;
2113  } // tj
2114 
2115  if(nok != nvxtj) {
2116  mf::LogVerbatim("TC")<<"StoreVertex: vertex "<<vx.ID<<" Topo "<<vx.Topo<<" has inconsistent CTP code "<<vx.CTP<<" with one or more Tjs\n";
2117  for(auto& tj : slc.tjs) {
2118  if(tj.AlgMod[kKilled]) continue;
2119  if(tj.VtxID[0] == vx.ID) tj.VtxID[0] = 0;
2120  if(tj.VtxID[1] == vx.ID) tj.VtxID[1] = 0;
2121  }
2122  return false;
2123  }
2124  vx.NTraj = nok;
2125  slc.vtxs.push_back(vx);
2126  return true;
2127 
2128  } // StoreVertex
2129 
2131  bool FitVertex(TCSlice& slc, VtxStore& vx, bool prt)
2132  {
2133  // A poor-mans fitting scheme. If the fitted vertex position error is within the supplied
2134  // value, the position and errors are updated and we return true, otherwise the vertex is
2135  // left unchanged and we return false
2136 
2137  // tcc.vtx2DCuts fcl input usage
2138  // 0 = maximum length of a short trajectory
2139  // 1 = max vertex - trajectory separation for short trajectories
2140  // 2 = max vertex - trajectory separation for long trajectories
2141  // 3 = max position pull for adding TJs to a vertex
2142  // 4 = max allowed vertex position error
2143  // 5 = min MCSMom
2144  // 6 = min Pts/Wire fraction
2145 
2146  if(vx.Stat[kFixed]) {
2147  if(prt) mf::LogVerbatim("TC")<<" vertex position fixed. No fit allowed";
2148  return true;
2149  }
2150 
2151  // Create a vector of trajectory points that will be used to fit the vertex position
2152  std::vector<TrajPoint> vxTp;
2153  for(auto& tj : slc.tjs) {
2154  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2155  if(tj.CTP != vx.CTP) continue;
2156  if(tj.AlgMod[kPhoton]) continue;
2157  if(tj.AlgMod[kNoFitToVx]) continue;
2158  bool added = false;
2159  if(tj.VtxID[0] == vx.ID) {
2160  vxTp.push_back(tj.Pts[tj.EndPt[0]]);
2161  added = true;
2162  }
2163  if(tj.VtxID[1] == vx.ID) {
2164  vxTp.push_back(tj.Pts[tj.EndPt[1]]);
2165  added = true;
2166  }
2167  // stash the ID in Step to help debugging
2168  if(added) {
2169  auto& tp = vxTp[vxTp.size()-1];
2170  if(tj.ID > 0) tp.Step = (int)tj.ID;
2171  // inflate the angle errors for Tjs with few fitted points
2172  if(tcc.useAlg[kNewVtxCuts] && tp.NTPsFit < 4) tp.AngErr *= 4;
2173  }
2174  } // tj
2175 
2176  bool success = FitVertex(slc, vx, vxTp, prt);
2177 
2178  if(!success) return false;
2179  return true;
2180 
2181  } // FitVertex
2182 
2184  bool FitVertex(TCSlice& slc, VtxStore& vx, std::vector<TrajPoint> vxTp, bool prt)
2185  {
2186  // Variant of FitVertex that fits the passed trajectory points to a vertex position but doesn't
2187  // require using information in TJStuff
2188 
2189  // tcc.vtx2DCuts fcl input usage
2190  // 0 = maximum length of a short trajectory
2191  // 1 = max vertex - trajectory separation for short trajectories
2192  // 2 = max vertex - trajectory separation for long trajectories
2193  // 3 = max position pull for adding TJs to a vertex
2194  // 4 = max allowed vertex position error
2195  // 5 = min MCSMom
2196  // 6 = min Pts/Wire fraction
2197  // 7 min Score
2198  // 8 Min charge fraction near a merge point (not a vertex)
2199  // 9 max MCSmom asymmetry for a merge
2200 
2201 
2202  vx.NTraj = vxTp.size();
2203 
2204  if(vxTp.size() < 2) return false;
2205 
2206  if(prt) {
2207  PrintHeader("FV");
2208  for(auto& tp : vxTp) PrintTrajPoint("FV", slc, 0, 1, 1, tp);
2209  }
2210 
2211  // Find trajectory intersections pair-wise tweaking the angle and position(?) within
2212  // +/- 1 sigma
2213  double sum0 = 0, sum02 = 0;
2214  double sum1 = 0, sum12 = 0;
2215  double sumw = 0;
2216  double wgt;
2217  double cnt = 0;
2218  // a temporary TP for tweaking the angle
2219  TrajPoint tmp;
2220  // another point to check for a signal at each intersection
2221  TrajPoint intTp;
2222  intTp.CTP = vxTp[0].CTP;
2223  for(unsigned short itj = 0; itj < vxTp.size() - 1; ++itj) {
2224  for(unsigned short jtj = itj + 1; jtj < vxTp.size(); ++jtj) {
2225  float p0, p1;
2226  TrajIntersection(vxTp[itj], vxTp[jtj], p0, p1);
2227  intTp.Pos[0] = p0; intTp.Pos[1] = p1;
2228  wgt = 1;
2229  // accumulate
2230  sum0 += wgt * p0; sum02 += wgt * p0 * p0;
2231  sum1 += wgt * p1; sum12 += wgt * p1 * p1; sumw += wgt; ++cnt;
2232  // tweak the itj angle +
2233  tmp = vxTp[itj];
2234  tmp.Ang += tmp.AngErr;
2235  tmp.Dir[0] = cos(tmp.Ang); tmp.Dir[1] = sin(tmp.Ang);
2236  TrajIntersection(tmp, vxTp[jtj], p0, p1);
2237  intTp.Pos[0] = p0; intTp.Pos[1] = p1;
2238  // adjust the weight for 4 points at +/1 1 sigma = 0.607 / 4
2239  wgt = 0.152;
2240  // accumulate
2241  sum0 += wgt * p0; sum02 += wgt * p0 * p0;
2242  sum1 += wgt * p1; sum12 += wgt * p1 * p1; sumw += wgt; ++cnt;
2243  // tweak the itj angle -
2244  tmp = vxTp[itj];
2245  tmp.Ang -= 2 * tmp.AngErr;
2246  tmp.Dir[0] = cos(tmp.Ang); tmp.Dir[1] = sin(tmp.Ang);
2247  TrajIntersection(tmp, vxTp[jtj], p0, p1);
2248  intTp.Pos[0] = p0; intTp.Pos[1] = p1;
2249  // accumulate
2250  sum0 += wgt * p0; sum02 += wgt * p0 * p0;
2251  sum1 += wgt * p1; sum12 += wgt * p1 * p1; sumw += wgt; ++cnt;
2252  // Repeat this process with jtj
2253  // tweak the jtj angle +
2254  tmp = vxTp[jtj];
2255  tmp.Ang += tmp.AngErr;
2256  tmp.Dir[0] = cos(tmp.Ang); tmp.Dir[1] = sin(tmp.Ang);
2257  TrajIntersection(vxTp[itj], tmp, p0, p1);
2258  intTp.Pos[0] = p0; intTp.Pos[1] = p1;
2259  // accumulate
2260  sum0 += wgt * p0; sum02 += wgt * p0 * p0;
2261  sum1 += wgt * p1; sum12 += wgt * p1 * p1; sumw += wgt; ++cnt;
2262  // tweak the itj angle -
2263  tmp = vxTp[itj];
2264  tmp.Ang -= 2 * tmp.AngErr;
2265  tmp.Dir[0] = cos(tmp.Ang); tmp.Dir[1] = sin(tmp.Ang);
2266  TrajIntersection(vxTp[itj], tmp, p0, p1);
2267  intTp.Pos[0] = p0; intTp.Pos[1] = p1;
2268  // accumulate
2269  sum0 += wgt * p0; sum02 += wgt * p0 * p0;
2270  sum1 += wgt * p1; sum12 += wgt * p1 * p1; sumw += wgt; ++cnt;
2271  } // jtj
2272  } // itj
2273 
2274  if(sumw == 0) return false;
2275 
2276  double vxP0 = sum0 / sumw;
2277  double vxP1 = sum1 / sumw;
2278  double vxP0rms = sqrt((sum02 - sumw * vxP0 * vxP0) / sumw);
2279  double vxP1rms = sqrt((sum12 - sumw * vxP1 * vxP1) / sumw);
2280  double rootN = sqrt(cnt);
2281  vxP0rms /= rootN;
2282  vxP1rms /= rootN;
2283  // don't let the errors get too small
2284  if(vxP0rms < 0.5) vxP0rms = 0.5;
2285  if(vxP1rms < 0.5) vxP1rms = 0.5;
2286 
2287  if(prt) mf::LogVerbatim("TC")<<"FitVertex 2V"<<vx.ID<<" CTP "<<vx.CTP<<" NTraj "<<vx.NTraj<<" in "<<std::fixed<<std::setprecision(1)<<vx.Pos[0]<<":"<<vx.Pos[1]/tcc.unitsPerTick<<" out wire "<<vxP0<<" +/- "<<vxP0rms<<" ticks "<<vxP1/tcc.unitsPerTick<<"+/-"<<vxP1rms/tcc.unitsPerTick;
2288 
2289  // apply Vertex2DCuts if this isn't a neutral vertex (which is expected to have very large
2290  // errors)
2291  if(vx.Topo != 11) {
2292  float inflate = 1;
2293  if(vx.Stat[kOnDeadWire]) inflate = 1.5;
2294  if(vxP0rms > inflate * tcc.vtx2DCuts[4] || vxP1rms > inflate * tcc.vtx2DCuts[4]) {
2295  if(prt) mf::LogVerbatim("TC")<<" fit failed. Max allowed position error "<<inflate * tcc.vtx2DCuts[4];
2296  return false;
2297  }
2298  } // not a neutral vertex
2299 
2300  vx.Pos[0] = vxP0;
2301  vx.PosErr[0] = vxP0rms;
2302  vx.Pos[1] = vxP1;
2303  vx.PosErr[1] = vxP1rms;
2304 
2305  // Calculate chisq
2306  vx.ChiDOF = 0;
2307  for(unsigned short itj = 0; itj < vxTp.size(); ++itj) {
2308  vx.ChiDOF += TrajPointVertexPull(slc, vxTp[itj], vx);
2309  } // itj
2310  vx.ChiDOF /= (float)vxTp.size();
2311 
2312  if(prt) {
2313  mf::LogVerbatim myprt("TC");
2314  myprt<<"Pull";
2315  for(unsigned short itj = 0; itj < vxTp.size(); ++itj) {
2316  float pull = TrajPointVertexPull(slc, vxTp[itj], vx);
2317  myprt<<" "<<PrintPos(slc, vxTp[itj])<<" - "<<std::fixed<<std::setprecision(2)<<pull;
2318  } // itj
2319  myprt<<" ChiDOF "<<vx.ChiDOF;
2320  }
2321  return true;
2322 
2323  } // FitVertex
2324 
2326  bool ChkVtxAssociations(TCSlice& slc, const CTP_t& inCTP)
2327  {
2328  // Check the associations
2329 
2330  // check the 2D -> 3D associations
2331  geo::PlaneID planeID = DecodeCTP(inCTP);
2332  unsigned short plane = planeID.Plane;
2333  for(auto& vx2 : slc.vtxs) {
2334  if(vx2.CTP != inCTP) continue;
2335  if(vx2.ID == 0) continue;
2336  if(vx2.Vx3ID == 0) continue;
2337  if(vx2.Vx3ID > int(slc.vtx3s.size())) {
2338  mf::LogVerbatim("TC")<<"ChkVtxAssociations: Invalid vx2.Vx3ID "<<vx2.Vx3ID<<" in 2D vtx "<<vx2.ID;
2339  return false;
2340  }
2341  auto& vx3 = slc.vtx3s[vx2.Vx3ID-1];
2342  if(vx3.ID == 0) {
2343  mf::LogVerbatim("TC")<<"ChkVtxAssociations: 2V"<<vx2.ID<<" thinks it is matched to 3V"<<vx3.ID<<" but vx3 is obsolete";
2344  return false;
2345  }
2346  if(vx3.Vx2ID[plane] != vx2.ID) {
2347  mf::LogVerbatim("TC")<<"ChkVtxAssociations: 2V"<<vx2.ID<<" thinks it is matched to 3V"<<vx3.ID<<" but vx3 says no!";
2348  return false;
2349  }
2350  } // vx2
2351  // check the 3D -> 2D associations
2352  for(auto& vx3 : slc.vtx3s) {
2353  if(vx3.ID == 0) continue;
2354  if(vx3.Vx2ID[plane] == 0) continue;
2355  if(vx3.Vx2ID[plane] > (int)slc.vtxs.size()) {
2356  mf::LogVerbatim("TC")<<"ChkVtxAssociations: Invalid vx3.Vx2ID "<<vx3.Vx2ID[plane]<<" in CTP "<<inCTP;
2357  return false;
2358  }
2359  auto& vx2 = slc.vtxs[vx3.Vx2ID[plane]-1];
2360  if(vx2.Vx3ID != vx3.ID) {
2361  mf::LogVerbatim("TC")<<"ChkVtxAssociations: 3V"<<vx3.ID<<" thinks it is matched to 2V"<<vx2.ID<<" but vx2 says no!";
2362  return false;
2363  }
2364  } // vx3
2365 
2366  // check the Tj -> 2D associations
2367  for(auto& tj : slc.tjs) {
2368  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2369  for(unsigned short end = 0; end < 2; ++end) {
2370  if(tj.VtxID[end] == 0) continue;
2371  if(tj.VtxID[end] > slc.vtxs.size()) {
2372  mf::LogVerbatim("TC")<<"ChkVtxAssociations: T"<<tj.ID<<" thinks it is matched to 2V"<<tj.VtxID[end]<<" on end "<<end<<" but no vertex exists. Recovering";
2373  tj.VtxID[end] = 0;
2374  return false;
2375  }
2376  unsigned short ivx = tj.VtxID[end] - 1;
2377  auto& vx2 = slc.vtxs[ivx];
2378  if(vx2.ID == 0) {
2379  mf::LogVerbatim("TC")<<"ChkVtxAssociations: T"<<tj.ID<<" thinks it is matched to 2V"<<tj.VtxID[end]<<" on end "<<end<<" but the vertex is killed. Fixing the Tj";
2380  tj.VtxID[end] = 0;
2381  return false;
2382  }
2383  } // end
2384  } // tj
2385 
2386  return true;
2387 
2388  } // ChkVtxAssociations
2389 
2391  void ScoreVertices(TCSlice& slc)
2392  {
2393  // reset all 3D vertex, 2D vertex and Tj high-score vertex bits in tpcid
2394 
2395  // reset the 2D vertex status bits
2396  for(auto& vx : slc.vtxs) {
2397  if(vx.ID == 0) continue;
2398  vx.Stat[kHiVx3Score] = false;
2399  } // vx
2400  // and the tj bits
2401  for(auto& tj : slc.tjs) {
2402  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2403  tj.AlgMod[kTjHiVx3Score] = false;
2404  } // tj
2405  // Score the 2D vertices
2406  for(auto& vx : slc.vtxs) {
2407  if(vx.ID == 0) continue;
2408  SetVx2Score(slc, vx);
2409  } // vx
2410  // Score the 3D vertices
2411  for(auto& vx3 : slc.vtx3s) {
2412  if(vx3.ID == 0) continue;
2413  SetVx3Score(slc, vx3);
2414  } // vx3
2415  } // ScoreVertices
2416 
2418  void KillPoorVertices(TCSlice& slc)
2419  {
2420  // kill 2D vertices that have low score and are not attached to a high-score 3D vertex
2421  if(slc.vtxs.empty()) return;
2422  for(auto& vx : slc.vtxs) {
2423  if(vx.ID == 0) continue;
2424  if(vx.Score > tcc.vtx2DCuts[7]) continue;
2425  if(vx.Vx3ID > 0) {
2426  auto& vx3 = slc.vtx3s[vx.Vx3ID - 1];
2427  if(vx3.Primary) continue;
2428  if(slc.vtx3s[vx.Vx3ID - 1].Score >= tcc.vtx2DCuts[7]) continue;
2429  }
2430  MakeVertexObsolete("KPV", slc, vx, false);
2431  } // vx
2432 
2433  } // KillPoorVertices
2434 
2436  void SetHighScoreBits(TCSlice& slc, Vtx3Store& vx3)
2437  {
2438  // Sets the tj and 2D vertex score bits to true
2439 
2440  if(vx3.ID == 0) return;
2441 
2442  for(unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
2443  if(vx3.Vx2ID[ipl] <= 0) continue;
2444  VtxStore& vx2 = slc.vtxs[vx3.Vx2ID[ipl] - 1];
2445  vx2.Stat[kHiVx3Score] = false;
2446  // transfer this to all attached tjs and vertices attached to those tjs
2447  std::vector<int> tjlist = GetVtxTjIDs(slc, vx2);
2448  std::vector<int> vxlist;
2449  while(true) {
2450  // tag Tjs and make a list of attached vertices whose high-score
2451  // bit needs to be set
2452  vxlist.clear();
2453  for(auto tjid : tjlist) {
2454  auto& tj = slc.tjs[tjid - 1];
2455  tj.AlgMod[kTjHiVx3Score] = true;
2456  for(unsigned short end = 0; end < 2; ++end) {
2457  if(tj.VtxID[end] == 0) continue;
2458  auto& vx2 = slc.vtxs[tj.VtxID[end] - 1];
2459  if(vx2.Stat[kHiVx3Score]) continue;
2460  vx2.Stat[kHiVx3Score] = true;
2461  vxlist.push_back(vx2.ID);
2462  } // end
2463  } // tjid
2464 
2465  if(vxlist.empty()) break;
2466  // re-build tjlist using vxlist
2467  std::vector<int> newtjlist;
2468  for(auto vxid : vxlist) {
2469  auto& vx2 = slc.vtxs[vxid - 1];
2470  auto tmp = GetVtxTjIDs(slc, vx2);
2471  for(auto tjid : tmp) {
2472  if(std::find(tjlist.begin(), tjlist.end(), tjid) == tjlist.end()) newtjlist.push_back(tjid);
2473  } // tjid
2474  } // vxid
2475  if(newtjlist.empty()) break;
2476  tjlist = newtjlist;
2477  } // true
2478  } // ipl
2479 
2480  } // SetHighScoreBits
2481 
2483  void SetVx3Score(TCSlice& slc, Vtx3Store& vx3)
2484  {
2485  // Calculate the 3D vertex score and flag Tjs that are attached to high score vertices as defined
2486  // by Vertex2DCuts
2487 
2488  if(vx3.ID == 0) return;
2489 
2490  vx3.Score = 0;
2491  for(unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
2492  if(vx3.Vx2ID[ipl] <= 0) continue;
2493  VtxStore& vx2 = slc.vtxs[vx3.Vx2ID[ipl] - 1];
2494  vx3.Score += vx2.Score;
2495  } // ipl
2496  vx3.Score /= (float)slc.nPlanes;
2497  // don't allow it to get too small or negative
2498  if(vx3.Score < 0.001) vx3.Score = 0.001;
2499  if(vx3.Score > tcc.vtx2DCuts[7]) SetHighScoreBits(slc, vx3);
2500 
2501  } // SetVx3Score
2502 
2504  void SetVx2Score(TCSlice& slc)
2505  {
2506  // A version that sets the score of the last added vertex
2507  if(slc.vtxs.empty()) return;
2508  auto& vx2 = slc.vtxs[slc.vtxs.size() - 1];
2509  SetVx2Score(slc, vx2);
2510  } // SetVx2Score
2511 
2513  void SetVx2Score(TCSlice& slc, VtxStore& vx2)
2514  {
2515  // Calculate the 2D vertex score
2516  if(vx2.ID == 0) return;
2517 
2518  // Don't score vertices from CheckTrajBeginChg, MakeJunkVertices or Neutral vertices. Set to the minimum
2519  if(vx2.Topo == 8 || vx2.Topo == 9 || vx2.Topo == 11) {
2520  vx2.Score = tcc.vtx2DCuts[7] + 0.1;
2521  auto vtxTjID = GetVtxTjIDs(slc, vx2);
2522  vx2.TjChgFrac = ChgFracNearPos(slc, vx2.Pos, vtxTjID);
2523  return;
2524  }
2525 
2526  // Cuts on Tjs attached to vertices
2527  constexpr float maxChgRMS = 0.25;
2528  constexpr float momBin = 50;
2529 
2530  vx2.Score = -1000;
2531  vx2.TjChgFrac = 0;
2532  if(vx2.ID == 0) return;
2533  if(tcc.vtxScoreWeights.size() < 4) return;
2534 
2535  auto vtxTjIDs = GetVtxTjIDs(slc, vx2);
2536  if(vtxTjIDs.empty()) return;
2537 
2538  // Vertex position error
2539  float vpeScore = -tcc.vtxScoreWeights[0] * (vx2.PosErr[0] + vx2.PosErr[1]);
2540 
2541  unsigned short m3Dcnt = 0;
2542  if(vx2.Vx3ID > 0) {
2543  m3Dcnt = 1;
2544  // Add another if the 3D vertex is complete
2545  unsigned short ivx3 = vx2.Vx3ID - 1;
2546  if(slc.vtx3s[ivx3].Wire < 0) m3Dcnt = 2;
2547  }
2548  float m3DScore = tcc.vtxScoreWeights[1] * m3Dcnt;
2549 
2550  vx2.TjChgFrac = ChgFracNearPos(slc, vx2.Pos, vtxTjIDs);
2551  float cfScore = tcc.vtxScoreWeights[2] * vx2.TjChgFrac;
2552 
2553  // Define a weight for each Tj
2554  std::vector<int> tjids;
2555  std::vector<float> tjwts;
2556  unsigned short cnt13 = 0;
2557  for(auto tjid : vtxTjIDs) {
2558  Trajectory& tj = slc.tjs[tjid - 1];
2559  // Feb 22 Ignore short Tjs and junk tjs
2560  if(tj.AlgMod[kJunkTj]) continue;
2561  unsigned short lenth = tj.EndPt[1] - tj.EndPt[0] + 1;
2562  if(lenth < 3) continue;
2563  float wght = (float)tj.MCSMom / momBin;
2564  if(!tcc.useAlg[kNewVtxCuts] && wght > 10) wght = 10;
2565  // weight by the first tagged muon
2566  if(tj.PDGCode == 13) {
2567  ++cnt13;
2568  if(tcc.useAlg[kNewVtxCuts] && cnt13 == 1) wght *= 2;
2569  }
2570  // weight by charge rms
2571  if(tj.ChgRMS < maxChgRMS) ++wght;
2572  // Shower Tj
2573  if(tj.AlgMod[kShowerTj]) ++wght;
2574  // ShowerLike
2575  if(tj.AlgMod[kShowerLike]) --wght;
2576  tjids.push_back(tjid);
2577  tjwts.push_back(wght);
2578  } // tjid
2579 
2580  if(tjids.empty()) return;
2581 
2582  float tjScore = 0;
2583  float sum = 0;
2584  float cnt = 0;
2585  for(unsigned short it1 = 0; it1 < tjids.size() - 1; ++it1) {
2586  Trajectory& tj1 = slc.tjs[tjids[it1] - 1];
2587  float wght1 = tjwts[it1];
2588  // the end that has a vertex
2589  unsigned short end1 = 0;
2590  if(tj1.VtxID[1] == vx2.ID) end1 = 1;
2591  unsigned short endPt1 = tj1.EndPt[end1];
2592  // bump up the weight if there is a Bragg peak at the other end
2593  unsigned short oend1 = 1 - end1;
2594  if(tj1.StopFlag[oend1][kBragg]) ++wght1;
2595  float ang1 = tj1.Pts[endPt1].Ang;
2596  float ang1Err2 = tj1.Pts[endPt1].AngErr * tj1.Pts[endPt1].AngErr;
2597  for(unsigned short it2 = it1 + 1; it2 < tjids.size(); ++it2) {
2598  Trajectory& tj2 = slc.tjs[tjids[it2] - 1];
2599  float wght2 = tjwts[it2];
2600  unsigned end2 = 0;
2601  if(tj2.VtxID[1] == vx2.ID) end2 = 1;
2602  // bump up the weight if there is a Bragg peak at the other end
2603  unsigned short oend2 = 1 - end2;
2604  if(tj2.StopFlag[oend2][kBragg]) ++wght2;
2605  unsigned short endPt2 = tj2.EndPt[end2];
2606  float ang2 = tj2.Pts[endPt2].Ang;
2607  float ang2Err2 = tj2.Pts[endPt2].AngErr * tj2.Pts[endPt2].AngErr;
2608  float dang = DeltaAngle(ang1, ang2);
2609  float dangErr = 0.5 * sqrt(ang1Err2 + ang2Err2);
2610  if((dang / dangErr) > 3 && wght1 > 0 && wght2 > 0) {
2611  sum += wght1 + wght2;
2612  ++cnt;
2613  }
2614  } // it2
2615  } // it1
2616  if(cnt > 0) {
2617  sum /= cnt;
2618  tjScore = tcc.vtxScoreWeights[3] * sum;
2619  }
2620  vx2.Score = vpeScore + m3DScore + cfScore + tjScore;
2621  if(tcc.dbg2V && tcc.dbgSlc && (int)DecodeCTP(vx2.CTP).Plane == debug.Plane) {
2622  // last call after vertices have been matched to the truth. Use to optimize vtxScoreWeights using
2623  // an ntuple
2624  mf::LogVerbatim myprt("TC");
2625  myprt<<" SVx2W 2v"<<vx2.ID;
2626 /*
2627  myprt<<" m3Dcnt "<<m3Dcnt;
2628  myprt<<" PosErr "<<std::fixed<<std::setprecision(2)<<(vx2.PosErr[0] + vx2.PosErr[1]);
2629  myprt<<" TjChgFrac "<<std::fixed<<std::setprecision(3)<<vx2.TjChgFrac;
2630  myprt<<" sum "<<std::fixed<<std::setprecision(1)<<sum;
2631  myprt<<" cnt "<<(int)cnt;
2632 */
2633  myprt<<std::fixed<<std::setprecision(1);
2634  myprt<<" vpeScore "<<vpeScore<<" m3DScore "<<m3DScore;
2635  myprt<<" cfScore "<<cfScore<<" tjScore "<<tjScore;
2636  myprt<<" Score "<<vx2.Score;
2637  }
2638  } // SetVx2Score
2639 
2641  unsigned short Vx3Topo(TCSlice& slc, Vtx3Store& vx3)
2642  {
2643  // Returns the most common value of Topo for the 2D vertices that are matched
2644  // to this 3D vertex. This **might** be a useful measure to identify neutrino interaction
2645  // vertices
2646 
2647  if(vx3.ID == 0) return USHRT_MAX;
2648  // Consider Topo values between 0 and 9
2649  std::array<short, 10> cnts;
2650  cnts.fill(0);
2651  for(auto vx2id : vx3.Vx2ID) {
2652  if(vx2id == 0) continue;
2653  auto& vx2 = slc.vtxs[vx2id - 1];
2654  if(vx2.Topo < 0 || vx2.Topo > 9) continue;
2655  ++cnts[vx2.Topo];
2656  } // vx2id
2657  short most = 0;
2658  unsigned short theMost = USHRT_MAX;
2659  for(unsigned short itp = 0; itp < 10; ++itp) {
2660  if(cnts[itp] > most) {
2661  most = cnts[itp];
2662  theMost = itp;
2663  }
2664  } // itp
2665  return theMost;
2666  } // Vx3Topo
2667 
2669  void CompleteIncomplete3DVerticesInGaps(TCSlice& slc)
2670  {
2671 
2672  if(!tcc.useAlg[kComp3DVxIG]) return;
2673  if(slc.nPlanes != 3) return;
2674 
2675  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kComp3DVxIG]);
2676  if(prt) mf::LogVerbatim("TC")<<"Inside CI3DVIG:";
2677 
2678  for(unsigned short iv3 = 0; iv3 < slc.vtx3s.size(); ++iv3) {
2679  Vtx3Store& vx3 = slc.vtx3s[iv3];
2680  // ignore obsolete vertices
2681  if(vx3.ID == 0) continue;
2682  // check for a completed 3D vertex
2683  if(vx3.Wire < 0) continue;
2684  unsigned short mPlane = USHRT_MAX;
2685  for(unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
2686  if(vx3.Vx2ID[ipl] > 0) continue;
2687  mPlane = ipl;
2688  break;
2689  } // ipl
2690  if(mPlane == USHRT_MAX) continue;
2691 // CTP_t mCTP = EncodeCTP(vx3.CStat, vx3.TPC, mPlane);
2692  CTP_t mCTP = EncodeCTP(vx3.TPCID.Cryostat, vx3.TPCID.TPC, mPlane);
2693  // require that the missing vertex be in a large block of dead wires
2694  float dwc = DeadWireCount(slc, vx3.Wire - 4, vx3.Wire + 4, mCTP);
2695  if(dwc < 5) continue;
2696  // X position of the purported missing vertex
2697  // A TP for the missing 2D vertex
2698  VtxStore aVtx;
2699  aVtx.ID = slc.vtxs.size() + 1;
2700  aVtx.Pos[0] = vx3.Wire;
2701  aVtx.Pos[1] = tcc.detprop->ConvertXToTicks(vx3.X, mPlane, vx3.TPCID.TPC, vx3.TPCID.Cryostat) * tcc.unitsPerTick;
2702  aVtx.CTP = mCTP;
2703  aVtx.Topo = 4;
2704  aVtx.NTraj = 0;
2705  // Give it a bogus pass to indicate it wasn't created while stepping
2706  aVtx.Pass = 9;
2707  if(prt) mf::LogVerbatim("TC")<<"CI3DVIG: Incomplete vertex "<<iv3<<" in plane "<<mPlane<<" wire "<<vx3.Wire<<" Made 2D vertex ";
2708  std::vector<int> tjIDs;
2709  std::vector<unsigned short> tjEnds;
2710  for(unsigned short itj = 0; itj < slc.tjs.size(); ++itj) {
2711  if(slc.tjs[itj].CTP != mCTP) continue;
2712  if(slc.tjs[itj].AlgMod[kKilled] || slc.tjs[itj].AlgMod[kHaloTj]) continue;
2713  for(unsigned short end = 0; end < 2; ++end) {
2714  unsigned short ept = slc.tjs[itj].EndPt[end];
2715  TrajPoint& tp = slc.tjs[itj].Pts[ept];
2716  unsigned short oept = slc.tjs[itj].EndPt[1 - end];
2717  TrajPoint& otp = slc.tjs[itj].Pts[oept];
2718  // ensure that this is the end closest to the vertex
2719  if(std::abs(tp.Pos[0] - aVtx.Pos[0]) > std::abs(otp.Pos[0] - aVtx.Pos[0])) continue;
2720  float doca = PointTrajDOCA(slc, aVtx.Pos[0], aVtx.Pos[1], tp);
2721  if(doca > 2) continue;
2722  float dwc = DeadWireCount(slc, aVtx.Pos[0], tp.Pos[0], tp.CTP);
2723  float ptSep;
2724  if(aVtx.Pos[0] > tp.Pos[0]) {
2725  ptSep = aVtx.Pos[0] - tp.Pos[0] - dwc;
2726  } else {
2727  ptSep = tp.Pos[0] - aVtx.Pos[0] - dwc;
2728  }
2729  if(prt) mf::LogVerbatim("TC")<<"CI3DVIG: tj ID "<<slc.tjs[itj].ID<<" doca "<<doca<<" ptSep "<<ptSep;
2730  if(ptSep < -2 || ptSep > 2) continue;
2731  // don't clobber an existing association
2732  if(slc.tjs[itj].VtxID[end] > 0) continue;
2733  tjIDs.push_back(slc.tjs[itj].ID);
2734  tjEnds.push_back(end);
2735  } // end
2736  } // itj
2737  if(!tjIDs.empty()) {
2738  // Determine how messy this region is
2739  aVtx.TjChgFrac = ChgFracNearPos(slc, aVtx.Pos, tjIDs);
2740  if(aVtx.TjChgFrac < 0.7) continue;
2741  aVtx.Vx3ID = vx3.ID;
2742  // Save the 2D vertex
2743  if(!StoreVertex(slc, aVtx)) continue;
2744  for(unsigned short ii = 0; ii < tjIDs.size(); ++ii) {
2745  unsigned short itj = tjIDs[ii] - 1;
2746  slc.tjs[itj].VtxID[tjEnds[ii]] = aVtx.ID;
2747  slc.tjs[itj].AlgMod[kComp3DVxIG] = true;
2748  } // ii
2749  SetVx2Score(slc);
2750  vx3.Vx2ID[mPlane] = aVtx.ID;
2751  vx3.Wire = -1;
2752  if(prt) mf::LogVerbatim("TC")<<"CI3DVIG: new vtx 2V"<<aVtx.ID<<" points to 3V"<<vx3.ID;
2753  }
2754  } // vx3
2755 
2756  } // CompleteIncomplete3DVerticesInGaps
2757 
2759  void CompleteIncomplete3DVertices(TCSlice& slc)
2760  {
2761  // Look for trajectories in a plane that lack a 2D vertex as listed in
2762  // 2DVtxID that are near the projected wire. This may trigger splitting trajectories,
2763  // assigning them to a new 2D vertex and completing 3D vertices
2764 
2765  if(!tcc.useAlg[kComp3DVx]) return;
2766  if(slc.nPlanes != 3) return;
2767 
2768  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kComp3DVx]);
2769 
2770  float maxdoca = 3;
2771  if(prt) mf::LogVerbatim("TC")<<"Inside CI3DV with maxdoca set to "<<maxdoca;
2772  unsigned short ivx3 = 0;
2773  for(auto& vx3 : slc.vtx3s) {
2774  // ignore obsolete vertices
2775  if(vx3.ID == 0) continue;
2776  // check for a completed 3D vertex
2777  if(vx3.Wire < 0) continue;
2778  unsigned short mPlane = USHRT_MAX;
2779  // look for vertices in the induction plane in which the charge requirement wasn't imposed
2780  bool indPlnNoChgVtx = false;
2781  for(unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
2782  if(vx3.Vx2ID[plane] > 0) {
2783  auto& vx2 = slc.vtxs[vx3.Vx2ID[plane] - 1];
2784  if(vx2.Stat[kVtxIndPlnNoChg]) indPlnNoChgVtx = true;
2785  continue;
2786  }
2787  mPlane = plane;
2788  } // ipl
2789  if(mPlane == USHRT_MAX) continue;
2790  if(indPlnNoChgVtx) continue;
2791  CTP_t mCTP = EncodeCTP(vx3.TPCID.Cryostat, vx3.TPCID.TPC, mPlane);
2792  // X position of the purported missing vertex
2793  // A TP for the missing 2D vertex
2794  TrajPoint vtp;
2795  vtp.Pos[0] = vx3.Wire;
2796  vtp.Pos[1] = tcc.detprop->ConvertXToTicks(vx3.X, mPlane, vx3.TPCID.TPC, vx3.TPCID.Cryostat) * tcc.unitsPerTick;
2797  if(prt) mf::LogVerbatim("TC")<<"CI3DV 3V"<<vx3.ID<<" Pos "<<mPlane<<":"<<PrintPos(slc, vtp.Pos);
2798  std::vector<int> tjIDs;
2799  std::vector<unsigned short> tjPts;
2800  for(auto& tj : slc.tjs) {
2801  if(tj.CTP != mCTP) continue;
2802  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2803  if(tj.Pts.size() < 6) continue;
2804  if(tj.AlgMod[kComp3DVx]) continue;
2805  float doca = maxdoca;
2806  // find the closest distance between the vertex and the trajectory
2807  unsigned short closePt = 0;
2808  TrajPointTrajDOCA(slc, vtp, tj, closePt, doca);
2809  if(closePt > tj.EndPt[1]) continue;
2810  // try to improve the location of the vertex by looking for a distinctive feature on the
2811  // trajectory, e.g. high multiplicity hits or larger than normal charge
2812 // if(RefineVtxPosition(slc, tj, closePt, 3, false)) vtp.Pos = tj.Pts[closePt].Pos;
2813  if(prt) mf::LogVerbatim("TC")<<"CI3DV 3V"<<vx3.ID<<" candidate itj ID "<<tj.ID<<" vtx pos "<<PrintPos(slc, vtp.Pos)<<" doca "<<doca;
2814  tjIDs.push_back(tj.ID);
2815  tjPts.push_back(closePt);
2816  } // itj
2817  if(tjIDs.empty()) continue;
2818  // compare the length of the Tjs used to make the vertex with the length of the
2819  // Tj that we want to split. Don't allow a vertex using very short Tjs to split a long
2820  // Tj in the 3rd plane
2821 // float score;
2822 // auto vxtjs = GetVtxTjIDs(slc, vx3, score);
2823  auto vxtjs = GetAssns(slc, "3V", vx3.ID, "T");
2824  unsigned short maxPts = 0;
2825  for(auto tjid : vxtjs) {
2826  auto& tj = slc.tjs[tjid - 1];
2827  unsigned short npwc = NumPtsWithCharge(slc, tj, false);
2828  if(npwc > maxPts) maxPts = npwc;
2829  } // tjid
2830  // skip this operation if any of the Tjs in the split list are > 3 * maxPts
2831  maxPts *= 3;
2832  bool skipit = false;
2833  for(auto tjid : tjIDs) {
2834  auto& tj = slc.tjs[tjid - 1];
2835  if(NumPtsWithCharge(slc, tj, false) > maxPts) skipit = true;
2836  } // tjid
2837  if(prt) mf::LogVerbatim("TC")<<" maxPts "<<maxPts<<" skipit? "<<skipit;
2838  if(skipit) continue;
2839  // 2D vertex
2840  VtxStore aVtx;
2841  unsigned short newVtxIndx = slc.vtxs.size();
2842  aVtx.ID = newVtxIndx + 1;
2843  aVtx.CTP = mCTP;
2844  aVtx.Topo = 3;
2845  aVtx.NTraj = 0;
2846  // Give it a bogus pass to indicate it wasn't created while stepping
2847  aVtx.Pass = 9;
2848  aVtx.Pos = vtp.Pos;
2849  // ensure this isn't in a messy region
2850  aVtx.TjChgFrac = ChgFracNearPos(slc, aVtx.Pos, tjIDs);
2851  if(prt) mf::LogVerbatim("TC")<<" charge fraction near position "<<aVtx.TjChgFrac;
2852  if(aVtx.TjChgFrac < 0.6) continue;
2853  if(!StoreVertex(slc, aVtx)) continue;
2854  // make a reference to the new vertex
2855  VtxStore& newVtx = slc.vtxs[slc.vtxs.size()-1];
2856  if(prt) mf::LogVerbatim("TC")<<" Stored 2D vertex "<<newVtx.ID;
2857  // make a temporary copy so we can nudge it a bit if there is only one Tj
2858  std::array<float, 2> vpos = aVtx.Pos;
2859  for(unsigned short ii = 0; ii < tjIDs.size(); ++ii) {
2860  unsigned short itj = tjIDs[ii] - 1;
2861  // Don't use a reference variable since it may get scrambled after SplitTraj
2862 // auto& tj = slc.tjs[itj];
2863  unsigned short closePt = tjPts[ii];
2864  // determine which end is the closest
2865  unsigned short end = 1;
2866  // closest to the beginning?
2867  if(fabs(closePt - slc.tjs[itj].EndPt[0]) < fabs(closePt - slc.tjs[itj].EndPt[1])) end = 0;
2868  short dpt = fabs(closePt - slc.tjs[itj].EndPt[end]);
2869  if(dpt < 3) {
2870  // close to an end
2871  if(slc.tjs[itj].VtxID[end] > 0) {
2872  // find the distance btw the existing vertex and the end of this tj
2873  auto& oldTj = slc.tjs[itj];
2874  auto& oldVx = slc.vtxs[oldTj.VtxID[end] - 1];
2875  short oldSep = fabs(oldVx.Pos[0] - oldTj.Pts[oldTj.EndPt[end]].Pos[0]);
2876  if(prt) mf::LogVerbatim("TC")<<" T"<<slc.tjs[itj].ID<<" has vertex 2V"<<slc.tjs[itj].VtxID[end]<<" at end "<<end<<". oldSep "<<oldSep;
2877  if(dpt < oldSep) {
2878  MakeVertexObsolete("CI3DV", slc, oldVx, true);
2879  } else {
2880  continue;
2881  }
2882  } // slc.tjs[itj].VtxID[end] > 0
2883  slc.tjs[itj].VtxID[end] = slc.vtxs[newVtxIndx].ID;
2884  ++newVtx.NTraj;
2885  if(prt) mf::LogVerbatim("TC")<<" attach Traj T"<<slc.tjs[itj].ID<<" at end "<<end;
2886  slc.tjs[itj].AlgMod[kComp3DVx] = true;
2887  vpos = slc.tjs[itj].Pts[slc.tjs[itj].EndPt[end]].Pos;
2888  } else {
2889  // closePt is not near an end, so split the trajectory
2890  if(SplitTraj(slc, itj, closePt, newVtxIndx, prt)) {
2891  if(prt) mf::LogVerbatim("TC")<<" SplitTraj success 2V"<<slc.vtxs[newVtxIndx].ID<<" at closePt "<<closePt;
2892  // successfully split the Tj
2893  newVtx.NTraj += 2;
2894  } else {
2895  // split failed. Give up
2896  if(prt) mf::LogVerbatim("TC")<<" SplitTraj failed";
2897  newVtx.NTraj = 0;
2898  break;
2899  }
2900  // Update the PDGCode for the chopped trajectory
2901  SetPDGCode(slc, itj, true);
2902  // and for the new trajectory
2903  SetPDGCode(slc, slc.tjs.size()-1, true);
2904  } // closePt is not near an end, so split the trajectory
2905  slc.tjs[itj].AlgMod[kComp3DVx] = true;
2906  unsigned short newtj = slc.tjs.size() - 1;
2907  slc.tjs[newtj].AlgMod[kComp3DVx] = true;
2908  } // ii
2909  if(newVtx.NTraj == 0) {
2910  // A failure occurred. Recover
2911  if(prt) mf::LogVerbatim("TC")<<" Failed. Recover and delete vertex "<<newVtx.ID;
2912  MakeVertexObsolete("CI3DV", slc, newVtx, true);
2913  } else {
2914  // success
2915  vx3.Vx2ID[mPlane] = newVtx.ID;
2916  newVtx.Vx3ID = vx3.ID;
2917  vx3.Wire = -1;
2918  // set the vertex position to the start of the Tj if there is only one and fix it
2919  if(newVtx.NTraj == 1) {
2920  newVtx.Pos = vpos;
2921  newVtx.Stat[kFixed] = true;
2922  }
2923  AttachAnyTrajToVertex(slc, newVtx.ID - 1, prt);
2924  SetVx2Score(slc);
2925  if(prt) {
2926  mf::LogVerbatim myprt("TC");
2927  myprt<<" Success: new 2V"<<newVtx.ID<<" at "<<(int)newVtx.Pos[0]<<":"<<(int)newVtx.Pos[1]/tcc.unitsPerTick;
2928  myprt<<" points to 3V"<<vx3.ID;
2929  myprt<<" TjIDs:";
2930  for(auto& tjID : tjIDs) myprt<<" T"<<std::to_string(tjID);
2931  } // prt
2932  } // success
2933  ++ivx3;
2934  } // vx3
2935 
2936  } // CompleteIncomplete3DVertices
2937 /*
2939  bool RefineVtxPosition(TCSlice& slc, const Trajectory& tj, unsigned short& nearPt, short nPtsToChk, bool prt)
2940  {
2941  // The tj has been slated to be split somewhere near point nearPt. This function will move
2942  // the near point a bit to the most likely point of a vertex
2943 
2944  float maxChg = tj.Pts[nearPt].Chg;
2945  short maxChgPt = nearPt;
2946 
2947  for(short ipt = nearPt - nPtsToChk; ipt < nearPt + nPtsToChk; ++ipt) {
2948  if(ipt < tj.EndPt[0] || ipt > tj.EndPt[1]) continue;
2949  auto& tp = tj.Pts[ipt];
2950  if(tp.Chg > maxChg) {
2951  maxChg = tp.Chg;
2952  maxChgPt = ipt;
2953  }
2954  if(prt) mf::LogVerbatim("TC")<<"RVP: ipt "<<ipt<<" Pos "<<tp.CTP<<":"<<PrintPos(slc, tp.Pos)<<" chg "<<(int)tp.Chg<<" nhits "<<tp.Hits.size();
2955  } // ipt
2956  if(nearPt == maxChgPt) return false;
2957  nearPt = maxChgPt;
2958  return true;
2959  } //RefineVtxPosition
2960 */
2962  void VtxHitsSwap(TCSlice& slc, const CTP_t inCTP){
2963 
2964  if(!tcc.useAlg[kVtxHitsSwap]) return;
2965 
2966  geo::PlaneID planeID = DecodeCTP(inCTP);
2967 
2968  bool prt = ((tcc.dbg2V && debug.Plane == (int)planeID.Plane) || tcc.dbgAlg[kVtxHitsSwap]);
2969  for (unsigned short iv = 0; iv < slc.vtxs.size(); ++iv){
2970  VtxStore& rvx = slc.vtxs[iv];
2971  if(rvx.CTP != inCTP) continue;
2972  // Only consider vertex with two trajectories
2973  if(rvx.NTraj != 2) continue;
2974  if (prt) mf::LogVerbatim("TC")<<"VtxHitsSwap: 2V"<<rvx.ID<<" Pos[0]: "<<rvx.Pos[0]<<" "<<rvx.Pos[1];
2975  std::array<unsigned short, 2> tjlist{{0,0}};
2976  for(unsigned short itj = 0; itj < slc.tjs.size(); ++itj) {
2977  if(slc.tjs[itj].AlgMod[kKilled] || slc.tjs[itj].AlgMod[kHaloTj]) continue;
2978  if(slc.tjs[itj].CTP != rvx.CTP) continue;
2979  Trajectory& tj = slc.tjs[itj];
2980  // ensure that the ID is OK so the code below doesn't choke
2981  for(unsigned short end = 0; end < 2; ++end) {
2982  if(tj.VtxID[end] == rvx.ID) {
2983  tjlist[end] = itj;
2984  }
2985  }
2986  }//all trajectories
2987 
2988  //Ignore short trajectories
2989  if (slc.tjs[tjlist[0]].EndPt[1]<5||
2990  slc.tjs[tjlist[1]].EndPt[1]<5) continue;
2991 
2992  for (unsigned short i = 0; i<2; ++i){
2993 
2994  //First check if first hit should be swapped
2995  Trajectory& tj0 = slc.tjs[tjlist[i]];
2996  Trajectory& tj1 = slc.tjs[tjlist[1-i]];
2997  unsigned short endPt0 = slc.tjs[tjlist[i]].EndPt[i];
2998  unsigned short endPt1 = slc.tjs[tjlist[1-i]].EndPt[1-i];
2999  //first TP on first trajectory
3000  float chg0 = TpSumHitChg(slc, tj0.Pts[endPt0]);
3001  float w0 = tj0.Pts[endPt0].Pos[0];
3002  //if (vtxPrt) mf::LogVerbatim("TC")<<PrintPos(slc, tj0.Pts[endPt0]);
3003  //second TP on first trajectory
3004  float chg1 = 0;
3005  float w1 = 0;
3006  unsigned short j = endPt0;
3007  while (j!=tj0.EndPt[1-i]){
3008  if (i==0) ++j;
3009  else --j;
3010  if (tj0.Pts[j].Chg){
3011  chg1 = TpSumHitChg(slc, tj0.Pts[j]);
3012  w1 = tj0.Pts[j].Pos[0];
3013  //if (vtxPrt) mf::LogVerbatim("TC")<<PrintPos(slc, tj0.Pts[j]);
3014  break;
3015  }
3016  }
3017  //first TP on second trajectory
3018  float chg2 = TpSumHitChg(slc,tj1.Pts[endPt1]);
3019  float w2 = tj1.Pts[endPt1].Pos[0];
3020  //DOCA between first TP on first TJ and first TP on second TJ
3021  float delta = 1000;
3022  for (size_t k = 0; k<tj0.Pts[endPt0].Hits.size(); ++k){
3023  if (!tj0.Pts[endPt0].UseHit[k]) continue;
3024  float this_delta = PointTrajDOCA(slc, tj0.Pts[endPt0].Hits[k], tj1.Pts[endPt1]);
3025  if (this_delta<delta) delta = this_delta;
3026  }
3027  // if (vtxPrt) mf::LogVerbatim("TC")<<PrintPos(slc, tj1.Pts[endPt1]);
3028  //if (vtxPrt) mf::LogVerbatim("TC")<<chg0<<" "<<chg1<<" "<<chg2<<" "<<delta;
3029  if (chg0>0&&std::abs((chg0-chg1)/chg0)-std::abs((chg0-chg2)/chg0)>0.2&&delta<1.5&&std::abs(w2-w1)<1.5){
3030  if (prt) {
3031  mf::LogVerbatim("TC")<<" chg0 = "<<chg0<<" chg1 = "<<chg1<<" chg2 = "<<chg2<<" delta = "<<delta<<" w0 = "<<w0<<" w1 = "<<w1<<" w2 = "<<w2;
3032  mf::LogVerbatim("TC")<<" Moving TP "<<PrintPos(slc, tj0.Pts[endPt0])<<" from TJ "<<tj0.ID<<" to TJ "<<tj1.ID;
3033  }
3034  //Add first TP of first trajectory to the second trajectory
3035  TrajPoint tp = tj0.Pts[endPt0];
3036  for (size_t j = 0; j<tp.Hits.size(); ++j){
3037  if (!tp.UseHit[j]) continue;
3038  slc.slHits[tp.Hits[j]].InTraj = tj1.ID;
3039  }
3040  if (i==0){
3041  //append to the end
3042  tj1.Pts.push_back(tp);
3043  }
3044  else{
3045  //insert at the beginning
3046  tj1.Pts.insert(tj1.Pts.begin(), tp);
3047  }
3048  SetEndPoints(tj1);
3049 
3050  //Remove first TP from first trajectory
3051  tj0.Pts[endPt0].Chg = 0;
3052  for (size_t j = 0; j<tj0.Pts[endPt0].Hits.size(); ++j){
3053  tj0.Pts[endPt0].UseHit[j] = false;
3054  }
3055  SetEndPoints(tj0);
3056  tj0.AlgMod[kVtxHitsSwap] = true;
3057  tj1.AlgMod[kVtxHitsSwap] = true;
3058  break;
3059  }
3060 
3061  //Now Check if the beginning of the first trajectory should be moved to the second trajectory.
3062  j = endPt0;
3063  std::vector<unsigned short> tplist;
3064  while (j!=tj0.EndPt[1-i]){
3065  if (tj0.Pts[j].Chg){
3066  float delta = 1000;
3067  for (size_t k = 0; k<tj0.Pts[j].Hits.size(); ++k){
3068  if (!tj0.Pts[j].UseHit[k]) continue;
3069  float this_delta = PointTrajDOCA(slc, tj0.Pts[j].Hits[k], tj1.Pts[endPt1]);
3070  if (this_delta<delta) delta = this_delta;
3071  //if (vtxPrt) mf::LogVerbatim("TC")<<j<<" "<<k<<" "<<PrintPos(slc, tj0.Pts[j])<<" "<<PointTrajDOCA(slc, tj0.Pts[j].Hits[k], tj1.Pts[endPt1]);
3072  }
3073  if (delta < 0.3 && tj0.Pts[j].Delta > 1.0 && (j==endPt0 || !tplist.empty())){
3074  tplist.push_back(j);
3075  }
3076  else break;
3077  }
3078  if (i==0) ++j;
3079  else --j;
3080  }
3081  //if (vtxPrt) mf::LogVerbatim("TC")<<tplist.size();
3082  //Need at least two TPs to make this change
3083  if (tplist.size()>1){
3084  if (prt) mf::LogVerbatim("TC")<<" Moving "<<tplist.size()<<" TPs from TJ "<<tj0.ID<<" to TJ "<<tj1.ID;
3085  //Append TPs to the second TJ
3086  for (unsigned short j = 0; j<tplist.size(); ++j){
3087  TrajPoint tp = tj0.Pts[tplist[j]];
3088  for (size_t k = 0; k<tp.Hits.size(); ++k){
3089  if (!tp.UseHit[k]) continue;
3090  slc.slHits[tp.Hits[k]].InTraj = tj1.ID;
3091  }
3092  if (i==0){
3093  //append to the end
3094  tj1.Pts.push_back(tp);
3095  }
3096  else{
3097  //insert at the beginning
3098  tj1.Pts.insert(tj1.Pts.begin(), tp);
3099  }
3100  }
3101  SetEndPoints(tj1);
3102 
3103  //Remove TPs from first trajectory
3104  for (unsigned short j = 0; j<tplist.size(); ++j){
3105  tj0.Pts[tplist[j]].Chg = 0;
3106  for (size_t k = 0; k<tj0.Pts[tplist[j]].Hits.size(); ++k){
3107  tj0.Pts[tplist[j]].UseHit[k] = false;
3108  }
3109  }
3110  SetEndPoints(tj0);
3111  tj0.AlgMod[kVtxHitsSwap] = true;
3112  tj1.AlgMod[kVtxHitsSwap] = true;
3113  break;
3114  }
3115  }//loop over two trajectories
3116  }//loop over vertices
3117  }
3118 
3120  bool MakeVertexObsolete(std::string fcnLabel, TCSlice& slc, VtxStore& vx2, bool forceKill)
3121  {
3122  // Makes a 2D vertex obsolete
3123 
3124  // check for a high-score 3D vertex
3125  bool hasHighScoreVx3 = (vx2.Vx3ID > 0);
3126  if(hasHighScoreVx3 && !forceKill && slc.vtx3s[vx2.Vx3ID - 1].Score >= tcc.vtx2DCuts[7]) return false;
3127 
3128  if(tcc.dbg2V || tcc.dbg3V) {
3129  mf::LogVerbatim("TC")<<fcnLabel<<" MVO: killing 2V"<<vx2.ID;
3130  }
3131 
3132  // Kill it
3133  int vx2id = vx2.ID;
3134  if(vx2.Vx3ID > 0) {
3135  auto& vx3 = slc.vtx3s[vx2.Vx3ID - 1];
3136  for(auto& v3v2id : vx3.Vx2ID) if(v3v2id == vx2.ID) v3v2id = 0;
3137  }
3138  vx2.ID = 0;
3139  for(auto& tj : slc.tjs) {
3140  if(tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
3141  for(unsigned short end = 0; end < 2; ++end) {
3142  if(tj.VtxID[end] != vx2id) continue;
3143  tj.VtxID[end] = 0;
3144  tj.AlgMod[kPhoton] = false;
3145  // clear the kEnvOverlap bits on the TPs
3146  for(unsigned short ii = 0; ii < tj.Pts.size(); ++ii) {
3147  if(end == 0) {
3148  unsigned short ipt = tj.EndPt[0] + ii;
3149  auto& tp = tj.Pts[ipt];
3150  if(!tp.Environment[kEnvOverlap]) break;
3151  if(ipt == tj.EndPt[1]) break;
3152  } else {
3153  unsigned short ipt = tj.EndPt[1] - ii;
3154  auto& tp = tj.Pts[ipt];
3155  if(!tp.Environment[kEnvOverlap]) break;
3156  if(ipt == tj.EndPt[0]) break;
3157  }
3158  } // ii
3159  if(tj.AlgMod[kTjHiVx3Score]) {
3160  // see if the vertex at the other end is high-score and if so, preserve the state
3161  unsigned short oend = 1 - end;
3162  if(tj.VtxID[oend] > 0) {
3163  auto& ovx2 = slc.vtxs[tj.VtxID[oend] - 1];
3164  if(!ovx2.Stat[kHiVx3Score]) tj.AlgMod[kTjHiVx3Score] = false;
3165  } // vertex at the other end
3166  } // tj.AlgMod[kTjHiVx3Score]
3167  } // end
3168  } // tj
3169 
3170  if(!hasHighScoreVx3) return true;
3171 
3172  // update the affected 3D vertex
3173  Vtx3Store& vx3 = slc.vtx3s[vx2.Vx3ID - 1];
3174  // make the 3D vertex incomplete
3175  geo::PlaneID planeID = DecodeCTP(vx2.CTP);
3176  unsigned short plane = planeID.Plane;
3177  if(vx3.Vx2ID[plane] != vx2id) return true;
3178  vx3.Vx2ID[plane] = 0;
3179  vx3.Wire = vx2.Pos[0];
3180  // Ensure that there are at least two 2D vertices left
3181  unsigned short n2D = 0;
3182  for(unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
3183  if(vx3.Vx2ID[plane] > 0) ++n2D;
3184  } // plane
3185 
3186  if(n2D > 1) {
3187  // 3D vertex is incomplete
3188  // correct the score
3189  SetVx3Score(slc, vx3);
3190  return true;
3191  }
3192 
3193  // 3D vertex is obsolete
3194  // Detach the all remaining 2D vertices from the 3D vertex
3195  for(auto& vx2 : slc.vtxs) {
3196  if(vx2.ID == 0) continue;
3197  if(vx2.Vx3ID == vx3.ID) vx2.Vx3ID = 0;
3198  } // vx2
3199  for(auto& pfp : slc.pfps) {
3200  for(unsigned short end = 0; end < 2; ++end) if(pfp.Vx3ID[end] == vx3.ID) pfp.Vx3ID[end] = 0;
3201  } // pfp
3202  vx3.ID = 0;
3203  return true;
3204 
3205  } // MakeVertexObsolete
3206 
3208  bool MakeVertexObsolete(TCSlice& slc, Vtx3Store& vx3)
3209  {
3210  // Deletes a 3D vertex and 2D vertices in all planes
3211  // The 2D and 3D vertices are NOT killed if forceKill is false and the 3D vertex
3212  // has a high score
3213  if(vx3.ID == 0) return true;
3214  if(vx3.ID > int(slc.vtx3s.size())) return false;
3215 
3216  for(auto vx2id : vx3.Vx2ID) {
3217  if(vx2id == 0 || vx2id > (int)slc.vtxs.size()) continue;
3218  auto& vx2 = slc.vtxs[vx2id - 1];
3219  MakeVertexObsolete("MVO3", slc, vx2, true);
3220  }
3221  vx3.ID = 0;
3222  return true;
3223  } // MakeVertexObsolete
3224 
3226  std::vector<int> GetVtxTjIDs(const TCSlice& slc, const VtxStore& vx2)
3227  {
3228  // returns a list of trajectory IDs that are attached to vx2
3229  std::vector<int> tmp;
3230  if(vx2.ID == 0) return tmp;
3231  for(auto& tj : slc.tjs) {
3232  if(tj.AlgMod[kKilled]) continue;
3233  if(tj.CTP != vx2.CTP) continue;
3234  for(unsigned short end = 0; end < 2; ++end) {
3235  if(tj.VtxID[end] == vx2.ID) tmp.push_back(tj.ID);
3236  } // end
3237  } // tj
3238  return tmp;
3239  } // GetVtxTjIDs
3240 
3241 
3243  std::vector<int> GetVtxTjIDs(const TCSlice& slc, const Vtx3Store& vx3, float& score)
3244  {
3245  // returns a list of Tjs in all planes that are attached to vx3
3246  std::vector<int> tmp;
3247  if(vx3.ID == 0) return tmp;
3248  float nvx2 = 0;
3249  score = 0;
3250  for(auto& vx2 : slc.vtxs) {
3251  if(vx2.ID == 0) continue;
3252  if(vx2.Vx3ID != vx3.ID) continue;
3253  auto vtxTjID2 = GetVtxTjIDs(slc, vx2);
3254  tmp.insert(tmp.end(), vtxTjID2.begin(), vtxTjID2.end());
3255  score += vx2.Score;
3256  ++nvx2;
3257  } // vx2
3258  if(nvx2 < 1) return tmp;
3259  // find the average score
3260  score /= nvx2;
3261  // sort by increasing ID
3262  std::sort(tmp.begin(), tmp.end());
3263  return tmp;
3264  } // GetVtxTjIDs
3265 
3267  std::vector<unsigned short> GetPFPVertices(const TCSlice& slc, const PFPStruct& pfp)
3268  {
3269  // returns a list of 3D vertices that are attached to Tjs in this pfp. No check is
3270  // made of the actual vertex attachment of the pfp.
3271  std::vector<unsigned short> tmp;
3272  if(pfp.TjIDs.empty()) return tmp;
3273  for(auto tjid : pfp.TjIDs) {
3274  auto& tj = slc.tjs[tjid - 1];
3275  for(unsigned short end = 0; end < 2; ++end) {
3276  if(tj.VtxID[end] == 0) continue;
3277  auto& vx2 = slc.vtxs[tj.VtxID[end] - 1];
3278  if(vx2.Vx3ID == 0) continue;
3279  if(std::find(tmp.begin(), tmp.end(), vx2.Vx3ID) != tmp.end()) continue;
3280  tmp.push_back(vx2.Vx3ID);
3281  } // end
3282  } // tjid
3283  return tmp;
3284  } // GetPFPVertices
3285 
3287  void PosInPlane(const TCSlice& slc, const Vtx3Store& vx3, unsigned short plane, Point2_t& pos)
3288  {
3289  // returns the 2D position of the vertex in the plane
3290  pos[0] = tcc.geom->WireCoordinate(vx3.Y, vx3.Z, plane, vx3.TPCID.TPC, vx3.TPCID.Cryostat);
3291  pos[1] = tcc.detprop->ConvertXToTicks(vx3.X, plane, vx3.TPCID.TPC, vx3.TPCID.Cryostat) * tcc.unitsPerTick;
3292 
3293  } // PosInPlane
3294 
3295 
3297  unsigned short IsCloseToVertex(TCSlice& slc, VtxStore& inVx2)
3298  {
3299  // Returns the ID of a 2D vertex having the minimum pull < user-specified cut
3300 
3301  float minPull = tcc.vtx2DCuts[3];
3302  unsigned short imBest = 0;
3303  for(auto& vx2 : slc.vtxs) {
3304  if(vx2.CTP != inVx2.CTP) continue;
3305  if(vx2.ID <= 0) continue;
3306  float pull = VertexVertexPull(slc, inVx2, vx2);
3307  if(pull < minPull) {
3308  minPull = pull;
3309  imBest = vx2.ID;
3310  }
3311  } // vx2
3312  return imBest;
3313  } // IsCloseToVertex
3314 
3316  unsigned short IsCloseToVertex(TCSlice& slc, Vtx3Store& vx3)
3317  {
3318  // Returns the ID of a 3D vertex having the minimum pull < user-specified cut
3319 
3320  float minPull = tcc.vtx3DCuts[1];
3321  unsigned short imBest = 0;
3322  for(auto& oldvx3 : slc.vtx3s) {
3323  if(oldvx3.ID == 0) continue;
3324  if(std::abs(oldvx3.X - vx3.X) > tcc.vtx3DCuts[0]) continue;
3325  float pull = VertexVertexPull(slc, vx3, oldvx3);
3326  if(pull < minPull) {
3327  minPull = pull;
3328  imBest = oldvx3.ID;
3329  }
3330  } // oldvx3
3331  return imBest;
3332 
3333  } // IsCloseToVertex
3334 
3335 } // namespace
tca::TrajPoint::UseHit
std::bitset< 16 > UseHit
Definition: DataStructs.h:137
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1189
tca::valDecreasing
bool valDecreasing(SortEntry c1, SortEntry c2)
Definition: TCVertex.cxx:10
tca::TrajPoint::Dir
Vector2_t Dir
Definition: DataStructs.h:122
tca::ReleaseHits
void ReleaseHits(TCSlice &slc, Trajectory &tj)
Definition: Utils.cxx:1150
tca::TrajPoint::Pos
Point2_t Pos
Definition: DataStructs.h:121
tca::VtxStore::PosErr
Point2_t PosErr
Definition: DataStructs.h:56
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:219
tca::TCSlice::tjs
std::vector< Trajectory > tjs
vector of all trajectories in each plane
Definition: DataStructs.h:570
tca::TCSlice::vtx3s
std::vector< Vtx3Store > vtx3s
3D vertices
Definition: DataStructs.h:577
tca::StepAway
void StepAway(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:15
tca::VtxHitsSwap
void VtxHitsSwap(TCSlice &slc, const CTP_t inCTP)
Definition: TCVertex.cxx:2962
tca::MergeWithVertex
bool MergeWithVertex(TCSlice &slc, VtxStore &vx, unsigned short oVxID)
Definition: TCVertex.cxx:701
tca::PDGCodeVote
int PDGCodeVote(TCSlice &slc, std::vector< int > &tjIDs, bool prt)
Definition: Utils.cxx:372
tca::SetEndPoints
void SetEndPoints(Trajectory &tj)
Definition: Utils.cxx:2864
tca::MakeVertexObsolete
bool MakeVertexObsolete(std::string fcnLabel, TCSlice &slc, VtxStore &vx2, bool forceKill)
Definition: TCVertex.cxx:3120
tca::TCConfig::kinkCuts
std::vector< float > kinkCuts
kink angle, nPts fit, (alternate) kink angle significance
Definition: DataStructs.h:475
tca::InTrajOK
bool InTrajOK(TCSlice &slc, std::string someText)
Definition: Utils.cxx:1333
tca::VtxStore
struct of temporary 2D vertices (end points)
Definition: DataStructs.h:54
tca::GetAssns
std::vector< int > GetAssns(TCSlice &slc, std::string type1Name, int id, std::string type2Name)
Definition: Utils.cxx:4279
tca::kVtxTrjTried
FindVtxTraj algorithm tried.
Definition: DataStructs.h:73
tca::TrajPointTrajDOCA
void TrajPointTrajDOCA(TCSlice &slc, TrajPoint const &tp, Trajectory const &tj, unsigned short &closePt, float &minSep)
Definition: Utils.cxx:2155
evd::kWire
static unsigned int kWire
Definition: TWQMultiTPCProjection.cxx:52
tca::TrajPoint::CTP
CTP_t CTP
Cryostat, TPC, Plane code.
Definition: DataStructs.h:119
tca::TCConfig::maxPos0
std::vector< float > maxPos0
Definition: DataStructs.h:489
tca::FitVertex
bool FitVertex(TCSlice &slc, VtxStore &vx, bool prt)
Definition: TCVertex.cxx:2131
tca::TCConfig::vtx3DCuts
std::vector< float > vtx3DCuts
2D vtx -> 3D vtx matching cuts
Definition: DataStructs.h:467
tca::MCSMom
short MCSMom(TCSlice &slc, const std::vector< int > &tjIDs)
Definition: Utils.cxx:2914
mf::MaybeLogger_
Definition: MessageLogger.h:114
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:6
tca::PrintPFP
void PrintPFP(std::string someText, TCSlice &slc, const PFPStruct &pfp, bool printHeader)
Definition: Utils.cxx:5587
tca::kFixed
vertex position fixed manually - no fitting done
Definition: DataStructs.h:74
y
Float_t y
Definition: compare.C:6
tca::Trajectory::StopFlag
std::array< std::bitset< 8 >, 2 > StopFlag
Definition: DataStructs.h:165
tca::SplitTrajCrossingVertices
void SplitTrajCrossingVertices(TCSlice &slc, CTP_t inCTP)
Definition: TCVertex.cxx:1268
z
Double_t z
Definition: plot.C:279
tca::Vtx3Store::Vx2ID
std::vector< int > Vx2ID
Definition: DataStructs.h:94
geo::kX
Planes which measure X direction.
Definition: geo_types.h:81
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:250
tca::valIncreasing
bool valIncreasing(SortEntry c1, SortEntry c2)
Definition: TCVertex.cxx:11
tca::kHiVx3Score
matched to a high-score 3D vertex
Definition: DataStructs.h:76
tca::TrajPoint::Hits
std::vector< unsigned int > Hits
Definition: DataStructs.h:136
tca::kStepDir
step from US -> DS (true) or DS -> US (false)
Definition: DataStructs.h:447
tca::StoreTraj
bool StoreTraj(TCSlice &slc, Trajectory &tj)
Definition: Utils.cxx:1175
tca::VtxStore::Pass
unsigned short Pass
Definition: DataStructs.h:58
tca::TrajPointVertexPull
float TrajPointVertexPull(TCSlice &slc, const TrajPoint &tp, const VtxStore &vx)
Definition: TCVertex.cxx:2029
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:130
tca::PrintAllTraj
void PrintAllTraj(std::string someText, TCSlice &slc, unsigned short itj, unsigned short ipt, bool prtVtx)
Definition: Utils.cxx:5197
tca::TrajClosestApproach
bool TrajClosestApproach(Trajectory const &tj, float x, float y, unsigned short &closePt, float &DOCA)
Definition: Utils.cxx:2377
tca::AttachAnyTrajToVertex
bool AttachAnyTrajToVertex(TCSlice &slc, unsigned short ivx, bool prt)
Definition: TCVertex.cxx:1868
tmp
Float_t tmp
Definition: plot.C:37
tca::DeltaAngle
double DeltaAngle(const Vector3_t v1, const Vector3_t v2)
Definition: PFPUtils.cxx:1636
tca::TCConfig::dbgSlc
bool dbgSlc
debug only in the user-defined slice? default is all slices
Definition: DataStructs.h:507
tca::StartTraj
bool StartTraj(TCSlice &slc, Trajectory &tj, unsigned int fromhit, unsigned int tohit, unsigned short pass)
Definition: Utils.cxx:4431
tca::PFPStruct::Vx3ID
std::array< int, 2 > Vx3ID
Definition: DataStructs.h:217
tca::SetVx3Score
void SetVx3Score(TCSlice &slc, Vtx3Store &vx3)
Definition: TCVertex.cxx:2483
tca::AttachPFPToVertex
bool AttachPFPToVertex(TCSlice &slc, PFPStruct &pfp, unsigned short end, unsigned short vx3ID, bool prt)
Definition: TCVertex.cxx:1830
tca::Trajectory::IsGood
bool IsGood
set false if there is a failure or the Tj fails quality cuts
Definition: DataStructs.h:168
tca::DeadWireCount
float DeadWireCount(TCSlice &slc, const TrajPoint &tp1, const TrajPoint &tp2)
Definition: Utils.cxx:1881
tca::StoreVertex
bool StoreVertex(TCSlice &slc, VtxStore &vx)
Definition: TCVertex.cxx:2093
tca::Split3DKink
bool Split3DKink(TCSlice &slc, PFPStruct &pfp, double sep, bool prt)
Definition: PFPUtils.cxx:1811
tca::SetIntersection
std::vector< T > SetIntersection(const std::vector< T > &set1, const std::vector< T > &set2)
Definition: Utils.h:231
tca::TrajTrajDOCA
bool TrajTrajDOCA(TCSlice &slc, const Trajectory &tj1, const Trajectory &tj2, unsigned short &ipt1, unsigned short &ipt2, float &minSep)
Definition: Utils.cxx:2175
detinfo::DetectorProperties::ConvertXToTicks
virtual double ConvertXToTicks(double X, int p, int t, int c) const =0
tca::Trajectory::Pass
unsigned short Pass
the pass on which it was created
Definition: DataStructs.h:161
tca::TCConfig::dbg3V
bool dbg3V
debug 3D vertex finding
Definition: DataStructs.h:514
tca::FindCloseHits
std::vector< unsigned int > FindCloseHits(TCSlice &slc, std::array< int, 2 > const &wireWindow, Point2_t const &timeWindow, const unsigned short plane, HitStatus_t hitRequest, bool usePeakTime, bool &hitsNear)
Definition: Utils.cxx:2496
tca::Vtx3Store
struct of temporary 3D vertices
Definition: DataStructs.h:84
c1
TCanvas * c1
Definition: plotHisto.C:7
tca::GetPFPVertices
std::vector< unsigned short > GetPFPVertices(const TCSlice &slc, const PFPStruct &pfp)
Definition: TCVertex.cxx:3267
c2
TCanvas * c2
Definition: plot_hist.C:75
tca::Point2_t
std::array< float, 2 > Point2_t
Definition: DataStructs.h:37
tca::TCConfig::maxPos1
std::vector< float > maxPos1
Definition: DataStructs.h:490
tca::TCConfig::unitsPerTick
float unitsPerTick
scale factor from Tick to WSE equivalent units
Definition: DataStructs.h:488
tca::PFPStruct::XYZ
std::array< Point3_t, 2 > XYZ
Definition: DataStructs.h:212
tca::TCConfig::detprop
const detinfo::DetectorProperties * detprop
Definition: DataStructs.h:494
tca::NearestPtWithChg
unsigned short NearestPtWithChg(TCSlice &slc, Trajectory &tj, unsigned short thePt)
Definition: Utils.cxx:2966
tca::Vtx3Store::TPCID
geo::TPCID TPCID
Definition: DataStructs.h:93
tca::debug
DebugStuff debug
Definition: DebugStruct.cxx:4
tca::Trajectory::CTP
CTP_t CTP
Cryostat, TPC, Plane code.
Definition: DataStructs.h:144
tca::SignalBetween
bool SignalBetween(TCSlice &slc, const TrajPoint &tp1, const TrajPoint &tp2, const float &MinWireSignalFraction)
Definition: Utils.cxx:1671
tca::TCConfig::dbg2V
bool dbg2V
debug 2D vertex finding
Definition: DataStructs.h:510
tca::Trajectory::Pts
std::vector< TrajPoint > Pts
Trajectory points.
Definition: DataStructs.h:143
tca::ChkVtxAssociations
bool ChkVtxAssociations(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:2326
tca::DebugStuff::Plane
int Plane
Select plane.
Definition: DebugStruct.h:21
tca::PFPStruct::Dir
std::array< Vector3_t, 2 > Dir
Definition: DataStructs.h:213
tca::DefinePFP
bool DefinePFP(std::string inFcnLabel, TCSlice &slc, PFPStruct &pfp, bool prt)
Definition: PFPUtils.cxx:2233
tca::FindCloseTjs
std::vector< int > FindCloseTjs(TCSlice &slc, const TrajPoint &fromTp, const TrajPoint &toTp, const float &maxDelta)
Definition: Utils.cxx:2604
tca::TrajLength
float TrajLength(Trajectory &tj)
Definition: Utils.cxx:2337
tca::TCConfig::neutralVxCuts
std::vector< float > neutralVxCuts
Definition: DataStructs.h:469
tca::Trajectory::EndPt
std::array< unsigned short, 2 > EndPt
First and last point in the trajectory that has charge.
Definition: DataStructs.h:156
tca::Trajectory::PDGCode
unsigned short PDGCode
shower-like or track-like {default is track-like}
Definition: DataStructs.h:160
tca::SignalAtTp
bool SignalAtTp(TCSlice &slc, const TrajPoint &tp)
Definition: Utils.cxx:1814
tca::ScoreVertices
void ScoreVertices(TCSlice &slc)
Definition: TCVertex.cxx:2391
tca::PosSep2
double PosSep2(const Point3_t &pos1, const Point3_t &pos2)
Definition: PFPUtils.cxx:1660
tca::CloseEnd
unsigned short CloseEnd(TCSlice &slc, const Trajectory &tj, const Point2_t &pos)
Definition: Utils.cxx:2251
tca::PrintHit
std::string PrintHit(const TCHit &tch)
Definition: Utils.cxx:5704
tca::SortEntry::index
unsigned int index
Definition: PFPUtils.cxx:8
tca::SetPDGCode
void SetPDGCode(TCSlice &slc, unsigned short itj, bool tjDone)
Definition: Utils.cxx:3795
tca::VtxStore::Pos
Point2_t Pos
Definition: DataStructs.h:55
tca::TCConfig::geom
const geo::GeometryCore * geom
Definition: DataStructs.h:493
tca::TCConfig::minPts
std::vector< unsigned short > minPts
min number of Pts required to make a trajectory
Definition: DataStructs.h:483
tca::TpSumHitChg
float TpSumHitChg(TCSlice &slc, TrajPoint const &tp)
Definition: Utils.cxx:1853
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:258
geo::GeometryCore::IntersectionPoint
bool IntersectionPoint(geo::WireID const &wid1, geo::WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Definition: GeometryCore.cxx:1749
tca::FarEnd
unsigned short FarEnd(TCSlice &slc, const PFPStruct &pfp, const Point3_t &pos)
Definition: PFPUtils.cxx:3032
tca::AnalyzePFP
bool AnalyzePFP(TCSlice &slc, PFPStruct &pfp, bool prt)
Definition: PFPUtils.cxx:2421
tca::kVtxIndPlnNoChg
vertex quality is suspect - No requirement made on chg btw it and the Tj
Definition: DataStructs.h:79
tca::Trajectory::ID
int ID
ID that is local to one slice.
Definition: DataStructs.h:157
tca::CompleteIncomplete3DVerticesInGaps
void CompleteIncomplete3DVerticesInGaps(TCSlice &slc)
Definition: TCVertex.cxx:2669
tca::Trajectory::VtxID
std::array< unsigned short, 2 > VtxID
ID of 2D vertex.
Definition: DataStructs.h:155
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
tca::NumPtsWithCharge
unsigned short NumPtsWithCharge(TCSlice &slc, const Trajectory &tj, bool includeDeadWires)
Definition: Utils.cxx:1863
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:524
tca::TCSlice::slHits
std::vector< TCHit > slHits
Definition: DataStructs.h:569
tca::MergeAndStore
bool MergeAndStore(TCSlice &slc, unsigned int itj1, unsigned int itj2, bool doPrt)
Definition: Utils.cxx:4095
tca::FindVtxTjs
void FindVtxTjs(TCSlice &slc, VtxStore &vx2)
Definition: TCVertex.cxx:409
tca::VtxStore::Stat
std::bitset< 16 > Stat
Vertex status bits using kVtxBit_t.
Definition: DataStructs.h:69
tca::IsCloseToVertex
unsigned short IsCloseToVertex(TCSlice &slc, VtxStore &inVx2)
Definition: TCVertex.cxx:3297
tca::VtxStore::ID
int ID
set to 0 if killed
Definition: DataStructs.h:64
detinfo::DetectorProperties::ConvertTicksToX
virtual double ConvertTicksToX(double ticks, int p, int t, int c) const =0
tca::ChgFracNearPos
float ChgFracNearPos(TCSlice &slc, const Point2_t &pos, const std::vector< int > &tjIDs)
Definition: Utils.cxx:2698
tca::TrajIntersection
void TrajIntersection(TrajPoint const &tp1, TrajPoint const &tp2, Point2_t &pos)
Definition: Utils.cxx:2299
tca::CheckTraj
void CheckTraj(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:911
tca::TCConfig::vtxScoreWeights
std::vector< float > vtxScoreWeights
Definition: DataStructs.h:468
tca::CTP_t
unsigned int CTP_t
Definition: DataStructs.h:41
tca::SetHighScoreBits
void SetHighScoreBits(TCSlice &slc, Vtx3Store &vx3)
Definition: TCVertex.cxx:2436
tca::TCSlice::TPCID
geo::TPCID TPCID
Definition: DataStructs.h:562
tca::AttachTrajToVertex
bool AttachTrajToVertex(TCSlice &slc, Trajectory &tj, VtxStore &vx, bool prt)
Definition: TCVertex.cxx:1893
tca::TCConfig::useAlg
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Definition: DataStructs.h:504
util::flags::to_string
std::string to_string(Flag_t< Storage > const flag)
Convert a flag into a stream (shows its index).
Definition: BitMask.h:187
close
in close()
dir
TDirectory * dir
Definition: macro.C:5
tca::CreatePFP
PFPStruct CreatePFP(TCSlice &slc)
Definition: PFPUtils.cxx:1973
tca::TCSlice::vtxs
std::vector< VtxStore > vtxs
2D vertices
Definition: DataStructs.h:576
tca::Trajectory::StepDir
short StepDir
-1 = going US (-> small wire#), 1 = going DS (-> large wire#)
Definition: DataStructs.h:162
tca::Trajectory::AlgMod
std::bitset< 128 > AlgMod
Bit set if algorithm AlgBit_t modifed the trajectory.
Definition: DataStructs.h:145
tca::VtxStore::NTraj
unsigned short NTraj
Definition: DataStructs.h:57
tca::DecodeCTP
geo::PlaneID DecodeCTP(CTP_t CTP)
Definition: DataStructs.cxx:108
tca::GetVtxTjIDs
std::vector< int > GetVtxTjIDs(const TCSlice &slc, const VtxStore &vx2)
Definition: TCVertex.cxx:3226
tca::PFPStruct::MatchVecIndex
unsigned short MatchVecIndex
Definition: DataStructs.h:226
tca::PosSep
double PosSep(const Point3_t &pos1, const Point3_t &pos2)
Definition: PFPUtils.cxx:1654
tca::PFPStruct::TjIDs
std::vector< int > TjIDs
Definition: DataStructs.h:207
tca::TCConfig::dbgAlg
std::bitset< 128 > dbgAlg
Allow user to turn on debug printing in algorithms (that print...)
Definition: DataStructs.h:505
tca::Vx3Topo
unsigned short Vx3Topo(TCSlice &slc, Vtx3Store &vx3)
Definition: TCVertex.cxx:2641
tca::TCEvent::allHits
std::vector< recob::Hit > const * allHits
Definition: DataStructs.h:536
tca::PrintPos
std::string PrintPos(TCSlice &slc, const TrajPoint &tp)
Definition: Utils.cxx:5712
tca::EncodeCTP
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
Definition: DataStructs.h:45
tca::Find2DVertices
void Find2DVertices(TCSlice &slc, const CTP_t &inCTP, unsigned short pass)
Definition: TCVertex.cxx:114
tca::TCSlice::nWires
std::vector< unsigned int > nWires
Definition: DataStructs.h:553
tca::PointTrajDOCA
float PointTrajDOCA(TCSlice &slc, unsigned int iht, TrajPoint const &tp)
Definition: Utils.cxx:2270
tca::KillPoorVertices
void KillPoorVertices(TCSlice &slc)
Definition: TCVertex.cxx:2418
tca::TCSlice::ID
int ID
ID of the recob::Slice (not the sub-slice)
Definition: DataStructs.h:564
tca::TPNearVertex
unsigned short TPNearVertex(TCSlice &slc, const TrajPoint &tp)
Definition: TCVertex.cxx:1816
tca::VtxStore::UID
int UID
unique global ID
Definition: DataStructs.h:65
tca::SetVx2Score
void SetVx2Score(TCSlice &slc)
Definition: TCVertex.cxx:2504
tca::SetDifference
std::vector< T > SetDifference(const std::vector< T > &set1, const std::vector< T > &set2)
Definition: Utils.h:253
tca::FindHammerVertices
void FindHammerVertices(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:1154
tca::TCConfig::vtx2DCuts
std::vector< float > vtx2DCuts
Max position pull, max Position error rms.
Definition: DataStructs.h:466
tca::TCSlice::nPlanes
unsigned short nPlanes
Definition: DataStructs.h:563
tca::TCSlice::matchVec
std::vector< MatchStruct > matchVec
3D matching vector
Definition: DataStructs.h:578
tca::CompleteIncomplete3DVertices
void CompleteIncomplete3DVertices(TCSlice &slc)
Definition: TCVertex.cxx:2759
evd::details::end
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
Definition: RawDataDrawer.cxx:317
tca::TCSlice::pfps
std::vector< PFPStruct > pfps
Definition: DataStructs.h:579
tca::FindNeutralVertices
void FindNeutralVertices(TCSlice &slc)
Definition: TCVertex.cxx:529
tca::MoveTPToWire
void MoveTPToWire(TrajPoint &tp, float wire)
Definition: Utils.cxx:2485
tca::PosInPlane
void PosInPlane(const TCSlice &slc, const Vtx3Store &vx3, unsigned short plane, Point2_t &pos)
Definition: TCVertex.cxx:3287
tca::VertexVertexPull
float VertexVertexPull(TCSlice &slc, const Vtx3Store &vx1, const Vtx3Store &vx2)
Definition: TCVertex.cxx:2064
tca::MakeJunkVertices
void MakeJunkVertices(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:14
tca::evt
TCEvent evt
Definition: DataStructs.cxx:5
tca::SplitTraj
bool SplitTraj(TCSlice &slc, unsigned short itj, float XPos, bool makeVx2, bool prt)
Definition: Utils.cxx:2004
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:203
tca::StorePFP
bool StorePFP(TCSlice &slc, PFPStruct &pfp)
Definition: PFPUtils.cxx:2749
tca::Find3DVertices
void Find3DVertices(TCSlice &slc)
Definition: TCVertex.cxx:1392
tca::PrintHeader
void PrintHeader(std::string someText)
Definition: Utils.cxx:5535
tca::FindHammerVertices2
void FindHammerVertices2(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:979
tca::VtxStore::TjChgFrac
float TjChgFrac
Fraction of charge near the vertex that is from hits on the vertex Tjs.
Definition: DataStructs.h:68
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:449
tca::PrintTrajPoint
void PrintTrajPoint(std::string someText, TCSlice &slc, unsigned short ipt, short dir, unsigned short pass, TrajPoint const &tp)
Definition: Utils.cxx:5541
tca::ChkVxTjs
void ChkVxTjs(TCSlice &slc, const CTP_t &inCTP, bool prt)
Definition: TCVertex.cxx:874
tca::TrajPoint::HitPosErr2
double HitPosErr2
Definition: DataStructs.h:123
tca::CompatibleMerge
bool CompatibleMerge(TCSlice &slc, std::vector< int > &tjIDs, bool prt)
Definition: Utils.cxx:534
geo::GeometryCore::HasWire
bool HasWire(geo::WireID const &wireid) const
Returns whether we have the specified wire.
Definition: GeometryCore.h:3334
tca::Match3DVtxTjs
void Match3DVtxTjs(TCSlice &slc, bool prt)
Definition: TCVertex.cxx:1693