LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PMAlgCosmicTagger.cxx
Go to the documentation of this file.
1 // Class: PMAlgCosmicTagger
3 // Author: L. Whitehead (leigh.howard.whitehead@cern.ch),
4 // R. Sulej (robert.sulej@cern.ch) March 2017
6 
15 
17 
19 
20 #include <numeric> // std::accumulate
21 
23  pma::TrkCandidateColl& tracks)
24 {
25  // Get the detector dimensions
26  GetDimensions();
27 
28  mf::LogInfo("pma::PMAlgCosmicTagger")
29  << "Passed " << tracks.size() << " tracks for tagging cosmics.";
30 
31  size_t n = 0;
32 
33  if (fTagOutOfDriftTracks) n += outOfDriftWindow(tracks);
34  if (fTagFullHeightTracks) n += fullHeightCrossing(tracks);
35  if (fTagFullWidthTracks) n += fullWidthCrossing(tracks);
36  if (fTagFullLengthTracks) n += fullLengthCrossing(tracks);
37  if (fTagNonBeamT0Tracks) n += nonBeamT0Tag(clockData, tracks);
38  if (fTagTopFrontBack) n += tagTopFrontBack(tracks);
39  if (fTagApparentStopper) n += tagApparentStopper(tracks);
40 
41  mf::LogInfo("pma::PMAlgCosmicTagger") << "...tagged " << n << " cosmic-like tracks.";
42 }
43 
45 {
46  mf::LogInfo("pma::PMAlgCosmicTagger") << " - tag tracks out of 1 drift window;";
47  size_t n = 0;
48 
49  auto const* geom = lar::providerFrom<geo::Geometry>();
50 
51  for (auto& t : tracks.tracks()) {
52 
53  pma::Track3D& trk = *(t.Track());
54 
55  double min, max, p;
56  bool node_out_of_drift_min = false;
57  bool node_out_of_drift_max = false;
58  for (size_t nidx = 0; nidx < trk.Nodes().size(); ++nidx) {
59  auto const& node = *(trk.Nodes()[nidx]);
60  auto const& tpcGeo = geom->TPC(geo::TPCID(node.Cryo(), node.TPC()));
61  // DetectDriftDirection returns a short int, but switch requires an int
62  int driftDir = abs(tpcGeo.DetectDriftDirection());
63  p = node.Point3D()[driftDir - 1];
64  switch (driftDir) {
65  case 1:
66  min = tpcGeo.MinX();
67  max = tpcGeo.MaxX();
68  break;
69  case 2:
70  min = tpcGeo.MinY();
71  max = tpcGeo.MaxY();
72  break;
73  case 3:
74  min = tpcGeo.MinZ();
75  max = tpcGeo.MaxZ();
76  break;
77  default:
78  throw cet::exception("PMAlgCosmicTagger")
79  << "Drift direction unknown: " << driftDir << std::endl;
80  }
81  if (p < min - fOutOfDriftMargin) { node_out_of_drift_min = true; }
82  if (p > max + fOutOfDriftMargin) { node_out_of_drift_max = true; }
83  }
84 
85  if (node_out_of_drift_min && node_out_of_drift_max) {
88  ++n;
89  }
90  else if (node_out_of_drift_min || node_out_of_drift_max) {
93  ++n;
94  }
95  }
96 
97  mf::LogInfo("pma::PMAlgCosmicTagger") << " - Tagged " << n << " tracks out of 1 drift window.";
98 
99  return n;
100 }
101 
102 // Leigh: Make use of the fact that our cathode and anode crossing tracks have a reconstructed T0.
103 // Check to see if this time is consistent with the beam
105  pma::TrkCandidateColl& tracks) const
106 {
107  size_t n = 0;
108 
109  // Search through all of the tracks
110  for (auto& t : tracks.tracks()) {
111 
112  // Non zero T0 means we reconstructed it
113  if (t.Track()->GetT0() != 0.0) {
114  mf::LogInfo("pma::PMAlgCosmicTagger") << " - track with T0 = " << t.Track()->GetT0();
115 
116  if (fabs(t.Track()->GetT0() - trigger_offset(clockData)) > fNonBeamT0Margin) {
117  ++n;
118  t.Track()->SetTagFlag(pma::Track3D::kCosmic);
119  t.Track()->SetTagFlag(pma::Track3D::kBeamIncompatible);
120  }
121  }
122  }
123 
124  mf::LogInfo("pma::PMAlgCosmicTagger") << " - Tagged " << n << " non-beam T0 tracks.";
125  return n;
126 }
127 
129 {
130 
131  size_t n = 0;
132 
133  auto const* geom = lar::providerFrom<geo::Geometry>();
134 
135  short int hIdx = ConvertDirToInt(geom->TPC().HeightDir());
136  short int lIdx = ConvertDirToInt(geom->TPC().LengthDir());
137 
138  // Loop over the tracks
139  for (auto& t : tracks.tracks()) {
140 
141  // Get the first and last positions from the track.
142  auto const& node0 = *(t.Track()->Nodes()[0]);
143  auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);
144 
145  // Check which end is the vertex (assume the largest height)
146  TVector3 vtx =
147  (node0.Point3D()[hIdx] > node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
148  TVector3 end =
149  (node0.Point3D()[hIdx] <= node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
150 
151  // Check we have a track starting at the top of the detector
152  bool top = isTopVertex(vtx, fTopFrontBackMargin, hIdx);
153 
154  // Check the track ends at the front or back of the detector
155  bool frontBack = isFrontBackVertex(end, fTopFrontBackMargin, lIdx);
156 
157  // Check we path both criteria but without letting either the start or end of the track fulfill both
158  if (top && frontBack) {
159  ++n;
160  t.Track()->SetTagFlag(pma::Track3D::kCosmic);
161  t.Track()->SetTagFlag(pma::Track3D::kGeometry_YZ);
162  }
163  }
164 
165  mf::LogInfo("pma::PMAlgCosmicTagger")
166  << " - Tagged " << n << " tracks crossing from top to front/back." << std::endl;
167 
168  return n;
169 }
170 
172 {
173 
174  size_t n = 0;
175 
176  // Tracks entering from the top of the detector that stop in the fiducial volume
177  // are very likely to be cosmics that leave through the APA, but have their
178  // drift coordinate incorrectly set due to lack of T0
179  auto const* geom = lar::providerFrom<geo::Geometry>();
180 
181  short int hIdx = ConvertDirToInt(geom->TPC().HeightDir());
182 
183  // Loop over the tracks
184  for (auto& t : tracks.tracks()) {
185 
186  // Get the first and last positions from the track.
187  auto const& node0 = *(t.Track()->Nodes()[0]);
188  auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);
189 
190  // Check which end is the vertex (assume the largest height)
191  TVector3 vtx =
192  (node0.Point3D()[hIdx] > node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
193  TVector3 end =
194  (node0.Point3D()[hIdx] <= node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
195 
196  if (fabs(vtx[hIdx] - fDimensionsMax[hIdx]) < fApparentStopperMargin) {
197  // Check the other element to see if it ends away from the bottom of the detector
198  if (fabs(end[hIdx] - fDimensionsMin[hIdx]) > 5. * fApparentStopperMargin) {
199 
200  // We now need to loop over all of the tracks to see if any start within fStopperBuffer of our end point.
201  bool foundTrack = false;
202  for (auto const& tt : tracks.tracks()) {
203  // Don't match with itself!
204  if ((&tt) == (&t)) continue;
205 
206  // Compare this track with our main track
207  TVector3 trkVtx = (tt.Track()->Nodes()[0])->Point3D();
208  TVector3 trkEnd = (tt.Track()->Nodes()[tt.Track()->Nodes().size() - 1])->Point3D();
209 
210  if ((end - trkVtx).Mag() < fStopperBuffer || (end - trkEnd).Mag() < fStopperBuffer) {
211  foundTrack = true;
212  break;
213  }
214  }
215  if (foundTrack) {
216  // This isn't really a stopping particle, so move on
217  continue;
218  }
219 
220  // If we don't mind about tagging all stopping particles then this satisfies our requirements
221  if (!fVetoActualStopper) {
222  ++n;
223  t.Track()->SetTagFlag(pma::Track3D::kCosmic);
224  t.Track()->SetTagFlag(pma::Track3D::kGeometry_Y);
225  continue;
226  }
227 
228  // If we want to actually ignore the stopping particles, use de/dx...
229  // Store the number of sigma from the mean for the final dedx point in each view
230  std::vector<float> nSigmaPerView;
231 
232  // Loop over the views
233  for (auto const view : geom->Views()) {
234 
235  // Get the dedx for this track and view
236  std::map<size_t, std::vector<double>> dedx;
237  t.Track()->GetRawdEdxSequence(dedx, view);
238 
239  std::vector<double> trk_dedx;
240 
241  for (int h = t.Track()->NextHit(-1, view); h != -1; h = t.Track()->NextHit(h, view)) {
242  // If this is the last hit then this method won't work
243  if (h > t.Track()->PrevHit(t.Track()->size(), view)) break;
244  // Make sure we have a reasonable value
245  if (dedx[h][5] / dedx[h][6] <= 0 || dedx[h][5] / dedx[h][6] > 1e6) continue;
246  trk_dedx.push_back(dedx[h][5] / dedx[h][6]);
247  }
248 
249  if (trk_dedx.size() == 0) {
250  mf::LogInfo("pma::PMAlgCosmicTagger")
251  << "View " << view << " has no hits." << std::endl;
252  continue;
253  }
254 
255  double sum = std::accumulate(std::begin(trk_dedx), std::end(trk_dedx), 0.0);
256  double mean = sum / static_cast<double>(trk_dedx.size());
257  double accum = 0.0;
258  std::for_each(std::begin(trk_dedx), std::end(trk_dedx), [&](const double d) {
259  accum += (d - mean) * (d - mean);
260  });
261  double stdev = sqrt(accum / static_cast<double>(trk_dedx.size() - 1));
262 
263  mf::LogInfo("pma::PMAlgCosmicTagger")
264  << " View " << view << " has average dedx " << mean << " +/- " << stdev
265  << " and final dedx " << trk_dedx[trk_dedx.size() - 1] << std::endl;
266 
267  nSigmaPerView.push_back(fabs((trk_dedx[trk_dedx.size() - 1] - mean) / stdev));
268  }
269 
270  bool notStopper = true;
271  short unsigned int n2Sigma = 0;
272  short unsigned int n3Sigma = 0;
273  for (auto const nSigma : nSigmaPerView) {
274  if (nSigma >= 2.0) ++n2Sigma;
275  if (nSigma >= 3.0) ++n3Sigma;
276  }
277 
278  if (n3Sigma > 0) notStopper = false;
279  if (n2Sigma == nSigmaPerView.size()) notStopper = false;
280 
281  if (notStopper) {
282  ++n;
283  t.Track()->SetTagFlag(pma::Track3D::kCosmic);
284  t.Track()->SetTagFlag(pma::Track3D::kGeometry_Y);
285  }
286  } // Check on bottom position
287  } // Check on top position
288  } // End loop over tracks
289 
290  mf::LogInfo("pma::PMAlgCosmicTagger")
291  << " - Tagged " << n << " tracks stopping in the detector after starting at the top."
292  << std::endl;
293 
294  return n;
295 }
296 
298 {
299 
300  // Just use the first tpc to get the height dimension (instead of assuming y).
301  auto const* geom = lar::providerFrom<geo::Geometry>();
302  auto const dir = geom->TPC().HeightDir();
303 
304  size_t n = fullCrossingTagger(tracks, ConvertDirToInt(dir));
305 
306  mf::LogInfo("pma::PMAlgCosmicTagger")
307  << " - Tagged " << n << " tracks crossing the full detector height";
308  return n;
309 }
310 
312 {
313 
314  // Just use the first tpc to get the width dimension (instead of assuming x).
315  auto const* geom = lar::providerFrom<geo::Geometry>();
316  auto const dir = geom->TPC().WidthDir();
317 
318  size_t n = fullCrossingTagger(tracks, ConvertDirToInt(dir));
319 
320  mf::LogInfo("pma::PMAlgCosmicTagger")
321  << " - Tagged " << n << " tracks crossing the full detector width";
322  return n;
323 }
324 
326 {
327 
328  // Just use the first tpc to get the length dimension (instead of assuming z).
329  auto const* geom = lar::providerFrom<geo::Geometry>();
330  auto const dir = geom->TPC().LengthDir();
331 
332  size_t n = fullCrossingTagger(tracks, ConvertDirToInt(dir));
333 
334  mf::LogInfo("pma::PMAlgCosmicTagger")
335  << " - Tagged " << n << " tracks crossing the full detector length";
336  return n;
337 }
338 
340  int direction) const
341 {
342 
343  if (direction == -1) {
344  mf::LogWarning("pma::PMAlgCosmicTagger")
345  << " - Could not recognise direction, not attempting to perform fullCrossingTagger.";
346  return 0;
347  }
348 
349  size_t n = 0;
350 
351  double detDim = fDimensionsMax[direction] - fDimensionsMin[direction];
352 
354  switch (direction) {
355  case 0: dirTag = pma::Track3D::kGeometry_XX; break;
356  case 1: dirTag = pma::Track3D::kGeometry_YY; break;
357  case 2: dirTag = pma::Track3D::kGeometry_ZZ; break;
358  default: dirTag = pma::Track3D::kNotTagged; break;
359  }
360 
361  // Loop over the tracks
362  for (auto& t : tracks.tracks()) {
363 
364  // Get the first and last positions from the track.
365  auto const& node0 = *(t.Track()->Nodes()[0]);
366  auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);
367 
368  // Get the length of the track in the requested direction
369  double trkDim = fabs(node0.Point3D()[direction] - node1.Point3D()[direction]);
370 
371  if ((detDim - trkDim) < fFullCrossingMargin) {
372  ++n;
373  t.Track()->SetTagFlag(pma::Track3D::kCosmic);
374  t.Track()->SetTagFlag(dirTag);
375  mf::LogInfo("pma::PMAlgCosmicTagger") << " -- track tagged in direction " << direction
376  << " with " << trkDim << " (c.f. " << detDim << ")";
377  }
378  }
379 
380  return n;
381 }
382 
383 bool pma::PMAlgCosmicTagger::isTopVertex(const TVector3& pos,
384  double tolerance,
385  short int dirIndx) const
386 {
387 
388  return (fabs(pos[dirIndx] - fDimensionsMax[dirIndx]) < tolerance);
389 }
390 
392  double tolerance,
393  short int dirIndx) const
394 {
395 
396  bool front = (fabs(pos[dirIndx] - fDimensionsMin[dirIndx]) < tolerance);
397  bool back = (fabs(pos[dirIndx] - fDimensionsMax[dirIndx]) < tolerance);
398 
399  return front || back;
400 }
401 
403 {
404 
405  // Need to find the minimum and maximum height values from the geometry.
406  double minX = 1.e6;
407  double maxX = -1.e6;
408  double minY = 1.e6;
409  double maxY = -1.e6;
410  double minZ = 1.e6;
411  double maxZ = -1.e6;
412 
413  auto const* geom = lar::providerFrom<geo::Geometry>();
414 
415  // Since we can stack TPCs, we can't just use geom::TPCGeom::Height()
416  for (geo::TPCGeo const& TPC : geom->Iterate<geo::TPCGeo>()) {
417  // We need to make sure we only include the real TPCs
418  // We have dummy TPCs in the protoDUNE and DUNE geometries
419  // The dummy ones have a drift distance of only ~13 cm.
420  if (TPC.DriftDistance() < 25.0) { continue; }
421 
422  // get center in world coordinates
423  auto const center = TPC.GetCenter();
424  double tpcDim[3] = {TPC.HalfWidth(), TPC.HalfHeight(), 0.5 * TPC.Length()};
425 
426  if (center.X() - tpcDim[0] < minX) minX = center.X() - tpcDim[0];
427  if (center.X() + tpcDim[0] > maxX) maxX = center.X() + tpcDim[0];
428  if (center.Y() - tpcDim[1] < minY) minY = center.Y() - tpcDim[1];
429  if (center.Y() + tpcDim[1] > maxY) maxY = center.Y() + tpcDim[1];
430  if (center.Z() - tpcDim[2] < minZ) minZ = center.Z() - tpcDim[2];
431  if (center.Z() + tpcDim[2] > maxZ) maxZ = center.Z() + tpcDim[2];
432  } // for all TPC
433 
434  fDimensionsMin.clear();
435  fDimensionsMax.clear();
436  fDimensionsMin.push_back(minX);
437  fDimensionsMin.push_back(minY);
438  fDimensionsMin.push_back(minZ);
439  fDimensionsMax.push_back(maxX);
440  fDimensionsMax.push_back(maxY);
441  fDimensionsMax.push_back(maxZ);
442 }
443 
445 {
446 
447  if (dir.X() > 0.99) return 0;
448  if (dir.Y() > 0.99) return 1;
449  if (dir.Z() > 0.99)
450  return 2;
451 
452  else
453  return -1;
454 }
size_t fullWidthCrossing(pma::TrkCandidateColl &tracks) const
size_t size() const
double minY
Definition: plot_hist.C:9
Utilities related to art service access.
void SetTagFlag(ETag value)
Definition: PmaTrack3D.h:68
Implementation of the Projection Matching Algorithm.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Geometry information for a single TPC.
Definition: TPCGeo.h:36
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
constexpr auto abs(T v)
Returns the absolute value of the argument.
Implementation of the Projection Matching Algorithm.
double maxY
Definition: plot_hist.C:10
short int ConvertDirToInt(const geo::Vector_t &dir) const
std::vector< double > fDimensionsMin
size_t tagApparentStopper(pma::TrkCandidateColl &tracks) const
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
Access the description of detector geometry.
Definition: type_traits.h:61
size_t outOfDriftWindow(pma::TrkCandidateColl &tracks) const
bool isFrontBackVertex(const TVector3 &pos, double tolerance, short int dirIndx) const
size_t nonBeamT0Tag(detinfo::DetectorClocksData const &clockData, pma::TrkCandidateColl &tracks) const
Float_t d
Definition: plot.C:235
std::vector< double > fDimensionsMax
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
void tag(detinfo::DetectorClocksData const &clockData, pma::TrkCandidateColl &tracks)
Track finding helper for the Projection Matching Algorithm.
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
Contains all timing reference information for the detector.
TDirectory * dir
Definition: macro.C:5
size_t tagTopFrontBack(pma::TrkCandidateColl &tracks) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
int trigger_offset(DetectorClocksData const &data)
size_t fullHeightCrossing(pma::TrkCandidateColl &tracks) const
Char_t n[5]
size_t fullCrossingTagger(pma::TrkCandidateColl &tracks, int direction) const
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
std::vector< pma::Node3D * > const & Nodes() const noexcept
Definition: PmaTrack3D.h:274
bool isTopVertex(const TVector3 &pos, double tolerance, short int dirIndx) const
Double_t sum
Definition: plot.C:31
size_t fullLengthCrossing(pma::TrkCandidateColl &tracks) const
art framework interface to geometry description
std::vector< TrkCandidate > const & tracks() const
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane.