LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Private Member Functions | Private Attributes | List of all members
tca::TrajClusterAlg Class Reference

#include "TrajClusterAlg.h"

Public Member Functions

 TrajClusterAlg (fhicl::ParameterSet const &pset)
 
bool SetInputHits (std::vector< recob::Hit > const &inputHits, unsigned int run, unsigned int event)
 
void SetInputSpts (std::vector< recob::SpacePoint > const &sptHandle)
 
void SetSourceHits (std::vector< recob::Hit > const &srcHits)
 
void ExpectSlicedHits ()
 
void RunTrajClusterAlg (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
 
bool CreateSlice (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
 
void FinishEvent ()
 
void DefineShTree (TTree *t)
 
unsigned short GetSlicesSize () const
 
TCSlice const & GetSlice (unsigned short sliceIndex) const
 
void MergeTPHits (std::vector< unsigned int > &tpHits, std::vector< recob::Hit > &newHitCol, std::vector< unsigned int > &newHitAssns) const
 
std::vector< unsigned int > const & GetAlgModCount () const
 
std::vector< std::string > const & GetAlgBitNames () const
 
void ClearResults ()
 Deletes all the results. More...
 

Private Member Functions

recob::Hit MergeTPHitsOnWire (std::vector< unsigned int > &tpHits) const
 
void ReconstructAllTraj (detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, CTP_t inCTP)
 
void FindJunkTraj (TCSlice &slc, CTP_t inCTP)
 
void ChkInTraj (std::string someText, TCSlice &slc)
 

Private Attributes

TTree * showertree
 
calo::CalorimetryAlg fCaloAlg
 
TMVA::Reader fMVAReader
 
std::vector< unsigned int > fAlgModCount
 

Detailed Description

Definition at line 37 of file TrajClusterAlg.h.

Constructor & Destructor Documentation

tca::TrajClusterAlg::TrajClusterAlg ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 35 of file TrajClusterAlg.cxx.

References tca::AlgBitNames, tca::TCConfig::angleRanges, tca::TCEvent::aveHitRMS, tca::TCEvent::aveHitRMSValid, tca::TCConfig::caloAlg, tca::TCConfig::chargeCuts, tca::TCConfig::chkStopCuts, art::errors::Configuration, tca::ConfigureMVA(), tca::TCConfig::dbgAlg, tca::DecodeDebugString(), tca::TCConfig::deltaRayTag, tca::TCConfig::doForecast, tca::TCConfig::electronTag, tca::EndFlagNames, tca::TCEvent::eventsProcessed, tca::evt, fAlgModCount, fCaloAlg, fMVAReader, for(), fhicl::ParameterSet::get(), fhicl::ParameterSet::get_if_present(), fhicl::ParameterSet::has_key(), tca::TCConfig::hitErrFac, tca::TCConfig::JTMaxHitSep2, tca::kAlgBitSize, tca::kChkInTraj, tca::kDebug, tca::kFlagBitSize, tca::TCConfig::kinkCuts, tca::kSaveCRTree, tca::kSaveShowerTree, tca::kStepDir, tca::kStudy1, tca::kStudy2, tca::kStudy3, tca::kStudy4, tca::kTagCosmics, tca::kTestBeam, tca::TCConfig::match3DCuts, tca::TCConfig::maxAngleCode, tca::TCConfig::maxChi, tca::TCConfig::maxWireSkipNoSignal, tca::TCConfig::maxWireSkipWithSignal, tca::TCConfig::minMCSMom, tca::TCConfig::minPts, tca::TCConfig::minPtsFit, tca::TCConfig::modes, tca::TCConfig::multHitSep, tca::TCConfig::muonTag, tca::TCConfig::neutralVxCuts, tca::TCConfig::nPtsAve, tca::TCConfig::pfpStitchCuts, tca::TCConfig::projectionErrFactor, tca::TCConfig::qualityCuts, tca::TCConfig::showerParentReader, tca::TCConfig::showerTag, tca::tcc, tca::TCConfig::useAlg, tca::TCConfig::useChannelStatus, tca::TCConfig::VLAStepSize, tca::TCConfig::vtx2DCuts, tca::TCConfig::vtx3DCuts, and tca::TCConfig::vtxScoreWeights.

36  : fCaloAlg(pset.get<fhicl::ParameterSet>("CaloAlg")), fMVAReader("Silent")
37  {
38  tcc.showerParentReader = &fMVAReader;
39 
40  bool badinput = false;
41  // set all configurable modes false
42  tcc.modes.reset();
43 
44  // default mode is stepping US -> DS
45  tcc.modes[kStepDir] = true;
46  short userMode = 1;
47  if (pset.has_key("Mode")) userMode = pset.get<short>("Mode");
48  if (userMode < 0) tcc.modes[kStepDir] = false;
49  if (pset.has_key("DoForecast")) tcc.doForecast = pset.get<bool>("DoForecast");
50  if (pset.has_key("UseChannelStatus")) tcc.useChannelStatus = pset.get<bool>("UseChannelStatus");
51  if (pset.has_key("StudyMode")) {
52  std::cout << "StudyMode is not valid anymore. Try Study: 1 or Study: 2, etc/n";
53  } // old StudyMode
54  if (pset.has_key("Study")) {
55  userMode = pset.get<short>("Study");
56  if (userMode == 1) tcc.modes[kStudy1] = true;
57  if (userMode == 2) tcc.modes[kStudy2] = true;
58  if (userMode == 3) tcc.modes[kStudy3] = true;
59  if (userMode == 4) tcc.modes[kStudy4] = true;
60  } // new Study mode
61  if (pset.has_key("SaveShowerTree"))
62  tcc.modes[kSaveShowerTree] = pset.get<bool>("SaveShowerTree");
63  if (pset.has_key("SaveCRTree")) tcc.modes[kSaveCRTree] = pset.get<bool>("SaveCRTree");
64  if (pset.has_key("TagCosmics")) tcc.modes[kTagCosmics] = pset.get<bool>("TagCosmics");
65  std::vector<std::string> skipAlgsVec;
66  if (pset.has_key("SkipAlgs")) skipAlgsVec = pset.get<std::vector<std::string>>("SkipAlgs");
67  std::vector<std::string> debugConfigVec;
68  if (pset.has_key("DebugConfig"))
69  debugConfigVec = pset.get<std::vector<std::string>>("DebugConfig");
70 
71  tcc.hitErrFac = pset.get<float>("HitErrFac", 0.4);
72  // Allow the user to specify the typical hit rms for small-angle tracks
73  std::vector<float> aveHitRMS;
74  if (pset.has_key("AveHitRMS")) aveHitRMS = pset.get<std::vector<float>>("AveHitRMS");
75  // Turn off the call to AnalyzeHits
76  if (!aveHitRMS.empty()) {
77  evt.aveHitRMSValid = true;
78  evt.aveHitRMS = aveHitRMS;
79  }
80  tcc.angleRanges = pset.get<std::vector<float>>("AngleRanges");
81  tcc.nPtsAve = pset.get<short>("NPtsAve", 20);
82  tcc.minPtsFit = pset.get<std::vector<unsigned short>>("MinPtsFit");
83  tcc.minPts = pset.get<std::vector<unsigned short>>("MinPts");
84  tcc.maxAngleCode = pset.get<std::vector<unsigned short>>("MaxAngleCode");
85  tcc.minMCSMom = pset.get<std::vector<short>>("MinMCSMom");
86  tcc.maxChi = pset.get<float>("MaxChi", 10);
87  tcc.chargeCuts = pset.get<std::vector<float>>("ChargeCuts", {3, 0.15, 0.25});
88  tcc.multHitSep = pset.get<float>("MultHitSep", 2.5);
89  tcc.kinkCuts = pset.get<std::vector<float>>("KinkCuts");
90  tcc.qualityCuts = pset.get<std::vector<float>>("QualityCuts", {0.8, 3});
91  tcc.maxWireSkipNoSignal = pset.get<float>("MaxWireSkipNoSignal", 1);
92  tcc.maxWireSkipWithSignal = pset.get<float>("MaxWireSkipWithSignal", 100);
93  tcc.projectionErrFactor = pset.get<float>("ProjectionErrFactor", 2);
94  tcc.VLAStepSize = pset.get<float>("VLAStepSize", 1.5);
95  tcc.JTMaxHitSep2 = pset.get<float>("JTMaxHitSep", 2);
96  tcc.deltaRayTag = pset.get<std::vector<short>>("DeltaRayTag", {-1, -1, -1});
97  tcc.muonTag = pset.get<std::vector<short>>("MuonTag", {-1, -1, -1, -1});
98  if (pset.has_key("ElectronTag")) tcc.electronTag = pset.get<std::vector<float>>("ElectronTag");
99  tcc.showerTag = pset.get<std::vector<float>>("ShowerTag", {-1, -1, -1, -1, -1, -1});
100  std::string fMVAShowerParentWeights = "NA";
101  pset.get_if_present<std::string>("MVAShowerParentWeights", fMVAShowerParentWeights);
102  tcc.chkStopCuts = pset.get<std::vector<float>>("ChkStopCuts", {-1, -1, -1});
103  if (pset.has_key("MatchTruth")) {
104  std::cout << "MatchTruth is not used. Use ClusterAnaV2 or DebugConfig to configure\n";
105  }
106  tcc.vtx2DCuts = pset.get<std::vector<float>>("Vertex2DCuts", {-1, -1, -1, -1, -1, -1, -1});
107  tcc.vtx3DCuts = pset.get<std::vector<float>>("Vertex3DCuts", {-1, -1});
108  tcc.vtxScoreWeights = pset.get<std::vector<float>>("VertexScoreWeights");
109  tcc.match3DCuts = pset.get<std::vector<float>>("Match3DCuts", {-1, -1, -1, -1, -1});
110  tcc.pfpStitchCuts = pset.get<std::vector<float>>("PFPStitchCuts", {-1});
111  // don't search for a neutrino vertex in test beam mode
112  tcc.modes[kTestBeam] = pset.get<bool>("TestBeam", false);
113  pset.get_if_present<std::vector<float>>("NeutralVxCuts", tcc.neutralVxCuts);
114  if (tcc.JTMaxHitSep2 > 0) tcc.JTMaxHitSep2 *= tcc.JTMaxHitSep2;
115 
116  // in the following section we ensure that the fcl vectors are appropriately sized so that later references are valid
117  if (tcc.minPtsFit.size() != tcc.minPts.size()) badinput = true;
118  if (tcc.maxAngleCode.size() != tcc.minPts.size()) badinput = true;
119  if (tcc.minMCSMom.size() != tcc.minPts.size()) badinput = true;
120  if (badinput)
121  throw art::Exception(art::errors::Configuration)
122  << "Bad input from fcl file. Vector lengths for MinPtsFit, MaxAngleRange and MinMCSMom "
123  "should be defined for each reconstruction pass";
124 
125  if (tcc.vtx2DCuts.size() < 10)
126  throw art::Exception(art::errors::Configuration)
127  << "Vertex2DCuts must be size 11\n 0 = Max length definition for short TJs\n 1 = Max "
128  "vtx-TJ sep short TJs\n 2 = Max vtx-TJ sep long TJs\n 3 = Max position pull for >2 "
129  "TJs\n 4 = Max vtx position error\n 5 = Min MCSMom for one of two TJs\n 6 = Min "
130  "fraction of wires hit btw vtx and Tjs\n 7 = Min Score\n 8 = min ChgFrac at a vtx or "
131  "merge point\n 9 = max MCSMom asymmetry, 10 = require chg on wires btw vtx and tjs in "
132  "induction planes?";
133  if (tcc.vtx2DCuts.size() == 10) {
134  // User didn't specify a requirement for the presence of charge between a vertex and the start of the
135  // vertex Tjs in induction planes. Assume that it is required
136  tcc.vtx2DCuts.resize(11, 1.);
137  }
138  if (tcc.vtx3DCuts.size() < 3)
139  throw art::Exception(art::errors::Configuration)
140  << "Vertex3DCuts must be size > 2\n 0 = 2D Vtx max dX (cm)\n 1 = 2D Vtx max pull\n 2 = max "
141  "3D separation (cm) btw PFP and vertex";
142  if (tcc.vtx3DCuts.size() == 2) {
143  std::cout << "WARNING: Vertex3DCuts is size 2 but should be size 3, where Vertex3DCuts[2] = "
144  "max 3D separation (cm) btw a PFP and a 3D vertex. Setting it to 3 cm\n";
145  tcc.vtx3DCuts.resize(3, 3.);
146  }
147  if (tcc.kinkCuts.size() < 4) {
148  throw art::Exception(art::errors::Configuration)
149  << "KinkCuts must be size 3\n 0 = Number of points to fit at the end of the trajectory\n 1 "
150  "= Minimum kink significance\n 2 = Use charge in significance calculation? (yes if > "
151  "0)\n 3 = 3D kink fit length (cm). \nYou are using an out-of-date specification?\n";
152  }
153  // throw an exception if the user appears to be using an old version of KinkCuts where
154  // KinkCuts[0] was a kink angle cut
155  if (tcc.kinkCuts[0] > 0 && tcc.kinkCuts[0] < 1.) {
156  throw art::Exception(art::errors::Configuration)
157  << "Are you using an out-of-date specification for KinkCuts? KinkCuts[0] is the number of "
158  "points to fit.\n";
159  }
160 
161  if (tcc.chargeCuts.size() != 3)
162  throw art::Exception(art::errors::Configuration)
163  << "ChargeCuts must be size 3\n 0 = Charge pull cut\n 1 = Min allowed fractional chg RMS\n "
164  "2 = Max allowed fractional chg RMS";
165  // dressed muons - change next line
166  if (tcc.muonTag.size() < 4)
167  throw art::Exception(art::errors::Configuration)
168  << "MuonTag must be size 4\n 0 = minPtsFit\n 1 = minMCSMom\n 2= maxWireSkipNoSignal\n 3 = "
169  "min delta ray length for tagging\n 4 = dress muon window size (optional)";
170  if (tcc.deltaRayTag.size() != 3)
171  throw art::Exception(art::errors::Configuration)
172  << "DeltaRayTag must be size 3\n 0 = Max endpoint sep\n 1 = min MCSMom\n 2 = max MCSMom";
173  if (tcc.chkStopCuts.size() < 3)
174  throw art::Exception(art::errors::Configuration)
175  << "ChkStopCuts must be size 3\n 0 = Min Charge ratio\n 1 = Charge slope pull cut\n 2 = "
176  "Charge fit chisq cut\n 3 = BraggSplit FOM (optional)";
177  if (tcc.showerTag.size() < 13) {
178  std::cout
179  << "ShowerTag must be size 13\n 0 = Mode\n 1 = max MCSMom\n 2 = max separation (WSE "
180  "units)\n 3 = Max angle diff\n 4 = Factor * rms width\n 5 = Min half width\n 6 = min "
181  "total Tps\n 7 = Min Tjs\n 8 = max parent FOM\n 9 = max direction FOM 10 = max "
182  "AspectRatio\n 11 = min Score to preserve a vertex\n 12 = Debug showers in CTP\n";
183  std::cout << " Fixing this problem...";
184  tcc.showerTag.resize(13);
185  // set the min score to 0
186  tcc.showerTag[11] = 0;
187  // turn off printing
188  tcc.showerTag[12] = -1;
189  }
190  if (tcc.match3DCuts.size() < 6)
191  throw art::Exception(art::errors::Configuration)
192  << "Match3DCuts must be size 5\n 0 = dx (cm) matching cut\n 1 = max number of 3D "
193  "combinations\n 2 = min length for 2-view match\n 3 = number of TP3Ds in each plane to "
194  "fit in each PFP section\n 4 = max pull for accepting TP3Ds in sections\n 5 = max "
195  "ChiDOF for a SectionFit";
196 
197  // check the angle ranges and convert from degrees to radians
198  if (tcc.angleRanges.back() < 90) {
199  mf::LogVerbatim("TC") << "Last element of AngleRange != 90 degrees. Fixing it\n";
200  tcc.angleRanges.back() = 90;
201  }
202 
203  // convert PFP stitch cuts
204  if (tcc.pfpStitchCuts.size() > 1 && tcc.pfpStitchCuts[0] > 0) {
205  // square the separation cut
206  tcc.pfpStitchCuts[0] *= tcc.pfpStitchCuts[0];
207  // convert angle to cos
208  tcc.pfpStitchCuts[1] = cos(tcc.pfpStitchCuts[1]);
209  }
210 
211  // configure algorithm debugging. Configuration for debugging standard stepping
212  // is done in Utils/AnalyzeHits when the input hit collection is passed to SetInputHits
213  tcc.modes[kDebug] = false;
214  tcc.dbgAlg.reset();
215  for (auto strng : debugConfigVec) {
216  // try to interpret this as a C:T:P:W:Tick specification or something similar
217  if (!DecodeDebugString(strng)) {
218  // try to set a dbgAlg bit
219  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib) {
220  if (strng == AlgBitNames[ib]) {
221  tcc.dbgAlg[ib] = true;
222  tcc.modes[kDebug] = true;
223  break;
224  }
225  } // ib
226  // print some instructions and quit if there was a failure
227  if (!tcc.modes[kDebug]) {
228  std::cout << "DecodeDebugString failed: " << strng << "\n";
229  DecodeDebugString("instruct");
230  exit(1);
231  }
232  } // DecodeDebugString failed
233  } // strng
234 
235  for (auto& range : tcc.angleRanges) {
236  if (range < 0 || range > 90)
237  throw art::Exception(art::errors::Configuration)
238  << "Invalid angle range " << range << " Must be 0 - 90 degrees";
239  range *= M_PI / 180;
240  } // range
241 
242  // Ensure that the size of the AlgBitNames vector is consistent with the AlgBit typedef enum
243  if (kAlgBitSize != AlgBitNames.size())
244  throw art::Exception(art::errors::Configuration)
245  << "kAlgBitSize " << kAlgBitSize << " != AlgBitNames size " << AlgBitNames.size();
246  if (kAlgBitSize > 128)
247  throw art::Exception(art::errors::Configuration)
248  << "Increase the size of UseAlgs to at least " << kAlgBitSize;
249  fAlgModCount.resize(kAlgBitSize);
250 
251  if (kFlagBitSize != EndFlagNames.size())
252  throw art::Exception(art::errors::Configuration)
253  << "kFlagBitSize " << kFlagBitSize << " != EndFlagNames size " << EndFlagNames.size();
254 
255  if (kFlagBitSize > 8)
256  throw art::Exception(art::errors::Configuration)
257  << "Increase the size of EndFlag to at least " << kFlagBitSize;
258 
259  bool printHelp = false;
260  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
261  tcc.useAlg[ib] = true;
262 
263  // turn off the special algs
264  // Do an exhaustive (and slow) check of the hit -> trajectory associations
265  tcc.useAlg[kChkInTraj] = false;
266 
267  for (auto strng : skipAlgsVec) {
268  bool gotit = false;
269  if (strng == "All") {
270  // turn everything off
271  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
272  tcc.useAlg[ib] = false;
273  gotit = true;
274  break;
275  } // All off
276  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib) {
277  if (strng == AlgBitNames[ib]) {
278  tcc.useAlg[ib] = false;
279  gotit = true;
280  break;
281  }
282  } // ib
283  if (!gotit) {
284  std::cout << "******* Unknown SkipAlgs input string '" << strng << "'\n";
285  printHelp = true;
286  }
287  } // strng
288  if (printHelp) {
289  std::cout << "Valid AlgNames:";
290  for (auto strng : AlgBitNames)
291  std::cout << " " << strng;
292  std::cout << "\n";
293  std::cout << "Or specify All to turn all algs off\n";
294  }
295  // Configure the TMVA reader for the shower parent BDT
296  if (fMVAShowerParentWeights != "NA" && tcc.showerTag[0] > 0)
297  ConfigureMVA(tcc, fMVAShowerParentWeights);
298 
299  evt.eventsProcessed = 0;
300 
301  tcc.caloAlg = &fCaloAlg;
302  }
tca::kTestBeam
Expect tracks entering from the front face. Don&#39;t create neutrino PFParticles.
Definition: DataStructs.h:525
tca::TCConfig::caloAlg
calo::CalorimetryAlg * caloAlg
Definition: DataStructs.h:570
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
tca::ConfigureMVA
void ConfigureMVA(TCConfig &tcc, std::string fMVAShowerParentWeights)
Definition: TCShower.cxx:34
tca::TCConfig::kinkCuts
std::vector< float > kinkCuts
kink finder algorithm
Definition: DataStructs.h:552
tca::AlgBitNames
const std::vector< std::string > AlgBitNames
Definition: DataStructs.cxx:16
tca::TCConfig::minPtsFit
std::vector< unsigned short > minPtsFit
Reconstruct in several passes.
Definition: DataStructs.h:558
tca::TCConfig::vtx3DCuts
std::vector< float > vtx3DCuts
2D vtx -> 3D vtx matching cuts
Definition: DataStructs.h:544
tca::TCConfig::qualityCuts
std::vector< float > qualityCuts
Min points/wire, min consecutive pts after a gap.
Definition: DataStructs.h:556
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::kStepDir
step from US -> DS (true) or DS -> US (false)
Definition: DataStructs.h:524
tca::TrajClusterAlg::fCaloAlg
calo::CalorimetryAlg fCaloAlg
Definition: TrajClusterAlg.h:85
tca::TCConfig::maxWireSkipWithSignal
float maxWireSkipWithSignal
max number of wires to skip with a signal on them
Definition: DataStructs.h:575
tca::kStudy3
call study functions to develop cuts, etc
Definition: DataStructs.h:529
for
for(Int_t i=0;i< nentries;i++)
Definition: comparison.C:30
tca::TCConfig::electronTag
std::vector< float > electronTag
Definition: DataStructs.h:549
tca::TCEvent::eventsProcessed
unsigned int eventsProcessed
Definition: DataStructs.h:630
tca::TCConfig::doForecast
bool doForecast
See TCMode_t above.
Definition: DataStructs.h:601
tca::TCConfig::angleRanges
std::vector< float > angleRanges
list of max angles for each angle range
Definition: DataStructs.h:562
tca::EndFlagNames
const std::vector< std::string > EndFlagNames
Definition: DataStructs.cxx:87
tca::TCConfig::showerTag
std::vector< float > showerTag
shower-like trajectory tagging + shower reconstruction
Definition: DataStructs.h:551
tca::TCConfig::maxAngleCode
std::vector< unsigned short > maxAngleCode
max allowed angle code for each pass
Definition: DataStructs.h:560
tca::TrajClusterAlg::fAlgModCount
std::vector< unsigned int > fAlgModCount
Definition: TrajClusterAlg.h:88
tca::TCConfig::projectionErrFactor
float projectionErrFactor
Definition: DataStructs.h:576
tca::kSaveShowerTree
save shower tree
Definition: DataStructs.h:533
tca::TCConfig::nPtsAve
short nPtsAve
dump trajectory points
Definition: DataStructs.h:599
tca::kFlagBitSize
don&#39;t mess with this line
Definition: DataStructs.h:508
tca::TCConfig::minMCSMom
std::vector< short > minMCSMom
Min MCSMom for each pass.
Definition: DataStructs.h:561
tca::TCEvent::aveHitRMS
std::vector< float > aveHitRMS
average RMS of an isolated hit
Definition: DataStructs.h:631
tca::TCConfig::showerParentReader
TMVA::Reader * showerParentReader
Definition: DataStructs.h:571
tca::TCConfig::neutralVxCuts
std::vector< float > neutralVxCuts
Definition: DataStructs.h:546
tca::TCConfig::match3DCuts
std::vector< float > match3DCuts
3D matching cuts
Definition: DataStructs.h:553
tca::kStudy4
call study functions to develop cuts, etc
Definition: DataStructs.h:530
tca::TrajClusterAlg::fMVAReader
TMVA::Reader fMVAReader
Definition: TrajClusterAlg.h:86
tca::TCConfig::minPts
std::vector< unsigned short > minPts
min number of Pts required to make a trajectory
Definition: DataStructs.h:559
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:600
tca::TCConfig::chargeCuts
std::vector< float > chargeCuts
Definition: DataStructs.h:555
tca::TCEvent::aveHitRMSValid
bool aveHitRMSValid
set true when the average hit RMS is well-known
Definition: DataStructs.h:640
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
tca::TCConfig::maxWireSkipNoSignal
float maxWireSkipNoSignal
max number of wires to skip w/o a signal on them
Definition: DataStructs.h:574
tca::TCConfig::vtxScoreWeights
std::vector< float > vtxScoreWeights
Definition: DataStructs.h:545
tca::DecodeDebugString
bool DecodeDebugString(std::string strng)
Definition: Utils.cxx:5065
tca::TCConfig::useAlg
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Definition: DataStructs.h:579
tca::TCConfig::pfpStitchCuts
std::vector< float > pfpStitchCuts
cuts for stitching between TPCs
Definition: DataStructs.h:557
tca::TCConfig::deltaRayTag
std::vector< short > deltaRayTag
Definition: DataStructs.h:547
tca::kSaveCRTree
save cosmic ray tree
Definition: DataStructs.h:531
tca::TCConfig::VLAStepSize
float VLAStepSize
Definition: DataStructs.h:577
tca::kStudy2
call study functions to develop cuts, etc
Definition: DataStructs.h:528
tca::TCConfig::muonTag
std::vector< short > muonTag
Definition: DataStructs.h:548
tca::TCConfig::dbgAlg
std::bitset< 128 > dbgAlg
Allow user to turn on debug printing in algorithms (that print...)
Definition: DataStructs.h:580
tca::TCConfig::chkStopCuts
std::vector< float > chkStopCuts
Bragg peak finder configuration.
Definition: DataStructs.h:550
tca::TCConfig::vtx2DCuts
std::vector< float > vtx2DCuts
Max position pull, max Position error rms.
Definition: DataStructs.h:543
tca::TCConfig::multHitSep
float multHitSep
preferentially "merge" hits with < this separation
Definition: DataStructs.h:567
tca::TCConfig::JTMaxHitSep2
float JTMaxHitSep2
Definition: DataStructs.h:578
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::kStudy1
call study functions to develop cuts, etc (see TCTruth.cxx)
Definition: DataStructs.h:527
tca::TCConfig::useChannelStatus
bool useChannelStatus
Definition: DataStructs.h:602
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:526
tca::kTagCosmics
tag cosmic rays
Definition: DataStructs.h:532
tca::kAlgBitSize
don&#39;t mess with this line
Definition: DataStructs.h:489

