LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TrackContainmentAlg.cxx
Go to the documentation of this file.
1 #include "TrackContainmentAlg.hh"
2 
5 
6 #include <iostream>
7 
8 #include "TTree.h"
9 
11 
13 {
14  fTrackTree = tfs_tree_trk;
15  fTrackTree->SetObject(fTrackTree->GetName(), "Track Tree");
16  fTrackTree->Branch("run", &fRun);
17  fTrackTree->Branch("event", &fEvent);
18  fTrackTree->Branch("trk", &fTrackTreeObj, fTrackTreeObj.Leaflist().c_str());
19  fTrackTree->Branch("trk_col", &fCollection);
20  fTrackTree->Branch("trk_id", &fTrkID);
21  fTrackTree->Branch("trk_mindist", &fDistance);
22  fTrackTree->Branch("trk_containment", &fContainment);
23 }
24 
26 {
27  fZBuffer = p.get<double>("ZBuffer");
28  fYBuffer = p.get<double>("YBuffer");
29  fXBuffer = p.get<double>("XBuffer");
30  fIsolation = p.get<double>("Isolation");
31 
32  fDebug = p.get<bool>("Debug", false);
33  fMakeCosmicTags = p.get<bool>("MakeCosmicTags", true);
34  fFillOutputTree = p.get<bool>("FillOutputTree", true);
35 }
36 
38 {
39  if (track.Vertex().Z() < (0 + fZBuffer) || track.Vertex().Z() > (geo.DetLength() - fZBuffer))
40  return false;
41  if (track.End().Z() < (0 + fZBuffer) || track.End().Z() > (geo.DetLength() - fZBuffer))
42  return false;
43  if (track.Vertex().Y() < (-1 * geo.DetHalfHeight() + fYBuffer) ||
44  track.Vertex().Y() > (geo.DetHalfHeight() - fYBuffer))
45  return false;
46  if (track.End().Y() < (-1 * geo.DetHalfHeight() + fYBuffer) ||
47  track.End().Y() > (geo.DetHalfHeight() - fYBuffer))
48  return false;
49  if (track.Vertex().X() < (0 + fXBuffer) ||
50  track.Vertex().X() > (2 * geo.DetHalfWidth() - fXBuffer))
51  return false;
52  if (track.End().X() < (0 + fXBuffer) || track.End().X() > (2 * geo.DetHalfWidth() - fXBuffer))
53  return false;
54 
55  return true;
56 }
57 
59  geo::GeometryCore const& geo)
60 {
61 
63 
64  if (track.Vertex().Z() < (0 + fZBuffer) || track.Vertex().Z() > (geo.DetLength() - fZBuffer))
66  else if (track.Vertex().Y() < (-1 * geo.DetHalfHeight() + fYBuffer) ||
67  track.Vertex().Y() > (geo.DetHalfHeight() - fYBuffer))
69  else if ((track.Vertex().X() > 0 && track.Vertex().X() < (0 + fXBuffer)) ||
70  (track.Vertex().X() < 2 * geo.DetHalfWidth() &&
71  track.Vertex().X() > (2 * geo.DetHalfWidth() - fXBuffer)))
73 
74  if (track.End().Z() < (0 + fZBuffer) || track.End().Z() > (geo.DetLength() - fZBuffer)) {
77  else if (id == anab::CosmicTagID_t::kGeometry_Z)
79  else if (id == anab::CosmicTagID_t::kGeometry_Y)
81  else if (id == anab::CosmicTagID_t::kGeometry_X)
83  }
84  else if (track.End().Y() < (-1 * geo.DetHalfHeight() + fYBuffer) ||
85  track.End().Y() > (geo.DetHalfHeight() - fYBuffer)) {
88  else if (id == anab::CosmicTagID_t::kGeometry_Z)
90  else if (id == anab::CosmicTagID_t::kGeometry_Y)
92  else if (id == anab::CosmicTagID_t::kGeometry_X)
94  }
95  else if ((track.End().X() > 0 && track.End().X() < (0 + fXBuffer)) ||
96  (track.End().X() < 2 * geo.DetHalfWidth() &&
97  track.End().X() > (2 * geo.DetHalfWidth() - fXBuffer))) {
100  else if (id == anab::CosmicTagID_t::kGeometry_Z)
102  else if (id == anab::CosmicTagID_t::kGeometry_Y)
104  else if (id == anab::CosmicTagID_t::kGeometry_X)
106  }
107 
108  if (track.Vertex().X() < 0 || track.Vertex().X() > 2 * geo.DetHalfWidth())
110  if (track.End().X() < 0 || track.End().X() > 2 * geo.DetHalfWidth()) {
113  else
115  }
116 
117  return id;
118 }
119 
121  recob::Track const& t_ref)
122 {
123  double min_distance = 9e12;
124  double tmp;
125  for (size_t i_p = 0; i_p < t_ref.NumberTrajectoryPoints(); ++i_p) {
126  tmp = (t_probe.Vertex().X() - t_ref.LocationAtPoint(i_p).X()) *
127  (t_probe.Vertex().X() - t_ref.LocationAtPoint(i_p).X()) +
128  (t_probe.Vertex().Y() - t_ref.LocationAtPoint(i_p).Y()) *
129  (t_probe.Vertex().Y() - t_ref.LocationAtPoint(i_p).Y()) +
130  (t_probe.Vertex().Z() - t_ref.LocationAtPoint(i_p).Z()) *
131  (t_probe.Vertex().Z() - t_ref.LocationAtPoint(i_p).Z());
132  //std::cout << "\t\t\ttmp=" << tmp << std::endl;
133  if (tmp < min_distance) min_distance = tmp;
134  }
135  return std::sqrt(min_distance);
136 }
137 
139  recob::Track const& t_ref)
140 {
141  double min_distance = 9e12;
142  double tmp;
143  for (size_t i_p = 0; i_p < t_ref.NumberTrajectoryPoints(); ++i_p) {
144  tmp = (t_probe.End().X() - t_ref.LocationAtPoint(i_p).X()) *
145  (t_probe.End().X() - t_ref.LocationAtPoint(i_p).X()) +
146  (t_probe.End().Y() - t_ref.LocationAtPoint(i_p).Y()) *
147  (t_probe.End().Y() - t_ref.LocationAtPoint(i_p).Y()) +
148  (t_probe.End().Z() - t_ref.LocationAtPoint(i_p).Z()) *
149  (t_probe.End().Z() - t_ref.LocationAtPoint(i_p).Z());
150  if (tmp < min_distance) min_distance = tmp;
151  }
152  return std::sqrt(min_distance);
153 }
154 
155 void trk::TrackContainmentAlg::SetRunEvent(unsigned int const& run, unsigned int const& event)
156 {
157  fRun = run;
158  fEvent = event;
159 }
160 
162  std::vector<std::vector<recob::Track>> const& tracksVec,
163  geo::GeometryCore const& geo)
164 {
165 
166  if (fDebug) {
167  std::cout << "Geometry:" << std::endl;
168  std::cout << "\t" << geo.DetHalfWidth() << " " << geo.DetHalfHeight() << " " << geo.DetLength()
169  << std::endl;
170  std::cout << "\t z:(" << fZBuffer << "," << geo.DetLength() - fZBuffer << ")"
171  << "\t y:(" << -1. * geo.DetHalfHeight() + fYBuffer << ","
172  << geo.DetHalfHeight() - fYBuffer << ")"
173  << "\t x:(" << 0 + fXBuffer << "," << 2 * geo.DetHalfWidth() - fXBuffer << ")"
174  << std::endl;
175  }
176 
177  int containment_level = 0;
178  bool track_linked = false;
179 
180  fTrackContainmentLevel.clear();
181  fTrackContainmentLevel.resize(tracksVec.size());
182  fMinDistances.clear();
183  fMinDistances.resize(tracksVec.size());
184 
185  fTrackContainmentIndices.clear();
186  fTrackContainmentIndices.push_back(std::vector<std::pair<int, int>>());
187 
188  fCosmicTags.clear();
189  fCosmicTags.resize(tracksVec.size());
190 
191  //first, loop through tracks and see what's not contained
192 
193  for (size_t i_tc = 0; i_tc < tracksVec.size(); ++i_tc) {
194  fTrackContainmentLevel[i_tc].resize(tracksVec[i_tc].size(), -1);
195  fMinDistances[i_tc].resize(tracksVec[i_tc].size(), 9e12);
196  fCosmicTags[i_tc].resize(tracksVec[i_tc].size(), anab::CosmicTag(-1));
197  for (size_t i_t = 0; i_t < tracksVec[i_tc].size(); ++i_t) {
198 
199  if (!IsContained(tracksVec[i_tc][i_t], geo)) {
200  if (!track_linked) track_linked = true;
201  fTrackContainmentLevel[i_tc][i_t] = 0;
202  fTrackContainmentIndices.back().emplace_back(i_tc, i_t);
203  if (fDebug) {
204  std::cout << "\tTrack (" << i_tc << "," << i_t << ")"
205  << " " << containment_level << std::endl;
206  }
207 
208  } //end if contained
209  } //end loop over tracks
210 
211  } //end loop over track collections
212 
213  //now, while we are still linking tracks, loop over all tracks and note anything
214  //close to an uncontained (or linked) track
215  while (track_linked) {
216 
217  track_linked = false;
218  ++containment_level;
219  fTrackContainmentIndices.push_back(std::vector<std::pair<int, int>>());
220 
221  for (size_t i_tc = 0; i_tc < tracksVec.size(); ++i_tc) {
222  for (size_t i_t = 0; i_t < tracksVec[i_tc].size(); ++i_t) {
223  if (fTrackContainmentLevel[i_tc][i_t] >= 0)
224  continue;
225  else {
226  for (auto const& i_tr : fTrackContainmentIndices[containment_level - 1]) {
227 
228  /*
229  if(fDebug){
230  std::cout << "\t\t" << MinDistanceStartPt(tracksVec[i_tc][i_t],tracksVec[i_tr.first][i_tr.second]) << std::endl;
231  std::cout << "\t\t" << MinDistanceEndPt(tracksVec[i_tc][i_t],tracksVec[i_tr.first][i_tr.second]) << std::endl;
232  }
233  */
234 
235  if (MinDistanceStartPt(tracksVec[i_tc][i_t], tracksVec[i_tr.first][i_tr.second]) <
236  fMinDistances[i_tc][i_t])
237  fMinDistances[i_tc][i_t] =
238  MinDistanceStartPt(tracksVec[i_tc][i_t], tracksVec[i_tr.first][i_tr.second]);
239  if (MinDistanceEndPt(tracksVec[i_tc][i_t], tracksVec[i_tr.first][i_tr.second]) <
240  fMinDistances[i_tc][i_t])
241  fMinDistances[i_tc][i_t] =
242  MinDistanceEndPt(tracksVec[i_tc][i_t], tracksVec[i_tr.first][i_tr.second]);
243 
244  if (MinDistanceStartPt(tracksVec[i_tc][i_t], tracksVec[i_tr.first][i_tr.second]) <
245  fIsolation ||
246  MinDistanceEndPt(tracksVec[i_tc][i_t], tracksVec[i_tr.first][i_tr.second]) <
247  fIsolation) {
248  if (!track_linked) track_linked = true;
249  fTrackContainmentLevel[i_tc][i_t] = containment_level;
250  fTrackContainmentIndices.back().emplace_back(i_tc, i_t);
251 
252  if (fDebug) {
253  std::cout << "\tTrackPair (" << i_tc << "," << i_t << ") and (" << i_tr.first << ","
254  << i_tr.second << ")"
255  << " " << containment_level << std::endl;
256  }
257 
258  } //end if track not isolated
259 
260  } //end loop over existing uncontained/linked tracks
261 
262  } //end if track not already linked
263 
264  } //end loops over tracks
265  } //end loop over track collections
266  } //end while linking tracks
267 
268  if (fDebug) std::cout << "All done! Now let's make the tree and tags!" << std::endl;
269 
270  //now we're going to will the tree and create tags if requested
271  for (size_t i_tc = 0; i_tc < tracksVec.size(); ++i_tc) {
272  for (size_t i_t = 0; i_t < tracksVec[i_tc].size(); ++i_t) {
273 
274  //fill ROOT Tree
275  if (fFillOutputTree) {
276  fTrackTreeObj = TrackTree_t(tracksVec[i_tc][i_t]);
277  fDistance = fMinDistances[i_tc][i_t];
278  fCollection = i_tc;
279  fTrkID = i_t;
281  fTrackTree->Fill();
282  }
283 
284  if (fMakeCosmicTags) {
285 
286  //default (if track looks contained and isolated)
287  float score = 0;
289 
290  //overwrite if track is not contained or not isolated
291  if (fTrackContainmentLevel[i_tc][i_t] >= 0) {
292  score = 1. / (1. + (float)fTrackContainmentLevel[i_tc][i_t]);
293  if (fTrackContainmentLevel[i_tc][i_t] == 0)
294  id = GetCosmicTagID(tracksVec[i_tc][i_t], geo);
295  else
297  }
298 
299  fCosmicTags[i_tc][i_t] =
300  anab::CosmicTag(std::vector<float>{(float)tracksVec[i_tc][i_t].Vertex().X(),
301  (float)tracksVec[i_tc][i_t].Vertex().Y(),
302  (float)tracksVec[i_tc][i_t].Vertex().Z()},
303  std::vector<float>{(float)tracksVec[i_tc][i_t].End().X(),
304  (float)tracksVec[i_tc][i_t].End().Y(),
305  (float)tracksVec[i_tc][i_t].End().Z()},
306  score,
307  id);
308  } //end cosmic tag making
309 
310  //some debug work
311  if (fTrackContainmentLevel[i_tc][i_t] < 0 && fDebug) {
312  std::cout << "Track (" << i_tc << "," << i_t << ")"
313  << " " << fTrackContainmentLevel[i_tc][i_t] << " " << fMinDistances[i_tc][i_t]
314  << std::endl;
315  auto const& vertex = tracksVec[i_tc][i_t].Vertex();
316  std::cout << "\tS_(X,Y,Z) = (" << vertex.X() << "," << vertex.Y() << "," << vertex.Z()
317  << ")\n";
318  std::cout << "\tNearest wire ..." << std::endl;
319  for (unsigned int i_p = 0; i_p < geo.Nplanes(); ++i_p)
320  std::cout << "\t\tPlane " << i_p << " "
321  << geo.NearestWireID(vertex, geo::PlaneID{0, 0, i_p}).Wire << std::endl;
322 
323  auto const& end = tracksVec[i_tc][i_t].End();
324  std::cout << "\tE_(X,Y,Z) = (" << end.X() << "," << end.Y() << "," << end.Z() << ")\n";
325  std::cout << "\tNearest wire ..." << std::endl;
326  for (unsigned int i_p = 0; i_p < geo.Nplanes(); ++i_p)
327  std::cout << "\t\tPlane " << i_p << " "
328  << geo.NearestWireID(end, geo::PlaneID{0, 0, i_p}).Wire << std::endl;
329  std::cout << "\tLength=" << tracksVec[i_tc][i_t].Length() << std::endl;
330  std::cout << "\tSimple_length=" << (end - vertex).R() << std::endl;
331  } //end debug statements if track contained
332 
333  } //end loops over tracks
334  } //end loop over track collections
335 
336 } //end ProcessTracks
337 
338 std::vector<std::vector<anab::CosmicTag>> const& trk::TrackContainmentAlg::GetTrackCosmicTags()
339 {
340  if (fMakeCosmicTags)
341  return fCosmicTags;
342  else
343  throw cet::exception("TrackContainmentAlg::GetTrackCosmicTags")
344  << "Cosmic tags not created. Set MakeCosmicTags to true in fcl paramters.";
345 }
double MinDistanceStartPt(recob::Track const &, recob::Track const &)
void ProcessTracks(std::vector< std::vector< recob::Track >> const &, geo::GeometryCore const &)
std::vector< std::vector< std::pair< int, int > > > fTrackContainmentIndices
enum anab::cosmic_tag_id CosmicTagID_t
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::vector< std::vector< anab::CosmicTag > > fCosmicTags
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
anab::CosmicTagID_t GetCosmicTagID(recob::Track const &, geo::GeometryCore const &)
Float_t tmp
Definition: plot.C:35
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Access the description of detector geometry.
struct trk::TrackTree TrackTree_t
std::vector< std::vector< int > > fTrackContainmentLevel
std::string Leaflist()
T get(std::string const &key) const
Definition: ParameterSet.h:314
Point_t const & Vertex() const
Access to track position at different points.
Definition: Track.h:158
void Configure(fhicl::ParameterSet const &)
void SetRunEvent(unsigned int const &, unsigned int const &)
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
double MinDistanceEndPt(recob::Track const &, recob::Track const &)
std::vector< std::vector< double > > fMinDistances
WireID NearestWireID(Point_t const &point, PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
std::vector< std::vector< anab::CosmicTag > > const & GetTrackCosmicTags()
TrackContainmentAlg()
Default constructor.
Point_t const & End() const
Access to track position at different points.
Definition: Track.h:159
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
Namespace collecting geometry-related classes utilities.
Float_t track
Definition: plot.C:35
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
bool IsContained(recob::Track const &, geo::GeometryCore const &)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Event finding and building.
vertex reconstruction