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