Member Function Documentation

void tca::TrajClusterAlg::ChkInTraj ( std::string  someText,
TCSlice &  slc 
)
private

Definition at line 888 of file TrajClusterAlg.cxx.

References tca::AlgBitNames, util::end(), fAlgModCount, tca::TCSlice::isValid, tca::kChkInTraj, tca::kKilled, tca::kUsedHits, tca::PrintHit(), tca::PrintTrajectory(), tca::PutTrajHitsInVector(), tca::TCSlice::slHits, tca::tcc, tca::TCSlice::tjs, tca::TCConfig::useAlg, and tca::TCSlice::vtxs.

889  {
890  // Check slc.tjs -> InTraj associations
891 
892  if (!tcc.useAlg[kChkInTraj]) return;
893 
894  ++fAlgModCount[kChkInTraj];
895 
896  int tID;
897  unsigned int iht;
898  unsigned short itj = 0;
899  std::vector<unsigned int> tHits;
900  std::vector<unsigned int> atHits;
901  for (auto& tj : slc.tjs) {
902  // ignore abandoned trajectories
903  if (tj.AlgMod[kKilled]) continue;
904  tID = tj.ID;
905  for (auto& tp : tj.Pts) {
906  if (tp.Hits.size() > 16) {
907  tj.AlgMod[kKilled] = true;
908  mf::LogVerbatim("TC")
909  << "ChkInTraj: More than 16 hits created a UseHit bitset overflow\n";
910  slc.isValid = false;
911  std::cout << "ChkInTraj major failure\n";
912  return;
913  }
914  } // tp
915  if (tj.AlgMod[kKilled]) {
916  std::cout << someText << " ChkInTraj hit size mis-match in tj ID " << tj.ID
917  << " AlgBitNames";
918  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
919  if (tj.AlgMod[ib]) std::cout << " " << AlgBitNames[ib];
920  std::cout << "\n";
921  continue;
922  }
923  tHits = PutTrajHitsInVector(tj, kUsedHits);
924  if (tHits.size() < 2) {
925  mf::LogVerbatim("TC") << someText << " ChkInTraj: Insufficient hits in traj " << tj.ID;
926  PrintTrajectory("CIT", slc, tj, USHRT_MAX);
927  tj.AlgMod[kKilled] = true;
928  continue;
929  }
930  std::sort(tHits.begin(), tHits.end());
931  atHits.clear();
932  for (iht = 0; iht < slc.slHits.size(); ++iht) {
933  if (slc.slHits[iht].InTraj == tID) atHits.push_back(iht);
934  } // iht
935  if (atHits.size() < 2) {
936  mf::LogVerbatim("TC") << someText << " ChkInTraj: Insufficient hits in atHits in traj "
937  << tj.ID << " Killing it";
938  tj.AlgMod[kKilled] = true;
939  continue;
940  }
941  if (!std::equal(tHits.begin(), tHits.end(), atHits.begin())) {
942  mf::LogVerbatim myprt("TC");
943  myprt << someText << " ChkInTraj failed: inTraj - UseHit mis-match for tj ID " << tID
944  << " tj.WorkID " << tj.WorkID << " atHits size " << atHits.size() << " tHits size "
945  << tHits.size() << " in CTP " << tj.CTP << "\n";
946  myprt << "AlgMods: ";
947  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
948  if (tj.AlgMod[ib]) myprt << " " << AlgBitNames[ib];
949  myprt << "\n";
950  myprt << "index inTraj UseHit \n";
951  for (iht = 0; iht < atHits.size(); ++iht) {
952  myprt << "iht " << iht << " " << PrintHit(slc.slHits[atHits[iht]]);
953  if (iht < tHits.size()) myprt << " " << PrintHit(slc.slHits[tHits[iht]]);
954  if (atHits[iht] != tHits[iht]) myprt << " <<< " << atHits[iht] << " != " << tHits[iht];
955  myprt << "\n";
956  slc.isValid = false;
957  } // iht
958  if (tHits.size() > atHits.size()) {
959  for (iht = atHits.size(); iht < atHits.size(); ++iht) {
960  myprt << "atHits " << iht << " " << PrintHit(slc.slHits[atHits[iht]]) << "\n";
961  } // iht
962  PrintTrajectory("CIT", slc, tj, USHRT_MAX);
963  } // tHit.size > atHits.size()
964  }
965  // check the VtxID
966  for (unsigned short end = 0; end < 2; ++end) {
967  if (tj.VtxID[end] > slc.vtxs.size()) {
968  mf::LogVerbatim("TC") << someText << " ChkInTraj: Bad VtxID " << tj.ID;
969  std::cout << someText << " ChkInTraj: Bad VtxID " << tj.ID << " vtx size "
970  << slc.vtxs.size() << "\n";
971  tj.AlgMod[kKilled] = true;
972  PrintTrajectory("CIT", slc, tj, USHRT_MAX);
973  slc.isValid = false;
974  return;
975  }
976  } // end
977  ++itj;
978  if (!slc.isValid) return;
979  } // tj
980 
981  } // ChkInTraj
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
tca::PutTrajHitsInVector
std::vector< unsigned int > PutTrajHitsInVector(const Trajectory &tj, HitStatus_t hitRequest)
Definition: Utils.cxx:2688
tca::AlgBitNames
const std::vector< std::string > AlgBitNames
Definition: DataStructs.cxx:16
mf::MaybeLogger_
Definition: MessageLogger.h:113
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::PrintTrajectory
void PrintTrajectory(std::string someText, const TCSlice &slc, const Trajectory &tj, unsigned short tPoint)
Definition: Utils.cxx:6032
util::end
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
tca::TrajClusterAlg::fAlgModCount
std::vector< unsigned int > fAlgModCount
Definition: TrajClusterAlg.h:88
tca::PrintHit
std::string PrintHit(const TCHit &tch)
Definition: Utils.cxx:6344
tca::TCConfig::useAlg
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Definition: DataStructs.h:579
void tca::TrajClusterAlg::ClearResults ( )
inline

