All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TCVertex.cxx
Go to the documentation of this file.
1 #include "larreco/RecoAlg/TCAlg/TCVertex.h"
2 
3 #include "larcorealg/Geometry/GeometryCore.h"
4 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
5 #include "lardataalg/DetectorInfo/DetectorPropertiesData.h"
6 #include "larreco/RecoAlg/TCAlg/DebugStruct.h"
7 #include "larreco/RecoAlg/TCAlg/PFPUtils.h"
8 #include "larreco/RecoAlg/TCAlg/Utils.h"
9 #include "messagefacility/MessageLogger/MessageLogger.h"
10 
11 #include <algorithm>
12 #include <array>
13 #include <bitset>
14 #include <cmath>
15 #include <iomanip>
16 #include <iostream>
17 #include <limits.h>
18 #include <stdlib.h>
19 #include <string>
20 #include <vector>
21 
22 #include "TMatrixD.h"
23 #include "TVectorD.h"
24 
25 namespace tca {
26 
27  using namespace detail; // SortEntry, valsDecreasing(), valsIncreasing();
28 
29  //////////////////////////////////////////
30  void
31  MakeJunkVertices(TCSlice& slc, const CTP_t& inCTP)
32  {
33  // Vertices between poorly reconstructed tjs (especially junk slc) and normal
34  // tjs can fail because the junk tj trajectory parameters are inaccurate. This function
35  // uses proximity and not pointing to make junk vertices
36  // Don't use this if standard vertex reconstruction is disabled
37  if (tcc.vtx2DCuts[0] <= 0) return;
38  if (!tcc.useAlg[kJunkVx]) return;
39  if (slc.tjs.size() < 2) return;
40 
41  // Look for tjs that are within maxSep of the end of a Tj
42  constexpr float maxSep = 4;
43 
44  geo::PlaneID planeID = DecodeCTP(inCTP);
45  bool prt = (tcc.dbgVxJunk && tcc.dbgSlc);
46  if (prt) {
47  mf::LogVerbatim("TC") << "MakeJunkVertices: prt set for plane " << planeID.Plane
48  << " maxSep btw tjs " << maxSep;
49  }
50 
51  // make a template vertex
52  VtxStore junkVx;
53  junkVx.CTP = inCTP;
54  junkVx.Topo = 9;
55  junkVx.Stat[kJunkVx] = true;
56  junkVx.Stat[kFixed] = true;
57  // set an invalid ID
58  junkVx.ID = USHRT_MAX;
59  // put in generous errors
60  junkVx.PosErr = {{2.0, 2.0}};
61  // define a minimal score so it won't get clobbered
62  junkVx.Score = tcc.vtx2DCuts[7] + 0.1;
63 
64  // look at both ends of long tjs
65  for (unsigned short it1 = 0; it1 < slc.tjs.size() - 1; ++it1) {
66  auto& tj1 = slc.tjs[it1];
67  if (tj1.AlgMod[kKilled] || tj1.AlgMod[kHaloTj]) continue;
68  if (tj1.SSID > 0 || tj1.AlgMod[kShowerLike]) continue;
69  if (tj1.CTP != inCTP) continue;
70  if (tj1.AlgMod[kJunkTj]) continue;
71  if (TrajLength(tj1) < 10) continue;
72  if (tj1.MCSMom < 100) continue;
73  for (unsigned short end1 = 0; end1 < 2; ++end1) {
74  // existing vertex?
75  if (tj1.VtxID[end1] > 0) continue;
76  auto& tp1 = tj1.Pts[tj1.EndPt[end1]];
77  // get a list of tjs in this vicinity
78  auto tjlist = FindCloseTjs(slc, tp1, tp1, maxSep);
79  if (tjlist.empty()) continue;
80  // set to an invalid ID
81  junkVx.ID = USHRT_MAX;
82  for (auto tj2id : tjlist) {
83  auto& tj2 = slc.tjs[tj2id - 1];
84  if (tj2.CTP != inCTP) continue;
85  if (tj2id == tj1.ID) continue;
86  if (tj2.SSID > 0 || tj2.AlgMod[kShowerLike]) continue;
87  float close = maxSep;
88  unsigned short closeEnd = USHRT_MAX;
89  for (unsigned short end2 = 0; end2 < 2; ++end2) {
90  auto& tp2 = tj2.Pts[tj2.EndPt[end2]];
91  float sep = PosSep(tp1.Pos, tp2.Pos);
92  if (sep < close) {
93  close = sep;
94  closeEnd = end2;
95  } // sep
96  } // end2
97  if (closeEnd > 1) continue;
98  auto& tp2 = tj2.Pts[tj2.EndPt[closeEnd]];
99  bool signalBetween = SignalBetween(slc, tp1, tp2, 0.8);
100  if (!signalBetween) continue;
101  if (junkVx.ID == USHRT_MAX) {
102  // define the new vertex
103  junkVx.ID = slc.vtxs.size() + 1;
104  junkVx.Pos = tp1.Pos;
105  } // new vertex
106  tj2.VtxID[closeEnd] = junkVx.ID;
107  tj1.VtxID[end1] = junkVx.ID;
108  } // tjid
109  if (junkVx.ID == USHRT_MAX) continue;
110  if (!StoreVertex(slc, junkVx)) {
111  mf::LogVerbatim("TC") << "MJV: StoreVertex failed";
112  for (auto& tj : slc.tjs) {
113  if (tj.AlgMod[kKilled]) continue;
114  if (tj.VtxID[0] == junkVx.ID) tj.VtxID[0] = 0;
115  if (tj.VtxID[1] == junkVx.ID) tj.VtxID[1] = 0;
116  } // tj
117  continue;
118  } // StoreVertex failed
119  if (prt) {
120  mf::LogVerbatim("TC") << " New junk 2V" << junkVx.ID << " at " << std::fixed
121  << std::setprecision(1) << junkVx.Pos[0] << ":"
122  << junkVx.Pos[1] / tcc.unitsPerTick;
123  } // prt
124  junkVx.ID = USHRT_MAX;
125  } // end1
126  } // it1
127 
128  } // MakeJunkVertices
129 
130  //////////////////////////////////////////
131  void
132  Find2DVertices(detinfo::DetectorPropertiesData const& detProp,
133  TCSlice& slc,
134  const CTP_t& inCTP,
135  unsigned short pass)
136  {
137  // Find 2D vertices between pairs of tjs that have a same-end topology. Using an example
138  // where StepDir = 1 (end 0 is at small wire number) vertices will be found with Topo = 0
139  // with a vertex US of the ends (<) or Topo = 2 with a vertex DS of the ends (>). This is reversed
140  // if StepDir = -1. Vertices with Topo = 1 (/\) and (\/) are found in EndMerge.
141 
142  // tcc.vtx2DCuts fcl input usage
143  // 0 = maximum length of a short trajectory
144  // 1 = max vertex - trajectory separation for short trajectories
145  // 2 = max vertex - trajectory separation for long trajectories
146  // 3 = max position pull for adding TJs to a vertex
147  // 4 = max allowed vertex position error
148  // 5 = min MCSMom
149  // 6 = min Pts/Wire fraction
150  // 7 min Score
151  // 8 Min charge fraction near a merge point (not a vertex)
152  // 9 max MCSmom asymmetry for a merge
153  // 10 Require charge on wires between a vtx and the start of the tjs in induction planes? (1 = yes)
154 
155  if (tcc.vtx2DCuts[0] <= 0) return;
156  if (slc.tjs.size() < 2) return;
157 
158  bool firstPassCuts = (pass == 0);
159 
160  geo::PlaneID planeID = DecodeCTP(inCTP);
161 
162  // require charge between the vertex and the tj start points?
163  bool requireVtxTjChg = true;
164  if (tcc.vtx2DCuts[10] == 0 && int(planeID.Plane) < slc.nPlanes - 1) requireVtxTjChg = false;
165 
166  bool prt = (tcc.dbg2V && tcc.dbgSlc && debug.Plane == (int)planeID.Plane);
167  if (prt) {
168  mf::LogVerbatim("TC") << "prt set for CTP " << inCTP << " in Find2DVertices. firstPassCuts? "
169  << firstPassCuts << " requireVtxTjChg " << requireVtxTjChg;
170  PrintAllTraj(detProp, "F2DVi", slc, USHRT_MAX, slc.tjs.size());
171  }
172 
173  unsigned short maxShortTjLen = tcc.vtx2DCuts[0];
174  for (unsigned short it1 = 0; it1 < slc.tjs.size() - 1; ++it1) {
175  auto& tj1 = slc.tjs[it1];
176  if (tj1.AlgMod[kKilled] || tj1.AlgMod[kHaloTj]) continue;
177  if (tj1.SSID > 0 || tj1.AlgMod[kShowerLike]) continue;
178  if (tj1.CTP != inCTP) continue;
179  bool tj1Short = (TrajLength(tj1) < maxShortTjLen);
180  for (unsigned short end1 = 0; end1 < 2; ++end1) {
181  // vertex assignment exists?
182  if (tj1.VtxID[end1] > 0) continue;
183  // wrong end of a high energy electron?
184  if (tj1.PDGCode == 111 && end1 != tj1.StartEnd) continue;
185  // default condition is to use the end point to define the trajectory and direction
186  // at the end
187  short endPt1 = tj1.EndPt[end1];
188  float wire1 = tj1.Pts[endPt1].Pos[0];
189  // unless there are few points fitted, indicating that the trajectory fit
190  // may have been biased by the presence of another trajectory at the vertex or by
191  // other close unresolved tracks
192  if (tj1.Pts.size() > 6 && tj1.Pts[endPt1].NTPsFit < 4) {
193  if (end1 == 0 && endPt1 < int(tj1.Pts.size()) - 3) { endPt1 += 3; }
194  else if (end1 == 1 && endPt1 >= 3) {
195  endPt1 -= 3;
196  }
197  if (tj1.Pts[endPt1].Chg == 0) endPt1 = NearestPtWithChg(slc, tj1, endPt1);
198  } // few points fit at end1
199  TrajPoint tp1 = tj1.Pts[endPt1];
200  MoveTPToWire(tp1, wire1);
201  // re-purpose endPt1 to reference the end point. This will be used the find the point on
202  // tj1 that is closest to the vertex position
203  endPt1 = tj1.EndPt[end1];
204  short oendPt1 = tj1.EndPt[1 - end1];
205  // reference to the other end of tj1
206  auto& otp1 = tj1.Pts[oendPt1];
207  for (unsigned short it2 = it1 + 1; it2 < slc.tjs.size(); ++it2) {
208  auto& tj2 = slc.tjs[it2];
209  if (tj2.AlgMod[kKilled] || tj2.AlgMod[kHaloTj]) continue;
210  if (tj2.SSID > 0 || tj2.AlgMod[kShowerLike]) continue;
211  if (tj2.CTP != inCTP) continue;
212  if (tj1.VtxID[end1] > 0) continue;
213  if (tj1.MCSMom < tcc.vtx2DCuts[5] && tj2.MCSMom < tcc.vtx2DCuts[5]) continue;
214  bool tj2Short = (TrajLength(tj2) < maxShortTjLen);
215  // find the end that is closer to tp1
216  unsigned short end2 = 0;
217  if (PosSep2(tj2.Pts[tj2.EndPt[1]].Pos, tp1.Pos) <
218  PosSep2(tj2.Pts[tj2.EndPt[0]].Pos, tp1.Pos))
219  end2 = 1;
220  if (tj2.VtxID[end2] > 0) continue;
221  // wrong end of a high energy electron?
222  if (tj2.PDGCode == 111 && end2 != tj2.StartEnd) continue;
223  // check for a vertex between these tjs at the other ends
224  if (tj1.VtxID[1 - end1] > 0 && tj1.VtxID[1 - end1] == tj2.VtxID[1 - end2]) continue;
225  // see if the other ends are closer
226  unsigned short oendPt2 = tj2.EndPt[1 - end2];
227  auto& otp2 = tj2.Pts[oendPt2];
228  if (PosSep2(otp1.Pos, otp2.Pos) < PosSep2(tp1.Pos, tj2.Pts[tj2.EndPt[end2]].Pos))
229  continue;
230  short endPt2 = tj2.EndPt[end2];
231  float wire2 = tj2.Pts[endPt2].Pos[0];
232  if (tj2.Pts.size() > 6 && tj2.Pts[endPt2].NTPsFit < 4) {
233  if (end2 == 0 && endPt2 < int(tj2.Pts.size()) - 3) { endPt2 += 3; }
234  else if (end2 == 1 && endPt2 >= 3) {
235  endPt2 -= 3;
236  }
237  if (tj2.Pts[endPt2].Chg == 0) endPt2 = NearestPtWithChg(slc, tj2, endPt2);
238  } // few points fit at end1
239  TrajPoint tp2 = tj2.Pts[endPt2];
240  MoveTPToWire(tp2, wire2);
241  // re-purpose endPt2
242  endPt2 = tj2.EndPt[end2];
243  // Rough first cut on the separation between the end points of the
244  // two trajectories
245  float sepCut = 100;
246  if (std::abs(tp1.Pos[0] - tp2.Pos[0]) > sepCut) continue;
247  if (std::abs(tp1.Pos[1] - tp2.Pos[1]) > sepCut) continue;
248  float wint, tint;
249  TrajIntersection(tp1, tp2, wint, tint);
250  // make sure this is inside the TPC.
251  if (wint < 0 || wint > tcc.maxPos0[planeID.Plane] - 3) continue;
252  if (tint < 0 || tint > tcc.maxPos1[planeID.Plane]) continue;
253  // Next cut on separation between the TPs and the intersection point
254  if (tj1Short || tj2Short) { sepCut = tcc.vtx2DCuts[1]; }
255  else {
256  sepCut = tcc.vtx2DCuts[2];
257  }
258  // NewVtxCuts: require close separation on the first pass
259  if (firstPassCuts) sepCut = tcc.vtx2DCuts[1];
260  Point2_t vPos{{wint, tint}};
261  float vt1Sep = PosSep(vPos, tp1.Pos);
262  float vt2Sep = PosSep(vPos, tp2.Pos);
263  float dwc1 = DeadWireCount(slc, wint, tp1.Pos[0], tp1.CTP);
264  float dwc2 = DeadWireCount(slc, wint, tp2.Pos[0], tp1.CTP);
265  vt1Sep -= dwc1;
266  vt2Sep -= dwc2;
267  bool vtxOnDeadWire = (DeadWireCount(slc, wint, wint, tp1.CTP) == 1);
268  if (prt && vt1Sep < 200 && vt2Sep < 200) {
269  mf::LogVerbatim myprt("TC");
270  myprt << "F2DV candidate T" << tj1.ID << "_" << end1 << "-T" << tj2.ID << "_" << end2;
271  myprt << " vtx pos " << (int)wint << ":" << (int)(tint / tcc.unitsPerTick) << " tp1 "
272  << PrintPos(slc, tp1) << " tp2 " << PrintPos(slc, tp2);
273  myprt << " dwc1 " << dwc1 << " dwc2 " << dwc2 << " on dead wire? " << vtxOnDeadWire;
274  myprt << " vt1Sep " << vt1Sep << " vt2Sep " << vt2Sep << " sepCut " << sepCut;
275  }
276  if (vt1Sep > sepCut || vt2Sep > sepCut) continue;
277  // make sure that the other end isn't closer
278  if (PosSep(vPos, slc.tjs[it1].Pts[oendPt1].Pos) < vt1Sep) {
279  if (prt)
280  mf::LogVerbatim("TC") << " tj1 other end " << PrintPos(slc, tj1.Pts[oendPt1])
281  << " is closer to the vertex";
282  continue;
283  }
284  if (PosSep(vPos, slc.tjs[it2].Pts[oendPt2].Pos) < vt2Sep) {
285  if (prt)
286  mf::LogVerbatim("TC") << " tj2 other end " << PrintPos(slc, tj2.Pts[oendPt2])
287  << " is closer to the vertex";
288  continue;
289  }
290  // Ensure that the vertex position is close to the end of each Tj
291  unsigned short closePt1;
292  float doca1 = sepCut;
293  if (!TrajClosestApproach(tj1, wint, tint, closePt1, doca1)) continue;
294  // dpt1 (and dpt2) will be 0 if the vertex is at the end
295  short stepDir = -1;
296  if (tcc.modes[kStepDir]) stepDir = 1;
297  short dpt1 = stepDir * (closePt1 - endPt1);
298  if (prt)
299  mf::LogVerbatim("TC") << " endPt1 " << endPt1 << " closePt1 " << closePt1 << " dpt1 "
300  << dpt1 << " doca1 " << doca1;
301  if (dpt1 < -1) continue;
302  if (slc.tjs[it1].EndPt[1] > 4) {
303  if (dpt1 > 3) continue;
304  }
305  else {
306  // tighter cut for short trajectories
307  if (dpt1 > 2) continue;
308  }
309  unsigned short closePt2;
310  float doca2 = sepCut;
311  if (!TrajClosestApproach(tj2, wint, tint, closePt2, doca2)) continue;
312  short dpt2 = stepDir * (closePt2 - endPt2);
313  if (prt)
314  mf::LogVerbatim("TC") << " endPt2 " << endPt2 << " closePt2 " << closePt2 << " dpt2 "
315  << dpt2 << " doca2 " << doca2;
316  if (dpt2 < -1) continue;
317  if (slc.tjs[it2].EndPt[1] > 4) {
318  if (dpt2 > 3) continue;
319  }
320  else {
321  // tighter cut for short trajectories
322  if (dpt2 > 2) continue;
323  }
324  bool fixVxPos = false;
325  // fix the vertex position if there is a charge kink here
326  if (tj1.EndFlag[end1][kAtKink]) fixVxPos = true;
327  if (prt)
328  mf::LogVerbatim("TC") << " wint:tint " << (int)wint << ":"
329  << (int)(tint / tcc.unitsPerTick) << " fixVxPos? " << fixVxPos;
330  if (requireVtxTjChg) {
331  // ensure that there is a signal between these TPs and the vertex on most of the wires
332  bool signalBetween = true;
333  short dpt = abs(wint - tp1.Pos[0]);
334  if (dpt > 2 && !SignalBetween(slc, tp1, wint, tcc.vtx2DCuts[6])) {
335  if (prt) mf::LogVerbatim("TC") << " Fails SignalBetween for tp1 " << dpt;
336  signalBetween = false;
337  }
338  dpt = abs(wint - tp2.Pos[0]);
339  if (dpt > 2 && !SignalBetween(slc, tp2, wint, tcc.vtx2DCuts[6])) {
340  if (prt) mf::LogVerbatim("TC") << " Fails SignalBetween for tp2 " << dpt;
341  signalBetween = false;
342  }
343  // consider the case where the intersection point is wrong because the
344  // end TP angles are screwed up but the Tjs are close to each other near the end
345  if (!signalBetween) {
346  unsigned short ipt1, ipt2;
347  float maxSep = 3;
348  bool isClose = TrajTrajDOCA(slc, tj1, tj2, ipt1, ipt2, maxSep, false);
349  // require that they are close at the correct end
350  if (isClose) isClose = (abs(ipt1 - endPt1) < 4 && abs(ipt2 - endPt2) < 4);
351  if (isClose) {
352  if (prt)
353  mf::LogVerbatim("TC")
354  << " TrajTrajDOCA are close with minSep " << maxSep << " near "
355  << PrintPos(slc, tj1.Pts[ipt1].Pos) << " " << PrintPos(slc, tj2.Pts[ipt2].Pos);
356  // put the vertex at the TP that is closest to the intersection point
357  Point2_t vpos = {{wint, tint}};
358  if (PosSep2(tp1.Pos, vpos) < PosSep2(tp2.Pos, vpos)) {
359  wint = tp1.Pos[0];
360  tint = tp1.Pos[1];
361  }
362  else {
363  wint = tp2.Pos[0];
364  tint = tp2.Pos[1];
365  }
366  fixVxPos = true;
367  if (prt)
368  mf::LogVerbatim("TC")
369  << " new wint:tint " << (int)wint << ":" << (int)(tint / tcc.unitsPerTick);
370  }
371  else {
372  // closest approach > 3
373  continue;
374  }
375  } // no signal between
376  } // requireVtxTjChg
377  // make a new temporary vertex
378  VtxStore aVtx;
379  aVtx.Pos[0] = wint;
380  aVtx.Pos[1] = tint;
381  aVtx.NTraj = 0;
382  aVtx.Pass = tj1.Pass;
383  // Topo 0 has this topology (<) and Topo 2 has this (>)
384  aVtx.Topo = 2 * end1;
385  aVtx.ChiDOF = 0;
386  aVtx.CTP = inCTP;
387  aVtx.Stat[kOnDeadWire] = vtxOnDeadWire;
388  // fix the vertex position if we needed to move it significantly, or if it is on a dead wire
389  aVtx.Stat[kFixed] = fixVxPos;
390  aVtx.Stat[kVxIndPlnNoChg] = !requireVtxTjChg;
391  // try to fit it. We need to give it an ID to do that. Take the next
392  // available ID
393  unsigned short newVtxID = slc.vtxs.size() + 1;
394  aVtx.ID = newVtxID;
395  tj1.VtxID[end1] = newVtxID;
396  tj2.VtxID[end2] = newVtxID;
397  if (!FitVertex(slc, aVtx, prt)) {
398  tj1.VtxID[end1] = 0;
399  tj2.VtxID[end2] = 0;
400  continue;
401  }
402  // check proximity to nearby vertices
403  unsigned short mergeMeWithVx = IsCloseToVertex(slc, aVtx);
404  if (mergeMeWithVx > 0 && MergeWithVertex(slc, aVtx, mergeMeWithVx)) {
405  if (prt) mf::LogVerbatim("TC") << " Merged with close vertex " << mergeMeWithVx;
406  continue;
407  }
408  // Save it
409  if (!StoreVertex(slc, aVtx)) continue;
410  if (prt) {
411  mf::LogVerbatim myprt("TC");
412  myprt << " New vtx 2V" << aVtx.ID;
413  myprt << " Tjs " << tj1.ID << "_" << end1 << "-" << tj2.ID << "_" << end2;
414  myprt << " at " << std::fixed << std::setprecision(1) << aVtx.Pos[0] << ":"
415  << aVtx.Pos[1] / tcc.unitsPerTick;
416  }
417  AttachAnyTrajToVertex(slc, slc.vtxs.size() - 1, prt);
418  SetVx2Score(slc);
419  } // it2
420  } // end1
421  } // it1
422 
423  // only call these on the last pass
424  if (pass == USHRT_MAX) {
425  FindHammerVertices(slc, inCTP);
426  FindHammerVertices2(slc, inCTP);
427  }
428 
429  if (prt) PrintAllTraj(detProp, "F2DVo", slc, USHRT_MAX, USHRT_MAX);
430 
431  } // Find2DVertices
432 
433  //////////////////////////////////////////
434  bool
435  MergeWithVertex(TCSlice& slc, VtxStore& vx, unsigned short oVxID)
436  {
437  // Attempts to merge the trajectories attached to vx with an existing 2D vertex
438  // referenced by existingVxID. This function doesn't use the existing end0/end1 vertex association.
439  // It returns true if the merging was successful in which case the calling function should
440  // not store vx. The calling function needs to have set VtxID to vx.ID for tjs that are currently attached
441  // to vx. It assumed that vx hasn't yet been pushed onto slc.vtxs
442 
443  if (!tcc.useAlg[kVxMerge]) return false;
444 
445  bool prt = tcc.dbgVxMerge && tcc.dbgSlc;
446 
447  if (oVxID > slc.vtxs.size()) return false;
448  auto& oVx = slc.vtxs[oVxID - 1];
449  if (vx.CTP != oVx.CTP) return false;
450 
451  // get a list of tjs attached to both vertices
452  std::vector<int> tjlist = GetVtxTjIDs(slc, vx);
453  if (tjlist.empty()) return false;
454  std::vector<int> tmp = GetVtxTjIDs(slc, oVx);
455  if (tmp.empty()) return false;
456  for (auto tjid : tmp) {
457  if (std::find(tjlist.begin(), tjlist.end(), tjid) == tjlist.end()) tjlist.push_back(tjid);
458  } // tjid
459  if (tjlist.size() < 2) return false;
460  // handle the simple case
461  if (tjlist.size() == 2) {
462  // Unset the fixed bit
463  vx.Stat[kFixed] = false;
464  oVx.Stat[kFixed] = false;
465  // assign the vx tjs to oVx
466  for (auto tjid : tjlist) {
467  auto& tj = slc.tjs[tjid - 1];
468  for (unsigned short end = 0; end < 2; ++end) {
469  if (tj.VtxID[end] == vx.ID) tj.VtxID[end] = oVx.ID;
470  } // end
471  } // tjid
472  if (!FitVertex(slc, oVx, prt)) {
473  if (prt)
474  mf::LogVerbatim("TC") << "MWV: merge failed " << vx.ID << " and existing " << oVx.ID;
475  return false;
476  }
477  return true;
478  } // size = 2
479 
480  // sort by decreasing length
481  std::vector<SortEntry> sortVec(tjlist.size());
482  for (unsigned int indx = 0; indx < sortVec.size(); ++indx) {
483  sortVec[indx].index = indx;
484  auto& tj = slc.tjs[tjlist[indx] - 1];
485  sortVec[indx].val = tj.Pts.size();
486  } // indx
487  std::sort(sortVec.begin(), sortVec.end(), valsDecreasing);
488  // re-order the list of Tjs
489  auto ttl = tjlist;
490  for (unsigned short ii = 0; ii < sortVec.size(); ++ii)
491  tjlist[ii] = ttl[sortVec[ii].index];
492  // Create a local vertex using the two longest slc, then add the shorter ones
493  // until the pull reaches the cut
494  VtxStore aVx;
495  aVx.CTP = vx.CTP;
496  std::vector<TrajPoint> tjpts(tjlist.size());
497  // determine which point on each Tj that will be used in the vertex fit and stash it in
498  // the traj point Step variable. This requires knowing the real position of the merged vertex
499  // which we estimate by averaging
500  std::array<float, 2> vpos;
501  vpos[0] = 0.5 * (vx.Pos[0] + oVx.Pos[0]);
502  vpos[1] = 0.5 * (vx.Pos[1] + oVx.Pos[1]);
503  for (unsigned short ii = 0; ii < tjpts.size(); ++ii) {
504  auto& tj = slc.tjs[tjlist[ii] - 1];
505  unsigned short npwc = NumPtsWithCharge(slc, tj, false);
506  unsigned short end = CloseEnd(slc, tj, vpos);
507  // assume that we will use the end point of the tj
508  unsigned short endPt = tj.EndPt[end];
509  if (npwc > 6 && tj.Pts[endPt].NTPsFit < 4) {
510  if (end == 0) { endPt += 3; }
511  else {
512  endPt -= 3;
513  }
514  endPt = NearestPtWithChg(slc, tj, endPt);
515  } // few points fit at the end
516  if (endPt < tj.EndPt[0]) endPt = tj.EndPt[0];
517  if (endPt > tj.EndPt[1]) endPt = tj.EndPt[1];
518  // define tjpts
519  tjpts[ii].CTP = tj.CTP;
520  tjpts[ii].Pos = tj.Pts[endPt].Pos;
521  tjpts[ii].Dir = tj.Pts[endPt].Dir;
522  tjpts[ii].Ang = tj.Pts[endPt].Ang;
523  tjpts[ii].AngErr = tj.Pts[endPt].AngErr;
524  // stash the point in Step
525  tjpts[ii].Step = endPt;
526  // and the end in AngleCode
527  tjpts[ii].AngleCode = end;
528  // stash the ID in Hits
529  tjpts[ii].Hits.resize(1, tj.ID);
530  } // tjid
531  if (prt) {
532  mf::LogVerbatim myprt("TC");
533  myprt << "MWV: " << oVxID;
534  myprt << " Fit TPs";
535  for (unsigned short ii = 0; ii < tjpts.size(); ++ii) {
536  auto& tjpt = tjpts[ii];
537  myprt << " " << tjlist[ii] << "_" << tjpt.Step << "_" << PrintPos(slc, tjpt.Pos);
538  }
539  } // prt
540  // create a subset of the first two for the first fit
541  auto fitpts = tjpts;
542  fitpts.resize(2);
543  if (!FitVertex(slc, aVx, fitpts, prt)) {
544  if (prt) mf::LogVerbatim("TC") << "MWV: first fit failed ";
545  return false;
546  }
547  // Fit and add tjs to the vertex
548  bool needsUpdate = false;
549  for (unsigned short ii = 2; ii < tjlist.size(); ++ii) {
550  fitpts.push_back(tjpts[ii]);
551  if (FitVertex(slc, aVx, fitpts, prt)) { needsUpdate = false; }
552  else {
553  // remove the last Tj point and keep going
554  fitpts.pop_back();
555  needsUpdate = true;
556  }
557  } // ii
558 
559  if (needsUpdate) FitVertex(slc, aVx, fitpts, prt);
560  if (prt) mf::LogVerbatim("TC") << "MWV: done " << vx.ID << " and existing " << oVx.ID;
561 
562  // update. Remove old associations
563  for (auto& tj : slc.tjs) {
564  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
565  if (tj.CTP != vx.CTP) continue;
566  for (unsigned short end = 0; end < 2; ++end) {
567  if (tj.VtxID[end] == vx.ID) tj.VtxID[end] = 0;
568  if (tj.VtxID[end] == oVxID) tj.VtxID[end] = 0;
569  }
570  } // tj
571  // set the new associations
572  for (unsigned short ii = 0; ii < fitpts.size(); ++ii) {
573  auto& tjpt = fitpts[ii];
574  unsigned short end = tjpt.AngleCode;
575  auto& tj = slc.tjs[tjpt.Hits[0] - 1];
576  if (tj.VtxID[end] != 0) return false;
577  tj.VtxID[end] = oVxID;
578  } // ii
579 
580  // Update oVx
581  oVx.Pos = aVx.Pos;
582  oVx.PosErr = aVx.PosErr;
583  oVx.ChiDOF = aVx.ChiDOF;
584  oVx.NTraj = fitpts.size();
585  // Update the score and the charge fraction
586  SetVx2Score(slc, oVx);
587  oVx.Stat[kVxMerged] = true;
588  oVx.Stat[kFixed] = false;
589  if (prt) {
590  mf::LogVerbatim myprt("TC");
591  myprt << "MWV: " << oVxID;
592  myprt << " Done TPs";
593  for (unsigned short ii = 0; ii < fitpts.size(); ++ii) {
594  auto& tjpt = fitpts[ii];
595  myprt << " " << tjpt.Hits[0] << "_" << tjpt.AngleCode << "_" << PrintPos(slc, tjpt.Pos);
596  }
597  } // prt
598 
599  return true;
600  } // MergeWithVertex
601 
602  //////////////////////////////////////////
603  void
604  FindHammerVertices2(TCSlice& slc, const CTP_t& inCTP)
605  {
606  // Variant of FindHammerVertices with slightly different requirements:
607  // 1) tj1 is a straight trajectory with most of the points fit
608  // 2) No angle requirement between tj1 and tj2
609  // 3) Large charge near the intersection point X on tj2
610  // tj2 ---X---
611  // tj1 /
612  // tj1 /
613  // tj1 /
614  // minimum^2 DOCA of tj1 endpoint with tj2
615 
616  if (!tcc.useAlg[kHamVx2]) return;
617  // This wasn't written for test beam events
618  if (tcc.modes[kTestBeam]) return;
619 
620  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kHamVx2]);
621  if (prt) mf::LogVerbatim("TC") << "Inside HamVx2";
622 
623  for (unsigned short it1 = 0; it1 < slc.tjs.size(); ++it1) {
624  if (slc.tjs[it1].CTP != inCTP) continue;
625  if (slc.tjs[it1].AlgMod[kKilled] || slc.tjs[it1].AlgMod[kHaloTj]) continue;
626  if (slc.tjs[it1].AlgMod[kHamVx]) continue;
627  if (slc.tjs[it1].AlgMod[kHamVx2]) continue;
628  if (slc.tjs[it1].AlgMod[kJunkTj]) continue;
629  if (slc.tjs[it1].PDGCode == 111) continue;
630  unsigned short numPtsWithCharge1 = NumPtsWithCharge(slc, slc.tjs[it1], false);
631  if (numPtsWithCharge1 < 6) continue;
632  // require high MCSMom
633  if (slc.tjs[it1].MCSMom < 200) continue;
634  // Check each end of tj1
635  bool didaSplit = false;
636  for (unsigned short end1 = 0; end1 < 2; ++end1) {
637  // vertex assignment exists?
638  if (slc.tjs[it1].VtxID[end1] > 0) continue;
639  unsigned short endPt1 = slc.tjs[it1].EndPt[end1];
640  for (unsigned short it2 = 0; it2 < slc.tjs.size(); ++it2) {
641  if (it1 == it2) continue;
642  if (slc.tjs[it2].AlgMod[kKilled] || slc.tjs[it2].AlgMod[kHaloTj]) continue;
643  if (slc.tjs[it2].AlgMod[kHamVx]) continue;
644  if (slc.tjs[it2].AlgMod[kHamVx2]) continue;
645  // require that both be in the same CTP
646  if (slc.tjs[it2].CTP != inCTP) continue;
647  if (slc.tjs[it2].AlgMod[kShowerLike]) continue;
648  if (slc.tjs[it2].AlgMod[kJunkTj]) continue;
649  if (slc.tjs[it2].PDGCode == 111) continue;
650  unsigned short numPtsWithCharge2 = NumPtsWithCharge(slc, slc.tjs[it2], true);
651  if (numPtsWithCharge2 < 6) continue;
652  // ignore muon-like tjs
653  if (numPtsWithCharge2 > 100 && slc.tjs[it2].MCSMom > 500) continue;
654  // Find the minimum separation between tj1 and tj2
655  float minDOCA = 5;
656  float doca = minDOCA;
657  unsigned short closePt2 = 0;
658  TrajPointTrajDOCA(slc, slc.tjs[it1].Pts[endPt1], slc.tjs[it2], closePt2, doca);
659  if (doca == minDOCA) continue;
660  if (prt) {
661  mf::LogVerbatim myprt("TC");
662  auto& tj1 = slc.tjs[it1];
663  auto& tj2 = slc.tjs[it2];
664  myprt << " FHV2 CTP" << tj1.CTP;
665  myprt << " t" << tj1.ID << "_" << end1 << " MCSMom " << tj1.MCSMom << " ChgRMS "
666  << tj1.ChgRMS;
667  myprt << " split t" << tj2.ID << "? MCSMom " << tj2.MCSMom << " ChgRMS " << tj2.ChgRMS;
668  myprt << " doca " << doca << " tj2.EndPt[0] " << tj2.EndPt[0] << " closePt2 "
669  << closePt2;
670  myprt << " tj2.EndPt[1] " << tj2.EndPt[1];
671  } // prt
672  // ensure that the closest point is not near an end
673  if (closePt2 < slc.tjs[it2].EndPt[0] + 3) continue;
674  if (closePt2 > slc.tjs[it2].EndPt[1] - 3) continue;
675  // Find the intersection point between the tj1 end and tj2 closest Pt
676  float wint, tint;
677  TrajIntersection(slc.tjs[it1].Pts[endPt1], slc.tjs[it2].Pts[closePt2], wint, tint);
678  // make an angle cut
679  float dang = DeltaAngle(slc.tjs[it1].Pts[endPt1].Ang, slc.tjs[it2].Pts[closePt2].Ang);
680  if (dang < 0.2) continue;
681  // ensure that tj1 doesn't cross tj2 but ensuring that the closest point on tj1 is at closePt1
682  doca = 5;
683  unsigned short closePt1 = 0;
684  TrajPointTrajDOCA(slc, slc.tjs[it2].Pts[closePt2], slc.tjs[it1], closePt1, doca);
685  if (closePt1 != endPt1) continue;
686  if (prt)
687  mf::LogVerbatim("TC") << " intersection W:T " << (int)wint << ":"
688  << (int)(tint / tcc.unitsPerTick) << " dang " << dang;
689  // Find the point on tj2 that is closest to this point, overwriting closePt
690  doca = minDOCA;
691  // the point on tj2 that is closest to the intersection point
692  unsigned short intPt2;
693  TrajClosestApproach(slc.tjs[it2], wint, tint, intPt2, doca);
694  if (prt)
695  mf::LogVerbatim("TC") << " intPt2 on tj2 " << intPt2 << " at Pos "
696  << PrintPos(slc, slc.tjs[it2].Pts[intPt2]) << " doca " << doca;
697  if (doca == minDOCA) continue;
698  // find the MCSMom for both sections of tj2
699  float mcsmom = slc.tjs[it2].MCSMom;
700  float mcsmom1 = MCSMom(slc, slc.tjs[it2], slc.tjs[it2].EndPt[0], intPt2);
701  float mcsmom2 = MCSMom(slc, slc.tjs[it2], intPt2, slc.tjs[it2].EndPt[1]);
702  // require that the both MCSMoms be greater than
703  if (prt)
704  mf::LogVerbatim("TC") << " Check MCSMom after split: mcsmom1 " << mcsmom1 << " mcsmom2 "
705  << mcsmom2;
706  if (mcsmom1 < mcsmom || mcsmom2 < mcsmom) continue;
707  // start scanning for the point on tj2 that has the best IP with the end of tj1 in the direction
708  // from closePt2 -> endPt2
709  short dir = 1;
710  if (intPt2 < closePt2) dir = -1;
711  unsigned short nit = 0;
712  unsigned short ipt = intPt2;
713  float mostChg = slc.tjs[it2].Pts[ipt].Chg;
714  if (prt)
715  mf::LogVerbatim("TC") << " ipt " << ipt << " at Pos "
716  << PrintPos(slc, slc.tjs[it2].Pts[ipt]) << " chg " << mostChg;
717  while (nit < 20) {
718  ipt += dir;
719  if (ipt < 3 || ipt > slc.tjs[it2].Pts.size() - 4) break;
720  float delta = PointTrajDOCA(slc,
721  slc.tjs[it2].Pts[ipt].Pos[0],
722  slc.tjs[it2].Pts[ipt].Pos[1],
723  slc.tjs[it1].Pts[endPt1]);
724  float sep = PosSep(slc.tjs[it2].Pts[ipt].Pos, slc.tjs[it1].Pts[endPt1].Pos);
725  float dang = delta / sep;
726  float pull = dang / slc.tjs[it1].Pts[endPt1].DeltaRMS;
727  if (pull < 2 && slc.tjs[it2].Pts[ipt].Chg > mostChg) {
728  mostChg = slc.tjs[it2].Pts[ipt].Chg;
729  intPt2 = ipt;
730  }
731  }
732  // require a lot of charge in tj2 in this vicinity compared with the average.
733  float chgPull = (mostChg - slc.tjs[it2].AveChg) / slc.tjs[it2].ChgRMS;
734  if (prt) mf::LogVerbatim("TC") << " chgPull at intersection point " << chgPull;
735  if (chgPull < 10) continue;
736  // require a signal on every wire between it1 and intPt2
737  TrajPoint ltp = slc.tjs[it1].Pts[endPt1];
738  if (slc.tjs[it1].Pts[endPt1].Pos[0] < -0.4) continue;
739  unsigned int fromWire = std::nearbyint(slc.tjs[it1].Pts[endPt1].Pos[0]);
740  if (slc.tjs[it2].Pts[intPt2].Pos[0] < -0.4) continue;
741  unsigned int toWire = std::nearbyint(slc.tjs[it2].Pts[intPt2].Pos[0]);
742  if (fromWire > toWire) {
743  unsigned int tmp = fromWire;
744  fromWire = toWire;
745  toWire = tmp;
746  }
747  bool skipIt = false;
748  for (unsigned int wire = fromWire + 1; wire < toWire; ++wire) {
749  MoveTPToWire(ltp, (float)wire);
750  if (!SignalAtTp(ltp)) {
751  skipIt = true;
752  break;
753  }
754  } // wire
755  if (skipIt) continue;
756  // we have a winner
757  // create a new vertex
758  VtxStore aVtx;
759  aVtx.Pos = slc.tjs[it2].Pts[intPt2].Pos;
760  aVtx.NTraj = 3;
761  aVtx.Pass = slc.tjs[it2].Pass;
762  aVtx.Topo = 6;
763  aVtx.ChiDOF = 0;
764  aVtx.CTP = inCTP;
765  aVtx.ID = slc.vtxs.size() + 1;
766  unsigned short ivx = slc.vtxs.size();
767  if (!StoreVertex(slc, aVtx)) continue;
768  if (!SplitTraj(slc, it2, intPt2, ivx, prt)) {
769  if (prt) mf::LogVerbatim("TC") << "FHV2: Failed to split trajectory";
770  MakeVertexObsolete("HamVx2", slc, slc.vtxs[ivx], true);
771  continue;
772  }
773  slc.tjs[it1].VtxID[end1] = slc.vtxs[ivx].ID;
774  slc.tjs[it1].AlgMod[kHamVx2] = true;
775  slc.tjs[it2].AlgMod[kHamVx2] = true;
776  unsigned short newTjIndex = slc.tjs.size() - 1;
777  slc.tjs[newTjIndex].AlgMod[kHamVx2] = true;
778  AttachAnyTrajToVertex(slc, ivx, prt);
779  SetVx2Score(slc);
780  // Update the PDGCode for the chopped trajectory
781  SetPDGCode(slc, it2);
782  // and for the new trajectory
783  SetPDGCode(slc, newTjIndex);
784  if (prt)
785  mf::LogVerbatim("TC") << " FHV2: New vtx 2V" << slc.vtxs[ivx].ID << " Score "
786  << slc.vtxs[ivx].Score;
787  didaSplit = true;
788  break;
789  } // it2
790  if (didaSplit) break;
791  } // end1
792  } // it1
793  } // FindHammerVertices2
794 
795  //////////////////////////////////////////
796  void
797  FindHammerVertices(TCSlice& slc, const CTP_t& inCTP)
798  {
799  // Look for a trajectory that intersects another. Split
800  // the trajectory and make a vertex. The convention used
801  // is shown pictorially here. Trajectory tj1 must be longer
802  // than tj2
803  // tj2 ------
804  // tj1 /
805  // tj1 /
806  // tj1 /
807 
808  if (!tcc.useAlg[kHamVx]) return;
809 
810  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kHamVx]);
811  if (prt) { mf::LogVerbatim("TC") << "Inside HamVx inCTP " << inCTP; }
812 
813  for (unsigned short it1 = 0; it1 < slc.tjs.size(); ++it1) {
814  if (slc.tjs[it1].CTP != inCTP) continue;
815  if (slc.tjs[it1].AlgMod[kKilled] || slc.tjs[it1].AlgMod[kHaloTj]) continue;
816  if (slc.tjs[it1].AlgMod[kShowerLike]) continue;
817  if (slc.tjs[it1].AlgMod[kJunkTj]) continue;
818  if (slc.tjs[it1].PDGCode == 111) continue;
819  // minimum length requirements
820  unsigned short tj1len = slc.tjs[it1].EndPt[1] - slc.tjs[it1].EndPt[0] + 1;
821  if (tj1len < 5) continue;
822  // require high MCSMom
823  if (slc.tjs[it1].MCSMom < 50) continue;
824  if (prt) mf::LogVerbatim("TC") << "FHV: intersection T" << slc.tjs[it1].ID << " candidate";
825  // Check each end of tj1
826  bool didaSplit = false;
827  for (unsigned short end1 = 0; end1 < 2; ++end1) {
828  // vertex assignment exists?
829  if (slc.tjs[it1].VtxID[end1] > 0) continue;
830  unsigned short endPt1 = slc.tjs[it1].EndPt[end1];
831  for (unsigned short it2 = 0; it2 < slc.tjs.size(); ++it2) {
832  if (slc.tjs[it2].CTP != inCTP) continue;
833  if (it1 == it2) continue;
834  if (slc.tjs[it2].AlgMod[kKilled] || slc.tjs[it2].AlgMod[kHaloTj]) continue;
835  if (slc.tjs[it2].AlgMod[kJunkTj]) continue;
836  if (slc.tjs[it2].PDGCode == 111) continue;
837  // length of tj2 cut
838  unsigned short tj2len = slc.tjs[it2].EndPt[1] - slc.tjs[it2].EndPt[0] + 1;
839  if (tj2len < 6) continue;
840  // ignore very long straight trajectories (probably a cosmic muon)
841  if (tj2len > 200 && slc.tjs[it2].PDGCode == 13 && slc.tjs[it2].MCSMom > 600) continue;
842  float minDOCA = 5;
843  minDOCA /= std::abs(slc.tjs[it1].Pts[endPt1].Dir[0]);
844  float doca = minDOCA;
845  unsigned short closePt2 = 0;
846  TrajPointTrajDOCA(slc, slc.tjs[it1].Pts[endPt1], slc.tjs[it2], closePt2, doca);
847  if (doca == minDOCA) continue;
848  // ensure that the closest point is not near an end
849  if (prt)
850  mf::LogVerbatim("TC") << "FHV: Candidate T" << slc.tjs[it1].ID << " T"
851  << slc.tjs[it2].ID << " doca " << doca << " tj2.EndPt[0] "
852  << slc.tjs[it2].EndPt[0] << " closePt2 " << closePt2
853  << " tj2.EndPt[1] " << slc.tjs[it2].EndPt[1];
854  if (closePt2 < slc.tjs[it2].EndPt[0] + 3) continue;
855  if (closePt2 > slc.tjs[it2].EndPt[1] - 3) continue;
856  // make an angle cut
857  float dang = DeltaAngle(slc.tjs[it1].Pts[endPt1].Ang, slc.tjs[it2].Pts[closePt2].Ang);
858  if (prt)
859  mf::LogVerbatim("TC") << " dang " << dang << " imposing a hard cut of 0.4 for now ";
860  if (dang < 0.4) continue;
861  // check the cleanliness in this area
862  std::vector<int> tjids(2);
863  tjids[0] = slc.tjs[it1].ID;
864  tjids[1] = slc.tjs[it2].ID;
865  float chgFrac = ChgFracNearPos(slc, slc.tjs[it2].Pts[closePt2].Pos, tjids);
866  if (prt) mf::LogVerbatim("TC") << " chgFrac " << chgFrac;
867  if (chgFrac < 0.9) continue;
868  Point2_t vxpos = slc.tjs[it2].Pts[closePt2].Pos;
869  // get a better estimate of the vertex position
870  TrajIntersection(
871  slc.tjs[it1].Pts[endPt1], slc.tjs[it2].Pts[closePt2], vxpos[0], vxpos[1]);
872  // and a better estimate of the point on tj2 where the split should be made
873  doca = minDOCA;
874  TrajClosestApproach(slc.tjs[it2], vxpos[0], vxpos[1], closePt2, doca);
875  if (prt)
876  mf::LogVerbatim("TC") << " better pos " << PrintPos(slc, vxpos) << " new closePt2 "
877  << closePt2;
878  // create a new vertex
879  VtxStore aVtx;
880  aVtx.Pos = vxpos;
881  aVtx.NTraj = 3;
882  aVtx.Pass = slc.tjs[it2].Pass;
883  aVtx.Topo = 5;
884  aVtx.ChiDOF = 0;
885  aVtx.CTP = inCTP;
886  aVtx.ID = slc.vtxs.size() + 1;
887  if (!StoreVertex(slc, aVtx)) continue;
888  unsigned short ivx = slc.vtxs.size() - 1;
889  if (!SplitTraj(slc, it2, closePt2, ivx, prt)) {
890  if (prt) mf::LogVerbatim("TC") << "FHV: Failed to split trajectory";
891  MakeVertexObsolete("HamVx", slc, slc.vtxs[slc.vtxs.size() - 1], true);
892  continue;
893  }
894  slc.tjs[it1].VtxID[end1] = slc.vtxs[ivx].ID;
895  slc.tjs[it1].AlgMod[kHamVx] = true;
896  slc.tjs[it2].AlgMod[kHamVx] = true;
897  unsigned short newTjIndex = slc.tjs.size() - 1;
898  slc.tjs[newTjIndex].AlgMod[kHamVx] = true;
899  SetVx2Score(slc);
900  // Update the PDGCode for the chopped trajectory
901  SetPDGCode(slc, it2);
902  // and for the new trajectory
903  SetPDGCode(slc, newTjIndex);
904  if (prt) {
905  auto& vx2 = slc.vtxs[ivx];
906  mf::LogVerbatim myprt("TC");
907  myprt << " new 2V" << vx2.ID << " Score " << vx2.Score << " Tjs";
908  auto tjlist = GetAssns(slc, "2V", vx2.ID, "T");
909  for (auto tid : tjlist)
910  myprt << " t" << tid;
911  }
912  didaSplit = true;
913  break;
914  } // tj2
915  if (didaSplit) break;
916  } // end1
917  } // tj1
918 
919  } // FindHammerVertices
920 
921  //////////////////////////////////////////
922  void
923  SplitTrajCrossingVertices(TCSlice& slc, CTP_t inCTP)
924  {
925  // This is kind of self-explanatory...
926 
927  if (!tcc.useAlg[kSplitTjCVx]) return;
928 
929  if (slc.vtxs.empty()) return;
930  if (slc.tjs.empty()) return;
931 
932  constexpr float docaCut = 4;
933 
934  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kSplitTjCVx]);
935  if (prt) mf::LogVerbatim("TC") << "Inside SplitTrajCrossingVertices inCTP " << inCTP;
936 
937  geo::PlaneID planeID = DecodeCTP(inCTP);
938 
939  unsigned short nTraj = slc.tjs.size();
940  for (unsigned short itj = 0; itj < nTraj; ++itj) {
941  // NOTE: Don't use a reference variable because it may get lost if the tj is split
942  if (slc.tjs[itj].CTP != inCTP) continue;
943  // obsolete trajectory
944  if (slc.tjs[itj].AlgMod[kKilled] || slc.tjs[itj].AlgMod[kHaloTj]) continue;
945  if (slc.tjs[itj].AlgMod[kJunkTj]) continue;
946  if (slc.tjs[itj].AlgMod[kSplitTjCVx]) continue;
947  // too short
948  if (slc.tjs[itj].EndPt[1] < 6) continue;
949  for (unsigned short iv = 0; iv < slc.vtxs.size(); ++iv) {
950  auto& vx2 = slc.vtxs[iv];
951  // obsolete vertex
952  if (vx2.NTraj == 0) continue;
953  // trajectory already associated with vertex
954  if (slc.tjs[itj].VtxID[0] == vx2.ID || slc.tjs[itj].VtxID[1] == vx2.ID) continue;
955  // not in the cryostat/tpc/plane
956  if (slc.tjs[itj].CTP != vx2.CTP) continue;
957  // poor quality
958  if (vx2.Score < tcc.vtx2DCuts[7]) continue;
959  float doca = docaCut;
960  // make the cut significantly larger if the vertex is in a dead
961  // wire gap to get the first TP that is just outside the gap.
962  if (vx2.Stat[kOnDeadWire]) doca = 100;
963  unsigned short closePt = 0;
964  if (!TrajClosestApproach(slc.tjs[itj], vx2.Pos[0], vx2.Pos[1], closePt, doca)) continue;
965  if (vx2.Stat[kOnDeadWire]) {
966  // special handling for vertices in dead wire regions. Find the IP between
967  // the closest point on the Tj and the vertex
968  doca = PointTrajDOCA(slc, vx2.Pos[0], vx2.Pos[1], slc.tjs[itj].Pts[closePt]);
969  }
970  if (doca > docaCut) continue;
971  if (prt)
972  mf::LogVerbatim("TC") << " doca " << doca << " btw T" << slc.tjs[itj].ID << " and 2V"
973  << slc.vtxs[iv].ID << " closePt " << closePt << " in plane "
974  << planeID.Plane;
975  // compare the length of the Tjs used to make the vertex with the length of the
976  // Tj that we want to split. Don't allow a vertex using very short Tjs to split a long
977  // Tj in the 3rd plane
978  auto vxtjs = GetVtxTjIDs(slc, vx2);
979  if (vxtjs.empty()) continue;
980  unsigned short maxPts = 0;
981  // ensure that there is a large angle between a Tj already attached to the vertex and the
982  // tj that we want to split. We might be considering a delta-ray here
983  float maxdang = 0.3;
984  float tjAng = slc.tjs[itj].Pts[closePt].Ang;
985  for (auto tjid : vxtjs) {
986  auto& vtj = slc.tjs[tjid - 1];
987  if (vtj.AlgMod[kDeltaRay]) continue;
988  unsigned short npwc = NumPtsWithCharge(slc, vtj, false);
989  if (npwc > maxPts) maxPts = npwc;
990  unsigned short end = 0;
991  if (vtj.VtxID[1] == slc.vtxs[iv].ID) end = 1;
992  auto& vtp = vtj.Pts[vtj.EndPt[end]];
993  float dang = DeltaAngle(vtp.Ang, tjAng);
994  if (dang > maxdang) maxdang = dang;
995  } // tjid
996  // skip this operation if any of the Tjs in the split list are > 3 * maxPts
997  maxPts *= 3;
998  bool skipit = false;
999  if (NumPtsWithCharge(slc, slc.tjs[itj], false) > maxPts && maxPts < 100) skipit = true;
1000  if (!skipit && maxdang < 0.3) skipit = true;
1001  if (prt)
1002  mf::LogVerbatim("TC") << " maxPts " << maxPts << " vxtjs[0] " << vxtjs[0] << " maxdang "
1003  << maxdang << " skipit? " << skipit;
1004  if (skipit) {
1005  // kill the vertex?
1006  if (doca < 1) MakeVertexObsolete("STCV", slc, vx2, true);
1007  continue;
1008  }
1009 
1010  // make some adjustments to closePt
1011  if (vx2.Stat[kOnDeadWire]) {
1012  // ensure that the tj will be split at the gap. The closePt point may be
1013  // on the wrong side of it
1014  auto& closeTP = slc.tjs[itj].Pts[closePt];
1015  if (slc.tjs[itj].StepDir > 0 && closePt > slc.tjs[itj].EndPt[0]) {
1016  if (closeTP.Pos[0] > vx2.Pos[0]) --closePt;
1017  }
1018  else if (slc.tjs[itj].StepDir < 0 && closePt < slc.tjs[itj].EndPt[1]) {
1019  if (closeTP.Pos[0] < vx2.Pos[0]) ++closePt;
1020  }
1021  }
1022  else {
1023  // improve closePt based on vertex position
1024  // check if closePt and EndPt[1] are the two sides of vertex
1025  // take dot product of closePt-vtx and EndPt[1]-vtx
1026  if ((slc.tjs[itj].Pts[closePt].Pos[0] - vx2.Pos[0]) *
1027  (slc.tjs[itj].Pts[slc.tjs[itj].EndPt[1]].Pos[0] - vx2.Pos[0]) +
1028  (slc.tjs[itj].Pts[closePt].Pos[1] - vx2.Pos[1]) *
1029  (slc.tjs[itj].Pts[slc.tjs[itj].EndPt[1]].Pos[1] - vx2.Pos[1]) <
1030  0 &&
1031  closePt < slc.tjs[itj].EndPt[1] - 1)
1032  ++closePt;
1033  else if ((slc.tjs[itj].Pts[closePt].Pos[0] - vx2.Pos[0]) *
1034  (slc.tjs[itj].Pts[slc.tjs[itj].EndPt[0]].Pos[0] - vx2.Pos[0]) +
1035  (slc.tjs[itj].Pts[closePt].Pos[1] - vx2.Pos[1]) *
1036  (slc.tjs[itj].Pts[slc.tjs[itj].EndPt[0]].Pos[1] - slc.vtxs[iv].Pos[1]) <
1037  0 &&
1038  closePt > slc.tjs[itj].EndPt[0] + 1)
1039  --closePt;
1040  }
1041 
1042  if (prt) {
1043  mf::LogVerbatim("TC") << "Good doca " << doca << " btw T" << slc.tjs[itj].ID << " and 2V"
1044  << vx2.ID << " closePt " << closePt << " in plane " << planeID.Plane
1045  << " CTP " << slc.vtxs[iv].CTP;
1046  PrintTP("STCV", slc, closePt, 1, slc.tjs[itj].Pass, slc.tjs[itj].Pts[closePt]);
1047  }
1048  // ensure that the closest point is not near an end
1049  if (closePt < slc.tjs[itj].EndPt[0] + 3) continue;
1050  if (closePt > slc.tjs[itj].EndPt[1] - 3) continue;
1051  if (!SplitTraj(slc, itj, closePt, iv, prt)) {
1052  if (prt) mf::LogVerbatim("TC") << "SplitTrajCrossingVertices: Failed to split trajectory";
1053  continue;
1054  }
1055  slc.tjs[itj].AlgMod[kSplitTjCVx] = true;
1056  unsigned short newTjIndex = slc.tjs.size() - 1;
1057  slc.tjs[newTjIndex].AlgMod[kSplitTjCVx] = true;
1058  // re-fit the vertex position
1059  FitVertex(slc, vx2, prt);
1060  } // iv
1061  } // itj
1062 
1063  } // SplitTrajCrossingVertices
1064 
1065  //////////////////////////////////////
1066  void
1067  Reconcile2Vs(TCSlice& slc)
1068  {
1069  // This function is called before Find3DVertices to identify (and possibly reconcile)
1070  // Tj and 2V inconsistencies using 2D and 3D(?) information
1071  if (!tcc.useAlg[kReconcile2Vs]) return;
1072  if (slc.vtxs.empty()) return;
1073 
1074  bool prt = (tcc.dbg2V && tcc.dbgSlc);
1075 
1076  // clusters of 2Vs
1077  std::vector<std::vector<int>> vx2Cls;
1078 
1079  // iterate over planes
1080  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
1081  // look for 2D vertices that are close to each other
1082  vx2Cls.clear();
1083  for (unsigned short ii = 0; ii < slc.vtxs.size() - 1; ++ii) {
1084  auto& i2v = slc.vtxs[ii];
1085  if (i2v.ID <= 0) continue;
1086  if (DecodeCTP(i2v.CTP).Plane != plane) continue;
1087  for (unsigned short jj = ii + 1; jj < slc.vtxs.size(); ++jj) {
1088  auto& j2v = slc.vtxs[jj];
1089  if (j2v.ID <= 0) continue;
1090  if (DecodeCTP(j2v.CTP).Plane != plane) continue;
1091  // make rough separation cuts
1092  float dp0 = std::abs(i2v.Pos[0] - j2v.Pos[0]);
1093  if (dp0 > 10) continue;
1094  float dp1 = std::abs(i2v.Pos[1] - j2v.Pos[1]);
1095  if (dp1 > 10) continue;
1096  // do a more careful look
1097  float err = i2v.PosErr[0];
1098  if (j2v.PosErr[0] > err) err = j2v.PosErr[0];
1099  float dp0Sig = dp0 / err;
1100  if (dp0Sig > 4) continue;
1101  err = i2v.PosErr[1];
1102  if (j2v.PosErr[1] > err) err = j2v.PosErr[1];
1103  float dp1Sig = dp1 / err;
1104  if (dp1Sig > 4) continue;
1105  // Look for one of the 2V IDs in a cluster
1106  bool gotit = false;
1107  for (auto& vx2cls : vx2Cls) {
1108  bool goti = (std::find(vx2cls.begin(), vx2cls.end(), i2v.ID) != vx2cls.end());
1109  bool gotj = (std::find(vx2cls.begin(), vx2cls.end(), j2v.ID) != vx2cls.end());
1110  if (goti && gotj) {
1111  gotit = true;
1112  break;
1113  }
1114  else if (goti) {
1115  vx2cls.push_back(j2v.ID);
1116  gotit = true;
1117  break;
1118  }
1119  else if (gotj) {
1120  gotit = true;
1121  vx2cls.push_back(i2v.ID);
1122  break;
1123  }
1124  } // vx2cls
1125  if (!gotit) {
1126  // start a new cluster with this pair
1127  std::vector<int> cls(2);
1128  cls[0] = i2v.ID;
1129  cls[1] = j2v.ID;
1130  vx2Cls.push_back(cls);
1131  } // !gotit
1132  } // jj
1133  } // ii
1134  if (vx2Cls.empty()) continue;
1135  if (prt) {
1136  mf::LogVerbatim myprt("TC");
1137  myprt << "2V clusters in plane " << plane;
1138  for (auto& vx2cls : vx2Cls) {
1139  myprt << "\n";
1140  for (auto vx2id : vx2cls)
1141  myprt << " 2V" << vx2id;
1142  } // vx2cls
1143  } // prt
1144  for (auto& vx2cls : vx2Cls) {
1145  Reconcile2VTs(slc, vx2cls, prt);
1146  } // vx2cls
1147  } // plane
1148 
1149  // See if any of the vertices have been altered. If so the environment near them,
1150  // specifically tagging overlapping trajectories, should be re-done
1151  bool VxEnvironmentNeedsUpdate = false;
1152  for (auto& vx : slc.vtxs) {
1153  if (vx.ID <= 0) continue;
1154  if (!vx.Stat[kVxEnvOK]) VxEnvironmentNeedsUpdate = true;
1155  } // vx
1156 
1157  if (VxEnvironmentNeedsUpdate) UpdateVxEnvironment(slc);
1158 
1159  } // Reconcile2Vs
1160 
1161  //////////////////////////////////////
1162  bool
1163  Reconcile2VTs(TCSlice& slc, std::vector<int>& vx2cls, bool prt)
1164  {
1165  // The 2D vertices IDs in vx2cls were clustered by the calling function. This function
1166  // checks the T -> 2V assns and possibly changes it. It returns true if an assn is changed
1167  // or if a vertex in vx2cls is made obsolete, necessitating a change to the list of 2V
1168  // clusters
1169  if (vx2cls.size() < 2) return false;
1170 
1171  // Form a list of all Tjs associated with this 2V cluster
1172  std::vector<int> t2vList;
1173 
1174  CTP_t inCTP;
1175  for (auto vx2id : vx2cls) {
1176  auto& vx2 = slc.vtxs[vx2id - 1];
1177  inCTP = vx2.CTP;
1178  // vertex clobbered? If so, vertex clustering needs to be re-done
1179  if (vx2.ID <= 0) return true;
1180  auto tlist = GetAssns(slc, "2V", vx2.ID, "T");
1181  for (auto tid : tlist)
1182  if (std::find(t2vList.begin(), t2vList.end(), tid) == t2vList.end()) t2vList.push_back(tid);
1183  } // vx2id
1184  if (t2vList.size() < 3) return false;
1185 
1186  // Sum the T -> 2V pulls
1187  float sumPulls = 0;
1188  float cnt = 0;
1189  for (auto tid : t2vList) {
1190  auto& tj = slc.tjs[tid - 1];
1191  for (unsigned short end = 0; end < 2; ++end) {
1192  if (tj.VtxID[end] <= 0) continue;
1193  if (std::find(vx2cls.begin(), vx2cls.end(), tj.VtxID[end]) == vx2cls.end()) continue;
1194  auto& vx = slc.vtxs[tj.VtxID[end] - 1];
1195  unsigned short nearEnd = 1 - FarEnd(slc, tj, vx.Pos);
1196  unsigned short fitPt = NearbyCleanPt(slc, tj, nearEnd);
1197  if (fitPt == USHRT_MAX) return false;
1198  auto& tp = tj.Pts[fitPt];
1199  sumPulls += TrajPointVertexPull(slc, tp, vx);
1200  ++cnt;
1201  } // end
1202  } // tid
1203 
1204  if (prt) {
1205  mf::LogVerbatim myprt("TC");
1206  myprt << "R2VTs: cluster:";
1207  for (auto vid : vx2cls)
1208  myprt << " 2V" << vid;
1209  myprt << " ->";
1210  for (auto tid : t2vList)
1211  myprt << " T" << tid;
1212  myprt << " sumPulls " << std::setprecision(2) << sumPulls << " cnt " << cnt;
1213  } // prt
1214 
1215  // try to fit all Tjs to one vertex. Find the average position of all of the
1216  // vertices. This will be used to find the end of the Tjs that are closest to the
1217  // presumed single vertex
1218  VtxStore oneVx;
1219  oneVx.CTP = inCTP;
1220  for (auto vid : vx2cls) {
1221  auto& vx = slc.vtxs[vid - 1];
1222  oneVx.Pos[0] += vx.Pos[0];
1223  oneVx.Pos[1] += vx.Pos[2];
1224  } // vid
1225  oneVx.Pos[0] /= vx2cls.size();
1226  oneVx.Pos[1] /= vx2cls.size();
1227  std::vector<TrajPoint> oneVxTPs(t2vList.size());
1228  for (unsigned short itj = 0; itj < t2vList.size(); ++itj) {
1229  auto& tj = slc.tjs[t2vList[itj] - 1];
1230  unsigned short nearEnd = 1 - FarEnd(slc, tj, oneVx.Pos);
1231  unsigned short fitPt = NearbyCleanPt(slc, tj, nearEnd);
1232  if (fitPt == USHRT_MAX) return false;
1233  oneVxTPs[itj] = tj.Pts[fitPt];
1234  // inflate the TP angle angle error if a TP without an overlap wasn't found
1235  if (oneVxTPs[itj].Environment[kEnvOverlap]) oneVxTPs[itj].AngErr *= 4;
1236  oneVxTPs[itj].Step = tj.ID;
1237  } // ii
1238  if (!FitVertex(slc, oneVx, oneVxTPs, prt)) return false;
1239 
1240  if (oneVx.ChiDOF < 3) {
1241  // Update the position of the first 2V in the list
1242  auto& vx = slc.vtxs[vx2cls[0] - 1];
1243  vx.Pos = oneVx.Pos;
1244  vx.PosErr = oneVx.PosErr;
1245  vx.NTraj = t2vList.size();
1246  vx.ChiDOF = oneVx.ChiDOF;
1247  vx.Topo = 14;
1248  // Set a flag that the environment near this vertex (and all other vertices in this slice)
1249  // should be revisited
1250  vx.Stat[kVxEnvOK] = false;
1251  for (unsigned short ivx = 1; ivx < vx2cls.size(); ++ivx) {
1252  auto& vx = slc.vtxs[vx2cls[ivx] - 1];
1253  MakeVertexObsolete("R2VTPs", slc, vx, true);
1254  } // ivx
1255  // now attach the trajectories
1256  for (auto tid : t2vList) {
1257  auto& tj = slc.tjs[tid - 1];
1258  unsigned short nearEnd = 1 - FarEnd(slc, tj, vx.Pos);
1259  tj.VtxID[nearEnd] = vx.ID;
1260  } // tid
1261  return true;
1262  } // oneVx.ChiDOF < 3
1263  return false;
1264  } // Reconcile2VTs
1265 
1266  //////////////////////////////////////
1267  void
1268  Find3DVertices(detinfo::DetectorPropertiesData const& detProp, TCSlice& slc)
1269  {
1270  // Create 3D vertices from 2D vertices. 3D vertices that are matched
1271  // in all three planes have Vtx2ID > 0 for all planes. This function re-scores all
1272  // 2D and 3D vertices and flags Tjs that have high-score 3D vertices
1273 
1274  if (tcc.vtx3DCuts[0] < 0) return;
1275  if (slc.vtxs.size() < 2) return;
1276 
1277  // create a array/vector of 2D vertex indices in each plane
1278  std::vector<std::vector<unsigned short>> vIndex(3);
1279  for (unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
1280  // obsolete vertex
1281  if (slc.vtxs[ivx].ID == 0) continue;
1282  geo::PlaneID planeID = DecodeCTP(slc.vtxs[ivx].CTP);
1283  if (planeID.TPC != slc.TPCID.TPC) continue;
1284  unsigned short plane = planeID.Plane;
1285  if (plane > 2) continue;
1286  vIndex[plane].push_back(ivx);
1287  }
1288 
1289  unsigned short vtxInPln = 0;
1290  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane)
1291  if (vIndex[plane].size() > 0) ++vtxInPln;
1292  if (vtxInPln < 2) return;
1293 
1294  float thirdPlanedXCut = 2 * tcc.vtx3DCuts[0];
1295  bool prt = (tcc.dbg3V && tcc.dbgSlc);
1296  if (prt) {
1297  mf::LogVerbatim("TC") << "Inside Find3DVertices. dX cut " << tcc.vtx3DCuts[0]
1298  << " thirdPlanedXCut " << thirdPlanedXCut;
1299  }
1300 
1301  size_t vsize = slc.vtxs.size();
1302  // vector of 2D vertices -> 3D vertices.
1303  std::vector<short> vPtr(vsize, -1);
1304  // fill temp vectors of 2D vertex X and X errors
1305  std::vector<float> vX(vsize, FLT_MAX);
1306 
1307  for (unsigned short ivx = 0; ivx < vsize; ++ivx) {
1308  if (slc.vtxs[ivx].ID <= 0) continue;
1309  if (slc.vtxs[ivx].Score < tcc.vtx2DCuts[7]) continue;
1310  if (slc.vtxs[ivx].Pos[0] < -0.4) continue;
1311  geo::PlaneID planeID = DecodeCTP(slc.vtxs[ivx].CTP);
1312  // Convert 2D vertex time error to X error
1313  double ticks = slc.vtxs[ivx].Pos[1] / tcc.unitsPerTick;
1314  vX[ivx] = detProp.ConvertTicksToX(ticks, planeID);
1315  } // ivx
1316 
1317  // temp vector of all 2D vertex matches
1318  std::vector<Vtx3Store> v3temp;
1319 
1320  unsigned int cstat = slc.TPCID.Cryostat;
1321  unsigned int tpc = slc.TPCID.TPC;
1322 
1323  TrajPoint tp;
1324  constexpr float maxSep = 4;
1325  float maxScore = 0;
1326  // i, j, k indicates 3 different wire planes
1327  // compare vertices in each view
1328  for (unsigned short ipl = 0; ipl < 2; ++ipl) {
1329  for (unsigned short ii = 0; ii < vIndex[ipl].size(); ++ii) {
1330  unsigned short ivx = vIndex[ipl][ii];
1331  if (vX[ivx] == FLT_MAX) continue;
1332  auto& ivx2 = slc.vtxs[ivx];
1333  if (ivx2.Pos[0] < -0.4) continue;
1334  unsigned int iWire = std::nearbyint(ivx2.Pos[0]);
1335  for (unsigned short jpl = ipl + 1; jpl < 3; ++jpl) {
1336  for (unsigned short jj = 0; jj < vIndex[jpl].size(); ++jj) {
1337  unsigned short jvx = vIndex[jpl][jj];
1338  if (vX[jvx] == FLT_MAX) continue;
1339  auto& jvx2 = slc.vtxs[jvx];
1340  if (jvx2.Pos[0] < -0.4) continue;
1341  unsigned int jWire = std::nearbyint(jvx2.Pos[0]);
1342  float dX = std::abs(vX[ivx] - vX[jvx]);
1343  if (dX > tcc.vtx3DCuts[0]) continue;
1344  if (prt) {
1345  mf::LogVerbatim("TC")
1346  << "F3DV: ipl " << ipl << " i2V" << ivx2.ID << " iX " << vX[ivx] << " jpl " << jpl
1347  << " j2V" << jvx2.ID << " jvX " << vX[jvx] << " W:T " << (int)jvx2.Pos[0] << ":"
1348  << (int)jvx2.Pos[1] << " dX " << dX;
1349  }
1350  double y = -1000, z = -1000;
1351  tcc.geom->IntersectionPoint(iWire, jWire, ipl, jpl, cstat, tpc, y, z);
1352  if (y < slc.yLo || y > slc.yHi || z < slc.zLo || z > slc.zHi) continue;
1353  unsigned short kpl = 3 - ipl - jpl;
1354  float kX = 0.5 * (vX[ivx] + vX[jvx]);
1355  int kWire = -1;
1356  if (slc.nPlanes > 2) {
1357  kWire = tcc.geom->WireCoordinate(y, z, kpl, tpc, cstat);
1358  std::array<int, 2> wireWindow;
1359  std::array<float, 2> timeWindow;
1360  wireWindow[0] = kWire - maxSep;
1361  wireWindow[1] = kWire + maxSep;
1362  float time =
1363  detProp.ConvertXToTicks(kX, kpl, (int)tpc, (int)cstat) * tcc.unitsPerTick;
1364  timeWindow[0] = time - maxSep;
1365  timeWindow[1] = time + maxSep;
1366  bool hitsNear;
1367  std::vector<unsigned int> closeHits =
1368  FindCloseHits(slc, wireWindow, timeWindow, kpl, kAllHits, true, hitsNear);
1369  if (prt) {
1370  mf::LogVerbatim myprt("TC");
1371  myprt << " Hits near " << kpl << ":" << kWire << ":"
1372  << (int)(time / tcc.unitsPerTick) << " = ";
1373  for (auto iht : closeHits)
1374  myprt << " " << PrintHit(slc.slHits[iht]);
1375  }
1376  if (!hitsNear) continue;
1377  } // 3-plane TPC
1378  // save this incomplete 3D vertex
1379  Vtx3Store v3d;
1380  // give it a non-zero ID so that SetVx3Score returns a valid score
1381  v3d.ID = 666;
1382  v3d.Vx2ID.resize(slc.nPlanes);
1383  v3d.Vx2ID[ipl] = ivx2.ID;
1384  v3d.Vx2ID[jpl] = jvx2.ID;
1385  if (slc.nPlanes == 2) v3d.Vx2ID[2] = -1;
1386  v3d.X = kX;
1387  // Use XErr to store dX
1388  v3d.XErr = dX;
1389  v3d.Y = y;
1390  v3d.Z = z;
1391  v3d.Wire = kWire;
1392  float posError = dX / tcc.vtx3DCuts[0];
1393  float vxScoreWght = 0;
1394  SetVx3Score(slc, v3d);
1395  vxScoreWght = tcc.vtx3DCuts[2] / v3d.Score;
1396  if (posError < 0.5) posError = 0;
1397  v3d.Score = posError + vxScoreWght;
1398  v3d.TPCID = slc.TPCID;
1399  // push the incomplete vertex onto the list
1400  v3temp.push_back(v3d);
1401 
1402  if (prt) {
1403  mf::LogVerbatim myprt("TC");
1404  myprt << " 2 Plane match i2V";
1405  myprt << slc.vtxs[ivx].ID << " P:W:T " << ipl << ":" << (int)slc.vtxs[ivx].Pos[0]
1406  << ":" << (int)slc.vtxs[ivx].Pos[1];
1407  myprt << " j2V" << slc.vtxs[jvx].ID << " P:W:T " << jpl << ":"
1408  << (int)slc.vtxs[jvx].Pos[0] << ":" << (int)slc.vtxs[jvx].Pos[1];
1409  myprt << std::fixed << std::setprecision(3);
1410  myprt << " dX " << dX << " posError " << posError << " vxScoreWght " << vxScoreWght
1411  << " Score " << v3d.Score;
1412  }
1413 
1414  if (slc.nPlanes == 2) continue;
1415 
1416  // look for a 3 plane match
1417  for (unsigned short kk = 0; kk < vIndex[kpl].size(); ++kk) {
1418  unsigned short kvx = vIndex[kpl][kk];
1419  float dX = std::abs(vX[kvx] - v3d.X);
1420  // Wire difference error
1421  float dW = tcc.wirePitch * std::abs(slc.vtxs[kvx].Pos[0] - kWire);
1422  if (dX > thirdPlanedXCut) continue;
1423  if (dW > tcc.vtx3DCuts[1]) continue;
1424  // put the Y,Z difference in YErr and ZErr
1425  double y = -1000, z = -1000;
1426  tcc.geom->IntersectionPoint(iWire, kWire, ipl, kpl, cstat, tpc, y, z);
1427  v3d.YErr = y - v3d.Y;
1428  v3d.ZErr = z - v3d.Z;
1429  v3d.Vx2ID[kpl] = slc.vtxs[kvx].ID;
1430  v3d.Wire = -1;
1431  // hijack the Score variable to hold the separation^2, weighted by the
1432  // vertex3DCuts
1433  dX = (vX[kvx] - v3d.X) / tcc.vtx3DCuts[0];
1434  float dY = v3d.YErr / tcc.vtx3DCuts[1];
1435  float dZ = v3d.ZErr / tcc.vtx3DCuts[1];
1436  posError = dX * dX + dY * dY + dZ * dZ;
1437  vxScoreWght = 0;
1438  SetVx3Score(slc, v3d);
1439  vxScoreWght = tcc.vtx3DCuts[2] / v3d.Score;
1440  if (posError < 0.5) posError = 0;
1441  v3d.Score = posError + vxScoreWght;
1442  if (v3d.Score > maxScore) maxScore = v3d.Score;
1443  if (prt)
1444  mf::LogVerbatim("TC") << " k2V" << kvx + 1 << " dX " << dX << " dW " << dW
1445  << " 3D score " << v3d.Score;
1446  v3temp.push_back(v3d);
1447  } // kk
1448  } // jj
1449  } // jpl
1450  } // ii
1451  } // ipl
1452 
1453  if (v3temp.empty()) return;
1454 
1455  // We will sort this list by increasing score. First add the maxScore for 2-plane matches so that
1456  // they are considered after the 3-plane matches
1457  maxScore += 1;
1458  for (auto& v3 : v3temp)
1459  if (v3.Wire >= 0) v3.Score += maxScore;
1460 
1461  if (prt) {
1462  mf::LogVerbatim("TC") << "v3temp list";
1463  for (auto& v3 : v3temp) {
1464  if (slc.nPlanes == 2) {
1465  mf::LogVerbatim("TC") << "2V" << v3.Vx2ID[0] << " 2V" << v3.Vx2ID[1] << " wire "
1466  << v3.Wire << " Score " << v3.Score;
1467  }
1468  else {
1469  mf::LogVerbatim("TC") << "2V" << v3.Vx2ID[0] << " 2V" << v3.Vx2ID[1] << " 2V"
1470  << v3.Vx2ID[2] << " wire " << v3.Wire << " Score " << v3.Score;
1471  }
1472  } // v3
1473  }
1474  SortEntry sEntry;
1475  std::vector<SortEntry> sortVec(v3temp.size());
1476  for (unsigned short ivx = 0; ivx < v3temp.size(); ++ivx) {
1477  sEntry.index = ivx;
1478  sEntry.val = v3temp[ivx].Score;
1479  sortVec[ivx] = sEntry;
1480  } // ivx
1481  if (sortVec.size() > 1) std::sort(sortVec.begin(), sortVec.end(), valsIncreasing);
1482  // create a new vector of selected 3D vertices
1483  std::vector<Vtx3Store> v3sel;
1484  for (unsigned short ii = 0; ii < sortVec.size(); ++ii) {
1485  unsigned short ivx = sortVec[ii].index;
1486  // ensure that all 2D vertices are unique
1487  bool skipit = false;
1488  for (auto& v3 : v3sel) {
1489  for (unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
1490  if (v3temp[ivx].Vx2ID[ipl] == 0) continue;
1491  if (v3temp[ivx].Vx2ID[ipl] == v3.Vx2ID[ipl]) {
1492  skipit = true;
1493  break;
1494  }
1495  } // ipl
1496  if (skipit) break;
1497  } // v3
1498  if (skipit) continue;
1499  v3sel.push_back(v3temp[ivx]);
1500  } // ii
1501  v3temp.clear();
1502 
1503  if (prt) {
1504  mf::LogVerbatim myprt("TC");
1505  myprt << "v3sel list\n";
1506  for (auto& v3d : v3sel) {
1507  for (auto vx2id : v3d.Vx2ID)
1508  if (vx2id > 0) myprt << " 2V" << vx2id;
1509  myprt << " wire " << v3d.Wire << " Score " << v3d.Score;
1510  myprt << "\n";
1511  } // v3d
1512  } // prt
1513 
1514  // Count the number of incomplete vertices and store
1515  unsigned short ninc = 0;
1516  for (auto& vx3 : v3sel) {
1517  if (slc.nPlanes == 3 && vx3.Wire >= 0) ++ninc;
1518  vx3.ID = slc.vtx3s.size() + 1;
1519  ++evt.global3V_UID;
1520  vx3.UID = evt.global3V_UID;
1521  if (prt) {
1522  mf::LogVerbatim myprt("TC");
1523  myprt << " 3V" << vx3.ID;
1524  for (auto v2id : vx3.Vx2ID)
1525  myprt << " 2V" << v2id;
1526  myprt << " wire " << vx3.Wire;
1527  } // prt
1528  slc.vtx3s.push_back(vx3);
1529  // make the 2D -> 3D associations
1530  for (unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
1531  if (vx3.Vx2ID[ipl] == 0) continue;
1532  VtxStore& vx2 = slc.vtxs[vx3.Vx2ID[ipl] - 1];
1533  vx2.Vx3ID = vx3.ID;
1534  } // ipl
1535  } // ivx
1536 
1537  // Try to complete incomplete vertices
1538  if (ninc > 0) {
1539  CompleteIncomplete3DVerticesInGaps(detProp, slc);
1540  CompleteIncomplete3DVertices(detProp, slc);
1541  }
1542 
1543  // Score and flag Tjs that are attached to high-score vertices
1544  // First remove Tj vertex flags
1545  for (auto& tj : slc.tjs) {
1546  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1547  geo::PlaneID planeID = DecodeCTP(tj.CTP);
1548  if (planeID.TPC != tpc || planeID.Cryostat != cstat) continue;
1549  tj.AlgMod[kTjHiVx3Score] = false;
1550  } // tj
1551  // Score the 2D vertices
1552  for (auto& vx2 : slc.vtxs) {
1553  if (vx2.ID == 0) continue;
1554  geo::PlaneID planeID = DecodeCTP(vx2.CTP);
1555  if (planeID.TPC != tpc || planeID.Cryostat != cstat) continue;
1556  SetVx2Score(slc, vx2);
1557  } // vx2
1558  // and the 3D vertices
1559  for (auto& vx3 : slc.vtx3s) {
1560  if (vx3.ID == 0) continue;
1561  SetVx3Score(slc, vx3);
1562  } // vx3
1563 
1564  } // Find3DVertices
1565 
1566  //////////////////////////////////////////
1567  unsigned short
1568  TPNearVertex(const TCSlice& slc, const TrajPoint& tp)
1569  {
1570  // Returns the index of a vertex if tp is nearby
1571  for (unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
1572  if (slc.vtxs[ivx].ID == 0) continue;
1573  if (slc.vtxs[ivx].CTP != tp.CTP) continue;
1574  if (std::abs(slc.vtxs[ivx].Pos[0] - tp.Pos[0]) > 1.2) continue;
1575  if (std::abs(slc.vtxs[ivx].Pos[1] - tp.Pos[1]) > 1.2) continue;
1576  return ivx;
1577  } // ivx
1578  return USHRT_MAX;
1579  } // TPNearVertex
1580 
1581  //////////////////////////////////////////
1582  bool
1583  AttachToAnyVertex(TCSlice& slc, PFPStruct& pfp, float maxSep, bool prt)
1584  {
1585  // Attaches to any 3D vertex but doesn't require consistency with
1586  // PFP -> Tj -> 2V -> 3V assns
1587  if (pfp.ID <= 0) return false;
1588 
1589  float pLen = Length(pfp);
1590  if (pLen == 0) return false;
1591 
1592  // save the old assignents and clear them
1593  // auto oldVx3ID = pfp.Vx3ID;
1594  for (unsigned short end = 0; end < 2; ++end)
1595  pfp.Vx3ID[end] = 0;
1596  std::array<Point3_t, 2> endPos;
1597  endPos[0] = PosAtEnd(pfp, 0);
1598  endPos[1] = PosAtEnd(pfp, 1);
1599 
1600  std::array<float, 2> foms{{100., 100.}};
1601  std::array<int, 2> vtxs{{0}};
1602  for (auto& vx3 : slc.vtx3s) {
1603  if (vx3.ID <= 0) continue;
1604  if (vx3.TPCID != pfp.TPCID) continue;
1605  std::array<float, 2> sep;
1606  Point3_t vpos = {{vx3.X, vx3.Y, vx3.Z}};
1607  sep[0] = PosSep(vpos, endPos[0]);
1608  sep[1] = PosSep(vpos, endPos[1]);
1609  unsigned short end = 0;
1610  if (sep[1] < sep[0]) end = 1;
1611  // ignore if separation is too large
1612  if (sep[end] > 100) continue;
1613  // find the direction vector between these points
1614  auto vpDir = PointDirection(vpos, endPos[end]);
1615  auto dir = DirAtEnd(pfp, end);
1616  double dotp = std::abs(DotProd(vpDir, dir));
1617  float fom = dotp * sep[end];
1618  if (prt)
1619  mf::LogVerbatim("TC") << "ATAV: P" << pfp.ID << " end " << end << " 3V" << vx3.ID << " sep "
1620  << sep[end] << " fom " << fom << " maxSep " << maxSep;
1621  // ignore if separation is too large
1622  if (sep[end] > maxSep) continue;
1623  if (fom < foms[end]) {
1624  foms[end] = fom;
1625  vtxs[end] = vx3.ID;
1626  }
1627  } // vx3
1628  bool bingo = false;
1629  for (unsigned short end = 0; end < 2; ++end) {
1630  if (vtxs[end] == 0) continue;
1631  if (prt)
1632  mf::LogVerbatim("TC") << "ATAV: set P" << pfp.ID << " end " << end << " -> 3V" << vtxs[end];
1633  pfp.Vx3ID[end] = vtxs[end];
1634  bingo = true;
1635  } // end
1636  return bingo;
1637  } // AttachToAnyVertex
1638 
1639  //////////////////////////////////////////
1640  bool
1641  AttachAnyVertexToTraj(TCSlice& slc, int tjID, bool prt)
1642  {
1643  // Try to attach a tj that is stored in slc.tjs with any vertex
1644  if (tjID <= 0 || tjID > (int)slc.tjs.size()) return false;
1645  if (slc.vtxs.empty()) return false;
1646  auto& tj = slc.tjs[tjID - 1];
1647  if (tj.AlgMod[kKilled]) return false;
1648  if (tcc.vtx2DCuts[0] <= 0) return false;
1649 
1650  unsigned short bestVx = USHRT_MAX;
1651  // Construct a FOM = (TP-Vtx pull) * (TP-Vtx sep + 1) * (Vtx Score).
1652  // The +1 keeps FOM from being 0
1653  float bestFOM = 2 * tcc.vtx2DCuts[3] * (tcc.vtx2DCuts[0] + 1) * tcc.vtx2DCuts[7];
1654  for (unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
1655  auto& vx = slc.vtxs[ivx];
1656  if (vx.ID == 0) continue;
1657  if (vx.CTP != tj.CTP) continue;
1658  // make some rough cuts
1659  std::array<float, 2> sep;
1660  sep[0] = PosSep(vx.Pos, tj.Pts[tj.EndPt[0]].Pos);
1661  sep[1] = PosSep(vx.Pos, tj.Pts[tj.EndPt[1]].Pos);
1662  unsigned short end = 0;
1663  if (sep[1] < sep[0]) end = 1;
1664  if (sep[end] > 100) continue;
1665  if (tj.VtxID[end] > 0) continue;
1666  auto& tp = tj.Pts[tj.EndPt[end]];
1667  // Pad the separation a bit so we don't get zero
1668  float fom = TrajPointVertexPull(slc, tp, vx) * (sep[end] + 1) * vx.Score;
1669  if (fom > bestFOM) continue;
1670  if (prt)
1671  mf::LogVerbatim("TC") << "AAVTT: T" << tjID << " 2V" << vx.ID << " FOM " << fom << " cut "
1672  << bestFOM;
1673  bestVx = ivx;
1674  bestFOM = fom;
1675  } // vx
1676  if (bestVx > slc.vtxs.size() - 1) return false;
1677  auto& vx = slc.vtxs[bestVx];
1678  return AttachTrajToVertex(slc, tj, vx, prt);
1679  } // AttachAnyVertexToTraj
1680 
1681  //////////////////////////////////////////
1682  bool
1683  AttachAnyTrajToVertex(TCSlice& slc, unsigned short ivx, bool prt)
1684  {
1685 
1686  if (ivx > slc.vtxs.size() - 1) return false;
1687  if (slc.vtxs[ivx].ID == 0) return false;
1688  if (tcc.vtx2DCuts[0] < 0) return false;
1689 
1690  VtxStore& vx = slc.vtxs[ivx];
1691  // Hammer vertices should be isolated and clean
1692  if (vx.Topo == 5 || vx.Topo == 6) return false;
1693 
1694  unsigned short bestTj = USHRT_MAX;
1695  // Construct a FOM = (TP-Vtx pull) * (TP-Vtx sep + 1).
1696  // The +1 keeps FOM from being 0
1697  float bestFOM = 2 * tcc.vtx2DCuts[3] * (tcc.vtx2DCuts[0] + 1);
1698  for (unsigned int itj = 0; itj < slc.tjs.size(); ++itj) {
1699  auto& tj = slc.tjs[itj];
1700  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1701  if (tj.CTP != vx.CTP) continue;
1702  // make some rough cuts
1703  std::array<float, 2> sep;
1704  sep[0] = PosSep(vx.Pos, tj.Pts[tj.EndPt[0]].Pos);
1705  sep[1] = PosSep(vx.Pos, tj.Pts[tj.EndPt[1]].Pos);
1706  unsigned short end = 0;
1707  if (sep[1] < sep[0]) end = 1;
1708  if (sep[end] > 100) continue;
1709  if (tj.VtxID[end] > 0) continue;
1710  auto& tp = tj.Pts[tj.EndPt[end]];
1711  // Pad the separation a bit so we don't get zero
1712  float fom = TrajPointVertexPull(slc, tp, vx) * (sep[end] + 1);
1713  if (fom > bestFOM) continue;
1714  if (prt) {
1715  mf::LogVerbatim("TC") << "AATTV: T" << tj.ID << " 2V" << vx.ID << " Topo " << vx.Topo
1716  << " FOM " << fom << " cut " << bestFOM;
1717  }
1718  bestTj = itj;
1719  bestFOM = fom;
1720  } // tj
1721  if (bestTj > slc.tjs.size() - 1) return false;
1722  auto& tj = slc.tjs[bestTj];
1723  return AttachTrajToVertex(slc, tj, vx, prt);
1724  } // AttachAnyTrajToVertex
1725 
1726  //////////////////////////////////////////
1727  bool
1728  AttachTrajToVertex(TCSlice& slc, Trajectory& tj, VtxStore& vx, bool prt)
1729  {
1730  // Note that this function does not require a signal between the end of the Tj and the vertex
1731 
1732  // tcc.vtx2DCuts fcl input usage
1733  // 0 = maximum length of a short trajectory
1734  // 1 = max vertex - trajectory separation for short trajectories
1735  // 2 = max vertex - trajectory separation for long trajectories
1736  // 3 = max position pull for adding TJs to a vertex
1737  // 4 = max allowed vertex position error
1738  // 5 = min MCSMom
1739  // 6 = min Pts/Wire fraction
1740 
1741  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) return false;
1742  if (tj.CTP != vx.CTP) return false;
1743  // already attached?
1744  if (tj.VtxID[0] == vx.ID || tj.VtxID[1] == vx.ID) return false;
1745 
1746  unsigned short maxShortTjLen = tcc.vtx2DCuts[0];
1747  // square the separation cut to simplify testing in the loop
1748  float maxSepCutShort2 = tcc.vtx2DCuts[1] * tcc.vtx2DCuts[1];
1749  float maxSepCutLong2 = tcc.vtx2DCuts[2] * tcc.vtx2DCuts[2];
1750 
1751  // assume that end 0 is closest to the vertex
1752  unsigned short end = 0;
1753  float vtxTjSep2 = PosSep2(vx.Pos, tj.Pts[tj.EndPt[0]].Pos);
1754  float sep1 = PosSep2(vx.Pos, tj.Pts[tj.EndPt[1]].Pos);
1755  if (sep1 < vtxTjSep2) {
1756  // End 1 is closer
1757  end = 1;
1758  vtxTjSep2 = sep1;
1759  }
1760  // is there a vertex already assigned to this end?
1761  if (tj.VtxID[end] > 0) return false;
1762 
1763  // is the trajectory short?
1764  bool tjShort = (tj.EndPt[1] - tj.EndPt[0] < maxShortTjLen);
1765  // use the short Tj cut if the trajectory looks like an electron
1766  if (!tjShort && tj.ChgRMS > 0.5) tjShort = true;
1767  float closestApproach;
1768  // ignore bad separation between the closest tj end and the vertex
1769  if (tjShort) {
1770  if (vtxTjSep2 > maxSepCutShort2) return false;
1771  closestApproach = tcc.vtx2DCuts[1];
1772  }
1773  else {
1774  closestApproach = tcc.vtx2DCuts[2];
1775  if (vtxTjSep2 > maxSepCutLong2) return false;
1776  }
1777 
1778  // Calculate the pull on the vertex
1779  TrajPoint& tp = tj.Pts[tj.EndPt[end]];
1780  float tpVxPull = TrajPointVertexPull(slc, tp, vx);
1781  bool signalBetween = SignalBetween(slc, tp, vx.Pos[0], 0.8);
1782 
1783  // See if the vertex position is close to an end
1784  unsigned short closePt;
1785  TrajClosestApproach(tj, vx.Pos[0], vx.Pos[1], closePt, closestApproach);
1786  // count the number of points between the end of the trajectory and the vertex.
1787  // tj ------------- tj ------------
1788  // vx * >> dpt = 0 vx * >> dpt = 2
1789  short dpt;
1790  if (end == 0) { dpt = closePt - tj.EndPt[end]; }
1791  else {
1792  dpt = tj.EndPt[end] - closePt;
1793  }
1794 
1795  float length = TrajLength(tj);
1796  // don't attach it if the tj length is shorter than the separation distance
1797  if (length > 4 && length < closestApproach) return false;
1798 
1799  float pullCut = tcc.vtx2DCuts[3];
1800  // Dec 21, 2017 Loosen up the pull cut for short close slc. These are likely to
1801  // be poorly reconstructed. It is better to have them associated with the vertex
1802  // than not.
1803  if (tjShort) pullCut *= 2;
1804 
1805  if (prt) {
1806  mf::LogVerbatim myprt("TC");
1807  myprt << "ATTV: 2V" << vx.ID;
1808  myprt << " NTraj " << vx.NTraj;
1809  myprt << " oldTJs";
1810  for (unsigned short itj = 0; itj < slc.tjs.size(); ++itj) {
1811  Trajectory& tj = slc.tjs[itj];
1812  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1813  if (tj.CTP != vx.CTP) continue;
1814  if (tj.VtxID[0] == vx.ID) myprt << " T" << tj.ID << "_0";
1815  if (tj.VtxID[1] == vx.ID) myprt << " T" << tj.ID << "_1";
1816  }
1817  myprt << " + T" << tj.ID << "_" << end << " vtxTjSep " << sqrt(vtxTjSep2) << " tpVxPull "
1818  << tpVxPull << " pullCut " << pullCut << " dpt " << dpt;
1819  }
1820  if (tpVxPull > pullCut) return false;
1821  if (dpt > 2) return true;
1822 
1823  // remove the fixed position flag if there are more than 2 tjs
1824  bool fixedBit = vx.Stat[kFixed];
1825  if (fixedBit && vx.NTraj < 2) vx.Stat[kFixed] = false;
1826 
1827  // Passed all the cuts. Attach it to the vertex and try a fit
1828  tj.VtxID[end] = vx.ID;
1829  // flag as a photon Tj so it isn't included in the fit
1830  tj.AlgMod[kPhoton] = !signalBetween;
1831  // make a copy of the vertex and fit it
1832  auto vxTmp = vx;
1833  if (FitVertex(slc, vxTmp, prt)) {
1834  SetVx2Score(slc, vxTmp);
1835  if (prt) mf::LogVerbatim("TC") << " Success";
1836  vx = vxTmp;
1837  return true;
1838  }
1839  // Keep the Tj -> Vx assn since we got this far, but don't include this end of the Tj in the fit
1840  tj.EndFlag[end][kNoFitVx] = true;
1841  if (prt)
1842  mf::LogVerbatim("TC") << " Poor fit. Keep T" << tj.ID << "-2V" << vx.ID
1843  << " assn with kNoFitVx";
1844  // restore the fixed flag
1845  vx.Stat[kFixed] = fixedBit;
1846  return true;
1847  } // AttachTrajToVertex
1848 
1849  /////////////////////////////////////////
1850  float
1851  TrajPointVertexPull(const TCSlice& slc, const TrajPoint& tp, const VtxStore& vx)
1852  {
1853  // Calculates the position pull between a trajectory point and a vertex
1854 
1855  // impact parameter between them
1856  double ip = PointTrajDOCA(slc, vx.Pos[0], vx.Pos[1], tp);
1857  // separation^2
1858  double sep2 = PosSep2(vx.Pos, tp.Pos);
1859 
1860  // Find the projection of the vertex error ellipse in a coordinate system
1861  // along the TP direction
1862  double vxErrW = vx.PosErr[0] * tp.Dir[1];
1863  double vxErrT = vx.PosErr[1] * tp.Dir[0];
1864  double vxErr2 = vxErrW * vxErrW + vxErrT * vxErrT;
1865  // add the TP position error^2
1866  vxErr2 += tp.HitPosErr2;
1867 
1868  // close together so ignore the TP projection error and return
1869  // the pull using the vertex error and TP position error
1870  if (sep2 < 1) return (float)(ip / sqrt(vxErr2));
1871 
1872  double dang = ip / sqrt(sep2);
1873 
1874  // calculate the angle error.
1875  // Start with the vertex error^2
1876  double angErr = vxErr2 / sep2;
1877  // Add the TP angle error^2
1878  angErr += tp.AngErr * tp.AngErr;
1879  if (angErr == 0) return 999;
1880  angErr = sqrt(angErr);
1881  return (float)(dang / angErr);
1882 
1883  } // TrajPointVertexPull
1884 
1885  /////////////////////////////////////////
1886  float
1887  VertexVertexPull(const TCSlice& slc, const Vtx3Store& vx1, const Vtx3Store& vx2)
1888  {
1889  // Calculates the position pull between two vertices
1890  double dx = vx1.X - vx2.X;
1891  double dy = vx1.Y - vx2.Y;
1892  double dz = vx1.Z - vx2.Z;
1893  double dxErr2 = (vx1.XErr * vx1.XErr + vx2.XErr * vx2.XErr) / 2;
1894  double dyErr2 = (vx1.YErr * vx1.YErr + vx2.YErr * vx2.YErr) / 2;
1895  double dzErr2 = (vx1.ZErr * vx1.ZErr + vx2.ZErr * vx2.ZErr) / 2;
1896  dx = dx * dx / dxErr2;
1897  dy = dy * dy / dyErr2;
1898  dz = dz * dz / dzErr2;
1899  return (float)(sqrt(dx + dy + dz) / 3);
1900  }
1901 
1902  /////////////////////////////////////////
1903  float
1904  VertexVertexPull(const TCSlice& slc, const VtxStore& vx1, const VtxStore& vx2)
1905  {
1906  // Calculates the position pull between two vertices
1907  double dw = vx1.Pos[0] - vx2.Pos[0];
1908  double dt = vx1.Pos[1] - vx2.Pos[1];
1909  double dwErr2 = (vx1.PosErr[0] * vx1.PosErr[0] + vx2.PosErr[0] * vx2.PosErr[0]) / 2;
1910  double dtErr2 = (vx1.PosErr[1] * vx1.PosErr[1] + vx2.PosErr[1] * vx2.PosErr[1]) / 2;
1911  dw = dw * dw / dwErr2;
1912  dt = dt * dt / dtErr2;
1913  return (float)sqrt(dw + dt);
1914  }
1915 
1916  ////////////////////////////////////////////////
1917  bool
1918  StoreVertex(TCSlice& slc, VtxStore& vx)
1919  {
1920  // jacket around the push to ensure that the Tj and vtx CTP is consistent.
1921  // The calling function should score the vertex after the trajectories are attached
1922 
1923  if (vx.ID != int(slc.vtxs.size() + 1)) return false;
1924 
1925  ++evt.global2V_UID;
1926  vx.UID = evt.global2V_UID;
1927 
1928  unsigned short nvxtj = 0;
1929  unsigned short nok = 0;
1930  for (auto& tj : slc.tjs) {
1931  if (tj.AlgMod[kKilled]) continue;
1932  if (vx.ID == tj.VtxID[0] || vx.ID == tj.VtxID[1]) ++nvxtj;
1933  if (vx.CTP != tj.CTP) continue;
1934  if (vx.ID == tj.VtxID[0] || vx.ID == tj.VtxID[1]) ++nok;
1935  } // tj
1936 
1937  if (nok != nvxtj) {
1938  mf::LogVerbatim("TC") << "StoreVertex: vertex " << vx.ID << " Topo " << vx.Topo
1939  << " has inconsistent CTP code " << vx.CTP << " with one or more Tjs\n";
1940  for (auto& tj : slc.tjs) {
1941  if (tj.AlgMod[kKilled]) continue;
1942  if (tj.VtxID[0] == vx.ID) tj.VtxID[0] = 0;
1943  if (tj.VtxID[1] == vx.ID) tj.VtxID[1] = 0;
1944  }
1945  return false;
1946  }
1947  vx.NTraj = nok;
1948  slc.vtxs.push_back(vx);
1949  return true;
1950  } // StoreVertex
1951 
1952  /////////////////////////////////////////
1953  bool
1954  FitVertex(TCSlice& slc, VtxStore& vx, bool prt)
1955  {
1956  // Fit the vertex using T -> 2V assns
1957 
1958  // tcc.vtx2DCuts fcl input usage
1959  // 0 = maximum length of a short trajectory
1960  // 1 = max vertex - trajectory separation for short trajectories
1961  // 2 = max vertex - trajectory separation for long trajectories
1962  // 3 = max position pull for adding TJs to a vertex
1963  // 4 = max allowed vertex position error
1964  // 5 = min MCSMom
1965  // 6 = min Pts/Wire fraction
1966 
1967  if (vx.Stat[kFixed]) {
1968  if (prt) mf::LogVerbatim("TC") << " vertex position fixed. No fit allowed";
1969  return true;
1970  }
1971 
1972  // Create a vector of trajectory points that will be used to fit the vertex position
1973  std::vector<TrajPoint> vxTp;
1974  for (auto& tj : slc.tjs) {
1975  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
1976  if (tj.CTP != vx.CTP) continue;
1977  if (tj.AlgMod[kPhoton]) continue;
1978  bool added = false;
1979  if (tj.VtxID[0] == vx.ID && !tj.EndFlag[0][kNoFitVx]) {
1980  vxTp.push_back(tj.Pts[tj.EndPt[0]]);
1981  added = true;
1982  }
1983  if (tj.VtxID[1] == vx.ID && !tj.EndFlag[1][kNoFitVx]) {
1984  vxTp.push_back(tj.Pts[tj.EndPt[1]]);
1985  added = true;
1986  }
1987  // stash the ID in Step to help debugging
1988  if (added) {
1989  auto& tp = vxTp[vxTp.size() - 1];
1990  if (tj.ID > 0) tp.Step = (int)tj.ID;
1991  // inflate the angle errors for Tjs with few fitted points
1992  if (tp.NTPsFit < 4) tp.AngErr *= 4;
1993  }
1994  } // tj
1995 
1996  bool success = FitVertex(slc, vx, vxTp, prt);
1997 
1998  if (!success) return false;
1999  return true;
2000  } // FitVertex
2001 
2002  /////////////////////////////////////////
2003  bool
2004  FitVertex(TCSlice& slc, VtxStore& vx, std::vector<TrajPoint>& vxTPs, bool prt)
2005  {
2006  // Version with LSQ fit. Each TP position (P0,P1) and slope S are fit to a vertex
2007  // at position (V0, V1), using the equation P1 = V1 + (P0 - V0) * S. This is put
2008  // in the form A * V = b. V is found using V = (A^T * A)^-1 * A^T * b. This is
2009  // usually done using the TDecompSVD Solve method but here we do it step-wise to
2010  // get access to the covariance matrix (A^T * A)^-1. The pull between the TP position
2011  // and the vertex position is stored in tp.Delta
2012 
2013  if (vxTPs.size() < 2) return false;
2014  if (vxTPs.size() == 2) {
2015  vx.ChiDOF = 0.;
2016  return true;
2017  }
2018 
2019  unsigned short npts = vxTPs.size();
2020  TMatrixD A(npts, 2);
2021  TVectorD b(npts);
2022  for (unsigned short itj = 0; itj < vxTPs.size(); ++itj) {
2023  auto& tp = vxTPs[itj];
2024  double dtdw = tp.Dir[1] / tp.Dir[0];
2025  double wt = 1 / (tp.AngErr * tp.AngErr);
2026  A(itj, 0) = -dtdw * wt;
2027  A(itj, 1) = 1. * wt;
2028  b(itj) = (tp.Pos[1] - tp.Pos[0] * dtdw) * wt;
2029  } // itj
2030 
2031  TMatrixD AT(2, npts);
2032  AT.Transpose(A);
2033  TMatrixD ATA = AT * A;
2034  double* det = 0;
2035  try {
2036  ATA.Invert(det);
2037  }
2038  catch (...) {
2039  return false;
2040  }
2041  if (det == NULL) return false;
2042  TVectorD vxPos = ATA * AT * b;
2043  vx.PosErr[0] = sqrt(ATA[0][0]);
2044  vx.PosErr[1] = sqrt(ATA[1][1]);
2045  vx.Pos[0] = vxPos[0];
2046  vx.Pos[1] = vxPos[1];
2047 
2048  // Calculate Chisq
2049  vx.ChiDOF = 0;
2050  if (vxTPs.size() > 2) {
2051  for (auto& tp : vxTPs) {
2052  // highjack TP Delta for the vertex pull
2053  tp.Delta = TrajPointVertexPull(slc, tp, vx);
2054  vx.ChiDOF += tp.Delta;
2055  } // itj
2056  vx.ChiDOF /= (float)(vxTPs.size() - 2);
2057  } // vxTPs.size() > 2
2058 
2059  if (prt) {
2060  mf::LogVerbatim("TC") << "Note: TP - 2V pull is stored in TP.Delta";
2061  PrintTPHeader("FV");
2062  for (auto& tp : vxTPs)
2063  PrintTP("FV", slc, 0, 1, 1, tp);
2064  }
2065 
2066  return true;
2067  } // FitVertex
2068 
2069  //////////////////////////////////////////
2070  bool
2071  ChkVtxAssociations(TCSlice& slc, const CTP_t& inCTP)
2072  {
2073  // Check the associations
2074 
2075  // check the 2D -> 3D associations
2076  geo::PlaneID planeID = DecodeCTP(inCTP);
2077  unsigned short plane = planeID.Plane;
2078  for (auto& vx2 : slc.vtxs) {
2079  if (vx2.CTP != inCTP) continue;
2080  if (vx2.ID == 0) continue;
2081  if (vx2.Vx3ID == 0) continue;
2082  if (vx2.Vx3ID > int(slc.vtx3s.size())) {
2083  mf::LogVerbatim("TC") << "ChkVtxAssociations: Invalid vx2.Vx3ID " << vx2.Vx3ID
2084  << " in 2D vtx " << vx2.ID;
2085  return false;
2086  }
2087  auto& vx3 = slc.vtx3s[vx2.Vx3ID - 1];
2088  if (vx3.ID == 0) {
2089  mf::LogVerbatim("TC") << "ChkVtxAssociations: 2V" << vx2.ID << " thinks it is matched to 3V"
2090  << vx3.ID << " but vx3 is obsolete";
2091  return false;
2092  }
2093  if (vx3.Vx2ID[plane] != vx2.ID) {
2094  mf::LogVerbatim("TC") << "ChkVtxAssociations: 2V" << vx2.ID << " thinks it is matched to 3V"
2095  << vx3.ID << " but vx3 says no!";
2096  return false;
2097  }
2098  } // vx2
2099  // check the 3D -> 2D associations
2100  for (auto& vx3 : slc.vtx3s) {
2101  if (vx3.ID == 0) continue;
2102  if (vx3.Vx2ID[plane] == 0) continue;
2103  if (vx3.Vx2ID[plane] > (int)slc.vtxs.size()) {
2104  mf::LogVerbatim("TC") << "ChkVtxAssociations: Invalid vx3.Vx2ID " << vx3.Vx2ID[plane]
2105  << " in CTP " << inCTP;
2106  return false;
2107  }
2108  auto& vx2 = slc.vtxs[vx3.Vx2ID[plane] - 1];
2109  if (vx2.Vx3ID != vx3.ID) {
2110  mf::LogVerbatim("TC") << "ChkVtxAssociations: 3V" << vx3.ID << " thinks it is matched to 2V"
2111  << vx2.ID << " but vx2 says no!";
2112  return false;
2113  }
2114  } // vx3
2115 
2116  // check the Tj -> 2D associations
2117  for (auto& tj : slc.tjs) {
2118  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2119  for (unsigned short end = 0; end < 2; ++end) {
2120  if (tj.VtxID[end] == 0) continue;
2121  if (tj.VtxID[end] > slc.vtxs.size()) {
2122  mf::LogVerbatim("TC") << "ChkVtxAssociations: T" << tj.ID << " thinks it is matched to 2V"
2123  << tj.VtxID[end] << " on end " << end
2124  << " but no vertex exists. Recovering";
2125  tj.VtxID[end] = 0;
2126  return false;
2127  }
2128  unsigned short ivx = tj.VtxID[end] - 1;
2129  auto& vx2 = slc.vtxs[ivx];
2130  if (vx2.ID == 0) {
2131  mf::LogVerbatim("TC") << "ChkVtxAssociations: T" << tj.ID << " thinks it is matched to 2V"
2132  << tj.VtxID[end] << " on end " << end
2133  << " but the vertex is killed. Fixing the Tj";
2134  tj.VtxID[end] = 0;
2135  return false;
2136  }
2137  } // end
2138  } // tj
2139 
2140  return true;
2141 
2142  } // ChkVtxAssociations
2143 
2144  //////////////////////////////////////////
2145  void
2146  ScoreVertices(TCSlice& slc)
2147  {
2148  // reset all 3D vertex, 2D vertex and Tj high-score vertex bits in tpcid
2149 
2150  // reset the 2D vertex status bits
2151  for (auto& vx : slc.vtxs) {
2152  if (vx.ID == 0) continue;
2153  vx.Stat[kHiVx3Score] = false;
2154  } // vx
2155  // and the tj bits
2156  for (auto& tj : slc.tjs) {
2157  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2158  tj.AlgMod[kTjHiVx3Score] = false;
2159  } // tj
2160  // Score the 2D vertices
2161  for (auto& vx : slc.vtxs) {
2162  if (vx.ID == 0) continue;
2163  SetVx2Score(slc, vx);
2164  } // vx
2165  // Score the 3D vertices
2166  for (auto& vx3 : slc.vtx3s) {
2167  if (vx3.ID == 0) continue;
2168  SetVx3Score(slc, vx3);
2169  } // vx3
2170  } // ScoreVertices
2171 
2172  //////////////////////////////////////////
2173  void
2174  KillPoorVertices(TCSlice& slc)
2175  {
2176  // kill 2D vertices that have low score and are not attached to a high-score 3D vertex
2177  if (slc.vtxs.empty()) return;
2178  for (auto& vx : slc.vtxs) {
2179  if (vx.ID == 0) continue;
2180  if (vx.Score > tcc.vtx2DCuts[7]) continue;
2181  if (vx.Vx3ID > 0) {
2182  auto& vx3 = slc.vtx3s[vx.Vx3ID - 1];
2183  if (vx3.Primary) continue;
2184  if (slc.vtx3s[vx.Vx3ID - 1].Score >= tcc.vtx2DCuts[7]) continue;
2185  }
2186  MakeVertexObsolete("KPV", slc, vx, false);
2187  } // vx
2188 
2189  } // KillPoorVertices
2190 
2191  //////////////////////////////////////////
2192  void
2193  SetHighScoreBits(TCSlice& slc, Vtx3Store& vx3)
2194  {
2195  // Sets the tj and 2D vertex score bits to true
2196 
2197  if (vx3.ID == 0) return;
2198 
2199  for (unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
2200  if (vx3.Vx2ID[ipl] <= 0) continue;
2201  VtxStore& vx2 = slc.vtxs[vx3.Vx2ID[ipl] - 1];
2202  vx2.Stat[kHiVx3Score] = false;
2203  // transfer this to all attached tjs and vertices attached to those tjs
2204  std::vector<int> tjlist = GetVtxTjIDs(slc, vx2);
2205  std::vector<int> vxlist;
2206  while (true) {
2207  // tag Tjs and make a list of attached vertices whose high-score
2208  // bit needs to be set
2209  vxlist.clear();
2210  for (auto tjid : tjlist) {
2211  auto& tj = slc.tjs[tjid - 1];
2212  tj.AlgMod[kTjHiVx3Score] = true;
2213  for (unsigned short end = 0; end < 2; ++end) {
2214  if (tj.VtxID[end] == 0) continue;
2215  auto& vx2 = slc.vtxs[tj.VtxID[end] - 1];
2216  if (vx2.Stat[kHiVx3Score]) continue;
2217  vx2.Stat[kHiVx3Score] = true;
2218  vxlist.push_back(vx2.ID);
2219  } // end
2220  } // tjid
2221 
2222  if (vxlist.empty()) break;
2223  // re-build tjlist using vxlist
2224  std::vector<int> newtjlist;
2225  for (auto vxid : vxlist) {
2226  auto& vx2 = slc.vtxs[vxid - 1];
2227  auto tmp = GetVtxTjIDs(slc, vx2);
2228  for (auto tjid : tmp) {
2229  if (std::find(tjlist.begin(), tjlist.end(), tjid) == tjlist.end())
2230  newtjlist.push_back(tjid);
2231  } // tjid
2232  } // vxid
2233  if (newtjlist.empty()) break;
2234  tjlist = newtjlist;
2235  } // true
2236  } // ipl
2237 
2238  } // SetHighScoreBits
2239 
2240  //////////////////////////////////////////
2241  void
2242  SetVx3Score(TCSlice& slc, Vtx3Store& vx3)
2243  {
2244  // Calculate the 3D vertex score and flag Tjs that are attached to high score vertices as defined
2245  // by vtx2DCuts
2246 
2247  if (vx3.ID == 0) return;
2248 
2249  vx3.Score = 0;
2250  for (unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
2251  if (vx3.Vx2ID[ipl] <= 0) continue;
2252  VtxStore& vx2 = slc.vtxs[vx3.Vx2ID[ipl] - 1];
2253  vx3.Score += vx2.Score;
2254  } // ipl
2255  vx3.Score /= (float)slc.nPlanes;
2256  // don't allow it to get too small or negative
2257  if (vx3.Score < 0.001) vx3.Score = 0.001;
2258  if (vx3.Score > tcc.vtx2DCuts[7]) SetHighScoreBits(slc, vx3);
2259 
2260  } // SetVx3Score
2261 
2262  //////////////////////////////////////////
2263  void
2264  SetVx2Score(TCSlice& slc)
2265  {
2266  // A version that sets the score of the last added vertex
2267  if (slc.vtxs.empty()) return;
2268  auto& vx2 = slc.vtxs[slc.vtxs.size() - 1];
2269  SetVx2Score(slc, vx2);
2270  } // SetVx2Score
2271 
2272  //////////////////////////////////////////
2273  void
2274  SetVx2Score(TCSlice& slc, VtxStore& vx2)
2275  {
2276  // Calculate the 2D vertex score
2277  if (vx2.ID == 0) return;
2278 
2279  // Don't score vertices from CheckTrajBeginChg, MakeJunkVertices or Neutral vertices. Set to the minimum
2280  if (vx2.Topo == 8 || vx2.Topo == 9 || vx2.Topo == 11 || vx2.Topo == 12) {
2281  vx2.Score = tcc.vtx2DCuts[7] + 0.1;
2282  auto vtxTjID = GetVtxTjIDs(slc, vx2);
2283  vx2.TjChgFrac = ChgFracNearPos(slc, vx2.Pos, vtxTjID);
2284  return;
2285  }
2286 
2287  // Cuts on Tjs attached to vertices
2288  constexpr float maxChgRMS = 0.25;
2289  constexpr float momBin = 50;
2290 
2291  vx2.Score = -1000;
2292  vx2.TjChgFrac = 0;
2293  if (vx2.ID == 0) return;
2294  if (tcc.vtxScoreWeights.size() < 4) return;
2295 
2296  auto vtxTjIDs = GetVtxTjIDs(slc, vx2);
2297  if (vtxTjIDs.empty()) return;
2298 
2299  // Vertex position error
2300  float vpeScore = -tcc.vtxScoreWeights[0] * (vx2.PosErr[0] + vx2.PosErr[1]);
2301 
2302  unsigned short m3Dcnt = 0;
2303  if (vx2.Vx3ID > 0) {
2304  m3Dcnt = 1;
2305  // Add another if the 3D vertex is complete
2306  unsigned short ivx3 = vx2.Vx3ID - 1;
2307  if (slc.vtx3s[ivx3].Wire < 0) m3Dcnt = 2;
2308  }
2309  float m3DScore = tcc.vtxScoreWeights[1] * m3Dcnt;
2310 
2311  vx2.TjChgFrac = ChgFracNearPos(slc, vx2.Pos, vtxTjIDs);
2312  float cfScore = tcc.vtxScoreWeights[2] * vx2.TjChgFrac;
2313 
2314  // Define a weight for each Tj
2315  std::vector<int> tjids;
2316  std::vector<float> tjwts;
2317  unsigned short cnt13 = 0;
2318  for (auto tjid : vtxTjIDs) {
2319  Trajectory& tj = slc.tjs[tjid - 1];
2320  // Feb 22 Ignore short Tjs and junk tjs
2321  if (tj.AlgMod[kJunkTj]) continue;
2322  unsigned short lenth = tj.EndPt[1] - tj.EndPt[0] + 1;
2323  if (lenth < 3) continue;
2324  float wght = (float)tj.MCSMom / momBin;
2325  // weight by the first tagged muon
2326  if (tj.PDGCode == 13) {
2327  ++cnt13;
2328  if (cnt13 == 1) wght *= 2;
2329  }
2330  // weight by charge rms
2331  if (tj.ChgRMS < maxChgRMS) ++wght;
2332  // Shower Tj
2333  if (tj.AlgMod[kShowerTj]) ++wght;
2334  // ShowerLike
2335  if (tj.AlgMod[kShowerLike]) --wght;
2336  tjids.push_back(tjid);
2337  tjwts.push_back(wght);
2338  } // tjid
2339 
2340  if (tjids.empty()) return;
2341 
2342  float tjScore = 0;
2343  float sum = 0;
2344  float cnt = 0;
2345  for (unsigned short it1 = 0; it1 < tjids.size() - 1; ++it1) {
2346  Trajectory& tj1 = slc.tjs[tjids[it1] - 1];
2347  float wght1 = tjwts[it1];
2348  // the end that has a vertex
2349  unsigned short end1 = 0;
2350  if (tj1.VtxID[1] == vx2.ID) end1 = 1;
2351  unsigned short endPt1 = tj1.EndPt[end1];
2352  // bump up the weight if there is a Bragg peak at the other end
2353  unsigned short oend1 = 1 - end1;
2354  if (tj1.EndFlag[oend1][kBragg]) ++wght1;
2355  float ang1 = tj1.Pts[endPt1].Ang;
2356  float ang1Err2 = tj1.Pts[endPt1].AngErr * tj1.Pts[endPt1].AngErr;
2357  for (unsigned short it2 = it1 + 1; it2 < tjids.size(); ++it2) {
2358  Trajectory& tj2 = slc.tjs[tjids[it2] - 1];
2359  float wght2 = tjwts[it2];
2360  unsigned end2 = 0;
2361  if (tj2.VtxID[1] == vx2.ID) end2 = 1;
2362  // bump up the weight if there is a Bragg peak at the other end
2363  unsigned short oend2 = 1 - end2;
2364  if (tj2.EndFlag[oend2][kBragg]) ++wght2;
2365  unsigned short endPt2 = tj2.EndPt[end2];
2366  float ang2 = tj2.Pts[endPt2].Ang;
2367  float ang2Err2 = tj2.Pts[endPt2].AngErr * tj2.Pts[endPt2].AngErr;
2368  float dang = DeltaAngle(ang1, ang2);
2369  float dangErr = 0.5 * sqrt(ang1Err2 + ang2Err2);
2370  if ((dang / dangErr) > 3 && wght1 > 0 && wght2 > 0) {
2371  sum += wght1 + wght2;
2372  ++cnt;
2373  }
2374  } // it2
2375  } // it1
2376  if (cnt > 0) {
2377  sum /= cnt;
2378  tjScore = tcc.vtxScoreWeights[3] * sum;
2379  }
2380  vx2.Score = vpeScore + m3DScore + cfScore + tjScore;
2381  if (tcc.dbg2V && tcc.dbgSlc && vx2.CTP == debug.CTP) {
2382  // last call after vertices have been matched to the truth. Use to optimize vtxScoreWeights using
2383  // an ntuple
2384  mf::LogVerbatim myprt("TC");
2385  bool printHeader = true;
2386  Print2V("SVx2S", myprt, vx2, printHeader);
2387  myprt << std::fixed << std::setprecision(1);
2388  myprt << " vpeScore " << vpeScore << " m3DScore " << m3DScore;
2389  myprt << " cfScore " << cfScore << " tjScore " << tjScore;
2390  myprt << " Score " << vx2.Score;
2391  }
2392  } // SetVx2Score
2393 
2394  //////////////////////////////////////////
2395  void
2396  CompleteIncomplete3DVerticesInGaps(detinfo::DetectorPropertiesData const& detProp, TCSlice& slc)
2397  {
2398 
2399  if (!tcc.useAlg[kComp3DVxIG]) return;
2400  if (slc.nPlanes != 3) return;
2401 
2402  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kComp3DVxIG]);
2403  if (prt) mf::LogVerbatim("TC") << "Inside CI3DVIG:";
2404 
2405  for (unsigned short iv3 = 0; iv3 < slc.vtx3s.size(); ++iv3) {
2406  Vtx3Store& vx3 = slc.vtx3s[iv3];
2407  // ignore obsolete vertices
2408  if (vx3.ID == 0) continue;
2409  // check for a completed 3D vertex
2410  if (vx3.Wire < 0) continue;
2411  unsigned short mPlane = USHRT_MAX;
2412  for (unsigned short ipl = 0; ipl < slc.nPlanes; ++ipl) {
2413  if (vx3.Vx2ID[ipl] > 0) continue;
2414  mPlane = ipl;
2415  break;
2416  } // ipl
2417  if (mPlane == USHRT_MAX) continue;
2418  CTP_t mCTP = EncodeCTP(vx3.TPCID.Cryostat, vx3.TPCID.TPC, mPlane);
2419  // require that the missing vertex be in a large block of dead wires
2420  float dwc = DeadWireCount(slc, vx3.Wire - 4, vx3.Wire + 4, mCTP);
2421  if (dwc < 5) continue;
2422  // X position of the purported missing vertex
2423  VtxStore aVtx;
2424  aVtx.ID = slc.vtxs.size() + 1;
2425  aVtx.Pos[0] = vx3.Wire;
2426  aVtx.Pos[1] = detProp.ConvertXToTicks(vx3.X, mPlane, vx3.TPCID.TPC, vx3.TPCID.Cryostat) *
2427  tcc.unitsPerTick;
2428  aVtx.CTP = mCTP;
2429  aVtx.Topo = 4;
2430  aVtx.NTraj = 0;
2431  // Give it a bogus pass to indicate it wasn't created while stepping
2432  aVtx.Pass = 9;
2433  if (prt)
2434  mf::LogVerbatim("TC") << "CI3DVIG: Incomplete vertex " << iv3 << " in plane " << mPlane
2435  << " wire " << vx3.Wire << " Made 2D vertex ";
2436  std::vector<int> tjIDs;
2437  std::vector<unsigned short> tjEnds;
2438  for (unsigned short itj = 0; itj < slc.tjs.size(); ++itj) {
2439  if (slc.tjs[itj].CTP != mCTP) continue;
2440  if (slc.tjs[itj].AlgMod[kKilled] || slc.tjs[itj].AlgMod[kHaloTj]) continue;
2441  for (unsigned short end = 0; end < 2; ++end) {
2442  unsigned short ept = slc.tjs[itj].EndPt[end];
2443  TrajPoint& tp = slc.tjs[itj].Pts[ept];
2444  unsigned short oept = slc.tjs[itj].EndPt[1 - end];
2445  TrajPoint& otp = slc.tjs[itj].Pts[oept];
2446  // ensure that this is the end closest to the vertex
2447  if (std::abs(tp.Pos[0] - aVtx.Pos[0]) > std::abs(otp.Pos[0] - aVtx.Pos[0])) continue;
2448  float doca = PointTrajDOCA(slc, aVtx.Pos[0], aVtx.Pos[1], tp);
2449  if (doca > 2) continue;
2450  float dwc = DeadWireCount(slc, aVtx.Pos[0], tp.Pos[0], tp.CTP);
2451  float ptSep;
2452  if (aVtx.Pos[0] > tp.Pos[0]) { ptSep = aVtx.Pos[0] - tp.Pos[0] - dwc; }
2453  else {
2454  ptSep = tp.Pos[0] - aVtx.Pos[0] - dwc;
2455  }
2456  if (prt)
2457  mf::LogVerbatim("TC") << "CI3DVIG: tj ID " << slc.tjs[itj].ID << " doca " << doca
2458  << " ptSep " << ptSep;
2459  if (ptSep < -2 || ptSep > 2) continue;
2460  // don't clobber an existing association
2461  if (slc.tjs[itj].VtxID[end] > 0) continue;
2462  tjIDs.push_back(slc.tjs[itj].ID);
2463  tjEnds.push_back(end);
2464  } // end
2465  } // itj
2466  if (!tjIDs.empty()) {
2467  // Determine how messy this region is
2468  aVtx.TjChgFrac = ChgFracNearPos(slc, aVtx.Pos, tjIDs);
2469  if (aVtx.TjChgFrac < 0.7) continue;
2470  aVtx.Vx3ID = vx3.ID;
2471  // Save the 2D vertex
2472  if (!StoreVertex(slc, aVtx)) continue;
2473  for (unsigned short ii = 0; ii < tjIDs.size(); ++ii) {
2474  unsigned short itj = tjIDs[ii] - 1;
2475  slc.tjs[itj].VtxID[tjEnds[ii]] = aVtx.ID;
2476  slc.tjs[itj].AlgMod[kComp3DVxIG] = true;
2477  } // ii
2478  SetVx2Score(slc);
2479  vx3.Vx2ID[mPlane] = aVtx.ID;
2480  vx3.Wire = -1;
2481  if (prt)
2482  mf::LogVerbatim("TC") << "CI3DVIG: new vtx 2V" << aVtx.ID << " points to 3V" << vx3.ID;
2483  }
2484  } // vx3
2485 
2486  } // CompleteIncomplete3DVerticesInGaps
2487 
2488  //////////////////////////////////////////
2489  void
2490  CompleteIncomplete3DVertices(detinfo::DetectorPropertiesData const& detProp, TCSlice& slc)
2491  {
2492  // Look for trajectories in a plane that lack a 2D vertex as listed in
2493  // 2DVtxID that are near the projected wire. This may trigger splitting trajectories,
2494  // assigning them to a new 2D vertex and completing 3D vertices
2495 
2496  if (!tcc.useAlg[kComp3DVx]) return;
2497  if (slc.nPlanes != 3) return;
2498 
2499  bool prt = (tcc.modes[kDebug] && tcc.dbgSlc && tcc.dbgAlg[kComp3DVx]);
2500 
2501  float maxdoca = 3;
2502  if (prt) mf::LogVerbatim("TC") << "Inside CI3DV with maxdoca set to " << maxdoca;
2503  unsigned short ivx3 = 0;
2504  for (auto& vx3 : slc.vtx3s) {
2505  // ignore obsolete vertices
2506  if (vx3.ID == 0) continue;
2507  // check for a completed 3D vertex
2508  if (vx3.Wire < 0) continue;
2509  unsigned short mPlane = USHRT_MAX;
2510  // look for vertices in the induction plane in which the charge requirement wasn't imposed
2511  bool indPlnNoChgVtx = false;
2512  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
2513  if (vx3.Vx2ID[plane] > 0) {
2514  auto& vx2 = slc.vtxs[vx3.Vx2ID[plane] - 1];
2515  if (vx2.Stat[kVxIndPlnNoChg]) indPlnNoChgVtx = true;
2516  continue;
2517  }
2518  mPlane = plane;
2519  } // ipl
2520  if (mPlane == USHRT_MAX) continue;
2521  if (indPlnNoChgVtx) continue;
2522  CTP_t mCTP = EncodeCTP(vx3.TPCID.Cryostat, vx3.TPCID.TPC, mPlane);
2523  // X position of the purported missing vertex
2524  // A TP for the missing 2D vertex
2525  TrajPoint vtp;
2526  vtp.Pos[0] = vx3.Wire;
2527  vtp.Pos[1] = detProp.ConvertXToTicks(vx3.X, mPlane, vx3.TPCID.TPC, vx3.TPCID.Cryostat) *
2528  tcc.unitsPerTick;
2529  if (prt)
2530  mf::LogVerbatim("TC") << "CI3DV 3V" << vx3.ID << " Pos " << mPlane << ":"
2531  << PrintPos(slc, vtp.Pos);
2532  std::vector<int> tjIDs;
2533  std::vector<unsigned short> tjPts;
2534  for (auto& tj : slc.tjs) {
2535  if (tj.CTP != mCTP) continue;
2536  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2537  if (tj.Pts.size() < 6) continue;
2538  if (tj.AlgMod[kComp3DVx]) continue;
2539  float doca = maxdoca;
2540  // find the closest distance between the vertex and the trajectory
2541  unsigned short closePt = 0;
2542  TrajPointTrajDOCA(slc, vtp, tj, closePt, doca);
2543  if (closePt > tj.EndPt[1]) continue;
2544  // try to improve the location of the vertex by looking for a distinctive feature on the
2545  // trajectory, e.g. high multiplicity hits or larger than normal charge
2546  if (RefineVtxPosition(slc, tj, closePt, 3, false)) vtp.Pos = tj.Pts[closePt].Pos;
2547  if (prt)
2548  mf::LogVerbatim("TC") << "CI3DV 3V" << vx3.ID << " candidate T" << tj.ID << " vtx pos "
2549  << PrintPos(slc, vtp.Pos) << " doca " << doca << " closePt "
2550  << closePt;
2551  tjIDs.push_back(tj.ID);
2552  tjPts.push_back(closePt);
2553  } // itj
2554  if (tjIDs.empty()) continue;
2555  // compare the length of the Tjs used to make the vertex with the length of the
2556  // Tj that we want to split. Don't allow a vertex using very short Tjs to split a long
2557  // Tj in the 3rd plane
2558  auto vxtjs = GetAssns(slc, "3V", vx3.ID, "T");
2559  unsigned short maxPts = 0;
2560  unsigned short minPts = USHRT_MAX;
2561  for (auto tjid : vxtjs) {
2562  auto& tj = slc.tjs[tjid - 1];
2563  unsigned short npwc = NumPtsWithCharge(slc, tj, false);
2564  if (npwc > maxPts) maxPts = npwc;
2565  if (npwc < minPts) minPts = npwc;
2566  } // tjid
2567  // skip this operation if any of the Tjs in the split list are > 3 * maxPts
2568  maxPts *= 3;
2569  bool skipit = false;
2570  for (auto tjid : tjIDs) {
2571  auto& tj = slc.tjs[tjid - 1];
2572  if (NumPtsWithCharge(slc, tj, false) > maxPts) skipit = true;
2573  } // tjid
2574  if (prt)
2575  mf::LogVerbatim("TC") << " maxPts " << maxPts << " skipit? " << skipit << " minPts "
2576  << minPts;
2577  if (skipit) continue;
2578  // 2D vertex
2579  VtxStore aVtx;
2580  unsigned short newVtxIndx = slc.vtxs.size();
2581  aVtx.ID = newVtxIndx + 1;
2582  aVtx.CTP = mCTP;
2583  aVtx.Topo = 3;
2584  aVtx.NTraj = 0;
2585  // Give it a bogus pass to indicate it wasn't created while stepping
2586  aVtx.Pass = 9;
2587  aVtx.Pos = vtp.Pos;
2588  // ensure this isn't in a messy region
2589  aVtx.TjChgFrac = ChgFracNearPos(slc, aVtx.Pos, tjIDs);
2590  if (prt)
2591  mf::LogVerbatim("TC") << " charge fraction near position " << aVtx.TjChgFrac
2592  << " cut if < 0.6";
2593  if (aVtx.TjChgFrac < 0.6) continue;
2594  if (!StoreVertex(slc, aVtx)) continue;
2595  // make a reference to the new vertex
2596  VtxStore& newVtx = slc.vtxs[slc.vtxs.size() - 1];
2597  if (prt) mf::LogVerbatim("TC") << " Stored new 2V" << newVtx.ID;
2598  // make a temporary copy so we can nudge it a bit if there is only one Tj
2599  std::array<float, 2> vpos = aVtx.Pos;
2600  for (unsigned short ii = 0; ii < tjIDs.size(); ++ii) {
2601  unsigned short itj = tjIDs[ii] - 1;
2602  unsigned short closePt = tjPts[ii];
2603  // determine which end is the closest
2604  unsigned short end = 1;
2605  // closest to the beginning?
2606  if (fabs(closePt - slc.tjs[itj].EndPt[0]) < fabs(closePt - slc.tjs[itj].EndPt[1])) end = 0;
2607  short dpt = fabs(closePt - slc.tjs[itj].EndPt[end]);
2608  if (prt) mf::LogVerbatim("TC") << " dpt " << dpt << " to end " << end;
2609  if (dpt < 2) {
2610  // close to an end
2611  if (slc.tjs[itj].VtxID[end] > 0) {
2612  // find the distance btw the existing vertex and the end of this tj
2613  auto& oldTj = slc.tjs[itj];
2614  auto& oldVx = slc.vtxs[oldTj.VtxID[end] - 1];
2615  short oldSep = fabs(oldVx.Pos[0] - oldTj.Pts[oldTj.EndPt[end]].Pos[0]);
2616  if (prt)
2617  mf::LogVerbatim("TC")
2618  << " T" << slc.tjs[itj].ID << " has vertex 2V" << slc.tjs[itj].VtxID[end]
2619  << " at end " << end << ". oldSep " << oldSep;
2620  if (dpt < oldSep) { MakeVertexObsolete("CI3DV", slc, oldVx, true); }
2621  else {
2622  continue;
2623  }
2624  } // slc.tjs[itj].VtxID[end] > 0
2625  slc.tjs[itj].VtxID[end] = slc.vtxs[newVtxIndx].ID;
2626  ++newVtx.NTraj;
2627  if (prt)
2628  mf::LogVerbatim("TC") << " attach Traj T" << slc.tjs[itj].ID << " at end " << end;
2629  slc.tjs[itj].AlgMod[kComp3DVx] = true;
2630  vpos = slc.tjs[itj].Pts[slc.tjs[itj].EndPt[end]].Pos;
2631  }
2632  else {
2633  // closePt is not near an end, so split the trajectory
2634  if (SplitTraj(slc, itj, closePt, newVtxIndx, prt)) {
2635  if (prt)
2636  mf::LogVerbatim("TC")
2637  << " SplitTraj success 2V" << slc.vtxs[newVtxIndx].ID << " at closePt " << closePt;
2638  // successfully split the Tj
2639  newVtx.NTraj += 2;
2640  }
2641  else {
2642  // split failed. Give up
2643  if (prt) mf::LogVerbatim("TC") << " SplitTraj failed";
2644  newVtx.NTraj = 0;
2645  break;
2646  }
2647  // Update the PDGCode for the chopped trajectory
2648  SetPDGCode(slc, itj);
2649  // and for the new trajectory
2650  SetPDGCode(slc, slc.tjs.size() - 1);
2651  } // closePt is not near an end, so split the trajectory
2652  slc.tjs[itj].AlgMod[kComp3DVx] = true;
2653  unsigned short newtj = slc.tjs.size() - 1;
2654  slc.tjs[newtj].AlgMod[kComp3DVx] = true;
2655  } // ii
2656  if (newVtx.NTraj == 0) {
2657  // A failure occurred. Recover
2658  if (prt) mf::LogVerbatim("TC") << " Failed. Recover and delete vertex " << newVtx.ID;
2659  MakeVertexObsolete("CI3DV", slc, newVtx, true);
2660  }
2661  else {
2662  // success
2663  vx3.Vx2ID[mPlane] = newVtx.ID;
2664  newVtx.Vx3ID = vx3.ID;
2665  vx3.Wire = -1;
2666  // set the vertex position to the start of the Tj if there is only one and fix it
2667  if (newVtx.NTraj == 1) {
2668  newVtx.Pos = vpos;
2669  newVtx.Stat[kFixed] = true;
2670  }
2671  AttachAnyTrajToVertex(slc, newVtx.ID - 1, prt);
2672  SetVx2Score(slc);
2673  if (prt) {
2674  mf::LogVerbatim myprt("TC");
2675  myprt << " Success: new 2V" << newVtx.ID << " at " << (int)newVtx.Pos[0] << ":"
2676  << (int)newVtx.Pos[1] / tcc.unitsPerTick;
2677  myprt << " points to 3V" << vx3.ID;
2678  myprt << " TjIDs:";
2679  for (auto& tjID : tjIDs)
2680  myprt << " T" << std::to_string(tjID);
2681  } // prt
2682  } // success
2683  ++ivx3;
2684  } // vx3
2685 
2686  } // CompleteIncomplete3DVertices
2687 
2688  ////////////////////////////////////////////////
2689  bool
2690  RefineVtxPosition(TCSlice& slc,
2691  const Trajectory& tj,
2692  unsigned short& nearPt,
2693  short nPtsToChk,
2694  bool prt)
2695  {
2696  // The tj has been slated to be split somewhere near point nearPt. This function will move
2697  // the near point a bit to the most likely point of a vertex
2698 
2699  float maxChg = tj.Pts[nearPt].Chg;
2700  short maxChgPt = nearPt;
2701  unsigned short fromPt = tj.EndPt[0];
2702  short spt = (short)nearPt - (short)nPtsToChk;
2703  if (spt > (short)fromPt) fromPt = nearPt - nPtsToChk;
2704  unsigned short toPt = nearPt + nPtsToChk;
2705  if (toPt > tj.EndPt[1]) toPt = tj.EndPt[1];
2706 
2707  for (short ipt = fromPt; ipt <= toPt; ++ipt) {
2708  if (ipt < tj.EndPt[0] || ipt > tj.EndPt[1]) continue;
2709  auto& tp = tj.Pts[ipt];
2710  if (tp.Chg > maxChg) {
2711  maxChg = tp.Chg;
2712  maxChgPt = ipt;
2713  }
2714  if (prt)
2715  mf::LogVerbatim("TC") << "RVP: ipt " << ipt << " Pos " << tp.CTP << ":"
2716  << PrintPos(slc, tp.Pos) << " chg " << (int)tp.Chg << " nhits "
2717  << tp.Hits.size();
2718  } // ipt
2719  if (nearPt == maxChgPt) return false;
2720  nearPt = maxChgPt;
2721  return true;
2722  } //RefineVtxPosition
2723 
2724  ////////////////////////////////////////////////
2725  bool
2726  MakeVertexObsolete(std::string fcnLabel, TCSlice& slc, VtxStore& vx2, bool forceKill)
2727  {
2728  // Makes a 2D vertex obsolete
2729 
2730  // check for a high-score 3D vertex
2731  bool hasHighScoreVx3 = (vx2.Vx3ID > 0);
2732  if (hasHighScoreVx3 && !forceKill && slc.vtx3s[vx2.Vx3ID - 1].Score >= tcc.vtx2DCuts[7])
2733  return false;
2734 
2735  if (tcc.dbg2V || tcc.dbg3V) {
2736  mf::LogVerbatim("TC") << fcnLabel << " MVO: killing 2V" << vx2.ID;
2737  }
2738 
2739  // Kill it
2740  int vx2id = vx2.ID;
2741  if (vx2.Vx3ID > 0) {
2742  auto& vx3 = slc.vtx3s[vx2.Vx3ID - 1];
2743  for (auto& v3v2id : vx3.Vx2ID)
2744  if (v3v2id == vx2.ID) v3v2id = 0;
2745  }
2746  vx2.ID = 0;
2747  for (auto& tj : slc.tjs) {
2748  if (tj.AlgMod[kKilled] || tj.AlgMod[kHaloTj]) continue;
2749  for (unsigned short end = 0; end < 2; ++end) {
2750  if (tj.VtxID[end] != vx2id) continue;
2751  tj.VtxID[end] = 0;
2752  tj.AlgMod[kPhoton] = false;
2753  // clear the kEnvOverlap bits on the TPs
2754  for (unsigned short ii = 0; ii < tj.Pts.size(); ++ii) {
2755  if (end == 0) {
2756  unsigned short ipt = tj.EndPt[0] + ii;
2757  auto& tp = tj.Pts[ipt];
2758  if (!tp.Environment[kEnvOverlap]) break;
2759  if (ipt == tj.EndPt[1]) break;
2760  }
2761  else {
2762  unsigned short ipt = tj.EndPt[1] - ii;
2763  auto& tp = tj.Pts[ipt];
2764  if (!tp.Environment[kEnvOverlap]) break;
2765  if (ipt == tj.EndPt[0]) break;
2766  }
2767  } // ii
2768  if (tj.AlgMod[kTjHiVx3Score]) {
2769  // see if the vertex at the other end is high-score and if so, preserve the state
2770  unsigned short oend = 1 - end;
2771  if (tj.VtxID[oend] > 0) {
2772  auto& ovx2 = slc.vtxs[tj.VtxID[oend] - 1];
2773  if (!ovx2.Stat[kHiVx3Score]) tj.AlgMod[kTjHiVx3Score] = false;
2774  } // vertex at the other end
2775  } // tj.AlgMod[kTjHiVx3Score]
2776  } // end
2777  } // tj
2778 
2779  if (!hasHighScoreVx3) return true;
2780 
2781  // update the affected 3D vertex
2782  Vtx3Store& vx3 = slc.vtx3s[vx2.Vx3ID - 1];
2783  // make the 3D vertex incomplete
2784  geo::PlaneID planeID = DecodeCTP(vx2.CTP);
2785  unsigned short plane = planeID.Plane;
2786  if (vx3.Vx2ID[plane] != vx2id) return true;
2787  vx3.Vx2ID[plane] = 0;
2788  vx3.Wire = vx2.Pos[0];
2789  // Ensure that there are at least two 2D vertices left
2790  unsigned short n2D = 0;
2791  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane)
2792  if (vx3.Vx2ID[plane] > 0) ++n2D;
2793 
2794  if (n2D > 1) {
2795  // 3D vertex is incomplete
2796  // correct the score
2797  SetVx3Score(slc, vx3);
2798  return true;
2799  }
2800 
2801  // 3D vertex is obsolete
2802  // Detach all remaining 2D vertices from the 3D vertex
2803  for (auto& vx2 : slc.vtxs) {
2804  if (vx2.ID == 0) continue;
2805  if (vx2.Vx3ID == vx3.ID) vx2.Vx3ID = 0;
2806  } // vx2
2807  for (auto& pfp : slc.pfps) {
2808  for (unsigned short end = 0; end < 2; ++end)
2809  if (pfp.Vx3ID[end] == vx3.ID) pfp.Vx3ID[end] = 0;
2810  } // pfp
2811  vx3.ID = 0;
2812  return true;
2813 
2814  } // MakeVertexObsolete
2815 
2816  ////////////////////////////////////////////////
2817  bool
2818  MakeVertexObsolete(TCSlice& slc, Vtx3Store& vx3)
2819  {
2820  // Deletes a 3D vertex and 2D vertices in all planes
2821  // The 2D and 3D vertices are NOT killed if forceKill is false and the 3D vertex
2822  // has a high score
2823  if (vx3.ID <= 0) return true;
2824  if (vx3.ID > int(slc.vtx3s.size())) return false;
2825 
2826  for (auto vx2id : vx3.Vx2ID) {
2827  if (vx2id == 0 || vx2id > (int)slc.vtxs.size()) continue;
2828  auto& vx2 = slc.vtxs[vx2id - 1];
2829  MakeVertexObsolete("MVO3", slc, vx2, true);
2830  }
2831  vx3.ID = 0;
2832  return true;
2833  } // MakeVertexObsolete
2834 
2835  //////////////////////////////////////////
2836  std::vector<int>
2837  GetVtxTjIDs(const TCSlice& slc, const VtxStore& vx2)
2838  {
2839  // returns a list of trajectory IDs that are attached to vx2
2840  std::vector<int> tmp;
2841  if (vx2.ID == 0) return tmp;
2842  for (auto& tj : slc.tjs) {
2843  if (tj.AlgMod[kKilled]) continue;
2844  if (tj.CTP != vx2.CTP) continue;
2845  for (unsigned short end = 0; end < 2; ++end) {
2846  if (tj.VtxID[end] == vx2.ID) tmp.push_back(tj.ID);
2847  } // end
2848  } // tj
2849  return tmp;
2850  } // GetVtxTjIDs
2851 
2852  //////////////////////////////////////////
2853  std::vector<int>
2854  GetVtxTjIDs(const TCSlice& slc, const Vtx3Store& vx3, float& score)
2855  {
2856  // returns a list of Tjs in all planes that are attached to vx3
2857  std::vector<int> tmp;
2858  if (vx3.ID == 0) return tmp;
2859  float nvx2 = 0;
2860  score = 0;
2861  for (auto& vx2 : slc.vtxs) {
2862  if (vx2.ID == 0) continue;
2863  if (vx2.Vx3ID != vx3.ID) continue;
2864  auto vtxTjID2 = GetVtxTjIDs(slc, vx2);
2865  tmp.insert(tmp.end(), vtxTjID2.begin(), vtxTjID2.end());
2866  score += vx2.Score;
2867  ++nvx2;
2868  } // vx2
2869  if (nvx2 < 1) return tmp;
2870  // find the average score
2871  score /= nvx2;
2872  // sort by increasing ID
2873  std::sort(tmp.begin(), tmp.end());
2874  return tmp;
2875  } // GetVtxTjIDs
2876 
2877  //////////////////////////////////////////
2878  void
2879  PosInPlane(detinfo::DetectorPropertiesData const& detProp,
2880  const TCSlice& slc,
2881  const Vtx3Store& vx3,
2882  unsigned short plane,
2883  Point2_t& pos)
2884  {
2885  // returns the 2D position of the vertex in the plane
2886  pos[0] = tcc.geom->WireCoordinate(vx3.Y, vx3.Z, plane, vx3.TPCID.TPC, vx3.TPCID.Cryostat);
2887  pos[1] =
2888  detProp.ConvertXToTicks(vx3.X, plane, vx3.TPCID.TPC, vx3.TPCID.Cryostat) * tcc.unitsPerTick;
2889 
2890  } // PosInPlane
2891 
2892  /////////////////////////////////////////
2893  unsigned short
2894  IsCloseToVertex(const TCSlice& slc, const VtxStore& inVx2)
2895  {
2896  // Returns the ID of a 2D vertex having the minimum pull < user-specified cut
2897 
2898  float minPull = tcc.vtx2DCuts[3];
2899  unsigned short imBest = 0;
2900  for (auto& vx2 : slc.vtxs) {
2901  if (vx2.CTP != inVx2.CTP) continue;
2902  if (vx2.ID <= 0) continue;
2903  float pull = VertexVertexPull(slc, inVx2, vx2);
2904  if (pull < minPull) {
2905  minPull = pull;
2906  imBest = vx2.ID;
2907  }
2908  } // vx2
2909  return imBest;
2910  } // IsCloseToVertex
2911 
2912  /////////////////////////////////////////
2913  unsigned short
2914  IsCloseToVertex(const TCSlice& slc, const Vtx3Store& vx3)
2915  {
2916  // Returns the ID of a 3D vertex having the minimum pull < user-specified cut
2917 
2918  float minPull = tcc.vtx3DCuts[1];
2919  unsigned short imBest = 0;
2920  for (auto& oldvx3 : slc.vtx3s) {
2921  if (oldvx3.ID == 0) continue;
2922  if (std::abs(oldvx3.X - vx3.X) > tcc.vtx3DCuts[0]) continue;
2923  float pull = VertexVertexPull(slc, vx3, oldvx3);
2924  if (pull < minPull) {
2925  minPull = pull;
2926  imBest = oldvx3.ID;
2927  }
2928  } // oldvx3
2929  return imBest;
2930 
2931  } // IsCloseToVertex
2932 
2933 } // namespace
tca::kTestBeam
Expect tracks entering from the front face. Don&#39;t create neutrino PFParticles.
Definition: DataStructs.h:532
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:1217
tca::TrajPoint::Dir
Vector2_t Dir
Definition: DataStructs.h:156
Environment
BEGIN_PROLOG or simple_flux see Environment
Definition: genie.fcl:8
tca::TrajPoint::Pos
Point2_t Pos
Definition: DataStructs.h:155
tca::MCSMom
short MCSMom(const TCSlice &slc, const std::vector< int > &tjIDs)
Definition: Utils.cxx:3468
z
process_name opflash particleana ie ie ie z
Definition: run_opflash_electron.fcl:92
tca::VtxStore::PosErr
Point2_t PosErr
Definition: DataStructs.h:74
tca::TCSlice::tjs
std::vector< Trajectory > tjs
vector of all trajectories in each plane
Definition: DataStructs.h:670
tca::FarEnd
unsigned short FarEnd(const TCSlice &slc, const PFPStruct &pfp, const Point3_t &pos)
Definition: PFPUtils.cxx:3339
tca::Length
float Length(const PFPStruct &pfp)
Definition: PFPUtils.cxx:3304
tca::MergeWithVertex
bool MergeWithVertex(TCSlice &slc, VtxStore &vx, unsigned short oVxID)
Definition: TCVertex.cxx:435
if
then if[["$THISISATEST"==1]]
Definition: neoSmazza.sh:95
tca::MakeVertexObsolete
bool MakeVertexObsolete(std::string fcnLabel, TCSlice &slc, VtxStore &vx2, bool forceKill)
Definition: TCVertex.cxx:2726
tca::CloseEnd
unsigned short CloseEnd(const TCSlice &slc, const Trajectory &tj, const Point2_t &pos)
Definition: Utils.cxx:2551
tca::VtxStore
struct of temporary 2D vertices (end points)
Definition: DataStructs.h:72
tca::GetAssns
std::vector< int > GetAssns(TCSlice &slc, std::string type1Name, int id, std::string type2Name)
Definition: Utils.cxx:4849
evd::kWire
static unsigned int kWire
Definition: TWQMultiTPCProjection.cxx:45
tca::TrajPointVertexPull
float TrajPointVertexPull(const TCSlice &slc, const TrajPoint &tp, const VtxStore &vx)
Definition: TCVertex.cxx:1851
tca::AttachAnyVertexToTraj
bool AttachAnyVertexToTraj(TCSlice &slc, int tjID, bool prt)
Definition: TCVertex.cxx:1641
tca::TrajPoint::CTP
CTP_t CTP
Cryostat, TPC, Plane code.
Definition: DataStructs.h:153
tca::TCConfig::maxPos0
std::vector< float > maxPos0
Definition: DataStructs.h:572
tca::CompleteIncomplete3DVerticesInGaps
void CompleteIncomplete3DVerticesInGaps(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: TCVertex.cxx:2396
tca::Point3_t
std::array< double, 3 > Point3_t
Definition: DataStructs.h:41
tca::SignalAtTp
bool SignalAtTp(TrajPoint &tp)
Definition: Utils.cxx:2004
tca::FitVertex
bool FitVertex(TCSlice &slc, VtxStore &vx, bool prt)
Definition: TCVertex.cxx:1954
tca::TCConfig::vtx3DCuts
std::vector< float > vtx3DCuts
2D vtx -&gt; 3D vtx matching cuts
Definition: DataStructs.h:551
tca::SetPDGCode
void SetPDGCode(TCSlice &slc, unsigned short itj)
Definition: Utils.cxx:4350
tca::Find3DVertices
void Find3DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: TCVertex.cxx:1268
tca::Print2V
void Print2V(std::string someText, mf::LogVerbatim &myprt, VtxStore &vx2, bool &printHeader)
Definition: Utils.cxx:5768
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::kFixed
vertex position fixed manually - no fitting done
Definition: DataStructs.h:93
evd::err
EResult err(const char *call)
Definition: WebEVDServer.cxx:140
tca::SplitTrajCrossingVertices
void SplitTrajCrossingVertices(TCSlice &slc, CTP_t inCTP)
Definition: TCVertex.cxx:923
tca::PrintTP
void PrintTP(std::string someText, const TCSlice &slc, unsigned short ipt, short dir, unsigned short pass, const TrajPoint &tp)
Definition: Utils.cxx:6296
tca::Vtx3Store::Vx2ID
std::vector< int > Vx2ID
Definition: DataStructs.h:114
tca::Reconcile2Vs
void Reconcile2Vs(TCSlice &slc)
Definition: TCVertex.cxx:1067
geo::kX
Planes which measure X direction.
Definition: geo_types.h:134
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
tca::kHiVx3Score
matched to a high-score 3D vertex
Definition: DataStructs.h:95
score
BEGIN_PROLOG or score(default)}sbnd_crttrackmatchingalg_crID
Definition: crttrackmatchingalg_sbnd.fcl:9
tca::kStepDir
step from US -&gt; DS (true) or DS -&gt; US (false)
Definition: DataStructs.h:531
tca::Trajectory::MCSMom
short MCSMom
Normalized RMS using ALL hits. Assume it is 50% to start.
Definition: DataStructs.h:200
tca::RefineVtxPosition
bool RefineVtxPosition(TCSlice &slc, const Trajectory &tj, unsigned short &nearPt, short nPtsToChk, bool prt)
Definition: TCVertex.cxx:2690
tca::VtxStore::Pass
unsigned short Pass
Definition: DataStructs.h:76
tca::PrintPos
std::string PrintPos(const TCSlice &slc, const TrajPoint &tp)
Definition: Utils.cxx:6526
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
tca::TrajClosestApproach
bool TrajClosestApproach(Trajectory const &tj, float x, float y, unsigned short &closePt, float &DOCA)
Definition: Utils.cxx:2690
tca::PosInPlane
void PosInPlane(detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, const Vtx3Store &vx3, unsigned short plane, Point2_t &pos)
Definition: TCVertex.cxx:2879
tca::AttachAnyTrajToVertex
bool AttachAnyTrajToVertex(TCSlice &slc, unsigned short ivx, bool prt)
Definition: TCVertex.cxx:1683
tca::TPNearVertex
unsigned short TPNearVertex(const TCSlice &slc, const TrajPoint &tp)
Definition: TCVertex.cxx:1568
tca::PosAtEnd
Point3_t PosAtEnd(const PFPStruct &pfp, unsigned short end)
Definition: PFPUtils.cxx:3293
tca::PointTrajDOCA
float PointTrajDOCA(const TCSlice &slc, unsigned int iht, TrajPoint const &tp)
Definition: Utils.cxx:2572
tca::DeltaAngle
double DeltaAngle(const Vector3_t v1, const Vector3_t v2)
Definition: PFPUtils.cxx:2540
tca::TCConfig::dbgSlc
bool dbgSlc
debug only in the user-defined slice? default is all slices
Definition: DataStructs.h:590
tca::PFPStruct::Vx3ID
std::array< int, 2 > Vx3ID
Definition: DataStructs.h:290
tca::SetVx3Score
void SetVx3Score(TCSlice &slc, Vtx3Store &vx3)
Definition: TCVertex.cxx:2242
util::quantities::ticks
tick ticks
Alias for common language habits.
Definition: electronics.h:78
tca::CompleteIncomplete3DVertices
void CompleteIncomplete3DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: TCVertex.cxx:2490
tca::StoreVertex
bool StoreVertex(TCSlice &slc, VtxStore &vx)
Definition: TCVertex.cxx:1918
tca::TCConfig::dbg3V
bool dbg3V
debug 3D vertex finding
Definition: DataStructs.h:597
GeometryCore.h
Access the description of detector geometry.
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
tca::Vtx3Store
struct of temporary 3D vertices
Definition: DataStructs.h:104
tca::PFPStruct::TPCID
geo::TPCID TPCID
Definition: DataStructs.h:295
tca::PrintAllTraj
void PrintAllTraj(detinfo::DetectorPropertiesData const &detProp, std::string someText, TCSlice &slc, unsigned short itj, unsigned short ipt, bool prtVtx)
Definition: Utils.cxx:5952
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
tca::Point2_t
std::array< float, 2 > Point2_t
Definition: DataStructs.h:43
tca::TCConfig::maxPos1
std::vector< float > maxPos1
Definition: DataStructs.h:573
tca::TCConfig::unitsPerTick
float unitsPerTick
scale factor from Tick to WSE equivalent units
Definition: DataStructs.h:571
tca::DeadWireCount
float DeadWireCount(const TCSlice &slc, const TrajPoint &tp1, const TrajPoint &tp2)
Definition: Utils.cxx:2141
tca::NearbyCleanPt
unsigned short NearbyCleanPt(const TCSlice &slc, const Trajectory &tj, unsigned short end)
Definition: Utils.cxx:2956
tca::Vtx3Store::TPCID
geo::TPCID TPCID
Definition: DataStructs.h:113
tca::debug
DebugStuff debug
Definition: DebugStruct.cxx:4
tca::PointDirection
Vector3_t PointDirection(const Point3_t p1, const Point3_t p2)
Definition: PFPUtils.cxx:2548
tca::Trajectory::CTP
CTP_t CTP
Cryostat, TPC, Plane code.
Definition: DataStructs.h:192
tca::kVxEnvOK
the environment near the vertex was checked - See UpdateVxEnvironment
Definition: DataStructs.h:99
tca::TCConfig::dbg2V
bool dbg2V
debug 2D vertex finding
Definition: DataStructs.h:593
tca::Trajectory::Pts
std::vector< TrajPoint > Pts
Trajectory points.
Definition: DataStructs.h:191
tca::ChkVtxAssociations
bool ChkVtxAssociations(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:2071
tca::Reconcile2VTs
bool Reconcile2VTs(TCSlice &slc, std::vector< int > &vx2cls, bool prt)
Definition: TCVertex.cxx:1163
tca::DebugStuff::Plane
int Plane
Select plane.
Definition: DebugStruct.h:22
for
for($it=0;$it< $RaceTrack_number;$it++)
Definition: generate_icarus_geometry_splitwire.pl:1671
tca::ChgFracNearPos
float ChgFracNearPos(const TCSlice &slc, const Point2_t &pos, const std::vector< int > &tjIDs)
Definition: Utils.cxx:3236
tca::FindCloseHits
std::vector< unsigned int > FindCloseHits(const 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:2843
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
tca::TCSlice::vtxs
std::vector< VtxStore > vtxs
2D vertices
Definition: DataStructs.h:676
tca::Trajectory::EndPt
std::array< unsigned short, 2 > EndPt
First and last point in the trajectory that has charge.
Definition: DataStructs.h:203
tca::Trajectory::PDGCode
unsigned short PDGCode
shower-like or track-like {default is track-like}
Definition: DataStructs.h:208
detinfo::DetectorPropertiesData::ConvertXToTicks
double ConvertXToTicks(double X, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:99
tca::ScoreVertices
void ScoreVertices(TCSlice &slc)
Definition: TCVertex.cxx:2146
tca::PosSep2
double PosSep2(const Point3_t &pos1, const Point3_t &pos2)
Definition: PFPUtils.cxx:2572
tca::PrintHit
std::string PrintHit(const TCHit &tch)
Definition: Utils.cxx:6516
tca::VtxStore::Pos
Point2_t Pos
Definition: DataStructs.h:73
tca::SplitTraj
bool SplitTraj(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, unsigned short itj, float XPos, bool makeVx2, bool prt)
Definition: Utils.cxx:2272
tca::TrajPointTrajDOCA
void TrajPointTrajDOCA(const TCSlice &slc, TrajPoint const &tp, Trajectory const &tj, unsigned short &closePt, float &minSep)
Definition: Utils.cxx:2435
tca::TCConfig::geom
const geo::GeometryCore * geom
Definition: DataStructs.h:576
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
tca::NearestPtWithChg
unsigned short NearestPtWithChg(const TCSlice &slc, const Trajectory &tj, unsigned short thePt)
Definition: Utils.cxx:3522
tca::detail::valsIncreasing
bool valsIncreasing(const SortEntry &c1, const SortEntry &c2)
Definition: Utils.cxx:39
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:1746
geo_types.h
Definition of data types for geometry description.
tca::Trajectory::ID
int ID
ID that is local to one slice.
Definition: DataStructs.h:205
tca::Trajectory::VtxID
std::array< unsigned short, 2 > VtxID
ID of 2D vertex.
Definition: DataStructs.h:202
tca::kVxIndPlnNoChg
vertex quality is suspect - No requirement made on chg btw it and the Tj
Definition: DataStructs.h:98
dropbox.dir
tuple dir
Definition: dropbox.py:28
detinfo::DetectorPropertiesData::ConvertTicksToX
double ConvertTicksToX(double ticks, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:115
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:607
tca::TCSlice::slHits
std::vector< TCHit > slHits
Definition: DataStructs.h:669
tca::IsCloseToVertex
unsigned short IsCloseToVertex(const TCSlice &slc, const VtxStore &inVx2)
Definition: TCVertex.cxx:2894
tca::VertexVertexPull
float VertexVertexPull(const TCSlice &slc, const Vtx3Store &vx1, const Vtx3Store &vx2)
Definition: TCVertex.cxx:1887
tca::VtxStore::Stat
std::bitset< 16 > Stat
Vertex status bits using kVtxBit_t.
Definition: DataStructs.h:88
tca::VtxStore::ID
int ID
set to 0 if killed
Definition: DataStructs.h:83
tca::TrajIntersection
void TrajIntersection(TrajPoint const &tp1, TrajPoint const &tp2, Point2_t &pos)
Definition: Utils.cxx:2604
tca::TCConfig::vtxScoreWeights
std::vector< float > vtxScoreWeights
Definition: DataStructs.h:552
tca::DotProd
double DotProd(const Vector3_t &v1, const Vector3_t &v2)
Definition: PFPUtils.h:128
tca::CTP_t
unsigned int CTP_t
Definition: DataStructs.h:47
tca::SetHighScoreBits
void SetHighScoreBits(TCSlice &slc, Vtx3Store &vx3)
Definition: TCVertex.cxx:2193
tca::TCSlice::vtx3s
std::vector< Vtx3Store > vtx3s
3D vertices
Definition: DataStructs.h:677
tca::TCSlice::TPCID
geo::TPCID TPCID
Definition: DataStructs.h:662
tca::AttachTrajToVertex
bool AttachTrajToVertex(TCSlice &slc, Trajectory &tj, VtxStore &vx, bool prt)
Definition: TCVertex.cxx:1728
tca::TCConfig::useAlg
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. &lt; 0 to turn off.
Definition: DataStructs.h:586
tca::TrajLength
float TrajLength(const Trajectory &tj)
Definition: Utils.cxx:2646
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
tca::Trajectory::AlgMod
std::bitset< 128 > AlgMod
Bit set if algorithm AlgBit_t modifed the trajectory.
Definition: DataStructs.h:193
tca::UpdateVxEnvironment
void UpdateVxEnvironment(TCSlice &slc)
Definition: Utils.cxx:3862
tca::VtxStore::NTraj
unsigned short NTraj
Definition: DataStructs.h:75
tca::DecodeCTP
geo::PlaneID DecodeCTP(CTP_t CTP)
Definition: DataStructs.cxx:116
tca::GetVtxTjIDs
std::vector< int > GetVtxTjIDs(const TCSlice &slc, const VtxStore &vx2)
Definition: TCVertex.cxx:2837
tca::PosSep
double PosSep(const Point3_t &pos1, const Point3_t &pos2)
Definition: PFPUtils.cxx:2565
tca::TCConfig::dbgAlg
std::bitset< 128 > dbgAlg
Allow user to turn on debug printing in algorithms (that print...)
Definition: DataStructs.h:587
tca::detail::valsDecreasing
bool valsDecreasing(const SortEntry &c1, const SortEntry &c2)
Definition: Utils.cxx:38
tca::EncodeCTP
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
Definition: DataStructs.h:52
tca::KillPoorVertices
void KillPoorVertices(TCSlice &slc)
Definition: TCVertex.cxx:2174
tca::Trajectory::EndFlag
std::array< std::bitset< 8 >, 2 > EndFlag
Definition: DataStructs.h:212
tca::TCSlice::ID
int ID
ID of the recob::Slice (not the sub-slice)
Definition: DataStructs.h:664
tca::VtxStore::UID
int UID
unique global ID
Definition: DataStructs.h:84
tca::SetVx2Score
void SetVx2Score(TCSlice &slc)
Definition: TCVertex.cxx:2264
tca::FindHammerVertices
void FindHammerVertices(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:797
tca::TCConfig::vtx2DCuts
std::vector< float > vtx2DCuts
Max position pull, max Position error rms.
Definition: DataStructs.h:550
tca::TCSlice::nPlanes
unsigned short nPlanes
Definition: DataStructs.h:663
tca::TCSlice::pfps
std::vector< PFPStruct > pfps
Definition: DataStructs.h:678
tca::FindCloseTjs
std::vector< int > FindCloseTjs(const TCSlice &slc, const TrajPoint &fromTp, const TrajPoint &toTp, const float &maxDelta)
Definition: Utils.cxx:2977
tca::SignalBetween
bool SignalBetween(const TCSlice &slc, const TrajPoint &tp1, const TrajPoint &tp2, const float &MinWireSignalFraction)
Definition: Utils.cxx:1807
tca::MoveTPToWire
void MoveTPToWire(TrajPoint &tp, float wire)
Definition: Utils.cxx:2831
tca::NumPtsWithCharge
unsigned short NumPtsWithCharge(const TCSlice &slc, const Trajectory &tj, bool includeDeadWires)
Definition: Utils.cxx:2116
tca::MakeJunkVertices
void MakeJunkVertices(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:31
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::Find2DVertices
void Find2DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, const CTP_t &inCTP, unsigned short pass)
Definition: TCVertex.cxx:132
close
print OUTPUT<< EOF;< setup name="Default"version="1.0">< worldref="volWorld"/></setup ></gdml > EOF close(OUTPUT)
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
lib.dedx.A
float A
Definition: dedx.py:137
tca::TrajTrajDOCA
bool TrajTrajDOCA(const TCSlice &slc, const Trajectory &tj1, const Trajectory &tj2, unsigned short &ipt1, unsigned short &ipt2, float &minSep)
Definition: Utils.cxx:2460
tca::DirAtEnd
Vector3_t DirAtEnd(const PFPStruct &pfp, unsigned short end)
Definition: PFPUtils.cxx:3284
tca::FindHammerVertices2
void FindHammerVertices2(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:604
tca::VtxStore::TjChgFrac
float TjChgFrac
Fraction of charge near the vertex that is from hits on the vertex Tjs.
Definition: DataStructs.h:87
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:533
tca::AttachToAnyVertex
bool AttachToAnyVertex(TCSlice &slc, PFPStruct &pfp, float maxSep, bool prt)
Definition: TCVertex.cxx:1583
tca::PrintTPHeader
void PrintTPHeader(std::string someText)
Definition: Utils.cxx:6287
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
tca::DebugStuff::CTP
CTP_t CTP
set to an invalid CTP
Definition: DebugStruct.h:23
tca::TrajPoint::HitPosErr2
double HitPosErr2
Definition: DataStructs.h:157
tca::detail::SortEntry::index
unsigned int index
Definition: DataStructs.h:36
Generated by   doxygen 1.8.5