LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
TrajClusterAlg.cxx
Go to the documentation of this file.
1 
24 
25 #include "fhiclcpp/ParameterSet.h"
27 
28 #include <iostream>
29 #include <string>
30 #include <vector>
31 
32 namespace tca {
33 
34  //------------------------------------------------------------------------------
35 
37  : fCaloAlg(pset.get<fhicl::ParameterSet>("CaloAlg")), fMVAReader("Silent")
38  {
40 
41  bool badinput = false;
42  // set all configurable modes false
43  tcc.modes.reset();
44 
45  // default mode is stepping US -> DS
46  tcc.modes[kStepDir] = true;
47  short userMode = 1;
48  if (pset.has_key("Mode")) userMode = pset.get<short>("Mode");
49  if (userMode < 0) tcc.modes[kStepDir] = false;
50  if (pset.has_key("DoForecast")) tcc.doForecast = pset.get<bool>("DoForecast");
51  if (pset.has_key("UseChannelStatus")) tcc.useChannelStatus = pset.get<bool>("UseChannelStatus");
52  if (pset.has_key("StudyMode")) {
53  std::cout << "StudyMode is not valid anymore. Try Study: 1 or Study: 2, etc/n";
54  } // old StudyMode
55  if (pset.has_key("Study")) {
56  userMode = pset.get<short>("Study");
57  if (userMode == 1) tcc.modes[kStudy1] = true;
58  if (userMode == 2) tcc.modes[kStudy2] = true;
59  if (userMode == 3) tcc.modes[kStudy3] = true;
60  if (userMode == 4) tcc.modes[kStudy4] = true;
61  } // new Study mode
62  if (pset.has_key("SaveShowerTree"))
63  tcc.modes[kSaveShowerTree] = pset.get<bool>("SaveShowerTree");
64  if (pset.has_key("SaveCRTree")) tcc.modes[kSaveCRTree] = pset.get<bool>("SaveCRTree");
65  if (pset.has_key("TagCosmics")) tcc.modes[kTagCosmics] = pset.get<bool>("TagCosmics");
66  std::vector<std::string> skipAlgsVec;
67  if (pset.has_key("SkipAlgs")) skipAlgsVec = pset.get<std::vector<std::string>>("SkipAlgs");
68  std::vector<std::string> debugConfigVec;
69  if (pset.has_key("DebugConfig"))
70  debugConfigVec = pset.get<std::vector<std::string>>("DebugConfig");
71 
72  tcc.hitErrFac = pset.get<float>("HitErrFac", 0.4);
73  // Allow the user to specify the typical hit rms for small-angle tracks
74  std::vector<float> aveHitRMS;
75  if (pset.has_key("AveHitRMS")) aveHitRMS = pset.get<std::vector<float>>("AveHitRMS");
76  // Turn off the call to AnalyzeHits
77  if (!aveHitRMS.empty()) {
78  evt.aveHitRMSValid = true;
79  evt.aveHitRMS = aveHitRMS;
80  }
81  tcc.angleRanges = pset.get<std::vector<float>>("AngleRanges");
82  tcc.nPtsAve = pset.get<short>("NPtsAve", 20);
83  tcc.minPtsFit = pset.get<std::vector<unsigned short>>("MinPtsFit");
84  tcc.minPts = pset.get<std::vector<unsigned short>>("MinPts");
85  tcc.maxAngleCode = pset.get<std::vector<unsigned short>>("MaxAngleCode");
86  tcc.minMCSMom = pset.get<std::vector<short>>("MinMCSMom");
87  tcc.maxChi = pset.get<float>("MaxChi", 10);
88  tcc.chargeCuts = pset.get<std::vector<float>>("ChargeCuts", {3, 0.15, 0.25});
89  tcc.multHitSep = pset.get<float>("MultHitSep", 2.5);
90  tcc.kinkCuts = pset.get<std::vector<float>>("KinkCuts");
91  tcc.qualityCuts = pset.get<std::vector<float>>("QualityCuts", {0.8, 3});
92  tcc.maxWireSkipNoSignal = pset.get<float>("MaxWireSkipNoSignal", 1);
93  tcc.maxWireSkipWithSignal = pset.get<float>("MaxWireSkipWithSignal", 100);
94  tcc.projectionErrFactor = pset.get<float>("ProjectionErrFactor", 2);
95  tcc.VLAStepSize = pset.get<float>("VLAStepSize", 1.5);
96  tcc.JTMaxHitSep2 = pset.get<float>("JTMaxHitSep", 2);
97  tcc.deltaRayTag = pset.get<std::vector<short>>("DeltaRayTag", {-1, -1, -1});
98  tcc.muonTag = pset.get<std::vector<short>>("MuonTag", {-1, -1, -1, -1});
99  if (pset.has_key("ElectronTag")) tcc.electronTag = pset.get<std::vector<float>>("ElectronTag");
100  tcc.showerTag = pset.get<std::vector<float>>("ShowerTag", {-1, -1, -1, -1, -1, -1});
101  std::string fMVAShowerParentWeights = "NA";
102  pset.get_if_present<std::string>("MVAShowerParentWeights", fMVAShowerParentWeights);
103  tcc.chkStopCuts = pset.get<std::vector<float>>("ChkStopCuts", {-1, -1, -1});
104  if (pset.has_key("MatchTruth")) {
105  std::cout << "MatchTruth is not used. Use ClusterAnaV2 or DebugConfig to configure\n";
106  }
107  tcc.vtx2DCuts = pset.get<std::vector<float>>("Vertex2DCuts", {-1, -1, -1, -1, -1, -1, -1});
108  tcc.vtx3DCuts = pset.get<std::vector<float>>("Vertex3DCuts", {-1, -1});
109  tcc.vtxScoreWeights = pset.get<std::vector<float>>("VertexScoreWeights");
110  tcc.match3DCuts = pset.get<std::vector<float>>("Match3DCuts", {-1, -1, -1, -1, -1});
111  tcc.pfpStitchCuts = pset.get<std::vector<float>>("PFPStitchCuts", {-1});
112  // don't search for a neutrino vertex in test beam mode
113  tcc.modes[kTestBeam] = pset.get<bool>("TestBeam", false);
114  pset.get_if_present<std::vector<float>>("NeutralVxCuts", tcc.neutralVxCuts);
116 
117  // in the following section we ensure that the fcl vectors are appropriately sized so that later references are valid
118  if (tcc.minPtsFit.size() != tcc.minPts.size()) badinput = true;
119  if (tcc.maxAngleCode.size() != tcc.minPts.size()) badinput = true;
120  if (tcc.minMCSMom.size() != tcc.minPts.size()) badinput = true;
121  if (badinput)
123  << "Bad input from fcl file. Vector lengths for MinPtsFit, MaxAngleRange and MinMCSMom "
124  "should be defined for each reconstruction pass";
125 
126  if (tcc.vtx2DCuts.size() < 10)
128  << "Vertex2DCuts must be size 11\n 0 = Max length definition for short TJs\n 1 = Max "
129  "vtx-TJ sep short TJs\n 2 = Max vtx-TJ sep long TJs\n 3 = Max position pull for >2 "
130  "TJs\n 4 = Max vtx position error\n 5 = Min MCSMom for one of two TJs\n 6 = Min "
131  "fraction of wires hit btw vtx and Tjs\n 7 = Min Score\n 8 = min ChgFrac at a vtx or "
132  "merge point\n 9 = max MCSMom asymmetry, 10 = require chg on wires btw vtx and tjs in "
133  "induction planes?";
134  if (tcc.vtx2DCuts.size() == 10) {
135  // User didn't specify a requirement for the presence of charge between a vertex and the start of the
136  // vertex Tjs in induction planes. Assume that it is required
137  tcc.vtx2DCuts.resize(11, 1.);
138  }
139  if (tcc.vtx3DCuts.size() < 3)
141  << "Vertex3DCuts must be size > 2\n 0 = 2D Vtx max dX (cm)\n 1 = 2D Vtx max pull\n 2 = max "
142  "3D separation (cm) btw PFP and vertex";
143  if (tcc.vtx3DCuts.size() == 2) {
144  std::cout << "WARNING: Vertex3DCuts is size 2 but should be size 3, where Vertex3DCuts[2] = "
145  "max 3D separation (cm) btw a PFP and a 3D vertex. Setting it to 3 cm\n";
146  tcc.vtx3DCuts.resize(3, 3.);
147  }
148  if (tcc.kinkCuts.size() < 4) {
150  << "KinkCuts must be size 3\n 0 = Number of points to fit at the end of the trajectory\n 1 "
151  "= Minimum kink significance\n 2 = Use charge in significance calculation? (yes if > "
152  "0)\n 3 = 3D kink fit length (cm). \nYou are using an out-of-date specification?\n";
153  }
154  // throw an exception if the user appears to be using an old version of KinkCuts where
155  // KinkCuts[0] was a kink angle cut
156  if (tcc.kinkCuts[0] > 0 && tcc.kinkCuts[0] < 1.) {
158  << "Are you using an out-of-date specification for KinkCuts? KinkCuts[0] is the number of "
159  "points to fit.\n";
160  }
161 
162  if (tcc.chargeCuts.size() != 3)
164  << "ChargeCuts must be size 3\n 0 = Charge pull cut\n 1 = Min allowed fractional chg RMS\n "
165  "2 = Max allowed fractional chg RMS";
166  // dressed muons - change next line
167  if (tcc.muonTag.size() < 4)
169  << "MuonTag must be size 4\n 0 = minPtsFit\n 1 = minMCSMom\n 2= maxWireSkipNoSignal\n 3 = "
170  "min delta ray length for tagging\n 4 = dress muon window size (optional)";
171  if (tcc.deltaRayTag.size() != 3)
173  << "DeltaRayTag must be size 3\n 0 = Max endpoint sep\n 1 = min MCSMom\n 2 = max MCSMom";
174  if (tcc.chkStopCuts.size() < 3)
176  << "ChkStopCuts must be size 3\n 0 = Min Charge ratio\n 1 = Charge slope pull cut\n 2 = "
177  "Charge fit chisq cut\n 3 = BraggSplit FOM (optional)";
178  if (tcc.showerTag.size() < 13) {
179  std::cout
180  << "ShowerTag must be size 13\n 0 = Mode\n 1 = max MCSMom\n 2 = max separation (WSE "
181  "units)\n 3 = Max angle diff\n 4 = Factor * rms width\n 5 = Min half width\n 6 = min "
182  "total Tps\n 7 = Min Tjs\n 8 = max parent FOM\n 9 = max direction FOM 10 = max "
183  "AspectRatio\n 11 = min Score to preserve a vertex\n 12 = Debug showers in CTP\n";
184  std::cout << " Fixing this problem...";
185  tcc.showerTag.resize(13);
186  // set the min score to 0
187  tcc.showerTag[11] = 0;
188  // turn off printing
189  tcc.showerTag[12] = -1;
190  }
191  if (tcc.match3DCuts.size() < 6)
193  << "Match3DCuts must be size 5\n 0 = dx (cm) matching cut\n 1 = max number of 3D "
194  "combinations\n 2 = min length for 2-view match\n 3 = number of TP3Ds in each plane to "
195  "fit in each PFP section\n 4 = max pull for accepting TP3Ds in sections\n 5 = max "
196  "ChiDOF for a SectionFit";
197 
198  // check the angle ranges and convert from degrees to radians
199  if (tcc.angleRanges.back() < 90) {
200  mf::LogVerbatim("TC") << "Last element of AngleRange != 90 degrees. Fixing it\n";
201  tcc.angleRanges.back() = 90;
202  }
203 
204  // convert PFP stitch cuts
205  if (tcc.pfpStitchCuts.size() > 1 && tcc.pfpStitchCuts[0] > 0) {
206  // square the separation cut
208  // convert angle to cos
209  tcc.pfpStitchCuts[1] = cos(tcc.pfpStitchCuts[1]);
210  }
211 
212  // configure algorithm debugging. Configuration for debugging standard stepping
213  // is done in Utils/AnalyzeHits when the input hit collection is passed to SetInputHits
214  tcc.modes[kDebug] = false;
215  tcc.dbgAlg.reset();
216  for (auto strng : debugConfigVec) {
217  // try to interpret this as a C:T:P:W:Tick specification or something similar
218  if (!DecodeDebugString(strng)) {
219  // try to set a dbgAlg bit
220  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib) {
221  if (strng == AlgBitNames[ib]) {
222  tcc.dbgAlg[ib] = true;
223  tcc.modes[kDebug] = true;
224  break;
225  }
226  } // ib
227  // print some instructions and quit if there was a failure
228  if (!tcc.modes[kDebug]) {
229  std::cout << "DecodeDebugString failed: " << strng << "\n";
230  DecodeDebugString("instruct");
231  exit(1);
232  }
233  } // DecodeDebugString failed
234  } // strng
235 
236  for (auto& range : tcc.angleRanges) {
237  if (range < 0 || range > 90)
239  << "Invalid angle range " << range << " Must be 0 - 90 degrees";
240  range *= M_PI / 180;
241  } // range
242 
243  // Ensure that the size of the AlgBitNames vector is consistent with the AlgBit typedef enum
244  if (kAlgBitSize != AlgBitNames.size())
246  << "kAlgBitSize " << kAlgBitSize << " != AlgBitNames size " << AlgBitNames.size();
247  if (kAlgBitSize > 128)
249  << "Increase the size of UseAlgs to at least " << kAlgBitSize;
250  fAlgModCount.resize(kAlgBitSize);
251 
252  if (kFlagBitSize != EndFlagNames.size())
254  << "kFlagBitSize " << kFlagBitSize << " != EndFlagNames size " << EndFlagNames.size();
255 
256  if (kFlagBitSize > 8)
258  << "Increase the size of EndFlag to at least " << kFlagBitSize;
259 
260  bool printHelp = false;
261  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
262  tcc.useAlg[ib] = true;
263 
264  // turn off the special algs
265  // Do an exhaustive (and slow) check of the hit -> trajectory associations
266  tcc.useAlg[kChkInTraj] = false;
267 
268  for (auto strng : skipAlgsVec) {
269  bool gotit = false;
270  if (strng == "All") {
271  // turn everything off
272  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
273  tcc.useAlg[ib] = false;
274  gotit = true;
275  break;
276  } // All off
277  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib) {
278  if (strng == AlgBitNames[ib]) {
279  tcc.useAlg[ib] = false;
280  gotit = true;
281  break;
282  }
283  } // ib
284  if (!gotit) {
285  std::cout << "******* Unknown SkipAlgs input string '" << strng << "'\n";
286  printHelp = true;
287  }
288  } // strng
289  if (printHelp) {
290  std::cout << "Valid AlgNames:";
291  for (auto strng : AlgBitNames)
292  std::cout << " " << strng;
293  std::cout << "\n";
294  std::cout << "Or specify All to turn all algs off\n";
295  }
296  // Configure the TMVA reader for the shower parent BDT
297  if (fMVAShowerParentWeights != "NA" && tcc.showerTag[0] > 0)
298  ConfigureMVA(tcc, fMVAShowerParentWeights);
299 
300  evt.eventsProcessed = 0;
301 
302  tcc.caloAlg = &fCaloAlg;
303  }
304 
306  bool TrajClusterAlg::SetInputHits(std::vector<recob::Hit> const& inputHits,
307  unsigned int run,
308  unsigned int event)
309  {
310  // defines the pointer to the input hit collection, analyzes them,
311  // initializes global counters and refreshes service references
312  ClearResults();
313  evt.allHits = &inputHits;
314  evt.run = run;
315  evt.event = event;
316  // refresh service references
317  tcc.geom = lar::providerFrom<geo::Geometry>();
319  evt.WorkID = 0;
320  evt.globalT_UID = 0;
321  evt.global2V_UID = 0;
322  evt.global3V_UID = 0;
323  evt.globalP_UID = 0;
324  evt.global2S_UID = 0;
325  evt.global3S_UID = 0;
326  // find the average hit RMS using the full hit collection and define the
327  // configuration for the current TPC
328 
330 
331  return AnalyzeHits();
332  } // SetInputHits
333 
335  void TrajClusterAlg::SetSourceHits(std::vector<recob::Hit> const& srcHits)
336  {
337  evt.srcHits = &srcHits;
338  evt.tpcSrcHitRange.resize(tcc.geom->NTPC());
339  for (auto& thr : evt.tpcSrcHitRange)
340  thr = {UINT_MAX, UINT_MAX};
341  for (unsigned int iht = 0; iht < (*evt.srcHits).size(); ++iht) {
342  auto& hit = (*evt.srcHits)[iht];
343  unsigned int tpc = hit.WireID().TPC;
344  if (tpc >= evt.tpcSrcHitRange.size()) return;
345  if (evt.tpcSrcHitRange[tpc].first == UINT_MAX) evt.tpcSrcHitRange[tpc].first = iht;
346  evt.tpcSrcHitRange[tpc].second = iht;
347  } // iht
348  } // SetSourceHits
349 
352  detinfo::DetectorPropertiesData const& detProp,
353  std::vector<unsigned int>& hitsInSlice,
354  int sliceID)
355  {
356  // Reconstruct everything using the hits in a slice
357 
358  if (slices.empty()) ++evt.eventsProcessed;
359  if (hitsInSlice.size() < 2) return;
360  if (tcc.recoSlice > 0 && sliceID != tcc.recoSlice) return;
361 
362  if (!CreateSlice(clockData, detProp, hitsInSlice, sliceID)) return;
363 
364  seeds.resize(0);
365  // get a reference to the stored slice
366  auto& slc = slices[slices.size() - 1];
367  // special debug mode reconstruction
368  if (tcc.recoTPC > 0 && (short)slc.TPCID.TPC != tcc.recoTPC) {
369  slices.pop_back();
370  return;
371  }
372 
373  if (evt.aveHitRMS.size() != slc.nPlanes)
375  << " AveHitRMS vector size != the number of planes ";
376  if (tcc.recoSlice)
377  std::cout << "Reconstruct " << hitsInSlice.size() << " hits in Slice " << sliceID
378  << " in TPC " << slc.TPCID.TPC << "\n";
379  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
380  CTP_t inCTP = EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
381  ReconstructAllTraj(detProp, slc, inCTP);
382  if (!slc.isValid) return;
383  } // plane
384  // Compare 2D vertices in each plane and try to reconcile T -> 2V attachments using
385  // 2D and 3D(?) information
386  Reconcile2Vs(slc);
387  Find3DVertices(detProp, slc);
388  ScoreVertices(slc);
389  // Define the ParentID of trajectories using the vertex score
390  DefineTjParents(slc, false);
391  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
392  CTP_t inCTP = EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
393  if (!ChkVtxAssociations(slc, inCTP)) {
394  std::cout << "RTC: ChkVtxAssociations found an error\n";
395  }
396  } // plane
397  if (tcc.match3DCuts[0] > 0) {
398  FindPFParticles(clockData, detProp, slc);
399  DefinePFPParents(slc);
400  } // 3D matching requested
401  KillPoorVertices(slc);
402  // Use 3D matching information to find showers in 2D. FindShowers3D returns
403  // true if the algorithm was successful indicating that the matching needs to be redone
404  if (tcc.showerTag[0] == 2 || tcc.showerTag[0] == 4) {
405  FindShowers3D(detProp, slc);
406  if (tcc.modes[kSaveShowerTree]) {
407  std::cout << "SHOWER TREE STAGE NUM SIZE: " << stv.StageNum.size() << std::endl;
408  showertree->Fill();
409  }
410  } // 3D shower code
411 
412  if (!slc.isValid) {
413  mf::LogVerbatim("TC") << "RunTrajCluster failed in MakeAllTrajClusters";
414  return;
415  }
416 
417  // dump a trajectory?
418  if (tcc.modes[kDebug] && tcc.dbgDump) DumpTj();
419 
420  Finish3DShowers(slc);
421 
422  // count algorithm usage
423  for (auto& tj : slc.tjs) {
424  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
425  if (tj.AlgMod[ib]) ++fAlgModCount[ib];
426  } // tj
427 
428  // clear vectors that are not needed later
429  slc.mallTraj.resize(0);
430 
431  } // RunTrajClusterAlg
432 
435  TCSlice& slc,
436  CTP_t inCTP)
437  {
438  // Reconstruct trajectories in inCTP and put them in allTraj
439 
440  unsigned short plane = DecodeCTP(inCTP).Plane;
441  if (slc.firstWire[plane] > slc.nWires[plane]) return;
442  unsigned int nwires = slc.lastWire[plane] - slc.firstWire[plane] - 1;
443  if (nwires > slc.nWires[plane]) return;
444 
445  // Make several passes through the hits with user-specified cuts for each
446  // pass. In general these are to not reconstruct large angle trajectories on
447  // the first pass
448  float maxHitsRMS = 4 * evt.aveHitRMS[plane];
449  for (unsigned short pass = 0; pass < tcc.minPtsFit.size(); ++pass) {
450  for (unsigned int ii = 0; ii < nwires; ++ii) {
451  // decide which way to step given the sign of StepDir
452  unsigned int iwire = 0;
453  unsigned int jwire = 0;
454  if (tcc.modes[kStepDir]) {
455  // step DS
456  iwire = slc.firstWire[plane] + ii;
457  jwire = iwire + 1;
458  }
459  else {
460  // step US
461  iwire = slc.lastWire[plane] - ii - 1;
462  jwire = iwire - 1;
463  }
464  if (iwire > slc.wireHitRange[plane].size() - 1) continue;
465  if (jwire > slc.wireHitRange[plane].size() - 1) continue;
466  // skip bad wires or no hits on the wire
467  if (slc.wireHitRange[plane][iwire].first == UINT_MAX) continue;
468  if (slc.wireHitRange[plane][jwire].first == UINT_MAX) continue;
469  unsigned int ifirsthit = slc.wireHitRange[plane][iwire].first;
470  unsigned int ilasthit = slc.wireHitRange[plane][iwire].second;
471  unsigned int jfirsthit = slc.wireHitRange[plane][jwire].first;
472  unsigned int jlasthit = slc.wireHitRange[plane][jwire].second;
473  if (ifirsthit > slc.slHits.size() || ilasthit > slc.slHits.size()) {
474  std::cout << "RAT: bad hit range " << ifirsthit << " " << ilasthit << " size "
475  << slc.slHits.size() << " inCTP " << inCTP << "\n";
476  return;
477  }
478  for (unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
479  tcc.dbgStp = (tcc.modes[kDebug] && (slc.slHits[iht].allHitsIndex == debug.Hit));
480  if (tcc.dbgStp) {
481  mf::LogVerbatim("TC") << "+++++++ Pass " << pass << " Found debug hit "
482  << slices.size() - 1 << ":" << PrintHit(slc.slHits[iht])
483  << " iht " << iht;
484  }
485  // Only consider hits that are available
486  if (slc.slHits[iht].InTraj != 0) continue;
487  // We hope to make a trajectory point at the hit position of iht in WSE units
488  // with a direction pointing to jht
489  if (slc.slHits[iht].allHitsIndex > (*evt.allHits).size() - 1) {
490  std::cout << "RAT: Bad allHitsIndex\n";
491  continue;
492  }
493  auto& iHit = (*evt.allHits)[slc.slHits[iht].allHitsIndex];
494  if (LongPulseHit(iHit)) continue;
495  unsigned int fromWire = iHit.WireID().Wire;
496  float fromTick = iHit.PeakTime();
497  float iqtot = iHit.Integral();
498  float hitsRMSTick = iHit.RMS();
499  std::vector<unsigned int> iHitsInMultiplet;
500  if (pass == 0) {
501  // only use the hit on the first pass
502  iHitsInMultiplet.resize(1);
503  iHitsInMultiplet[0] = iht;
504  }
505  else {
506  GetHitMultiplet(slc, iht, iHitsInMultiplet, false);
507  if (iHitsInMultiplet.empty()) continue;
508  // ignore very high multiplicities
509  if (iHitsInMultiplet.size() > 4) continue;
510  }
511  if (iHitsInMultiplet.size() > 1) {
512  fromTick = HitsPosTick(slc, iHitsInMultiplet, iqtot, kUnusedHits);
513  hitsRMSTick = HitsRMSTick(slc, iHitsInMultiplet, kUnusedHits);
514  }
515  if (hitsRMSTick == 0) continue;
516  bool fatIHit = (hitsRMSTick > maxHitsRMS);
517  if (tcc.dbgStp)
518  mf::LogVerbatim("TC") << " hit RMS " << iHit.RMS() << " BB Multiplicity "
519  << iHitsInMultiplet.size() << " AveHitRMS[" << plane << "] "
520  << evt.aveHitRMS[plane] << " HitsRMSTick " << hitsRMSTick
521  << " fatIHit " << fatIHit;
522  for (unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
523  // Only consider hits that are available
524  if (slc.slHits[iht].InTraj != 0) break;
525  if (slc.slHits[jht].InTraj != 0) continue;
526  // clear out any leftover work inTraj's that weren't cleaned up properly
527  for (auto& slHit : slc.slHits) {
528  if (slHit.InTraj < 0) {
529  std::cout << "RAT: Dirty hit " << PrintHit(slHit) << " EventsProcessed "
530  << evt.eventsProcessed << " WorkID " << evt.WorkID << "\n";
531  slHit.InTraj = 0;
532  }
533  }
534  unsigned int toWire = jwire;
535  auto& jHit = (*evt.allHits)[slc.slHits[jht].allHitsIndex];
536  if (LongPulseHit(jHit)) continue;
537  float toTick = jHit.PeakTime();
538  float jqtot = jHit.Integral();
539  std::vector<unsigned int> jHitsInMultiplet;
540  if (pass == 0) {
541  // only use the hit on the first pass
542  jHitsInMultiplet.resize(1);
543  jHitsInMultiplet[0] = jht;
544  }
545  else {
546  GetHitMultiplet(slc, jht, jHitsInMultiplet, false);
547  if (jHitsInMultiplet.empty()) continue;
548  // ignore very high multiplicities
549  if (jHitsInMultiplet.size() > 4) continue;
550  }
551  hitsRMSTick = HitsRMSTick(slc, jHitsInMultiplet, kUnusedHits);
552  if (hitsRMSTick == 0) continue;
553  bool fatJHit = (hitsRMSTick > maxHitsRMS);
554  if (pass == 0) {
555  // require both hits to be consistent
556  if ((fatIHit && !fatJHit) || (!fatIHit && fatJHit)) { continue; }
557  }
558  else {
559  // pass > 0
560  if (jHitsInMultiplet.size() > 1)
561  toTick = HitsPosTick(slc, jHitsInMultiplet, jqtot, kUnusedHits);
562  }
563  bool hitsOK = TrajHitsOK(slc, iHitsInMultiplet, jHitsInMultiplet);
564  // Ensure that the hits StartTick and EndTick have the proper overlap
565  if (!hitsOK) continue;
566  // start a trajectory with direction from iht -> jht
567  Trajectory work;
568  if (!StartTraj(work, fromWire, fromTick, toWire, toTick, inCTP, pass)) continue;
569  // check for a major failure
570  if (!slc.isValid) {
571  std::cout << "RAT: StartTraj major failure\n";
572  return;
573  }
574  if (work.Pts.empty()) {
575  if (tcc.dbgStp) mf::LogVerbatim("TC") << "ReconstructAllTraj: StartTraj failed";
576  continue;
577  }
578  work.Pts[0].DeltaRMS = tcc.hitErrFac * tcc.unitsPerTick * hitsRMSTick;
579  // don't include the charge of iht since it will be low if this
580  // is a starting/ending track
581  work.AveChg = jqtot;
582  // try to add close hits
583  bool sigOK;
584  AddHits(slc, work, 0, sigOK);
585  // check for a major failure
586  if (!slc.isValid) {
587  std::cout << "RAT: AddHits major failure\n";
588  return;
589  }
590  if (!sigOK || work.Pts[0].Chg == 0) {
591  if (tcc.dbgStp) mf::LogVerbatim("TC") << " No hits at initial trajectory point ";
592  ReleaseHits(slc, work);
593  continue;
594  }
595  // move the TP position to the hit position but don't mess with the direction
596  work.Pts[0].Pos = work.Pts[0].HitPos;
597  // print the header and the first TP
598  if (tcc.dbgStp) PrintTrajectory("RAT", slc, work, USHRT_MAX);
599  // We can't update the trajectory yet because there is only one TP.
600  work.EndPt[0] = 0;
601  // now try stepping away
602  StepAway(slc, work);
603  // check for a major failure
604  if (!slc.isValid) {
605  std::cout << "RAT: StepAway major failure\n";
606  return;
607  }
608  if (tcc.dbgStp)
609  mf::LogVerbatim("TC") << " After first StepAway. IsGood " << work.IsGood;
610  // Check the quality of the work trajectory
611  CheckTraj(slc, work);
612  // check for a major failure
613  if (!slc.isValid) {
614  std::cout << "RAT: CheckTraj major failure\n";
615  return;
616  }
617  if (tcc.dbgStp)
618  mf::LogVerbatim("TC") << "ReconstructAllTraj: After CheckTraj EndPt " << work.EndPt[0]
619  << "-" << work.EndPt[1] << " IsGood " << work.IsGood;
620  if (tcc.dbgStp)
621  mf::LogVerbatim("TC")
622  << "StepAway done: IsGood " << work.IsGood << " NumPtsWithCharge "
623  << NumPtsWithCharge(slc, work, true) << " cut " << tcc.minPts[work.Pass];
624  // decide if the trajectory is long enough for this pass
625  if (!work.IsGood || NumPtsWithCharge(slc, work, true) < tcc.minPts[work.Pass]) {
626  if (tcc.dbgStp)
627  mf::LogVerbatim("TC")
628  << " xxxxxxx Not enough points " << NumPtsWithCharge(slc, work, true)
629  << " minimum " << tcc.minPts[work.Pass] << " or !IsGood";
630  ReleaseHits(slc, work);
631  continue;
632  }
633  if (!StoreTraj(slc, work)) continue;
634  if (tcc.dbgStp) {
635  auto& tj = slc.tjs[slc.tjs.size() - 1];
636  PrintTrajectory("RAT", slc, tj, USHRT_MAX);
637  if (!InTrajOK(slc, "RAT")) {
638  std::cout << "RAT: InTrajOK major failure. " << tj.ID << "\n";
639  return;
640  }
641  } // dbgStp
642  CheckTrajBeginChg(slc, slc.tjs.size() - 1);
643  BraggSplit(slc, slc.tjs.size() - 1);
644  break;
645  } // jht
646  } // iht
647  } // iwire
648 
649  // See if there are any seed trajectory points that were saved before reverse
650  // propagation and try to make Tjs from them
651  for (auto tp : seeds) {
652  unsigned short nAvailable = 0;
653  for (unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
654  if (!tp.UseHit[ii]) continue;
655  unsigned int iht = tp.Hits[ii];
656  if (slc.slHits[iht].InTraj == 0) ++nAvailable;
657  tcc.dbgStp = (tcc.modes[kDebug] && (slc.slHits[iht].allHitsIndex == debug.Hit));
658  if (tcc.dbgStp) {
659  mf::LogVerbatim("TC") << "+++++++ Seed debug hit " << slices.size() - 1 << ":"
660  << PrintHit(slc.slHits[iht]) << " iht " << iht;
661  }
662  } // ii
663  if (nAvailable == 0) continue;
664  Trajectory work;
665  work.ID = evt.WorkID;
666  for (unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
667  if (!tp.UseHit[ii]) continue;
668  unsigned int iht = tp.Hits[ii];
669  if (slc.slHits[iht].InTraj == 0) slc.slHits[iht].InTraj = work.ID;
670  } // ii
671  work.Pass = pass;
672  work.StepDir = 1;
673  if (tp.Dir[0] < 0) work.StepDir = -1;
674  work.CTP = tp.CTP;
675  work.ParentID = -1;
676  work.Strategy.reset();
677  work.Strategy[kNormal] = true;
678  // don't allow yet another reverse propagation
679  work.AlgMod[kRvPrp] = true;
680  work.Pts.push_back(tp);
681  StepAway(slc, work);
682  CheckTraj(slc, work);
683  // check for a major failure
684  if (!slc.isValid) {
685  std::cout << "RAT: CheckTraj major failure\n";
686  return;
687  }
688  // decide if the trajectory is long enough for this pass
689  if (!work.IsGood || NumPtsWithCharge(slc, work, true) < tcc.minPts[work.Pass]) {
690  if (tcc.dbgStp)
691  mf::LogVerbatim("TC") << " xxxxxxx Not enough points "
692  << NumPtsWithCharge(slc, work, true) << " minimum "
693  << tcc.minPts[work.Pass] << " or !IsGood";
694  ReleaseHits(slc, work);
695  continue;
696  }
697  if (!StoreTraj(slc, work)) continue;
698  if (tcc.dbgStp) {
699  auto& tj = slc.tjs[slc.tjs.size() - 1];
700  mf::LogVerbatim("TC") << "TRP RAT Stored T" << tj.ID << " using seed TP " << PrintPos(tp);
701  }
702  BraggSplit(slc, slc.tjs.size() - 1);
703  } // seed
704 
705  seeds.resize(0);
706 
707  bool lastPass = (pass == tcc.minPtsFit.size() - 1);
708  // don't use lastPass cuts if we will use LastEndMerge
709  if (tcc.useAlg[kLastEndMerge]) lastPass = false;
710  EndMerge(slc, inCTP, lastPass);
711  if (!slc.isValid) return;
712 
713  Find2DVertices(detProp, slc, inCTP, pass);
714 
715  } // pass
716 
717  // Last attempt to merge long straight Tjs that failed the EndMerge cuts
718  LastEndMerge(slc, inCTP);
719 
720  // make junk trajectories using nearby un-assigned hits
721  FindJunkTraj(slc, inCTP);
722  // dressed muons with halo trajectories
723  if (tcc.muonTag.size() > 4 && tcc.muonTag[4] > 0) {
724  for (auto& tj : slc.tjs) {
725  if (tj.AlgMod[kKilled]) continue;
726  if (tj.PDGCode != 13) continue;
727  MakeHaloTj(slc, tj, tcc.dbgSlc);
728  } // tj
729  } // dressed muons
730 
731  // Tag ShowerLike Tjs
732  if (tcc.showerTag[0] > 0) TagShowerLike(slc, inCTP);
733  // Set TP Environment bits
734  SetTPEnvironment(slc, inCTP);
735  Find2DVertices(detProp, slc, inCTP, USHRT_MAX);
736  SplitTrajCrossingVertices(slc, inCTP);
737  // Make vertices between long Tjs and junk Tjs
738  MakeJunkVertices(slc, inCTP);
739  // check for a major failure
740  if (!slc.isValid) {
741  std::cout << "RAT: MakeJunkVertices major failure\n";
742  return;
743  }
744 
745  // last attempt to attach Tjs to vertices
746  for (unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
747  auto& vx2 = slc.vtxs[ivx];
748  if (vx2.ID == 0) continue;
749  if (vx2.CTP != inCTP) continue;
750  AttachAnyTrajToVertex(slc, ivx, tcc.dbgStp || tcc.dbg2V);
751  } // ivx
752 
753  // Set the kEnvOverlap bit true for all TPs that are close to other
754  // trajectories that are close to vertices. The positions of TPs that
755  // overlap are biased and shouldn't be used in a vertex fit. Also, these
756  // TPs shouldn't be used to calculate dE/dx
757  UpdateVxEnvironment(slc);
758 
759  // Check the Tj <-> vtx associations and define the vertex quality
760  if (!ChkVtxAssociations(slc, inCTP)) {
761  std::cout << "RAT: ChkVtxAssociations found an error. Events processed "
762  << evt.eventsProcessed << " WorkID " << evt.WorkID << "\n";
763  }
764 
765  } // ReconstructAllTraj
766 
769  {
770  // Makes junk trajectories using unassigned hits
771 
772  if (tcc.JTMaxHitSep2 <= 0) return;
773  if (!tcc.useAlg[kJunkTj]) return;
774  unsigned short plane = DecodeCTP(inCTP).Plane;
775  if ((int)slc.lastWire[plane] - 3 < (int)slc.firstWire[plane]) return;
776 
777  // shouldn't have to do this but...
778  for (auto& slHit : slc.slHits)
779  if (slHit.InTraj < 0) slHit.InTraj = 0;
780 
781  bool prt = false;
782 
783  std::vector<unsigned int> tHits;
784  // Stay well away from the last wire in the plane
785  for (unsigned int iwire = slc.firstWire[plane]; iwire < slc.lastWire[plane] - 3; ++iwire) {
786  // skip bad wires or no hits on the wire
787  if (slc.wireHitRange[plane][iwire].first > slc.slHits.size()) continue;
788  if (slc.wireHitRange[plane][iwire].second > slc.slHits.size()) continue;
789  unsigned int jwire = iwire + 1;
790  if (slc.wireHitRange[plane][jwire].first > slc.slHits.size()) continue;
791  if (slc.wireHitRange[plane][jwire].second > slc.slHits.size()) continue;
792  unsigned int ifirsthit = slc.wireHitRange[plane][iwire].first;
793  unsigned int ilasthit = slc.wireHitRange[plane][iwire].second;
794  unsigned int jfirsthit = slc.wireHitRange[plane][jwire].first;
795  unsigned int jlasthit = slc.wireHitRange[plane][jwire].second;
796  for (unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
797  if (iht >= slc.slHits.size()) break;
798  auto& islHit = slc.slHits[iht];
799  if (islHit.InTraj != 0) continue;
800  std::vector<unsigned int> iHits;
801  GetHitMultiplet(slc, iht, iHits, true);
802  prt =
803  (tcc.modes[kDebug] && std::find(iHits.begin(), iHits.end(), debug.Hit) != iHits.end());
804  if (prt) mf::LogVerbatim("TC") << "FJT: debug iht multiplet size " << iHits.size();
805  if (iHits.empty()) continue;
806  for (unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
807  if (jht >= slc.slHits.size()) break;
808  auto& jslHit = slc.slHits[jht];
809  if (jslHit.InTraj != 0) continue;
810  if (prt && HitSep2(slc, iht, jht) < 100)
811  mf::LogVerbatim("TC") << " use " << PrintHit(jslHit) << " hitSep2 "
812  << HitSep2(slc, iht, jht);
813  if (HitSep2(slc, iht, jht) > tcc.JTMaxHitSep2) continue;
814  std::vector<unsigned int> jHits;
815  GetHitMultiplet(slc, jht, jHits, true);
816  if (jHits.empty()) continue;
817  // check for hit overlap consistency
818  if (!TrajHitsOK(slc, iHits, jHits)) continue;
819  tHits.clear();
820  // add the available hits and flag them
821  for (auto iht : iHits)
822  if (slc.slHits[iht].InTraj == 0) tHits.push_back(iht);
823  for (auto jht : jHits)
824  if (slc.slHits[jht].InTraj == 0) tHits.push_back(jht);
825  for (auto tht : tHits)
826  slc.slHits[tht].InTraj = -4;
827  unsigned int loWire;
828  if (iwire != 0) { loWire = iwire - 1; }
829  else {
830  loWire = 0;
831  }
832  unsigned int hiWire = jwire + 1;
833  if (hiWire > slc.nWires[plane]) break;
834  unsigned short nit = 0;
835  while (nit < 100) {
836  bool hitsAdded = false;
837  for (unsigned int kwire = loWire; kwire <= hiWire; ++kwire) {
838  if (slc.wireHitRange[plane][kwire].first == UINT_MAX) continue;
839  if (slc.wireHitRange[plane][kwire].second == UINT_MAX) continue;
840  unsigned int kfirsthit = slc.wireHitRange[plane][kwire].first;
841  unsigned int klasthit = slc.wireHitRange[plane][kwire].second;
842  for (unsigned int kht = kfirsthit; kht <= klasthit; ++kht) {
843  if (kht >= slc.slHits.size()) continue;
844  if (slc.slHits[kht].InTraj != 0) continue;
845  // this shouldn't be needed but do it anyway
846  if (std::find(tHits.begin(), tHits.end(), kht) != tHits.end()) continue;
847  // re-purpose jHits and check for consistency
848  GetHitMultiplet(slc, kht, jHits, true);
849  if (!TrajHitsOK(slc, tHits, jHits)) continue;
850  // add them all and update the wire range
851  for (auto jht : jHits) {
852  if (slc.slHits[jht].InTraj != 0) continue;
853  tHits.push_back(jht);
854  slc.slHits[jht].InTraj = -4;
855  if (kwire > hiWire) hiWire = kwire;
856  if (kwire < loWire) loWire = kwire;
857  hitsAdded = true;
858  } // jht
859  } // kht
860  } // kwire
861  if (!hitsAdded) break;
862  ++nit;
863  ++hiWire;
864  if (hiWire >= slc.nWires[plane]) break;
865  } // nit < 100
866  // clear InTraj
867  for (auto iht : tHits)
868  slc.slHits[iht].InTraj = 0;
869  if (tHits.size() < 3) continue;
870  if (prt) {
871  mf::LogVerbatim myprt("TC");
872  myprt << "FJT: tHits";
873  for (auto tht : tHits)
874  myprt << " " << PrintHit(slc.slHits[tht]);
875  myprt << "\n";
876  } // prt
877  // See if this is a ghost trajectory
878  if (IsGhost(slc, tHits)) break;
879  if (!MakeJunkTraj(slc, tHits)) {
880  if (prt) mf::LogVerbatim() << "FJT: MakeJunkTraj failed";
881  break;
882  }
883  if (slc.slHits[jht].InTraj > 0) break;
884  } // jht
885  } // iht
886  } // iwire
887  } // FindJunkTraj
888 
890  void TrajClusterAlg::ChkInTraj(std::string someText, TCSlice& slc)
891  {
892  // Check slc.tjs -> InTraj associations
893 
894  if (!tcc.useAlg[kChkInTraj]) return;
895 
897 
898  int tID;
899  unsigned int iht;
900  unsigned short itj = 0;
901  std::vector<unsigned int> tHits;
902  std::vector<unsigned int> atHits;
903  for (auto& tj : slc.tjs) {
904  // ignore abandoned trajectories
905  if (tj.AlgMod[kKilled]) continue;
906  tID = tj.ID;
907  for (auto& tp : tj.Pts) {
908  if (tp.Hits.size() > 16) {
909  tj.AlgMod[kKilled] = true;
910  mf::LogVerbatim("TC")
911  << "ChkInTraj: More than 16 hits created a UseHit bitset overflow\n";
912  slc.isValid = false;
913  std::cout << "ChkInTraj major failure\n";
914  return;
915  }
916  } // tp
917  if (tj.AlgMod[kKilled]) {
918  std::cout << someText << " ChkInTraj hit size mis-match in tj ID " << tj.ID
919  << " AlgBitNames";
920  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
921  if (tj.AlgMod[ib]) std::cout << " " << AlgBitNames[ib];
922  std::cout << "\n";
923  continue;
924  }
925  tHits = PutTrajHitsInVector(tj, kUsedHits);
926  if (tHits.size() < 2) {
927  mf::LogVerbatim("TC") << someText << " ChkInTraj: Insufficient hits in traj " << tj.ID;
928  PrintTrajectory("CIT", slc, tj, USHRT_MAX);
929  tj.AlgMod[kKilled] = true;
930  continue;
931  }
932  std::sort(tHits.begin(), tHits.end());
933  atHits.clear();
934  for (iht = 0; iht < slc.slHits.size(); ++iht) {
935  if (slc.slHits[iht].InTraj == tID) atHits.push_back(iht);
936  } // iht
937  if (atHits.size() < 2) {
938  mf::LogVerbatim("TC") << someText << " ChkInTraj: Insufficient hits in atHits in traj "
939  << tj.ID << " Killing it";
940  tj.AlgMod[kKilled] = true;
941  continue;
942  }
943  if (!std::equal(tHits.begin(), tHits.end(), atHits.begin())) {
944  mf::LogVerbatim myprt("TC");
945  myprt << someText << " ChkInTraj failed: inTraj - UseHit mis-match for tj ID " << tID
946  << " tj.WorkID " << tj.WorkID << " atHits size " << atHits.size() << " tHits size "
947  << tHits.size() << " in CTP " << tj.CTP << "\n";
948  myprt << "AlgMods: ";
949  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
950  if (tj.AlgMod[ib]) myprt << " " << AlgBitNames[ib];
951  myprt << "\n";
952  myprt << "index inTraj UseHit \n";
953  for (iht = 0; iht < atHits.size(); ++iht) {
954  myprt << "iht " << iht << " " << PrintHit(slc.slHits[atHits[iht]]);
955  if (iht < tHits.size()) myprt << " " << PrintHit(slc.slHits[tHits[iht]]);
956  if (atHits[iht] != tHits[iht]) myprt << " <<< " << atHits[iht] << " != " << tHits[iht];
957  myprt << "\n";
958  slc.isValid = false;
959  } // iht
960  if (tHits.size() > atHits.size()) {
961  for (iht = atHits.size(); iht < atHits.size(); ++iht) {
962  myprt << "atHits " << iht << " " << PrintHit(slc.slHits[atHits[iht]]) << "\n";
963  } // iht
964  PrintTrajectory("CIT", slc, tj, USHRT_MAX);
965  } // tHit.size > atHits.size()
966  }
967  // check the VtxID
968  for (unsigned short end = 0; end < 2; ++end) {
969  if (tj.VtxID[end] > slc.vtxs.size()) {
970  mf::LogVerbatim("TC") << someText << " ChkInTraj: Bad VtxID " << tj.ID;
971  std::cout << someText << " ChkInTraj: Bad VtxID " << tj.ID << " vtx size "
972  << slc.vtxs.size() << "\n";
973  tj.AlgMod[kKilled] = true;
974  PrintTrajectory("CIT", slc, tj, USHRT_MAX);
975  slc.isValid = false;
976  return;
977  }
978  } // end
979  ++itj;
980  if (!slc.isValid) return;
981  } // tj
982 
983  } // ChkInTraj
984 
986  void TrajClusterAlg::MergeTPHits(std::vector<unsigned int>& tpHits,
987  std::vector<recob::Hit>& newHitCol,
988  std::vector<unsigned int>& newHitAssns) const
989  {
990  // merge the hits indexed by tpHits into one or more hits with the requirement that the hits
991  // are on different wires
992 
993  if (tpHits.empty()) return;
994 
995  // no merge required. Just put a close copy of the single hit in the output hit collection
996  if (tpHits.size() == 1) {
997  if (tpHits[0] > (*evt.allHits).size() - 1) {
998  std::cout << "MergeTPHits Bad input hit index " << tpHits[0] << " allHits size "
999  << (*evt.allHits).size() << "\n";
1000  return;
1001  }
1002  newHitCol.push_back(MergeTPHitsOnWire(tpHits));
1003  newHitAssns[tpHits[0]] = newHitCol.size() - 1;
1004  return;
1005  } // tpHits.size() == 1
1006 
1007  // split the hit list into sub-lists of hits on a single wire
1008  std::vector<unsigned int> wires;
1009  std::vector<std::vector<unsigned int>> wireHits;
1010  auto& firstHit = (*evt.allHits)[tpHits[0]];
1011  wires.push_back(firstHit.WireID().Wire);
1012  std::vector<unsigned int> tmp(1, tpHits[0]);
1013  wireHits.push_back(tmp);
1014  for (unsigned short ii = 1; ii < tpHits.size(); ++ii) {
1015  auto& hit = (*evt.allHits)[tpHits[ii]];
1016  unsigned int wire = hit.WireID().Wire;
1017  unsigned short indx = 0;
1018  for (indx = 0; indx < wires.size(); ++indx)
1019  if (wires[indx] == wire) break;
1020  if (indx == wires.size()) {
1021  wires.push_back(wire);
1022  wireHits.resize(wireHits.size() + 1);
1023  }
1024  wireHits[indx].push_back(tpHits[ii]);
1025  } // ii
1026 
1027  // now merge hits in each sub-list.
1028  for (unsigned short indx = 0; indx < wireHits.size(); ++indx) {
1029  auto& hitsOnWire = wireHits[indx];
1030  newHitCol.push_back(MergeTPHitsOnWire(hitsOnWire));
1031  for (unsigned short ii = 0; ii < hitsOnWire.size(); ++ii) {
1032  newHitAssns[hitsOnWire[ii]] = newHitCol.size() - 1;
1033  }
1034  } // hitsOnWire
1035 
1036  return;
1037 
1038  } // MergeTPHits
1039 
1041  recob::Hit TrajClusterAlg::MergeTPHitsOnWire(std::vector<unsigned int>& tpHits) const
1042  {
1043  // merge the hits indexed by tpHits into one hit
1044 
1045  if (tpHits.empty()) return recob::Hit();
1046 
1047  // no merge required. Just return a slightly modified copy of the single hit
1048  if (tpHits.size() == 1) {
1049  if (tpHits[0] > (*evt.allHits).size() - 1) {
1050  std::cout << "MergeTPHits Bad input hit index " << tpHits[0] << " allHits size "
1051  << (*evt.allHits).size() << "\n";
1052  return recob::Hit();
1053  }
1054  auto& oldHit = (*evt.allHits)[tpHits[0]];
1055  raw::TDCtick_t startTick = oldHit.PeakTime() - 3 * oldHit.RMS();
1056  raw::TDCtick_t endTick = oldHit.PeakTime() + 3 * oldHit.RMS();
1057 
1058  return recob::Hit(oldHit.Channel(),
1059  startTick,
1060  endTick,
1061  oldHit.PeakTime(),
1062  oldHit.SigmaPeakTime(),
1063  oldHit.RMS(),
1064  oldHit.PeakAmplitude(),
1065  oldHit.SigmaPeakAmplitude(),
1066  oldHit.ROISummedADC(),
1067  oldHit.HitSummedADC(),
1068  oldHit.Integral(),
1069  oldHit.SigmaIntegral(),
1070  1,
1071  0, // Multiplicity, LocalIndex
1072  1,
1073  0, // GoodnessOfFit, DOF
1074  oldHit.View(),
1075  oldHit.SignalType(),
1076  oldHit.WireID());
1077  } // tpHits.size() == 1
1078 
1079  double integral = 0;
1080  double sIntegral = 0;
1081  double peakTime = 0;
1082  double sPeakTime = 0;
1083  double peakAmp = 0;
1084  double sPeakAmp = 0;
1085  float ROIsumADC = 0;
1086  float HitsumADC = 0;
1087  raw::TDCtick_t startTick = INT_MAX;
1088  raw::TDCtick_t endTick = 0;
1089  for (auto allHitsIndex : tpHits) {
1090  if (allHitsIndex > (*evt.allHits).size() - 1) return recob::Hit();
1091  auto& hit = (*evt.allHits)[allHitsIndex];
1092  if (hit.StartTick() < startTick) startTick = hit.StartTick();
1093  if (hit.EndTick() > endTick) endTick = hit.EndTick();
1094  double intgrl = hit.Integral();
1095  sPeakTime += intgrl * hit.SigmaPeakTime();
1096  sPeakAmp += intgrl * hit.SigmaPeakAmplitude();
1097  ROIsumADC += hit.ROISummedADC();
1098  HitsumADC += hit.HitSummedADC();
1099  integral += intgrl;
1100  sIntegral += intgrl * hit.SigmaIntegral();
1101  // Get the charge normalization from an input hit
1102  } // tpHit
1103  if (integral <= 0) {
1104  std::cout << "MergeTPHits found bad hit integral " << integral << "\n";
1105  return recob::Hit();
1106  }
1107 
1108  // Create a signal shape vector to find the rms and peak tick
1109  std::vector<double> shape(endTick - startTick + 1, 0.);
1110  for (auto allHitsIndex : tpHits) {
1111  auto& hit = (*evt.allHits)[allHitsIndex];
1112  double peakTick = hit.PeakTime();
1113  double rms = hit.RMS();
1114  double peakAmp = hit.PeakAmplitude();
1115  for (raw::TDCtick_t tick = startTick; tick <= endTick; ++tick) {
1116  double arg = ((double)tick - peakTick) / rms;
1117  unsigned short indx = tick - startTick;
1118  shape[indx] += peakAmp * exp(-0.5 * arg * arg);
1119  } // tick
1120  } // allHitsIndex
1121 
1122  // find the peak tick
1123  double sigsum = 0;
1124  double sigsumt = 0;
1125  for (raw::TDCtick_t tick = startTick; tick <= endTick; ++tick) {
1126  unsigned short indx = tick - startTick;
1127  sigsum += shape[indx];
1128  sigsumt += shape[indx] * tick;
1129  } // tick
1130 
1131  peakTime = sigsumt / sigsum;
1132  // Use the sigma peak time calculated in the first loop
1133  sPeakTime /= integral;
1134 
1135  // find the RMS
1136  sigsumt = 0.;
1137  for (raw::TDCtick_t tick = startTick; tick <= endTick; ++tick) {
1138  double dTick = tick - peakTime;
1139  unsigned short indx = tick - startTick;
1140  sigsumt += shape[indx] * dTick * dTick;
1141  }
1142  double rms = std::sqrt(sigsumt / sigsum);
1143  // get a reference to the first hit to get the charge normalization, channel, view, etc
1144  auto& firstHit = (*evt.allHits)[tpHits[0]];
1145  double chgNorm = 2.507 * firstHit.PeakAmplitude() * firstHit.RMS() / firstHit.Integral();
1146  // find the amplitude from the integrated charge and the RMS
1147  peakAmp = (float)(integral * chgNorm / (2.507 * rms));
1148  // Use the sigma integral calculated in the first loop
1149  sPeakAmp /= integral;
1150  // reset the start and end tick
1151  startTick = peakTime - 3 * rms;
1152  endTick = peakTime + 3 * rms;
1153 
1154  // construct the hit
1155  return recob::Hit(firstHit.Channel(),
1156  startTick,
1157  endTick,
1158  peakTime,
1159  sPeakTime,
1160  rms,
1161  peakAmp,
1162  sPeakAmp,
1163  ROIsumADC,
1164  HitsumADC,
1165  integral,
1166  sIntegral,
1167  1,
1168  0, // Multiplicity, LocalIndex
1169  1,
1170  0, // GoodnessOfFit, DOF
1171  firstHit.View(),
1172  firstHit.SignalType(),
1173  firstHit.WireID());
1174 
1175  } // MergeTPHits
1176 
1179  {
1180  showertree = t;
1181 
1182  showertree->Branch("run", &evt.run, "run/I");
1183  showertree->Branch("subrun", &evt.subRun, "subrun/I");
1184  showertree->Branch("event", &evt.event, "event/I");
1185 
1186  showertree->Branch("BeginWir", &stv.BeginWir);
1187  showertree->Branch("BeginTim", &stv.BeginTim);
1188  showertree->Branch("BeginAng", &stv.BeginAng);
1189  showertree->Branch("BeginChg", &stv.BeginChg);
1190  showertree->Branch("BeginVtx", &stv.BeginVtx);
1191 
1192  showertree->Branch("EndWir", &stv.EndWir);
1193  showertree->Branch("EndTim", &stv.EndTim);
1194  showertree->Branch("EndAng", &stv.EndAng);
1195  showertree->Branch("EndChg", &stv.EndChg);
1196  showertree->Branch("EndVtx", &stv.EndVtx);
1197 
1198  showertree->Branch("MCSMom", &stv.MCSMom);
1199 
1200  showertree->Branch("PlaneNum", &stv.PlaneNum);
1201  showertree->Branch("TjID", &stv.TjID);
1202  showertree->Branch("IsShowerTj", &stv.IsShowerTj);
1203  showertree->Branch("ShowerID", &stv.ShowerID);
1204  showertree->Branch("IsShowerParent", &stv.IsShowerParent);
1205  showertree->Branch("StageNum", &stv.StageNum);
1206  showertree->Branch("StageName", &stv.StageName);
1207 
1208  showertree->Branch("Envelope", &stv.Envelope);
1209  showertree->Branch("EnvPlane", &stv.EnvPlane);
1210  showertree->Branch("EnvStage", &stv.EnvStage);
1211  showertree->Branch("EnvShowerID", &stv.EnvShowerID);
1212 
1213  showertree->Branch("nStages", &stv.nStages);
1214  showertree->Branch("nPlanes", &stv.nPlanes);
1215 
1216  } // end DefineShTree
1217 
1220  detinfo::DetectorPropertiesData const& detProp,
1221  std::vector<unsigned int>& hitsInSlice,
1222  int sliceID)
1223  {
1224  // Defines a TCSlice struct and pushes the slice onto slices.
1225  // Sets the isValid flag true if successful.
1226  if ((*evt.allHits).empty()) return false;
1227  if (hitsInSlice.size() < 2) return false;
1228 
1229  TCSlice slc;
1230  slc.ID = sliceID;
1231  slc.slHits.resize(hitsInSlice.size());
1232  bool first = true;
1233  unsigned int cstat = 0;
1234  unsigned int tpc = UINT_MAX;
1235  unsigned int cnt = 0;
1236  std::vector<unsigned int> nHitsInPln;
1237  for (auto iht : hitsInSlice) {
1238  if (iht >= (*evt.allHits).size()) return false;
1239  auto& hit = (*evt.allHits)[iht];
1240  if (first) {
1241  cstat = hit.WireID().Cryostat;
1242  tpc = hit.WireID().TPC;
1243  slc.TPCID = geo::TPCID(cstat, tpc);
1244  nHitsInPln.resize(tcc.wireReadoutGeom->Nplanes(slc.TPCID));
1245  first = false;
1246  }
1247  if (hit.WireID().Cryostat != cstat || hit.WireID().TPC != tpc) return false;
1248  slc.slHits[cnt].allHitsIndex = iht;
1249  ++nHitsInPln[hit.WireID().Plane];
1250  ++cnt;
1251  } // iht
1252  // require at least two hits in each plane
1253  for (auto hip : nHitsInPln)
1254  if (hip < 2) return false;
1255  // Define the TCEvent wire hit range vector for this new TPC for ALL hits
1256  FillWireHitRange(slc.TPCID);
1257  // next define the Slice wire hit range vectors, UnitsPerTick, etc for this
1258  // slice
1259  if (!FillWireHitRange(clockData, detProp, slc)) return false;
1260  slc.isValid = true;
1261  slices.push_back(slc);
1262  if (tcc.modes[kDebug] && debug.Slice >= 0 && !tcc.dbgSlc) {
1263  tcc.dbgSlc = ((int)(slices.size() - 1) == debug.Slice);
1264  if (tcc.dbgSlc) std::cout << "Enabled debugging in sub-slice " << slices.size() - 1 << "\n";
1265  if (tcc.modes[kDebug] && (unsigned int)debug.Cryostat == cstat &&
1266  (unsigned int)debug.TPC == tpc && debug.Plane >= 0) {
1267  debug.CTP = EncodeCTP(
1268  (unsigned int)debug.Cryostat, (unsigned int)debug.TPC, (unsigned int)debug.Plane);
1269  }
1270  }
1271  // do a sanity check
1272  for (unsigned int iht = 0; iht < slc.slHits.size(); ++iht) {
1273  unsigned int ahi = slc.slHits[iht].allHitsIndex;
1274  if (ahi >= (*evt.allHits).size()) {
1275  std::cout << "CreateSlice: Bad allHitsIndex " << ahi << " " << (*evt.allHits).size()
1276  << "\n";
1277  return false;
1278  }
1279  } // iht
1280  return true;
1281  } // CreateSlice
1282 
1285  {
1286  // final steps that involve correlations between slices
1287  // Stitch PFParticles between TPCs
1288 
1289  // define the PFP TjUIDs vector before calling StitchPFPs
1290  for (auto& slc : slices) {
1291  if (!slc.isValid) continue;
1292  MakePFPTjs(slc);
1293  for (auto& pfp : slc.pfps) {
1294  if (pfp.ID <= 0) continue;
1295  pfp.TjUIDs.resize(pfp.TjIDs.size());
1296  for (unsigned short ii = 0; ii < pfp.TjIDs.size(); ++ii) {
1297  // do a sanity check while we are here
1298  if (pfp.TjIDs[ii] <= 0 || pfp.TjIDs[ii] > (int)slc.tjs.size()) {
1299  std::cout << "FinishEvent found an invalid T" << pfp.TjIDs[ii] << " in P" << pfp.UID
1300  << "\n";
1301  slc.isValid = false;
1302  continue;
1303  } // sanity check
1304  auto& tj = slc.tjs[pfp.TjIDs[ii] - 1];
1305  pfp.TjUIDs[ii] = tj.UID;
1306  } // ii
1307  } // pfp
1308  } // slc
1309 
1310  StitchPFPs();
1311  // Ensure that all PFParticles have a start vertex
1312  for (auto& slc : slices)
1313  PFPVertexCheck(slc);
1314  } // FinishEvent
1315 
1316 } // namespace cluster
Expect tracks entering from the front face. Don&#39;t create neutrino PFParticles.
Definition: DataStructs.h:523
calo::CalorimetryAlg * caloAlg
Definition: DataStructs.h:569
void ClearResults()
Deletes all the results.
void CheckTrajBeginChg(TCSlice &slc, unsigned short itj)
Definition: Utils.cxx:1317
float AveChg
Calculated using ALL hits.
Definition: DataStructs.h:192
std::vector< int > EnvStage
Definition: DataStructs.h:402
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< int > IsShowerParent
Definition: DataStructs.h:395
std::vector< Trajectory > tjs
vector of all trajectories in each plane
Definition: DataStructs.h:662
void StepAway(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:30
void ConfigureMVA(TCConfig &tcc, std::string fMVAShowerParentWeights)
Definition: TCShower.cxx:35
bool TrajHitsOK(TCSlice const &slc, const std::vector< unsigned int > &iHitsInMultiplet, const std::vector< unsigned int > &jHitsInMultiplet)
Definition: Utils.cxx:1839
std::vector< float > kinkCuts
kink finder algorithm
Definition: DataStructs.h:550
bool IsGhost(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:3084
unsigned int event
Definition: DataStructs.h:626
Utilities related to art service access.
bool InTrajOK(TCSlice &slc, std::string someText)
Definition: Utils.cxx:1257
std::vector< float > EndWir
Definition: DataStructs.h:382
std::vector< float > EndAng
Definition: DataStructs.h:384
std::vector< unsigned int > PutTrajHitsInVector(const Trajectory &tj, HitStatus_t hitRequest)
Definition: Utils.cxx:2689
const std::vector< std::string > AlgBitNames
Definition: DataStructs.cxx:16
unsigned int run
Definition: DataStructs.h:627
short recoTPC
only reconstruct in the seleted TPC
Definition: DataStructs.h:581
std::vector< unsigned short > minPtsFit
Reconstruct in several passes.
Definition: DataStructs.h:556
const geo::WireReadoutGeom * wireReadoutGeom
Definition: DataStructs.h:568
std::vector< float > vtx3DCuts
2D vtx -> 3D vtx matching cuts
Definition: DataStructs.h:542
std::vector< float > qualityCuts
Min points/wire, min consecutive pts after a gap.
Definition: DataStructs.h:554
void Find3DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: TCVertex.cxx:1258
std::vector< float > BeginTim
Definition: DataStructs.h:378
void AddHits(TCSlice &slc, Trajectory &tj, unsigned short ipt, bool &sigOK)
Definition: StepUtils.cxx:1113
TCConfig tcc
Definition: DataStructs.cxx:9
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:416
std::vector< float > BeginAng
Definition: DataStructs.h:379
void PrintTrajectory(std::string someText, const TCSlice &slc, const Trajectory &tj, unsigned short tPoint)
Definition: Utils.cxx:6035
void SplitTrajCrossingVertices(TCSlice &slc, CTP_t inCTP)
Definition: TCVertex.cxx:916
bool StartTraj(TCSlice const &slc, Trajectory &tj, unsigned int fromhit, unsigned int tohit, unsigned short pass)
Definition: Utils.cxx:4858
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > wireHitRange
Definition: DataStructs.h:667
void Reconcile2Vs(TCSlice &slc)
Definition: TCVertex.cxx:1059
int ParentID
ID of the parent, or the ID of the Tj this one was merged with if it is killed.
Definition: DataStructs.h:190
step from US -> DS (true) or DS -> US (false)
Definition: DataStructs.h:522
bool StoreTraj(TCSlice &slc, Trajectory &tj)
Definition: Utils.cxx:1074
void RunTrajClusterAlg(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
calo::CalorimetryAlg fCaloAlg
float maxWireSkipWithSignal
max number of wires to skip with a signal on them
Definition: DataStructs.h:574
call study functions to develop cuts, etc
Definition: DataStructs.h:527
for(Int_t i=0;i< nentries;i++)
Definition: comparison.C:30
void FillWireHitRange(geo::TPCID inTPCID)
Definition: Utils.cxx:4325
bool AttachAnyTrajToVertex(TCSlice &slc, unsigned short ivx, bool prt)
Definition: TCVertex.cxx:1678
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
std::vector< std::pair< unsigned int, unsigned int > > tpcSrcHitRange
Definition: DataStructs.h:620
Float_t tmp
Definition: plot.C:35
bool MakeJunkTraj(TCSlice &slc, std::vector< unsigned int > tHits)
Definition: StepUtils.cxx:4146
std::vector< float > electronTag
Definition: DataStructs.h:547
bool FindShowers3D(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: TCShower.cxx:286
bool dbgSlc
debug only in the user-defined slice? default is all slices
Definition: DataStructs.h:582
void TagShowerLike(TCSlice &slc, const CTP_t &inCTP)
Definition: TCShower.cxx:3232
std::vector< unsigned int > lastWire
the last wire with a hit
Definition: DataStructs.h:647
bool CreateSlice(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
bool LongPulseHit(const recob::Hit &hit)
Definition: Utils.cxx:4317
unsigned int eventsProcessed
Definition: DataStructs.h:629
bool IsGood
set false if there is a failure or the Tj fails quality cuts
Definition: DataStructs.h:210
std::string PrintPos(const TrajPoint &tp)
Definition: Utils.cxx:6356
bool doForecast
See TCMode_t above.
Definition: DataStructs.h:600
void MakePFPTjs(TCSlice &slc)
Definition: PFPUtils.cxx:508
std::vector< float > angleRanges
list of max angles for each angle range
Definition: DataStructs.h:560
const std::vector< std::string > EndFlagNames
Definition: DataStructs.cxx:87
std::vector< float > showerTag
shower-like trajectory tagging + shower reconstruction
Definition: DataStructs.h:549
unsigned short Pass
the pass on which it was created
Definition: DataStructs.h:204
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
std::vector< float > EndTim
Definition: DataStructs.h:383
std::vector< int > ShowerID
Definition: DataStructs.h:394
ShowerTreeVars stv
Definition: DataStructs.cxx:11
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::vector< unsigned short > maxAngleCode
max allowed angle code for each pass
Definition: DataStructs.h:558
bool SetInputHits(std::vector< recob::Hit > const &inputHits, unsigned int run, unsigned int event)
void DefineShTree(TTree *t)
std::vector< unsigned int > fAlgModCount
float projectionErrFactor
Definition: DataStructs.h:575
unsigned int Hit
set to the hit index in evt.allHits if a Plane:Wire:Tick match is found
Definition: DebugStruct.h:26
save shower tree
Definition: DataStructs.h:531
std::vector< std::string > StageName
Definition: DataStructs.h:397
float HitSep2(const TCSlice &slc, unsigned int iht, unsigned int jht)
Definition: Utils.cxx:2483
int Cryostat
Select Cryostat.
Definition: DebugStruct.h:20
float HitsRMSTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, HitStatus_t hitRequest)
Definition: Utils.cxx:4133
short nPtsAve
dump trajectory points
Definition: DataStructs.h:598
if(nlines<=0)
std::vector< int > TjID
Definition: DataStructs.h:392
void PFPVertexCheck(TCSlice &slc)
Definition: PFPUtils.cxx:2822
bool dbgStp
debug stepping using debug.Cryostat, debug.TPC, etc
Definition: DataStructs.h:583
recob::Hit MergeTPHitsOnWire(std::vector< unsigned int > &tpHits) const
don&#39;t mess with this line
Definition: DataStructs.h:506
float unitsPerTick
scale factor from Tick to WSE equivalent units
Definition: DataStructs.h:562
std::vector< short > minMCSMom
Min MCSMom for each pass.
Definition: DataStructs.h:559
bool BraggSplit(TCSlice &slc, unsigned short itj)
Definition: Utils.cxx:1423
std::vector< short > BeginVtx
Definition: DataStructs.h:381
DebugStuff debug
Definition: DebugStruct.cxx:4
CTP_t CTP
Cryostat, TPC, Plane code.
Definition: DataStructs.h:187
void DumpTj()
Definition: Utils.cxx:5251
bool dbg2V
debug 2D vertex finding
Definition: DataStructs.h:585
parameter set interface
void GetHitMultiplet(const TCSlice &slc, unsigned int theHit, std::vector< unsigned int > &hitsInMultiplet, bool useLongPulseHits)
Definition: StepUtils.cxx:1521
std::vector< float > aveHitRMS
average RMS of an isolated hit
Definition: DataStructs.h:630
std::vector< TrajPoint > Pts
Trajectory points.
Definition: DataStructs.h:186
TMVA::Reader * showerParentReader
Definition: DataStructs.h:570
bool ChkVtxAssociations(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:2058
T get(std::string const &key) const
Definition: ParameterSet.h:314
int Plane
Select plane.
Definition: DebugStruct.h:22
std::vector< float > neutralVxCuts
Definition: DataStructs.h:544
void MakeHaloTj(TCSlice &slc, Trajectory &muTj, bool prt)
Definition: Utils.cxx:53
std::vector< short > EndVtx
Definition: DataStructs.h:386
std::vector< VtxStore > vtxs
2D vertices
Definition: DataStructs.h:668
std::array< unsigned short, 2 > EndPt
First and last point in the trajectory that has charge.
Definition: DataStructs.h:198
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:73
std::vector< float > match3DCuts
3D matching cuts
Definition: DataStructs.h:551
void ScoreVertices(TCSlice &slc)
Definition: TCVertex.cxx:2132
std::string PrintHit(const TCHit &tch)
Definition: Utils.cxx:6347
std::vector< float > Envelope
Definition: DataStructs.h:400
call study functions to develop cuts, etc
Definition: DataStructs.h:528
const geo::GeometryCore * geom
Definition: DataStructs.h:567
TMVA::Reader fMVAReader
std::vector< unsigned short > minPts
min number of Pts required to make a trajectory
Definition: DataStructs.h:557
bool has_key(std::string const &key) const
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
void ReleaseHits(TCSlice &slc, Trajectory const &tj)
Definition: Utils.cxx:1049
void DefineTjParents(TCSlice &slc, bool prt)
Definition: Utils.cxx:173
float HitsPosTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, float &sum, HitStatus_t hitRequest)
Definition: Utils.cxx:4174
std::vector< unsigned int > firstWire
the first wire with a hit
Definition: DataStructs.h:646
std::vector< float > BeginChg
Definition: DataStructs.h:380
int ID
ID that is local to one slice.
Definition: DataStructs.h:200
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
Detector simulation of raw signals on wires.
unsigned int Nplanes(TPCID const &tpcid=details::tpc_zero) const
Returns the total number of planes in the specified TPC.
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:599
std::vector< TCHit > slHits
Definition: DataStructs.h:661
std::vector< float > chargeCuts
Definition: DataStructs.h:553
bool aveHitRMSValid
set true when the average hit RMS is well-known
Definition: DataStructs.h:639
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::vector< int > EnvPlane
Definition: DataStructs.h:401
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:14
float maxWireSkipNoSignal
max number of wires to skip w/o a signal on them
Definition: DataStructs.h:573
std::vector< short > MCSMom
Definition: DataStructs.h:388
int TPC
Select TPC.
Definition: DebugStruct.h:21
void CheckTraj(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:997
std::vector< float > vtxScoreWeights
Definition: DataStructs.h:543
unsigned int CTP_t
Definition: DataStructs.h:45
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
geo::TPCID TPCID
Definition: DataStructs.h:654
bool DecodeDebugString(std::string strng)
Definition: Utils.cxx:5068
std::vector< int > StageNum
Definition: DataStructs.h:396
std::vector< recob::Hit > const * srcHits
Definition: DataStructs.h:615
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Definition: DataStructs.h:578
bool AnalyzeHits()
Definition: Utils.cxx:4261
Contains all timing reference information for the detector.
std::vector< float > pfpStitchCuts
cuts for stitching between TPCs
Definition: DataStructs.h:555
std::vector< short > deltaRayTag
Definition: DataStructs.h:545
save cosmic ray tree
Definition: DataStructs.h:529
unsigned short nPlanes
Definition: DataStructs.h:406
float VLAStepSize
Definition: DataStructs.h:576
short StepDir
-1 = going US (-> small wire#), 1 = going DS (-> large wire#)
Definition: DataStructs.h:205
std::bitset< 128 > AlgMod
Bit set if algorithm AlgBit_t modifed the trajectory.
Definition: DataStructs.h:188
call study functions to develop cuts, etc
Definition: DataStructs.h:526
void LastEndMerge(TCSlice &slc, CTP_t inCTP)
Definition: StepUtils.cxx:3339
std::vector< short > muonTag
Definition: DataStructs.h:546
void UpdateVxEnvironment(TCSlice &slc)
Definition: Utils.cxx:3765
std::vector< float > BeginWir
Definition: DataStructs.h:377
geo::PlaneID DecodeCTP(CTP_t CTP)
std::optional< T > get_if_present(std::string const &key) const
Definition: ParameterSet.h:267
std::vector< float > EndChg
Definition: DataStructs.h:385
TrajClusterAlg(fhicl::ParameterSet const &pset)
std::bitset< 128 > dbgAlg
Allow user to turn on debug printing in algorithms (that print...)
Definition: DataStructs.h:579
std::vector< recob::Hit > const * allHits
Definition: DataStructs.h:614
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
Definition: DataStructs.h:49
std::vector< unsigned int > nWires
Definition: DataStructs.h:645
void KillPoorVertices(TCSlice &slc)
Definition: TCVertex.cxx:2159
std::vector< float > chkStopCuts
Bragg peak finder configuration.
Definition: DataStructs.h:548
int ID
ID of the recob::Slice (not the sub-slice)
Definition: DataStructs.h:656
void PrintDebugMode()
Definition: Utils.cxx:5300
void EndMerge(TCSlice &slc, CTP_t inCTP, bool lastPass)
Definition: StepUtils.cxx:3486
unsigned int subRun
Definition: DataStructs.h:628
void MergeTPHits(std::vector< unsigned int > &tpHits, std::vector< recob::Hit > &newHitCol, std::vector< unsigned int > &newHitAssns) const
std::vector< float > vtx2DCuts
Max position pull, max Position error rms.
Definition: DataStructs.h:541
std::bitset< 8 > Strategy
Definition: DataStructs.h:208
void Finish3DShowers(TCSlice &slc)
Definition: TCShower.cxx:154
IDparameter< geo::TPCID > TPCID
Member type of validated geo::TPCID parameter.
void SetSourceHits(std::vector< recob::Hit > const &srcHits)
void FindPFParticles(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: PFPUtils.cxx:190
float multHitSep
preferentially "merge" hits with < this separation
Definition: DataStructs.h:565
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
void ChkInTraj(std::string someText, TCSlice &slc)
void ReconstructAllTraj(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, CTP_t inCTP)
float JTMaxHitSep2
Definition: DataStructs.h:577
unsigned short NumPtsWithCharge(const TCSlice &slc, const Trajectory &tj, bool includeDeadWires)
Definition: Utils.cxx:2072
std::vector< int > EnvShowerID
Definition: DataStructs.h:403
void MakeJunkVertices(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:31
TCEvent evt
Definition: DataStructs.cxx:8
call study functions to develop cuts, etc (see TCTruth.cxx)
Definition: DataStructs.h:525
void Find2DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, const CTP_t &inCTP, unsigned short pass)
Definition: TCVertex.cxx:131
std::vector< int > IsShowerTj
Definition: DataStructs.h:393
std::vector< short > PlaneNum
Definition: DataStructs.h:390
short recoSlice
only reconstruct the slice with ID (0 = all)
Definition: DataStructs.h:580
void StitchPFPs()
Definition: PFPUtils.cxx:42
bool useChannelStatus
Definition: DataStructs.h:601
master switch for turning on debug mode
Definition: DataStructs.h:524
tag cosmic rays
Definition: DataStructs.h:530
void DefinePFPParents(TCSlice &slc)
Definition: PFPUtils.cxx:2860
don&#39;t mess with this line
Definition: DataStructs.h:487
art framework interface to geometry description
CTP_t CTP
set to an invalid CTP
Definition: DebugStruct.h:23
Event finding and building.
void SetTPEnvironment(TCSlice &slc, CTP_t inCTP)
Definition: Utils.cxx:3531
void FindJunkTraj(TCSlice &slc, CTP_t inCTP)