Deletes all the results.

Definition at line 72 of file TrajClusterAlg.h.

References tca::evt, tca::slices, tca::TCEvent::sptHits, and tca::TCEvent::wireHitRange.

Referenced by cluster::TrajCluster::produce(), and SetInputHits().

73  {
74  slices.resize(0);
75  evt.sptHits.resize(0);
76  evt.wireHitRange.resize(0);
77  }
tca::TCEvent::sptHits
std::vector< std::array< unsigned int, 3 > > sptHits
SpacePoint -> Hits assns by plane.
Definition: DataStructs.h:626
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::TCEvent::wireHitRange
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > wireHitRange
Definition: DataStructs.h:619
bool tca::TrajClusterAlg::CreateSlice ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< unsigned int > &  hitsInSlice,
int  sliceID 
)

Definition at line 1213 of file TrajClusterAlg.cxx.

References tca::TCEvent::allHits, tca::DebugStuff::Cryostat, tca::DebugStuff::CTP, tca::TCConfig::dbgSlc, tca::debug, tca::EncodeCTP(), tca::evt, tca::FillWireHitRange(), tca::TCConfig::geom, tca::TCSlice::ID, tca::TCSlice::isValid, tca::kDebug, tca::TCConfig::modes, geo::GeometryCore::Nplanes(), tca::DebugStuff::Plane, util::size(), tca::TCSlice::slHits, tca::DebugStuff::Slice, tca::slices, tca::tcc, tca::DebugStuff::TPC, and tca::TCSlice::TPCID.

Referenced by RunTrajClusterAlg().

1217  {
1218  // Defines a TCSlice struct and pushes the slice onto slices.
1219  // Sets the isValid flag true if successful.
1220  if ((*evt.allHits).empty()) return false;
1221  if (hitsInSlice.size() < 2) return false;
1222 
1223  TCSlice slc;
1224  slc.ID = sliceID;
1225  slc.slHits.resize(hitsInSlice.size());
1226  bool first = true;
1227  unsigned int cstat = 0;
1228  unsigned int tpc = UINT_MAX;
1229  unsigned int cnt = 0;
1230  std::vector<unsigned int> nHitsInPln;
1231  for (auto iht : hitsInSlice) {
1232  if (iht >= (*evt.allHits).size()) return false;
1233  auto& hit = (*evt.allHits)[iht];
1234  if (first) {
1235  cstat = hit.WireID().Cryostat;
1236  tpc = hit.WireID().TPC;
1237  slc.TPCID = geo::TPCID(cstat, tpc);
1238  nHitsInPln.resize(tcc.geom->Nplanes(slc.TPCID));
1239  first = false;
1240  }
1241  if (hit.WireID().Cryostat != cstat || hit.WireID().TPC != tpc) return false;
1242  slc.slHits[cnt].allHitsIndex = iht;
1243  ++nHitsInPln[hit.WireID().Plane];
1244  ++cnt;
1245  } // iht
1246  // require at least two hits in each plane
1247  for (auto hip : nHitsInPln)
1248  if (hip < 2) return false;
1249  // Define the TCEvent wire hit range vector for this new TPC for ALL hits
1250  FillWireHitRange(slc.TPCID);
1251  // next define the Slice wire hit range vectors, UnitsPerTick, etc for this
1252  // slice
1253  if (!FillWireHitRange(clockData, detProp, slc)) return false;
1254  slc.isValid = true;
1255  slices.push_back(slc);
1256  if (tcc.modes[kDebug] && debug.Slice >= 0 && !tcc.dbgSlc) {
1257  tcc.dbgSlc = ((int)(slices.size() - 1) == debug.Slice);
1258  if (tcc.dbgSlc) std::cout << "Enabled debugging in sub-slice " << slices.size() - 1 << "\n";
1259  if (tcc.modes[kDebug] && (unsigned int)debug.Cryostat == cstat &&
1260  (unsigned int)debug.TPC == tpc && debug.Plane >= 0) {
1261  debug.CTP = EncodeCTP(
1262  (unsigned int)debug.Cryostat, (unsigned int)debug.TPC, (unsigned int)debug.Plane);
1263  }
1264  }
1265  // do a sanity check
1266  for (unsigned int iht = 0; iht < slc.slHits.size(); ++iht) {
1267  unsigned int ahi = slc.slHits[iht].allHitsIndex;
1268  if (ahi >= (*evt.allHits).size()) {
1269  std::cout << "CreateSlice: Bad allHitsIndex " << ahi << " " << (*evt.allHits).size()
1270  << "\n";
1271  return false;
1272  }
1273  } // iht
1274  return true;
1275  } // CreateSlice
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::FillWireHitRange
void FillWireHitRange(geo::TPCID inTPCID)
Definition: Utils.cxx:4323
tca::TCConfig::dbgSlc
bool dbgSlc
debug only in the user-defined slice? default is all slices
Definition: DataStructs.h:583
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
tca::DebugStuff::Cryostat
int Cryostat
Select Cryostat.
Definition: DebugStruct.h:20
tca::debug
DebugStuff debug
Definition: DebugStruct.cxx:4
tca::DebugStuff::Plane
int Plane
Select plane.
Definition: DebugStruct.h:22
tca::TCConfig::geom
const geo::GeometryCore * geom
Definition: DataStructs.h:569
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:600
tca::DebugStuff::TPC
int TPC
Select TPC.
Definition: DebugStruct.h:21
tca::TCEvent::allHits
std::vector< recob::Hit > const * allHits
Definition: DataStructs.h:615
tca::EncodeCTP
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
Definition: DataStructs.h:51
geo::GeometryCore::Nplanes
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
geo::fhicl::TPCID
IDparameter< geo::TPCID > TPCID
Member type of validated geo::TPCID parameter.
Definition: geo_types_fhicl.h:298
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:526
tca::DebugStuff::CTP
CTP_t CTP
set to an invalid CTP
Definition: DebugStruct.h:23
void tca::TrajClusterAlg::DefineShTree ( TTree *  t)

Definition at line 1172 of file TrajClusterAlg.cxx.

References tca::ShowerTreeVars::BeginAng, tca::ShowerTreeVars::BeginChg, tca::ShowerTreeVars::BeginTim, tca::ShowerTreeVars::BeginVtx, tca::ShowerTreeVars::BeginWir, tca::ShowerTreeVars::EndAng, tca::ShowerTreeVars::EndChg, tca::ShowerTreeVars::EndTim, tca::ShowerTreeVars::EndVtx, tca::ShowerTreeVars::EndWir, tca::ShowerTreeVars::Envelope, tca::ShowerTreeVars::EnvPlane, tca::ShowerTreeVars::EnvShowerID, tca::ShowerTreeVars::EnvStage, tca::TCEvent::event, tca::evt, tca::ShowerTreeVars::IsShowerParent, tca::ShowerTreeVars::IsShowerTj, tca::ShowerTreeVars::MCSMom, tca::ShowerTreeVars::nPlanes, tca::ShowerTreeVars::nStages, tca::ShowerTreeVars::PlaneNum, tca::TCEvent::run, tca::ShowerTreeVars::ShowerID, showertree, tca::ShowerTreeVars::StageName, tca::ShowerTreeVars::StageNum, tca::stv, tca::TCEvent::subRun, and tca::ShowerTreeVars::TjID.

Referenced by cluster::TrajCluster::beginJob().

1173  {
1174  showertree = t;
1175 
1176  showertree->Branch("run", &evt.run, "run/I");
1177  showertree->Branch("subrun", &evt.subRun, "subrun/I");
1178  showertree->Branch("event", &evt.event, "event/I");
1179 
1180  showertree->Branch("BeginWir", &stv.BeginWir);
1181  showertree->Branch("BeginTim", &stv.BeginTim);
1182  showertree->Branch("BeginAng", &stv.BeginAng);
1183  showertree->Branch("BeginChg", &stv.BeginChg);
1184  showertree->Branch("BeginVtx", &stv.BeginVtx);
1185 
1186  showertree->Branch("EndWir", &stv.EndWir);
1187  showertree->Branch("EndTim", &stv.EndTim);
1188  showertree->Branch("EndAng", &stv.EndAng);
1189  showertree->Branch("EndChg", &stv.EndChg);
1190  showertree->Branch("EndVtx", &stv.EndVtx);
1191 
1192  showertree->Branch("MCSMom", &stv.MCSMom);
1193 
1194  showertree->Branch("PlaneNum", &stv.PlaneNum);
1195  showertree->Branch("TjID", &stv.TjID);
1196  showertree->Branch("IsShowerTj", &stv.IsShowerTj);
1197  showertree->Branch("ShowerID", &stv.ShowerID);
1198  showertree->Branch("IsShowerParent", &stv.IsShowerParent);
1199  showertree->Branch("StageNum", &stv.StageNum);
1200  showertree->Branch("StageName", &stv.StageName);
1201 
1202  showertree->Branch("Envelope", &stv.Envelope);
1203  showertree->Branch("EnvPlane", &stv.EnvPlane);
1204  showertree->Branch("EnvStage", &stv.EnvStage);
1205  showertree->Branch("EnvShowerID", &stv.EnvShowerID);
1206 
1207  showertree->Branch("nStages", &stv.nStages);
1208  showertree->Branch("nPlanes", &stv.nPlanes);
1209 
1210  } // end DefineShTree
tca::ShowerTreeVars::EnvStage
std::vector< int > EnvStage
Definition: DataStructs.h:404
tca::ShowerTreeVars::IsShowerParent
std::vector< int > IsShowerParent
Definition: DataStructs.h:397
tca::TCEvent::event
unsigned int event
Definition: DataStructs.h:627
tca::ShowerTreeVars::EndWir
std::vector< float > EndWir
Definition: DataStructs.h:384
tca::ShowerTreeVars::EndAng
std::vector< float > EndAng
Definition: DataStructs.h:386
tca::TCEvent::run
unsigned int run
Definition: DataStructs.h:628
tca::ShowerTreeVars::BeginTim
std::vector< float > BeginTim
Definition: DataStructs.h:380
tca::ShowerTreeVars::BeginAng
std::vector< float > BeginAng
Definition: DataStructs.h:381
tca::ShowerTreeVars::EndTim
std::vector< float > EndTim
Definition: DataStructs.h:385
tca::ShowerTreeVars::ShowerID
std::vector< int > ShowerID
Definition: DataStructs.h:396
tca::stv
ShowerTreeVars stv
Definition: DataStructs.cxx:11
tca::ShowerTreeVars::StageName
std::vector< std::string > StageName
Definition: DataStructs.h:399
tca::ShowerTreeVars::TjID
std::vector< int > TjID
Definition: DataStructs.h:394
tca::ShowerTreeVars::BeginVtx
std::vector< short > BeginVtx
Definition: DataStructs.h:383
tca::ShowerTreeVars::EndVtx
std::vector< short > EndVtx
Definition: DataStructs.h:388
tca::ShowerTreeVars::Envelope
std::vector< float > Envelope
Definition: DataStructs.h:402
tca::ShowerTreeVars::BeginChg
std::vector< float > BeginChg
Definition: DataStructs.h:382
tca::ShowerTreeVars::EnvPlane
std::vector< int > EnvPlane
Definition: DataStructs.h:403
tca::ShowerTreeVars::MCSMom
std::vector< short > MCSMom
Definition: DataStructs.h:390
tca::ShowerTreeVars::StageNum
std::vector< int > StageNum
Definition: DataStructs.h:398
tca::ShowerTreeVars::nPlanes
unsigned short nPlanes
Definition: DataStructs.h:408
tca::ShowerTreeVars::BeginWir
std::vector< float > BeginWir
Definition: DataStructs.h:379
tca::ShowerTreeVars::EndChg
std::vector< float > EndChg
Definition: DataStructs.h:387
tca::TCEvent::subRun
unsigned int subRun
Definition: DataStructs.h:629
tca::ShowerTreeVars::EnvShowerID
std::vector< int > EnvShowerID
Definition: DataStructs.h:405
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::ShowerTreeVars::IsShowerTj
std::vector< int > IsShowerTj
Definition: DataStructs.h:395
tca::ShowerTreeVars::PlaneNum
std::vector< short > PlaneNum
Definition: DataStructs.h:392
void tca::TrajClusterAlg::ExpectSlicedHits ( )
inline

Definition at line 49 of file TrajClusterAlg.h.

References tca::evt, and tca::TCEvent::expectSlicedHits.

Referenced by cluster::TrajCluster::produce().

49 { evt.expectSlicedHits = true; }
tca::TCEvent::expectSlicedHits
bool expectSlicedHits
info passed from the module - used to (not) define wireHitRange
Definition: DataStructs.h:641
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
void tca::TrajClusterAlg::FindJunkTraj ( TCSlice &  slc,
CTP_t  inCTP 
)
private

Definition at line 766 of file TrajClusterAlg.cxx.

References tca::debug, tca::DecodeCTP(), tca::TCSlice::firstWire, tca::GetHitMultiplet(), tca::DebugStuff::Hit, tca::HitSep2(), tca::IsGhost(), tca::TCConfig::JTMaxHitSep2, tca::kDebug, tca::kJunkTj, tca::TCSlice::lastWire, tca::MakeJunkTraj(), tca::TCConfig::modes, tca::TCSlice::nWires, geo::PlaneID::Plane, tca::PrintHit(), tca::TCSlice::slHits, tca::tcc, tca::TrajHitsOK(), tca::TCConfig::useAlg, and tca::TCSlice::wireHitRange.

Referenced by ReconstructAllTraj().

767  {
768  // Makes junk trajectories using unassigned hits
769 
770  if (tcc.JTMaxHitSep2 <= 0) return;
771  if (!tcc.useAlg[kJunkTj]) return;
772  unsigned short plane = DecodeCTP(inCTP).Plane;
773  if ((int)slc.lastWire[plane] - 3 < (int)slc.firstWire[plane]) return;
774 
775  // shouldn't have to do this but...
776  for (auto& slHit : slc.slHits)
777  if (slHit.InTraj < 0) slHit.InTraj = 0;
778 
779  bool prt = false;
780 
781  std::vector<unsigned int> tHits;
782  // Stay well away from the last wire in the plane
783  for (unsigned int iwire = slc.firstWire[plane]; iwire < slc.lastWire[plane] - 3; ++iwire) {
784  // skip bad wires or no hits on the wire
785  if (slc.wireHitRange[plane][iwire].first > slc.slHits.size()) continue;
786  if (slc.wireHitRange[plane][iwire].second > slc.slHits.size()) continue;
787  unsigned int jwire = iwire + 1;
788  if (slc.wireHitRange[plane][jwire].first > slc.slHits.size()) continue;
789  if (slc.wireHitRange[plane][jwire].second > slc.slHits.size()) continue;
790  unsigned int ifirsthit = slc.wireHitRange[plane][iwire].first;
791  unsigned int ilasthit = slc.wireHitRange[plane][iwire].second;
792  unsigned int jfirsthit = slc.wireHitRange[plane][jwire].first;
793  unsigned int jlasthit = slc.wireHitRange[plane][jwire].second;
794  for (unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
795  if (iht >= slc.slHits.size()) break;
796  auto& islHit = slc.slHits[iht];
797  if (islHit.InTraj != 0) continue;
798  std::vector<unsigned int> iHits;
799  GetHitMultiplet(slc, iht, iHits, true);
800  prt =
801  (tcc.modes[kDebug] && std::find(iHits.begin(), iHits.end(), debug.Hit) != iHits.end());
802  if (prt) mf::LogVerbatim("TC") << "FJT: debug iht multiplet size " << iHits.size();
803  if (iHits.empty()) continue;
804  for (unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
805  if (jht >= slc.slHits.size()) break;
806  auto& jslHit = slc.slHits[jht];
807  if (jslHit.InTraj != 0) continue;
808  if (prt && HitSep2(slc, iht, jht) < 100)
809  mf::LogVerbatim("TC") << " use " << PrintHit(jslHit) << " hitSep2 "
810  << HitSep2(slc, iht, jht);
811  if (HitSep2(slc, iht, jht) > tcc.JTMaxHitSep2) continue;
812  std::vector<unsigned int> jHits;
813  GetHitMultiplet(slc, jht, jHits, true);
814  if (jHits.empty()) continue;
815  // check for hit overlap consistency
816  if (!TrajHitsOK(slc, iHits, jHits)) continue;
817  tHits.clear();
818  // add the available hits and flag them
819  for (auto iht : iHits)
820  if (slc.slHits[iht].InTraj == 0) tHits.push_back(iht);
821  for (auto jht : jHits)
822  if (slc.slHits[jht].InTraj == 0) tHits.push_back(jht);
823  for (auto tht : tHits)
824  slc.slHits[tht].InTraj = -4;
825  unsigned int loWire;
826  if (iwire != 0) { loWire = iwire - 1; }
827  else {
828  loWire = 0;
829  }
830  unsigned int hiWire = jwire + 1;
831  if (hiWire > slc.nWires[plane]) break;
832  unsigned short nit = 0;
833  while (nit < 100) {
834  bool hitsAdded = false;
835  for (unsigned int kwire = loWire; kwire <= hiWire; ++kwire) {
836  if (slc.wireHitRange[plane][kwire].first == UINT_MAX) continue;
837  if (slc.wireHitRange[plane][kwire].second == UINT_MAX) continue;
838  unsigned int kfirsthit = slc.wireHitRange[plane][kwire].first;
839  unsigned int klasthit = slc.wireHitRange[plane][kwire].second;
840  for (unsigned int kht = kfirsthit; kht <= klasthit; ++kht) {
841  if (kht >= slc.slHits.size()) continue;
842  if (slc.slHits[kht].InTraj != 0) continue;
843  // this shouldn't be needed but do it anyway
844  if (std::find(tHits.begin(), tHits.end(), kht) != tHits.end()) continue;
845  // re-purpose jHits and check for consistency
846  GetHitMultiplet(slc, kht, jHits, true);
847  if (!TrajHitsOK(slc, tHits, jHits)) continue;
848  // add them all and update the wire range
849  for (auto jht : jHits) {
850  if (slc.slHits[jht].InTraj != 0) continue;
851  tHits.push_back(jht);
852  slc.slHits[jht].InTraj = -4;
853  if (kwire > hiWire) hiWire = kwire;
854  if (kwire < loWire) loWire = kwire;
855  hitsAdded = true;
856  } // jht
857  } // kht
858  } // kwire
859  if (!hitsAdded) break;
860  ++nit;
861  ++hiWire;
862  if (hiWire >= slc.nWires[plane]) break;
863  } // nit < 100
864  // clear InTraj
865  for (auto iht : tHits)
866  slc.slHits[iht].InTraj = 0;
867  if (tHits.size() < 3) continue;
868  if (prt) {
869  mf::LogVerbatim myprt("TC");
870  myprt << "FJT: tHits";
871  for (auto tht : tHits)
872  myprt << " " << PrintHit(slc.slHits[tht]);
873  myprt << "\n";
874  } // prt
875  // See if this is a ghost trajectory
876  if (IsGhost(slc, tHits)) break;
877  if (!MakeJunkTraj(slc, tHits)) {
878  if (prt) mf::LogVerbatim() << "FJT: MakeJunkTraj failed";
879  break;
880  }
881  if (slc.slHits[jht].InTraj > 0) break;
882  } // jht
883  } // iht
884  } // iwire
885  } // FindJunkTraj
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
tca::TrajHitsOK
bool TrajHitsOK(TCSlice const &slc, const std::vector< unsigned int > &iHitsInMultiplet, const std::vector< unsigned int > &jHitsInMultiplet)
Definition: Utils.cxx:1838
tca::IsGhost
bool IsGhost(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:3084
mf::MaybeLogger_
Definition: MessageLogger.h:113
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::MakeJunkTraj
bool MakeJunkTraj(TCSlice &slc, std::vector< unsigned int > tHits)
Definition: StepUtils.cxx:4146
tca::DebugStuff::Hit
unsigned int Hit
set to the hit index in evt.allHits if a Plane:Wire:Tick match is found
Definition: DebugStruct.h:26
tca::HitSep2
float HitSep2(const TCSlice &slc, unsigned int iht, unsigned int jht)
Definition: Utils.cxx:2482
tca::debug
DebugStuff debug
Definition: DebugStruct.cxx:4
tca::GetHitMultiplet
void GetHitMultiplet(const TCSlice &slc, unsigned int theHit, std::vector< unsigned int > &hitsInMultiplet, bool useLongPulseHits)
Definition: StepUtils.cxx:1521
tca::PrintHit
std::string PrintHit(const TCHit &tch)
Definition: Utils.cxx:6344
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:600
tca::TCConfig::useAlg
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Definition: DataStructs.h:579
tca::DecodeCTP
geo::PlaneID DecodeCTP(CTP_t CTP)
Definition: DataStructs.cxx:100
tca::TCConfig::JTMaxHitSep2
float JTMaxHitSep2
Definition: DataStructs.h:578
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:526
void tca::TrajClusterAlg::FinishEvent ( )

Definition at line 1278 of file TrajClusterAlg.cxx.

References tca::MakePFPTjs(), tca::PFPVertexCheck(), tca::slices, and tca::StitchPFPs().

Referenced by cluster::TrajCluster::produce().

1279  {
1280  // final steps that involve correlations between slices
1281  // Stitch PFParticles between TPCs
1282 
1283  // define the PFP TjUIDs vector before calling StitchPFPs
1284  for (auto& slc : slices) {
1285  if (!slc.isValid) continue;
1286  MakePFPTjs(slc);
1287  for (auto& pfp : slc.pfps) {
1288  if (pfp.ID <= 0) continue;
1289  pfp.TjUIDs.resize(pfp.TjIDs.size());
1290  for (unsigned short ii = 0; ii < pfp.TjIDs.size(); ++ii) {
1291  // do a sanity check while we are here
1292  if (pfp.TjIDs[ii] <= 0 || pfp.TjIDs[ii] > (int)slc.tjs.size()) {
1293  std::cout << "FinishEvent found an invalid T" << pfp.TjIDs[ii] << " in P" << pfp.UID
1294  << "\n";
1295  slc.isValid = false;
1296  continue;
1297  } // sanity check
1298  auto& tj = slc.tjs[pfp.TjIDs[ii] - 1];
1299  pfp.TjUIDs[ii] = tj.UID;
1300  } // ii
1301  } // pfp
1302  } // slc
1303 
1304  StitchPFPs();
1305  // Ensure that all PFParticles have a start vertex
1306  for (auto& slc : slices)
1307  PFPVertexCheck(slc);
1308  } // FinishEvent
tca::MakePFPTjs
void MakePFPTjs(TCSlice &slc)
Definition: PFPUtils.cxx:507
tca::PFPVertexCheck
void PFPVertexCheck(TCSlice &slc)
Definition: PFPUtils.cxx:2825
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
tca::StitchPFPs
void StitchPFPs()
Definition: PFPUtils.cxx:41
std::vector<std::string> const& tca::TrajClusterAlg::GetAlgBitNames ( ) const
inline

Definition at line 69 of file TrajClusterAlg.h.

References tca::AlgBitNames.

Referenced by cluster::TrajCluster::endJob().

69 { return AlgBitNames; }
tca::AlgBitNames
const std::vector< std::string > AlgBitNames
Definition: DataStructs.cxx:16
std::vector<unsigned int> const& tca::TrajClusterAlg::GetAlgModCount ( ) const
inline

Definition at line 68 of file TrajClusterAlg.h.

Referenced by cluster::TrajCluster::endJob().

68 { return fAlgModCount; }
tca::TrajClusterAlg::fAlgModCount
std::vector< unsigned int > fAlgModCount
Definition: TrajClusterAlg.h:88
TCSlice const& tca::TrajClusterAlg::GetSlice ( unsigned short  sliceIndex) const
inline

Definition at line 63 of file TrajClusterAlg.h.

References tca::slices.

Referenced by cluster::TrajCluster::produce().

63 { return slices[sliceIndex]; }
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
unsigned short tca::TrajClusterAlg::GetSlicesSize ( ) const
inline

Definition at line 62 of file TrajClusterAlg.h.

References tca::slices.

Referenced by cluster::TrajCluster::produce().

62 { return slices.size(); }
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
void tca::TrajClusterAlg::MergeTPHits ( std::vector< unsigned int > &  tpHits,
std::vector< recob::Hit > &  newHitCol,
std::vector< unsigned int > &  newHitAssns 
) const

Definition at line 984 of file TrajClusterAlg.cxx.

References tca::TCEvent::allHits, tca::evt, MergeTPHitsOnWire(), util::size(), and tmp.

Referenced by cluster::TrajCluster::produce().

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
tmp
Float_t tmp
Definition: plot.C:35
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
tca::TrajClusterAlg::MergeTPHitsOnWire
recob::Hit MergeTPHitsOnWire(std::vector< unsigned int > &tpHits) const
Definition: TrajClusterAlg.cxx:1039
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
tca::TCEvent::allHits
std::vector< recob::Hit > const * allHits
Definition: DataStructs.h:615
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
recob::Hit tca::TrajClusterAlg::MergeTPHitsOnWire ( std::vector< unsigned int > &  tpHits) const
private

Definition at line 1039 of file TrajClusterAlg.cxx.

References tca::TCEvent::allHits, tca::evt, and util::size().

Referenced by MergeTPHits().

1040  {
1041  // merge the hits indexed by tpHits into one hit
1042 
1043  if (tpHits.empty()) return recob::Hit();
1044 
1045  // no merge required. Just return a slightly modified copy of the single hit
1046  if (tpHits.size() == 1) {
1047  if (tpHits[0] > (*evt.allHits).size() - 1) {
1048  std::cout << "MergeTPHits Bad input hit index " << tpHits[0] << " allHits size "
1049  << (*evt.allHits).size() << "\n";
1050  return recob::Hit();
1051  }
1052  auto& oldHit = (*evt.allHits)[tpHits[0]];
1053  raw::TDCtick_t startTick = oldHit.PeakTime() - 3 * oldHit.RMS();
1054  raw::TDCtick_t endTick = oldHit.PeakTime() + 3 * oldHit.RMS();
1055 
1056  return recob::Hit(oldHit.Channel(),
1057  startTick,
1058  endTick,
1059  oldHit.PeakTime(),
1060  oldHit.SigmaPeakTime(),
1061  oldHit.RMS(),
1062  oldHit.PeakAmplitude(),
1063  oldHit.SigmaPeakAmplitude(),
1064  oldHit.SummedADC(),
1065  oldHit.Integral(),
1066  oldHit.SigmaIntegral(),
1067  1,
1068  0, // Multiplicity, LocalIndex
1069  1,
1070  0, // GoodnessOfFit, DOF
1071  oldHit.View(),
1072  oldHit.SignalType(),
1073  oldHit.WireID());
1074  } // tpHits.size() == 1
1075 
1076  double integral = 0;
1077  double sIntegral = 0;
1078  double peakTime = 0;
1079  double sPeakTime = 0;
1080  double peakAmp = 0;
1081  double sPeakAmp = 0;
1082  float sumADC = 0;
1083  raw::TDCtick_t startTick = INT_MAX;
1084  raw::TDCtick_t endTick = 0;
1085  for (auto allHitsIndex : tpHits) {
1086  if (allHitsIndex > (*evt.allHits).size() - 1) return recob::Hit();
1087  auto& hit = (*evt.allHits)[allHitsIndex];
1088  if (hit.StartTick() < startTick) startTick = hit.StartTick();
1089  if (hit.EndTick() > endTick) endTick = hit.EndTick();
1090  double intgrl = hit.Integral();
1091  sPeakTime += intgrl * hit.SigmaPeakTime();
1092  sPeakAmp += intgrl * hit.SigmaPeakAmplitude();
1093  sumADC += hit.SummedADC();
1094  integral += intgrl;
1095  sIntegral += intgrl * hit.SigmaIntegral();
1096  // Get the charge normalization from an input hit
1097  } // tpHit
1098  if (integral <= 0) {
1099  std::cout << "MergeTPHits found bad hit integral " << integral << "\n";
1100  return recob::Hit();
1101  }
1102 
1103  // Create a signal shape vector to find the rms and peak tick
1104  std::vector<double> shape(endTick - startTick + 1, 0.);
1105  for (auto allHitsIndex : tpHits) {
1106  auto& hit = (*evt.allHits)[allHitsIndex];
1107  double peakTick = hit.PeakTime();
1108  double rms = hit.RMS();
1109  double peakAmp = hit.PeakAmplitude();
1110  for (raw::TDCtick_t tick = startTick; tick <= endTick; ++tick) {
1111  double arg = ((double)tick - peakTick) / rms;
1112  unsigned short indx = tick - startTick;
1113  shape[indx] += peakAmp * exp(-0.5 * arg * arg);
1114  } // tick
1115  } // allHitsIndex
1116 
1117  // find the peak tick
1118  double sigsum = 0;
1119  double sigsumt = 0;
1120  for (raw::TDCtick_t tick = startTick; tick <= endTick; ++tick) {
1121  unsigned short indx = tick - startTick;
1122  sigsum += shape[indx];
1123  sigsumt += shape[indx] * tick;
1124  } // tick
1125 
1126  peakTime = sigsumt / sigsum;
1127  // Use the sigma peak time calculated in the first loop
1128  sPeakTime /= integral;
1129 
1130  // find the RMS
1131  sigsumt = 0.;
1132  for (raw::TDCtick_t tick = startTick; tick <= endTick; ++tick) {
1133  double dTick = tick - peakTime;
1134  unsigned short indx = tick - startTick;
1135  sigsumt += shape[indx] * dTick * dTick;
1136  }
1137  double rms = std::sqrt(sigsumt / sigsum);
1138  // get a reference to the first hit to get the charge normalization, channel, view, etc
1139  auto& firstHit = (*evt.allHits)[tpHits[0]];
1140  double chgNorm = 2.507 * firstHit.PeakAmplitude() * firstHit.RMS() / firstHit.Integral();
1141  // find the amplitude from the integrated charge and the RMS
1142  peakAmp = (float)(integral * chgNorm / (2.507 * rms));
1143  // Use the sigma integral calculated in the first loop
1144  sPeakAmp /= integral;
1145  // reset the start and end tick
1146  startTick = peakTime - 3 * rms;
1147  endTick = peakTime + 3 * rms;
1148 
1149  // construct the hit
1150  return recob::Hit(firstHit.Channel(),
1151  startTick,
1152  endTick,
1153  peakTime,
1154  sPeakTime,
1155  rms,
1156  peakAmp,
1157  sPeakAmp,
1158  sumADC,
1159  integral,
1160  sIntegral,
1161  1,
1162  0, // Multiplicity, LocalIndex
1163  1,
1164  0, // GoodnessOfFit, DOF
1165  firstHit.View(),
1166  firstHit.SignalType(),
1167  firstHit.WireID());
1168 
1169  } // MergeTPHits
raw::TDCtick_t
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
util::quantities::tick
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:73
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
tca::TCEvent::allHits
std::vector< recob::Hit > const * allHits
Definition: DataStructs.h:615
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
void tca::TrajClusterAlg::ReconstructAllTraj ( detinfo::DetectorPropertiesData const &  detProp,
TCSlice &  slc,
CTP_t  inCTP 
)
private

Definition at line 432 of file TrajClusterAlg.cxx.

References tca::AddHits(), tca::Trajectory::AlgMod, tca::TCEvent::allHits, tca::AttachAnyTrajToVertex(), tca::Trajectory::AveChg, tca::TCEvent::aveHitRMS, tca::BraggSplit(), tca::CheckTraj(), tca::CheckTrajBeginChg(), tca::ChkVtxAssociations(), tca::Trajectory::CTP, tca::TCConfig::dbg2V, tca::TCConfig::dbgSlc, tca::TCConfig::dbgStp, tca::debug, tca::DecodeCTP(), tca::EndMerge(), tca::Trajectory::EndPt, tca::TCEvent::eventsProcessed, tca::evt, tca::Find2DVertices(), FindJunkTraj(), tca::TCSlice::firstWire, tca::GetHitMultiplet(), tca::DebugStuff::Hit, tca::TCConfig::hitErrFac, tca::HitsPosTick(), tca::HitsRMSTick(), tca::Trajectory::ID, if(), tca::InTrajOK(), tca::Trajectory::IsGood, tca::TCSlice::isValid, tca::kDebug, tca::kKilled, tca::kLastEndMerge, tca::kNormal, tca::kRvPrp, tca::kStepDir, tca::kUnusedHits, tca::LastEndMerge(), tca::TCSlice::lastWire, tca::LongPulseHit(), tca::MakeHaloTj(), tca::MakeJunkVertices(), tca::TCConfig::minPts, tca::TCConfig::minPtsFit, tca::TCConfig::modes, tca::TCConfig::muonTag, tca::NumPtsWithCharge(), tca::TCSlice::nWires, tca::Trajectory::ParentID, tca::Trajectory::Pass, geo::PlaneID::Plane, tca::PrintHit(), tca::PrintPos(), tca::PrintTrajectory(), tca::Trajectory::Pts, tca::ReleaseHits(), tca::seeds, tca::SetTPEnvironment(), tca::TCConfig::showerTag, tca::TCSlice::slHits, tca::slices, tca::SplitTrajCrossingVertices(), tca::StartTraj(), tca::StepAway(), tca::Trajectory::StepDir, tca::StoreTraj(), tca::Trajectory::Strategy, tca::TagShowerLike(), tca::tcc, tca::TCSlice::tjs, tca::TrajHitsOK(), tca::TCConfig::unitsPerTick, tca::UpdateVxEnvironment(), tca::TCConfig::useAlg, tca::TCSlice::vtxs, tca::TCSlice::wireHitRange, and tca::TCEvent::WorkID.

Referenced by RunTrajClusterAlg().

435  {
436  // Reconstruct trajectories in inCTP and put them in allTraj
437 
438  unsigned short plane = DecodeCTP(inCTP).Plane;
439  if (slc.firstWire[plane] > slc.nWires[plane]) return;
440  unsigned int nwires = slc.lastWire[plane] - slc.firstWire[plane] - 1;
441  if (nwires > slc.nWires[plane]) return;
442 
443  // Make several passes through the hits with user-specified cuts for each
444  // pass. In general these are to not reconstruct large angle trajectories on
445  // the first pass
446  float maxHitsRMS = 4 * evt.aveHitRMS[plane];
447  for (unsigned short pass = 0; pass < tcc.minPtsFit.size(); ++pass) {
448  for (unsigned int ii = 0; ii < nwires; ++ii) {
449  // decide which way to step given the sign of StepDir
450  unsigned int iwire = 0;
451  unsigned int jwire = 0;
452  if (tcc.modes[kStepDir]) {
453  // step DS
454  iwire = slc.firstWire[plane] + ii;
455  jwire = iwire + 1;
456  }
457  else {
458  // step US
459  iwire = slc.lastWire[plane] - ii - 1;
460  jwire = iwire - 1;
461  }
462  if (iwire > slc.wireHitRange[plane].size() - 1) continue;
463  if (jwire > slc.wireHitRange[plane].size() - 1) continue;
464  // skip bad wires or no hits on the wire
465  if (slc.wireHitRange[plane][iwire].first == UINT_MAX) continue;
466  if (slc.wireHitRange[plane][jwire].first == UINT_MAX) continue;
467  unsigned int ifirsthit = slc.wireHitRange[plane][iwire].first;
468  unsigned int ilasthit = slc.wireHitRange[plane][iwire].second;
469  unsigned int jfirsthit = slc.wireHitRange[plane][jwire].first;
470  unsigned int jlasthit = slc.wireHitRange[plane][jwire].second;
471  if (ifirsthit > slc.slHits.size() || ilasthit > slc.slHits.size()) {
472  std::cout << "RAT: bad hit range " << ifirsthit << " " << ilasthit << " size "
473  << slc.slHits.size() << " inCTP " << inCTP << "\n";
474  return;
475  }
476  for (unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
477  tcc.dbgStp = (tcc.modes[kDebug] && (slc.slHits[iht].allHitsIndex == debug.Hit));
478  if (tcc.dbgStp) {
479  mf::LogVerbatim("TC") << "+++++++ Pass " << pass << " Found debug hit "
480  << slices.size() - 1 << ":" << PrintHit(slc.slHits[iht])
481  << " iht " << iht;
482  }
483  // Only consider hits that are available
484  if (slc.slHits[iht].InTraj != 0) continue;
485  // We hope to make a trajectory point at the hit position of iht in WSE units
486  // with a direction pointing to jht
487  if (slc.slHits[iht].allHitsIndex > (*evt.allHits).size() - 1) {
488  std::cout << "RAT: Bad allHitsIndex\n";
489  continue;
490  }
491  auto& iHit = (*evt.allHits)[slc.slHits[iht].allHitsIndex];
492  if (LongPulseHit(iHit)) continue;
493  unsigned int fromWire = iHit.WireID().Wire;
494  float fromTick = iHit.PeakTime();
495  float iqtot = iHit.Integral();
496  float hitsRMSTick = iHit.RMS();
497  std::vector<unsigned int> iHitsInMultiplet;
498  if (pass == 0) {
499  // only use the hit on the first pass
500  iHitsInMultiplet.resize(1);
501  iHitsInMultiplet[0] = iht;
502  }
503  else {
504  GetHitMultiplet(slc, iht, iHitsInMultiplet, false);
505  if (iHitsInMultiplet.empty()) continue;
506  // ignore very high multiplicities
507  if (iHitsInMultiplet.size() > 4) continue;
508  }
509  if (iHitsInMultiplet.size() > 1) {
510  fromTick = HitsPosTick(slc, iHitsInMultiplet, iqtot, kUnusedHits);
511  hitsRMSTick = HitsRMSTick(slc, iHitsInMultiplet, kUnusedHits);
512  }
513  if (hitsRMSTick == 0) continue;
514  bool fatIHit = (hitsRMSTick > maxHitsRMS);
515  if (tcc.dbgStp)
516  mf::LogVerbatim("TC") << " hit RMS " << iHit.RMS() << " BB Multiplicity "
517  << iHitsInMultiplet.size() << " AveHitRMS[" << plane << "] "
518  << evt.aveHitRMS[plane] << " HitsRMSTick " << hitsRMSTick
519  << " fatIHit " << fatIHit;
520  for (unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
521  // Only consider hits that are available
522  if (slc.slHits[iht].InTraj != 0) break;
523  if (slc.slHits[jht].InTraj != 0) continue;
524  // clear out any leftover work inTraj's that weren't cleaned up properly
525  for (auto& slHit : slc.slHits) {
526  if (slHit.InTraj < 0) {
527  std::cout << "RAT: Dirty hit " << PrintHit(slHit) << " EventsProcessed "
528  << evt.eventsProcessed << " WorkID " << evt.WorkID << "\n";
529  slHit.InTraj = 0;
530  }
531  }
532  unsigned int toWire = jwire;
533  auto& jHit = (*evt.allHits)[slc.slHits[jht].allHitsIndex];
534  if (LongPulseHit(jHit)) continue;
535  float toTick = jHit.PeakTime();
536  float jqtot = jHit.Integral();
537  std::vector<unsigned int> jHitsInMultiplet;
538  if (pass == 0) {
539  // only use the hit on the first pass
540  jHitsInMultiplet.resize(1);
541  jHitsInMultiplet[0] = jht;
542  }
543  else {
544  GetHitMultiplet(slc, jht, jHitsInMultiplet, false);
545  if (jHitsInMultiplet.empty()) continue;
546  // ignore very high multiplicities
547  if (jHitsInMultiplet.size() > 4) continue;
548  }
549  hitsRMSTick = HitsRMSTick(slc, jHitsInMultiplet, kUnusedHits);
550  if (hitsRMSTick == 0) continue;
551  bool fatJHit = (hitsRMSTick > maxHitsRMS);
552  if (pass == 0) {
553  // require both hits to be consistent
554  if ((fatIHit && !fatJHit) || (!fatIHit && fatJHit)) { continue; }
555  }
556  else {
557  // pass > 0
558  if (jHitsInMultiplet.size() > 1)
559  toTick = HitsPosTick(slc, jHitsInMultiplet, jqtot, kUnusedHits);
560  }
561  bool hitsOK = TrajHitsOK(slc, iHitsInMultiplet, jHitsInMultiplet);
562  // Ensure that the hits StartTick and EndTick have the proper overlap
563  if (!hitsOK) continue;
564  // start a trajectory with direction from iht -> jht
565  Trajectory work;
566  if (!StartTraj(work, fromWire, fromTick, toWire, toTick, inCTP, pass)) continue;
567  // check for a major failure
568  if (!slc.isValid) {
569  std::cout << "RAT: StartTraj major failure\n";
570  return;
571  }
572  if (work.Pts.empty()) {
573  if (tcc.dbgStp) mf::LogVerbatim("TC") << "ReconstructAllTraj: StartTraj failed";
574  continue;
575  }
576  work.Pts[0].DeltaRMS = tcc.hitErrFac * tcc.unitsPerTick * hitsRMSTick;
577  // don't include the charge of iht since it will be low if this
578  // is a starting/ending track
579  work.AveChg = jqtot;
580  // try to add close hits
581  bool sigOK;
582  AddHits(slc, work, 0, sigOK);
583  // check for a major failure
584  if (!slc.isValid) {
585  std::cout << "RAT: AddHits major failure\n";
586  return;
587  }
588  if (!sigOK || work.Pts[0].Chg == 0) {
589  if (tcc.dbgStp) mf::LogVerbatim("TC") << " No hits at initial trajectory point ";
590  ReleaseHits(slc, work);
591  continue;
592  }
593  // move the TP position to the hit position but don't mess with the direction
594  work.Pts[0].Pos = work.Pts[0].HitPos;
595  // print the header and the first TP
596  if (tcc.dbgStp) PrintTrajectory("RAT", slc, work, USHRT_MAX);
597  // We can't update the trajectory yet because there is only one TP.
598  work.EndPt[0] = 0;
599  // now try stepping away
600  StepAway(slc, work);
601  // check for a major failure
602  if (!slc.isValid) {
603  std::cout << "RAT: StepAway major failure\n";
604  return;
605  }
606  if (tcc.dbgStp)
607  mf::LogVerbatim("TC") << " After first StepAway. IsGood " << work.IsGood;
608  // Check the quality of the work trajectory
609  CheckTraj(slc, work);
610  // check for a major failure
611  if (!slc.isValid) {
612  std::cout << "RAT: CheckTraj major failure\n";
613  return;
614  }
615  if (tcc.dbgStp)
616  mf::LogVerbatim("TC") << "ReconstructAllTraj: After CheckTraj EndPt " << work.EndPt[0]
617  << "-" << work.EndPt[1] << " IsGood " << work.IsGood;
618  if (tcc.dbgStp)
619  mf::LogVerbatim("TC")
620  << "StepAway done: IsGood " << work.IsGood << " NumPtsWithCharge "
621  << NumPtsWithCharge(slc, work, true) << " cut " << tcc.minPts[work.Pass];
622  // decide if the trajectory is long enough for this pass
623  if (!work.IsGood || NumPtsWithCharge(slc, work, true) < tcc.minPts[work.Pass]) {
624  if (tcc.dbgStp)
625  mf::LogVerbatim("TC")
626  << " xxxxxxx Not enough points " << NumPtsWithCharge(slc, work, true)
627  << " minimum " << tcc.minPts[work.Pass] << " or !IsGood";
628  ReleaseHits(slc, work);
629  continue;
630  }
631  if (!StoreTraj(slc, work)) continue;
632  if (tcc.dbgStp) {
633  auto& tj = slc.tjs[slc.tjs.size() - 1];
634  PrintTrajectory("RAT", slc, tj, USHRT_MAX);
635  if (!InTrajOK(slc, "RAT")) {
636  std::cout << "RAT: InTrajOK major failure. " << tj.ID << "\n";
637  return;
638  }
639  } // dbgStp
640  CheckTrajBeginChg(slc, slc.tjs.size() - 1);
641  BraggSplit(slc, slc.tjs.size() - 1);
642  break;
643  } // jht
644  } // iht
645  } // iwire
646 
647  // See if there are any seed trajectory points that were saved before reverse
648  // propagation and try to make Tjs from them
649  for (auto tp : seeds) {
650  unsigned short nAvailable = 0;
651  for (unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
652  if (!tp.UseHit[ii]) continue;
653  unsigned int iht = tp.Hits[ii];
654  if (slc.slHits[iht].InTraj == 0) ++nAvailable;
655  tcc.dbgStp = (tcc.modes[kDebug] && (slc.slHits[iht].allHitsIndex == debug.Hit));
656  if (tcc.dbgStp) {
657  mf::LogVerbatim("TC") << "+++++++ Seed debug hit " << slices.size() - 1 << ":"
658  << PrintHit(slc.slHits[iht]) << " iht " << iht;
659  }
660  } // ii
661  if (nAvailable == 0) continue;
662  Trajectory work;
663  work.ID = evt.WorkID;
664  for (unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
665  if (!tp.UseHit[ii]) continue;
666  unsigned int iht = tp.Hits[ii];
667  if (slc.slHits[iht].InTraj == 0) slc.slHits[iht].InTraj = work.ID;
668  } // ii
669  work.Pass = pass;
670  work.StepDir = 1;
671  if (tp.Dir[0] < 0) work.StepDir = -1;
672  work.CTP = tp.CTP;
673  work.ParentID = -1;
674  work.Strategy.reset();
675  work.Strategy[kNormal] = true;
676  // don't allow yet another reverse propagation
677  work.AlgMod[kRvPrp] = true;
678  work.Pts.push_back(tp);
679  StepAway(slc, work);
680  CheckTraj(slc, work);
681  // check for a major failure
682  if (!slc.isValid) {
683  std::cout << "RAT: CheckTraj major failure\n";
684  return;
685  }
686  // decide if the trajectory is long enough for this pass
687  if (!work.IsGood || NumPtsWithCharge(slc, work, true) < tcc.minPts[work.Pass]) {
688  if (tcc.dbgStp)
689  mf::LogVerbatim("TC") << " xxxxxxx Not enough points "
690  << NumPtsWithCharge(slc, work, true) << " minimum "
691  << tcc.minPts[work.Pass] << " or !IsGood";
692  ReleaseHits(slc, work);
693  continue;
694  }
695  if (!StoreTraj(slc, work)) continue;
696  if (tcc.dbgStp) {
697  auto& tj = slc.tjs[slc.tjs.size() - 1];
698  mf::LogVerbatim("TC") << "TRP RAT Stored T" << tj.ID << " using seed TP " << PrintPos(tp);
699  }
700  BraggSplit(slc, slc.tjs.size() - 1);
701  } // seed
702 
703  seeds.resize(0);
704 
705  bool lastPass = (pass == tcc.minPtsFit.size() - 1);
706  // don't use lastPass cuts if we will use LastEndMerge
707  if (tcc.useAlg[kLastEndMerge]) lastPass = false;
708  EndMerge(slc, inCTP, lastPass);
709  if (!slc.isValid) return;
710 
711  Find2DVertices(detProp, slc, inCTP, pass);
712 
713  } // pass
714 
715  // Last attempt to merge long straight Tjs that failed the EndMerge cuts
716  LastEndMerge(slc, inCTP);
717 
718  // make junk trajectories using nearby un-assigned hits
719  FindJunkTraj(slc, inCTP);
720  // dressed muons with halo trajectories
721  if (tcc.muonTag.size() > 4 && tcc.muonTag[4] > 0) {
722  for (auto& tj : slc.tjs) {
723  if (tj.AlgMod[kKilled]) continue;
724  if (tj.PDGCode != 13) continue;
725  MakeHaloTj(slc, tj, tcc.dbgSlc);
726  } // tj
727  } // dressed muons
728 
729  // Tag ShowerLike Tjs
730  if (tcc.showerTag[0] > 0) TagShowerLike(slc, inCTP);
731  // Set TP Environment bits
732  SetTPEnvironment(slc, inCTP);
733  Find2DVertices(detProp, slc, inCTP, USHRT_MAX);
734  SplitTrajCrossingVertices(slc, inCTP);
735  // Make vertices between long Tjs and junk Tjs
736  MakeJunkVertices(slc, inCTP);
737  // check for a major failure
738  if (!slc.isValid) {
739  std::cout << "RAT: MakeJunkVertices major failure\n";
740  return;
741  }
742 
743  // last attempt to attach Tjs to vertices
744  for (unsigned short ivx = 0; ivx < slc.vtxs.size(); ++ivx) {
745  auto& vx2 = slc.vtxs[ivx];
746  if (vx2.ID == 0) continue;
747  if (vx2.CTP != inCTP) continue;
748  AttachAnyTrajToVertex(slc, ivx, tcc.dbgStp || tcc.dbg2V);
749  } // ivx
750 
751  // Set the kEnvOverlap bit true for all TPs that are close to other
752  // trajectories that are close to vertices. The positions of TPs that
753  // overlap are biased and shouldn't be used in a vertex fit. Also, these
754  // TPs shouldn't be used to calculate dE/dx
755  UpdateVxEnvironment(slc);
756 
757  // Check the Tj <-> vtx associations and define the vertex quality
758  if (!ChkVtxAssociations(slc, inCTP)) {
759  std::cout << "RAT: ChkVtxAssociations found an error. Events processed "
760  << evt.eventsProcessed << " WorkID " << evt.WorkID << "\n";
761  }
762 
763  } // ReconstructAllTraj
tca::CheckTrajBeginChg
void CheckTrajBeginChg(TCSlice &slc, unsigned short itj)
Definition: Utils.cxx:1316
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
tca::StepAway
void StepAway(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:30
tca::TrajHitsOK
bool TrajHitsOK(TCSlice const &slc, const std::vector< unsigned int > &iHitsInMultiplet, const std::vector< unsigned int > &jHitsInMultiplet)
Definition: Utils.cxx:1838
tca::InTrajOK
bool InTrajOK(TCSlice &slc, std::string someText)
Definition: Utils.cxx:1256
tca::TCConfig::minPtsFit
std::vector< unsigned short > minPtsFit
Reconstruct in several passes.
Definition: DataStructs.h:558
tca::AddHits
void AddHits(TCSlice &slc, Trajectory &tj, unsigned short ipt, bool &sigOK)
Definition: StepUtils.cxx:1113
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::PrintTrajectory
void PrintTrajectory(std::string someText, const TCSlice &slc, const Trajectory &tj, unsigned short tPoint)
Definition: Utils.cxx:6032
tca::SplitTrajCrossingVertices
void SplitTrajCrossingVertices(TCSlice &slc, CTP_t inCTP)
Definition: TCVertex.cxx:915
tca::StartTraj
bool StartTraj(TCSlice const &slc, Trajectory &tj, unsigned int fromhit, unsigned int tohit, unsigned short pass)
Definition: Utils.cxx:4855
tca::kStepDir
step from US -> DS (true) or DS -> US (false)
Definition: DataStructs.h:524
tca::StoreTraj
bool StoreTraj(TCSlice &slc, Trajectory &tj)
Definition: Utils.cxx:1073
tca::AttachAnyTrajToVertex
bool AttachAnyTrajToVertex(TCSlice &slc, unsigned short ivx, bool prt)
Definition: TCVertex.cxx:1670
tca::TCConfig::dbgSlc
bool dbgSlc
debug only in the user-defined slice? default is all slices
Definition: DataStructs.h:583
tca::TagShowerLike
void TagShowerLike(TCSlice &slc, const CTP_t &inCTP)
Definition: TCShower.cxx:3232
tca::LongPulseHit
bool LongPulseHit(const recob::Hit &hit)
Definition: Utils.cxx:4315
tca::TCEvent::eventsProcessed
unsigned int eventsProcessed
Definition: DataStructs.h:630
tca::PrintPos
std::string PrintPos(const TrajPoint &tp)
Definition: Utils.cxx:6353
tca::TCConfig::showerTag
std::vector< float > showerTag
shower-like trajectory tagging + shower reconstruction
Definition: DataStructs.h:551
tca::DebugStuff::Hit
unsigned int Hit
set to the hit index in evt.allHits if a Plane:Wire:Tick match is found
Definition: DebugStruct.h:26
tca::HitsRMSTick
float HitsRMSTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, HitStatus_t hitRequest)
Definition: Utils.cxx:4131
if
if(nlines<=0)
Definition: comparison_ascii.C:19
tca::TCConfig::dbgStp
bool dbgStp
debug stepping using debug.Cryostat, debug.TPC, etc
Definition: DataStructs.h:584
tca::TCConfig::unitsPerTick
float unitsPerTick
scale factor from Tick to WSE equivalent units
Definition: DataStructs.h:564
tca::BraggSplit
bool BraggSplit(TCSlice &slc, unsigned short itj)
Definition: Utils.cxx:1422
tca::debug
DebugStuff debug
Definition: DebugStruct.cxx:4
tca::TCConfig::dbg2V
bool dbg2V
debug 2D vertex finding
Definition: DataStructs.h:586
tca::GetHitMultiplet
void GetHitMultiplet(const TCSlice &slc, unsigned int theHit, std::vector< unsigned int > &hitsInMultiplet, bool useLongPulseHits)
Definition: StepUtils.cxx:1521
tca::TCEvent::aveHitRMS
std::vector< float > aveHitRMS
average RMS of an isolated hit
Definition: DataStructs.h:631
tca::ChkVtxAssociations
bool ChkVtxAssociations(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:2050
tca::MakeHaloTj
void MakeHaloTj(TCSlice &slc, Trajectory &muTj, bool prt)
Definition: Utils.cxx:52
tca::PrintHit
std::string PrintHit(const TCHit &tch)
Definition: Utils.cxx:6344
tca::TCConfig::minPts
std::vector< unsigned short > minPts
min number of Pts required to make a trajectory
Definition: DataStructs.h:559
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
tca::ReleaseHits
void ReleaseHits(TCSlice &slc, Trajectory const &tj)
Definition: Utils.cxx:1048
tca::HitsPosTick
float HitsPosTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, float &sum, HitStatus_t hitRequest)
Definition: Utils.cxx:4172
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:600
tca::seeds
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:14
tca::CheckTraj
void CheckTraj(TCSlice &slc, Trajectory &tj)
Definition: StepUtils.cxx:997
tca::TCConfig::useAlg
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Definition: DataStructs.h:579
tca::LastEndMerge
void LastEndMerge(TCSlice &slc, CTP_t inCTP)
Definition: StepUtils.cxx:3339
tca::TCConfig::muonTag
std::vector< short > muonTag
Definition: DataStructs.h:548
tca::UpdateVxEnvironment
void UpdateVxEnvironment(TCSlice &slc)
Definition: Utils.cxx:3764
tca::DecodeCTP
geo::PlaneID DecodeCTP(CTP_t CTP)
Definition: DataStructs.cxx:100
tca::TCEvent::allHits
std::vector< recob::Hit > const * allHits
Definition: DataStructs.h:615
tca::EndMerge
void EndMerge(TCSlice &slc, CTP_t inCTP, bool lastPass)
Definition: StepUtils.cxx:3486
tca::NumPtsWithCharge
unsigned short NumPtsWithCharge(const TCSlice &slc, const Trajectory &tj, bool includeDeadWires)
Definition: Utils.cxx:2071
tca::MakeJunkVertices
void MakeJunkVertices(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:30
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::Find2DVertices
void Find2DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, const CTP_t &inCTP, unsigned short pass)
Definition: TCVertex.cxx:130
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:526
tca::SetTPEnvironment
void SetTPEnvironment(TCSlice &slc, CTP_t inCTP)
Definition: Utils.cxx:3530
tca::TrajClusterAlg::FindJunkTraj
void FindJunkTraj(TCSlice &slc, CTP_t inCTP)
Definition: TrajClusterAlg.cxx:766
void tca::TrajClusterAlg::RunTrajClusterAlg ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< unsigned int > &  hitsInSlice,
int  sliceID 
)

Definition at line 349 of file TrajClusterAlg.cxx.

References tca::AlgBitNames, tca::TCEvent::aveHitRMS, tca::ChkVtxAssociations(), art::errors::Configuration, CreateSlice(), tca::TCConfig::dbgDump, tca::DefinePFPParents(), tca::DefineTjParents(), tca::DumpTj(), tca::EncodeCTP(), tca::TCEvent::eventsProcessed, tca::evt, fAlgModCount, tca::Find3DVertices(), tca::FindPFParticles(), tca::FindShowers3D(), tca::Finish3DShowers(), if(), tca::kDebug, tca::KillPoorVertices(), tca::kSaveShowerTree, tca::TCConfig::match3DCuts, tca::TCConfig::modes, tca::Reconcile2Vs(), ReconstructAllTraj(), tca::TCConfig::recoSlice, tca::TCConfig::recoTPC, tca::ScoreVertices(), tca::seeds, tca::TCConfig::showerTag, showertree, tca::slices, tca::ShowerTreeVars::StageNum, tca::stv, and tca::tcc.

Referenced by cluster::TrajCluster::produce().

353  {
354  // Reconstruct everything using the hits in a slice
355 
356  if (slices.empty()) ++evt.eventsProcessed;
357  if (hitsInSlice.size() < 2) return;
358  if (tcc.recoSlice > 0 && sliceID != tcc.recoSlice) return;
359 
360  if (!CreateSlice(clockData, detProp, hitsInSlice, sliceID)) return;
361 
362  seeds.resize(0);
363  // get a reference to the stored slice
364  auto& slc = slices[slices.size() - 1];
365  // special debug mode reconstruction
366  if (tcc.recoTPC > 0 && (short)slc.TPCID.TPC != tcc.recoTPC) {
367  slices.pop_back();
368  return;
369  }
370 
371  if (evt.aveHitRMS.size() != slc.nPlanes)
372  throw art::Exception(art::errors::Configuration)
373  << " AveHitRMS vector size != the number of planes ";
374  if (tcc.recoSlice)
375  std::cout << "Reconstruct " << hitsInSlice.size() << " hits in Slice " << sliceID
376  << " in TPC " << slc.TPCID.TPC << "\n";
377  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
378  CTP_t inCTP = EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
379  ReconstructAllTraj(detProp, slc, inCTP);
380  if (!slc.isValid) return;
381  } // plane
382  // Compare 2D vertices in each plane and try to reconcile T -> 2V attachments using
383  // 2D and 3D(?) information
384  Reconcile2Vs(slc);
385  Find3DVertices(detProp, slc);
386  ScoreVertices(slc);
387  // Define the ParentID of trajectories using the vertex score
388  DefineTjParents(slc, false);
389  for (unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
390  CTP_t inCTP = EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
391  if (!ChkVtxAssociations(slc, inCTP)) {
392  std::cout << "RTC: ChkVtxAssociations found an error\n";
393  }
394  } // plane
395  if (tcc.match3DCuts[0] > 0) {
396  FindPFParticles(clockData, detProp, slc);
397  DefinePFPParents(slc);
398  } // 3D matching requested
399  KillPoorVertices(slc);
400  // Use 3D matching information to find showers in 2D. FindShowers3D returns
401  // true if the algorithm was successful indicating that the matching needs to be redone
402  if (tcc.showerTag[0] == 2 || tcc.showerTag[0] == 4) {
403  FindShowers3D(detProp, slc);
404  if (tcc.modes[kSaveShowerTree]) {
405  std::cout << "SHOWER TREE STAGE NUM SIZE: " << stv.StageNum.size() << std::endl;
406  showertree->Fill();
407  }
408  } // 3D shower code
409 
410  if (!slc.isValid) {
411  mf::LogVerbatim("TC") << "RunTrajCluster failed in MakeAllTrajClusters";
412  return;
413  }
414 
415  // dump a trajectory?
416  if (tcc.modes[kDebug] && tcc.dbgDump) DumpTj();
417 
418  Finish3DShowers(slc);
419 
420  // count algorithm usage
421  for (auto& tj : slc.tjs) {
422  for (unsigned short ib = 0; ib < AlgBitNames.size(); ++ib)
423  if (tj.AlgMod[ib]) ++fAlgModCount[ib];
424  } // tj
425 
426  // clear vectors that are not needed later
427  slc.mallTraj.resize(0);
428 
429  } // RunTrajClusterAlg
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
tca::AlgBitNames
const std::vector< std::string > AlgBitNames
Definition: DataStructs.cxx:16
tca::TCConfig::recoTPC
short recoTPC
only reconstruct in the seleted TPC
Definition: DataStructs.h:582
tca::Find3DVertices
void Find3DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: TCVertex.cxx:1257
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::Reconcile2Vs
void Reconcile2Vs(TCSlice &slc)
Definition: TCVertex.cxx:1058
tca::FindShowers3D
bool FindShowers3D(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: TCShower.cxx:285
tca::TrajClusterAlg::CreateSlice
bool CreateSlice(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
Definition: TrajClusterAlg.cxx:1213
tca::TCEvent::eventsProcessed
unsigned int eventsProcessed
Definition: DataStructs.h:630
tca::TCConfig::showerTag
std::vector< float > showerTag
shower-like trajectory tagging + shower reconstruction
Definition: DataStructs.h:551
tca::stv
ShowerTreeVars stv
Definition: DataStructs.cxx:11
tca::TrajClusterAlg::fAlgModCount
std::vector< unsigned int > fAlgModCount
Definition: TrajClusterAlg.h:88
tca::kSaveShowerTree
save shower tree
Definition: DataStructs.h:533
if
if(nlines<=0)
Definition: comparison_ascii.C:19
tca::DumpTj
void DumpTj()
Definition: Utils.cxx:5248
tca::TCEvent::aveHitRMS
std::vector< float > aveHitRMS
average RMS of an isolated hit
Definition: DataStructs.h:631
tca::ChkVtxAssociations
bool ChkVtxAssociations(TCSlice &slc, const CTP_t &inCTP)
Definition: TCVertex.cxx:2050
tca::TCConfig::match3DCuts
std::vector< float > match3DCuts
3D matching cuts
Definition: DataStructs.h:553
tca::ScoreVertices
void ScoreVertices(TCSlice &slc)
Definition: TCVertex.cxx:2124
tca::DefineTjParents
void DefineTjParents(TCSlice &slc, bool prt)
Definition: Utils.cxx:172
tca::slices
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:600
art::Exception
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
tca::seeds
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:14
tca::CTP_t
unsigned int CTP_t
Definition: DataStructs.h:47
tca::ShowerTreeVars::StageNum
std::vector< int > StageNum
Definition: DataStructs.h:398
tca::EncodeCTP
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
Definition: DataStructs.h:51
tca::KillPoorVertices
void KillPoorVertices(TCSlice &slc)
Definition: TCVertex.cxx:2151
tca::Finish3DShowers
void Finish3DShowers(TCSlice &slc)
Definition: TCShower.cxx:153
tca::FindPFParticles
void FindPFParticles(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
Definition: PFPUtils.cxx:189
tca::TrajClusterAlg::ReconstructAllTraj
void ReconstructAllTraj(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, CTP_t inCTP)
Definition: TrajClusterAlg.cxx:432
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::TCConfig::recoSlice
short recoSlice
only reconstruct the slice with ID (0 = all)
Definition: DataStructs.h:581
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:526
tca::DefinePFPParents
void DefinePFPParents(TCSlice &slc)
Definition: PFPUtils.cxx:2863
bool tca::TrajClusterAlg::SetInputHits ( std::vector< recob::Hit > const &  inputHits,
unsigned int  run,
unsigned int  event 
)

Definition at line 305 of file TrajClusterAlg.cxx.

References tca::TCEvent::allHits, tca::AnalyzeHits(), ClearResults(), tca::TCEvent::event, tca::TCEvent::eventsProcessed, tca::evt, tca::TCConfig::geom, tca::TCEvent::global2S_UID, tca::TCEvent::global2V_UID, tca::TCEvent::global3S_UID, tca::TCEvent::global3V_UID, tca::TCEvent::globalP_UID, tca::TCEvent::globalT_UID, tca::kDebug, tca::TCConfig::modes, tca::PrintDebugMode(), tca::TCEvent::run, tca::tcc, and tca::TCEvent::WorkID.

Referenced by cluster::TrajCluster::produce().

308  {
309  // defines the pointer to the input hit collection, analyzes them,
310  // initializes global counters and refreshes service references
311  ClearResults();
312  evt.allHits = &inputHits;
313  evt.run = run;
314  evt.event = event;
315  // refresh service references
316  tcc.geom = lar::providerFrom<geo::Geometry>();
317  evt.WorkID = 0;
318  evt.globalT_UID = 0;
319  evt.global2V_UID = 0;
320  evt.global3V_UID = 0;
321  evt.globalP_UID = 0;
322  evt.global2S_UID = 0;
323  evt.global3S_UID = 0;
324  // find the average hit RMS using the full hit collection and define the
325  // configuration for the current TPC
326 
327  if (tcc.modes[kDebug] && evt.eventsProcessed == 0) PrintDebugMode();
328 
329  return AnalyzeHits();
330  } // SetInputHits
tca::TrajClusterAlg::ClearResults
void ClearResults()
Deletes all the results.
Definition: TrajClusterAlg.h:72
tca::TCEvent::event
unsigned int event
Definition: DataStructs.h:627
tca::TCEvent::run
unsigned int run
Definition: DataStructs.h:628
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::TCEvent::eventsProcessed
unsigned int eventsProcessed
Definition: DataStructs.h:630
tca::TCConfig::geom
const geo::GeometryCore * geom
Definition: DataStructs.h:569
tca::TCConfig::modes
std::bitset< 16 > modes
number of points to find AveChg
Definition: DataStructs.h:600
tca::AnalyzeHits
bool AnalyzeHits()
Definition: Utils.cxx:4259
tca::TCEvent::allHits
std::vector< recob::Hit > const * allHits
Definition: DataStructs.h:615
tca::PrintDebugMode
void PrintDebugMode()
Definition: Utils.cxx:5297
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
tca::kDebug
master switch for turning on debug mode
Definition: DataStructs.h:526
event
Event finding and building.
Definition: EventCheater_module.cc:32
void tca::TrajClusterAlg::SetInputSpts ( std::vector< recob::SpacePoint > const &  sptHandle)
inline

Definition at line 44 of file TrajClusterAlg.h.

References tca::evt, and tca::TCEvent::sptHandle.

Referenced by cluster::TrajCluster::produce().

45  {
46  evt.sptHandle = &sptHandle;
47  }
tca::TCEvent::sptHandle
std::vector< recob::SpacePoint > const * sptHandle
handle to SpacePoints in the event
Definition: DataStructs.h:624
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
void tca::TrajClusterAlg::SetSourceHits ( std::vector< recob::Hit > const &  srcHits)

Definition at line 333 of file TrajClusterAlg.cxx.

References tca::evt, tca::TCConfig::geom, geo::GeometryCore::NTPC(), util::size(), tca::TCEvent::srcHits, tca::tcc, and tca::TCEvent::tpcSrcHitRange.

Referenced by cluster::TrajCluster::produce().

334  {
335  evt.srcHits = &srcHits;
336  evt.tpcSrcHitRange.resize(tcc.geom->NTPC());
337  for (auto& thr : evt.tpcSrcHitRange)
338  thr = {UINT_MAX, UINT_MAX};
339  for (unsigned int iht = 0; iht < (*evt.srcHits).size(); ++iht) {
340  auto& hit = (*evt.srcHits)[iht];
341  unsigned int tpc = hit.WireID().TPC;
342  if (tpc >= evt.tpcSrcHitRange.size()) return;
343  if (evt.tpcSrcHitRange[tpc].first == UINT_MAX) evt.tpcSrcHitRange[tpc].first = iht;
344  evt.tpcSrcHitRange[tpc].second = iht;
345  } // iht
346  } // SetSourceHits
geo::GeometryCore::NTPC
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:686
tca::tcc
TCConfig tcc
Definition: DataStructs.cxx:9
tca::TCEvent::tpcSrcHitRange
std::vector< std::pair< unsigned int, unsigned int > > tpcSrcHitRange
Definition: DataStructs.h:621
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
tca::TCConfig::geom
const geo::GeometryCore * geom
Definition: DataStructs.h:569
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
tca::TCEvent::srcHits
std::vector< recob::Hit > const * srcHits
Definition: DataStructs.h:616
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8

Member Data Documentation

std::vector<unsigned int> tca::TrajClusterAlg::fAlgModCount
private

Definition at line 88 of file TrajClusterAlg.h.

Referenced by ChkInTraj(), RunTrajClusterAlg(), and TrajClusterAlg().

calo::CalorimetryAlg tca::TrajClusterAlg::fCaloAlg
private

Definition at line 85 of file TrajClusterAlg.h.

Referenced by TrajClusterAlg().

TMVA::Reader tca::TrajClusterAlg::fMVAReader
private

Definition at line 86 of file TrajClusterAlg.h.

Referenced by TrajClusterAlg().

TTree* tca::TrajClusterAlg::showertree
private

Definition at line 83 of file TrajClusterAlg.h.

Referenced by DefineShTree(), and RunTrajClusterAlg().

The documentation for this class was generated from the following files: