LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
EMShowerAlg.cxx
Go to the documentation of this file.
1 // Implementation of the EMShower algorithm
3 //
4 // Forms EM showers from clusters and associated tracks.
5 // Also provides methods for finding the vertex and further
6 // properties of the shower.
7 //
8 // Mike Wallbank (m.wallbank@sheffield.ac.uk), September 2015
10 
11 #include "cetlib/container_algorithms.h"
12 #include "cetlib/pow.h"
14 
19 
20 #include "TAxis.h"
21 #include "TCanvas.h"
22 #include "TFile.h"
23 #include "TGraph.h"
24 #include "TH1.h"
25 #include "TLine.h"
26 #include "TMath.h"
27 #include "TMathBase.h"
28 #include "TMultiGraph.h"
29 #include "TProfile.h"
30 
31 #include "range/v3/numeric.hpp"
32 #include "range/v3/view.hpp"
33 
34 using lar::to_element;
35 using namespace ranges;
36 
38  : fDebug{debug}
39  , fMinTrackLength{pset.get<double>("MinTrackLength")}
40  , fdEdxTrackLength{pset.get<double>("dEdxTrackLength")}
41  , fSpacePointSize{pset.get<double>("SpacePointSize")}
42  , fNfitpass{pset.get<unsigned int>("Nfitpass")}
43  , fNfithits{pset.get<std::vector<unsigned int>>("Nfithits")}
44  , fToler{pset.get<std::vector<double>>("Toler")}
45  , fShowerEnergyAlg(pset.get<fhicl::ParameterSet>("ShowerEnergyAlg"))
46  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
47  , fProjectionMatchingAlg(pset.get<fhicl::ParameterSet>("ProjectionMatchingAlg"))
48  , fDetector{pset.get<std::string>("Detector", "dune35t")} // tmp
49  , fMakeGradientPlot{pset.get<bool>("MakeGradientPlot", false)}
50  , fMakeRMSGradientPlot{pset.get<bool>("MakeRMSGradientPlot", false)}
51  , fNumShowerSegments{pset.get<int>("NumShowerSegments", 4)}
52 {
53  if (fNfitpass != fNfithits.size() || fNfitpass != fToler.size()) {
55  << "EMShowerAlg: fNfithits and fToler need to have size fNfitpass";
56  }
57 }
58 
60  std::vector<art::Ptr<recob::Cluster>> const& clusters,
61  art::FindManyP<recob::Hit> const& fmh,
63  std::map<int, std::vector<int>>& clusterToTracks,
64  std::map<int, std::vector<int>>& trackToClusters) const
65 {
66  std::vector<int> clustersToIgnore = {-999};
68  clusters, fmh, fmt, clustersToIgnore, clusterToTracks, trackToClusters);
69 }
70 
72  std::vector<art::Ptr<recob::Cluster>> const& clusters,
73  art::FindManyP<recob::Hit> const& fmh,
75  std::vector<int> const& clustersToIgnore,
76  std::map<int, std::vector<int>>& clusterToTracks,
77  std::map<int, std::vector<int>>& trackToClusters) const
78 {
79  // Look through all the clusters
80  for (auto const& clusterPtr : clusters) {
81 
82  // Get the hits in this cluster
83  auto const& clusterHits = fmh.at(clusterPtr.key());
84 
85  // Look at all these hits and find the associated tracks
86  for (auto const& clusterHitPtr : clusterHits) {
87  // Get the tracks associated with this hit
88  auto const& clusterHitTracks = fmt.at(clusterHitPtr.key());
89  if (clusterHitTracks.size() > 1) {
90  std::cout << "More than one track associated with this hit!\n";
91  continue;
92  }
93 
94  if (clusterHitTracks.size() < 1) continue;
95 
96  auto const& clusterHitTrackPtr = clusterHitTracks[0];
97  if (clusterHitTrackPtr->Length() < fMinTrackLength) {
98  if (fDebug > 2)
99  std::cout << "Track " << clusterHitTrackPtr->ID() << " is too short! ("
100  << clusterHitTrackPtr->Length() << ")\n";
101  continue;
102  }
103 
104  // Add this cluster to the track map
105  int track = clusterHitTrackPtr.key();
106  int cluster = clusterPtr.key();
107  if (cet::search_all(clustersToIgnore, cluster)) continue;
108  if (not cet::search_all(trackToClusters[track], cluster))
109  trackToClusters[track].push_back(cluster);
110  if (not cet::search_all(clusterToTracks[cluster], track))
111  clusterToTracks[cluster].push_back(track);
112  }
113  }
114 }
115 
117  std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
118 {
119  std::map<int, std::vector<int>> firstTPC;
120  for (auto const& [plane, hits] : showerHitsMap)
121  firstTPC[hits.at(0)->WireID().TPC].push_back(plane);
122 
123  // If all in the same TPC then that's great!
124  if (firstTPC.size() != 2) return;
125 
126  // If they are in more than two TPCs, not much we can do
127  else if (firstTPC.size() > 2)
128  return;
129 
130  // If we get to this point, there should be something we can do!
131 
132  // Find the problem plane
133  int problemPlane = -1;
134  for (auto const& planes : firstTPC | views::values)
135  if (planes.size() == 1) problemPlane = planes[0];
136 
137  // Require three hits
138  if (showerHitsMap.at(problemPlane).size() < 3) return;
139 
140  // and get the other planes with at least three hits
141  std::vector<int> otherPlanes;
142  for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane)
143  if (plane != problemPlane and showerHitsMap.count(plane) and
144  showerHitsMap.at(plane).size() >= 3)
145  otherPlanes.push_back(plane);
146 
147  if (otherPlanes.size() == 0) return;
148 
149  // Look at the hits after the first one
150  unsigned int wrongTPC = showerHitsMap.at(problemPlane).at(0)->WireID().TPC;
151  unsigned int nHits = 0;
152  for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsMap.at(problemPlane).begin();
153  hitIt != showerHitsMap.at(problemPlane).end() and (*hitIt)->WireID().TPC == wrongTPC;
154  ++hitIt)
155  ++nHits;
156 
157  // If there are more than two hits in the 'wrong TPC', we can't be sure it is indeed wrong
158  if (nHits > 2) return;
159 
160  // See if at least the next four times as many hits are in a different TPC
161  std::map<int, int> otherTPCs;
163  std::next(showerHitsMap.at(problemPlane).begin(), nHits);
164  hitIt != showerHitsMap.at(problemPlane).end() and
165  std::distance(std::next(showerHitsMap.at(problemPlane).begin(), nHits), hitIt) < 4 * nHits;
166  ++hitIt)
167  ++otherTPCs[(*hitIt)->WireID().TPC];
168 
169  if (otherTPCs.size() > 1) return;
170 
171  // If we get this far, we can move the problem hits from the front of the shower to the back
172  std::map<int, int> tpcCount;
173  for (int const otherPlane : otherPlanes)
175  std::next(showerHitsMap.at(otherPlane).begin());
176  hitIt != showerHitsMap.at(otherPlane).end() and
177  hitIt != std::next(showerHitsMap.at(otherPlane).begin(), 2);
178  ++hitIt)
179  ++tpcCount[(*hitIt)->WireID().TPC];
180 
181  // Remove the first hit if it is in the wrong TPC
182  if (tpcCount.size() == 1 and
183  tpcCount.begin()->first ==
184  (int)(*std::next(showerHitsMap.at(problemPlane).begin(), nHits))->WireID().TPC) {
185  std::vector<art::Ptr<recob::Hit>> naughty_hits;
186  for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsMap.at(problemPlane).begin();
187  hitIt != std::next(showerHitsMap.at(problemPlane).begin(), nHits);
188  ++hitIt) {
189  naughty_hits.push_back(*hitIt);
190  showerHitsMap.at(problemPlane).erase(hitIt);
191  }
192  for (auto const& naughty_hit : naughty_hits)
193  showerHitsMap.at(problemPlane).push_back(naughty_hit);
194  }
195 }
196 
198  detinfo::DetectorPropertiesData const& detProp,
199  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& showerHitsMap) const
200 {
201  bool consistencyCheck = true;
202 
203  if (showerHitsMap.size() < 2) { consistencyCheck = true; }
204  else if (showerHitsMap.size() == 2) {
205 
206  // With two views, we can check:
207  // -- timing between views is consistent
208  // -- the 3D start point makes sense when projected back onto the individual planes
209 
210  std::vector<art::Ptr<recob::Hit>> startHits;
211  std::vector<int> planes;
212  for (auto const& [plane, hits] : showerHitsMap) {
213  startHits.push_back(hits.front());
214  planes.push_back(plane);
215  }
216 
217  auto const showerStartPos = Construct3DPoint_(detProp, startHits.at(0), startHits.at(1));
218  TVector2 proj1 = Project3DPointOntoPlane_(detProp, showerStartPos, planes.at(0));
219  TVector2 proj2 = Project3DPointOntoPlane_(detProp, showerStartPos, planes.at(1));
220 
221  double timingDifference = std::abs(startHits.at(0)->PeakTime() - startHits.at(1)->PeakTime());
222  double projectionDifference = ((HitPosition_(detProp, *startHits.at(0)) - proj1).Mod() +
223  (HitPosition_(detProp, *startHits.at(1)) - proj2).Mod()) /
224  (double)2;
225 
226  if (timingDifference > 40 or projectionDifference > 5 or showerStartPos.X() == -9999 or
227  showerStartPos.Y() == -9999 or showerStartPos.Z() == -9999)
228  consistencyCheck = false;
229 
230  if (fDebug > 1)
231  std::cout << "Timing difference is " << timingDifference << " and projection distance is "
232  << projectionDifference << " (start is (" << showerStartPos.X() << ", "
233  << showerStartPos.Y() << ", " << showerStartPos.Z() << ")\n";
234  }
235  else if (showerHitsMap.size() == 3) {
236 
237  // With three views, we can check:
238  // -- the timing between views is consistent
239  // -- the 3D start point formed by two views and projected back into the third is close to the start point in that view
240 
241  std::map<int, art::Ptr<recob::Hit>> start2DMap;
242  for (auto const& [plane, hits] : showerHitsMap) {
243  start2DMap[plane] = hits.front();
244  }
245 
246  std::map<int, double> projDiff;
247  std::map<int, double> timingDiff;
248 
249  for (int plane = 0; plane < 3; ++plane) {
250 
251  std::vector<int> otherPlanes;
252  for (int otherPlane = 0; otherPlane < 3; ++otherPlane)
253  if (otherPlane != plane) otherPlanes.push_back(otherPlane);
254 
255  auto const showerStartPos = Construct3DPoint_(
256  detProp, start2DMap.at(otherPlanes.at(0)), start2DMap.at(otherPlanes.at(1)));
257  TVector2 showerStartProj = Project3DPointOntoPlane_(detProp, showerStartPos, plane);
258 
259  if (fDebug > 2) {
260  std::cout << "Plane... " << plane;
261  std::cout << "\nStart position in this plane is "
262  << HitPosition_(detProp, *start2DMap.at(plane)).X() << ", "
263  << HitPosition_(detProp, *start2DMap.at(plane)).Y() << ")\n";
264  std::cout << "Shower start from other two planes is (" << showerStartPos.X() << ", "
265  << showerStartPos.Y() << ", " << showerStartPos.Z() << ")\n";
266  std::cout << "Projecting the other two planes gives position (" << showerStartProj.X()
267  << ", " << showerStartProj.Y() << ")\n";
268  }
269 
270  double projDiff =
271  std::abs((showerStartProj - HitPosition_(detProp, *start2DMap.at(plane))).Mod());
272  double timeDiff = TMath::Max(
273  std::abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(0))->PeakTime()),
274  std::abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(1))->PeakTime()));
275 
276  if (fDebug > 1)
277  std::cout << "Plane " << plane << " has projDiff " << projDiff << " and timeDiff "
278  << timeDiff << '\n';
279 
280  if (projDiff > 5 or timeDiff > 40) consistencyCheck = false;
281  }
282  }
283 
284  if (fDebug > 1) std::cout << "Consistency check is " << consistencyCheck << '\n';
285 
286  return consistencyCheck;
287 }
288 
290  std::vector<std::vector<int>> const& initialShowers,
291  std::vector<art::Ptr<recob::Cluster>> const& clusters,
292  art::FindManyP<recob::Hit> const& fmh) const
293 {
294  std::vector<int> clustersToIgnore;
295 
296  // Look at each shower
297  for (auto initialShowerIt = initialShowers.cbegin(); initialShowerIt != initialShowers.cend();
298  ++initialShowerIt) {
299 
300  if (std::distance(initialShowers.cbegin(), initialShowerIt) > 0) continue;
301 
302  // Map the clusters and cluster hits to each view
303  std::map<int, std::vector<art::Ptr<recob::Cluster>>> planeClusters;
304  std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHits;
305  for (int const clusterIndex : *initialShowerIt) {
306  art::Ptr<recob::Cluster> const& cluster = clusters.at(clusterIndex);
307  planeClusters[cluster->Plane().Plane].push_back(cluster);
308  for (auto const& hit : fmh.at(cluster.key()))
309  planeHits[hit->WireID().Plane].push_back(hit);
310  }
311 
312  TFile* outFile = new TFile("chargeDistributions.root", "RECREATE");
313  std::map<int, TH1D*> chargeDist;
314  for (auto const& [plane, clusterPtrs] : planeClusters) {
315  for (auto const& clusterPtr : clusterPtrs) {
316  chargeDist[plane] = new TH1D(std::string("chargeDist_Plane" + std::to_string(plane) +
317  "_Cluster" + std::to_string(clusterPtr.key()))
318  .c_str(),
319  "",
320  150,
321  0,
322  1000);
323  auto const& hits = fmh.at(clusterPtr.key());
324  for (auto const& hit : hits | views::transform(to_element)) {
325  chargeDist[plane]->Fill(hit.Integral());
326  }
327  outFile->cd();
328  chargeDist[plane]->Write();
329  }
330  }
331  outFile->Close();
332  delete outFile;
333 
334  // Can't do much with fewer than three views
335  if (planeClusters.size() < 3) continue;
336 
337  // Look at how many clusters each plane has, and the proportion of hits each one uses
338  std::map<int, std::vector<double>> planeClusterSizes;
339  for (std::map<int, std::vector<art::Ptr<recob::Cluster>>>::iterator planeClustersIt =
340  planeClusters.begin();
341  planeClustersIt != planeClusters.end();
342  ++planeClustersIt) {
343  for (std::vector<art::Ptr<recob::Cluster>>::iterator planeClusterIt =
344  planeClustersIt->second.begin();
345  planeClusterIt != planeClustersIt->second.end();
346  ++planeClusterIt) {
347  std::vector<art::Ptr<recob::Hit>> hits = fmh.at(planeClusterIt->key());
348  planeClusterSizes[planeClustersIt->first].push_back(
349  (double)hits.size() / (double)planeHits.at(planeClustersIt->first).size());
350  }
351  }
352 
353  // Find the average hit fraction across all clusters in the plane
354  std::map<int, double> planeClustersAvSizes;
355  for (auto const& [plane, cluster_sizes] : planeClusterSizes) {
356  double const average = accumulate(cluster_sizes, 0.) / cluster_sizes.size();
357  planeClustersAvSizes[plane] = average;
358  }
359 
360  // Now decide if there is one plane which is ruining the reconstruction
361  // If two planes have a low average cluster fraction and one high, this plane likely merges two particle deposits together
362  int badPlane = -1;
363  for (auto const [plane, avg] : planeClustersAvSizes) {
364 
365  // Get averages from other planes and add in quadrature
366  std::vector<double> otherAverages;
367  for (auto const [other_plane, other_avg] : planeClustersAvSizes)
368  if (plane != other_plane) otherAverages.push_back(other_avg);
369 
370  double const sumSquareAvgsOtherPlanes = accumulate(
371  otherAverages, 0., [](double sum, double avg) { return sum + cet::square(avg); });
372  double const quadOtherPlanes = std::sqrt(sumSquareAvgsOtherPlanes);
373 
374  // If this plane has an average higher than the quadratic sum of the
375  // others, it may be bad
376  if (avg >= quadOtherPlanes) badPlane = plane;
377  }
378 
379  if (badPlane != -1) {
380  if (fDebug > 1) std::cout << "Bad plane is " << badPlane << '\n';
381  for (auto const& cluster : planeClusters.at(badPlane))
382  clustersToIgnore.push_back(cluster.key());
383  }
384  }
385 
386  return clustersToIgnore;
387 }
388 
390  art::Ptr<recob::Hit> const& hit1,
391  art::Ptr<recob::Hit> const& hit2) const
392 {
393 
394  // x is average of the two x's
395  double x = (detProp.ConvertTicksToX(hit1->PeakTime(), hit1->WireID().planeID()) +
396  detProp.ConvertTicksToX(hit2->PeakTime(), hit2->WireID().planeID())) /
397  (double)2;
398 
399  // y and z got from the wire interections
400  geo::WireIDIntersection intersection;
401  fGeom->WireIDsIntersect(hit1->WireID(), hit2->WireID(), intersection);
402 
403  return {x, intersection.y, intersection.z};
404 }
405 
406 std::unique_ptr<recob::Track> shower::EMShowerAlg::ConstructTrack(
407  detinfo::DetectorPropertiesData const& detProp,
408  std::vector<art::Ptr<recob::Hit>> const& hits1,
409  std::vector<art::Ptr<recob::Hit>> const& hits2,
410  std::map<int, TVector2> const& showerCentreMap) const
411 {
412 
413  std::unique_ptr<recob::Track> track;
414 
415  std::vector<art::Ptr<recob::Hit>> track1, track2;
416 
417  // Check the TPCs
418  if ((*hits1.begin())->WireID().TPC != (*hits2.begin())->WireID().TPC) {
419  mf::LogWarning("EMShowerAlg") << "Warning: attempting to construct a track from two different "
420  "TPCs. Returning a null track.";
421  return track;
422  }
423 
424  // Check for tracks crossing TPC boundaries
425  std::map<int, int> tpcMap;
426  for (auto const& hit : hits1)
427  ++tpcMap[hit->WireID().TPC];
428  for (auto const& hit : hits2)
429  ++tpcMap[hit->WireID().TPC];
430  if (tpcMap.size() > 1) {
431  mf::LogWarning("EMShowerAlg")
432  << "Warning: attempting to construct a track which crosses more than one TPC -- PMTrack "
433  "can't handle this right now. Returning a track made just from hits in the first TPC it "
434  "traverses.";
435  unsigned int firstTPC1 = hits1.at(0)->WireID().TPC, firstTPC2 = hits2.at(0)->WireID().TPC;
436  for (auto const& hit : hits1)
437  if (hit->WireID().TPC == firstTPC1) track1.push_back(hit);
438  for (auto const& hit : hits2)
439  if (hit->WireID().TPC == firstTPC2) track2.push_back(hit);
440  }
441  else {
442  track1 = hits1;
443  track2 = hits2;
444  }
445 
446  if (fDebug > 1) {
447  std::cout << "About to make a track from these hits:\n";
448  auto print_hits = [this](auto const& track) {
449  for (auto const& hit : track | views::transform(to_element)) {
450  std::cout << "Hit (" << HitCoordinates_(hit).X() << ", " << HitCoordinates_(hit).Y()
451  << ") (real wire " << hit.WireID().Wire << ") in TPC " << hit.WireID().TPC
452  << '\n';
453  }
454  };
455  print_hits(track1);
456  print_hits(track2);
457  }
458 
459  auto const trackStart = Construct3DPoint_(detProp, track1.at(0), track2.at(0));
460  pma::Track3D* pmatrack = fProjectionMatchingAlg.buildSegment(detProp, track1, track2, trackStart);
461 
462  if (!pmatrack) {
463  mf::LogInfo("EMShowerAlg") << "Skipping this event because not enough hits in two views";
464  return track;
465  }
466 
467  std::vector<TVector3> xyz, dircos;
468 
469  for (unsigned int i = 0; i < pmatrack->size(); i++) {
470 
471  xyz.push_back((*pmatrack)[i]->Point3D());
472 
473  if (i < pmatrack->size() - 1) {
474  size_t j = i + 1;
475  double mag = 0.0;
476  TVector3 dc(0., 0., 0.);
477  while ((mag == 0.0) and (j < pmatrack->size())) {
478  dc = (*pmatrack)[j]->Point3D();
479  dc -= (*pmatrack)[i]->Point3D();
480  mag = dc.Mag();
481  ++j;
482  }
483  if (mag > 0.0)
484  dc *= 1.0 / mag;
485  else if (!dircos.empty())
486  dc = dircos.back();
487  dircos.push_back(dc);
488  }
489  else
490  dircos.push_back(dircos.back());
491  }
492 
493  // Orient the track correctly
494  std::map<int, double> distanceToVertex, distanceToEnd;
495  using geo::vect::toPoint;
496  geo::Point_t const vertex = toPoint(*xyz.begin());
497  geo::Point_t const end = toPoint(*xyz.rbegin());
498 
499  // Loop over all the planes and find the distance from the vertex and end
500  // projections to the centre in each plane
501  for (std::map<int, TVector2>::const_iterator showerCentreIt = showerCentreMap.begin();
502  showerCentreIt != showerCentreMap.end();
503  ++showerCentreIt) {
504 
505  // Project the vertex and the end point onto this plane
506  TVector2 vertexProj = Project3DPointOntoPlane_(detProp, vertex, showerCentreIt->first);
507  TVector2 endProj = Project3DPointOntoPlane_(detProp, end, showerCentreIt->first);
508 
509  // Find the distance of each to the centre of the cluster
510  distanceToVertex[showerCentreIt->first] = (vertexProj - showerCentreIt->second).Mod();
511  distanceToEnd[showerCentreIt->first] = (endProj - showerCentreIt->second).Mod();
512  }
513 
514  // Find the average distance to the vertex and the end across the planes
515  double avDistanceToVertex = 0, avDistanceToEnd = 0;
516  for (std::map<int, double>::iterator distanceToVertexIt = distanceToVertex.begin();
517  distanceToVertexIt != distanceToVertex.end();
518  ++distanceToVertexIt)
519  avDistanceToVertex += distanceToVertexIt->second;
520  avDistanceToVertex /= distanceToVertex.size();
521 
522  for (std::map<int, double>::iterator distanceToEndIt = distanceToEnd.begin();
523  distanceToEndIt != distanceToEnd.end();
524  ++distanceToEndIt)
525  avDistanceToEnd += distanceToEndIt->second;
526  if (distanceToEnd.size() != 0) avDistanceToEnd /= distanceToEnd.size();
527 
528  if (fDebug > 2)
529  std::cout << "Distance to vertex is " << avDistanceToVertex << " and distance to end is "
530  << avDistanceToEnd << '\n';
531 
532  // Change order if necessary
533  if (avDistanceToEnd > avDistanceToVertex) {
534  std::reverse(xyz.begin(), xyz.end());
535  std::transform(
536  dircos.begin(), dircos.end(), dircos.begin(), [](TVector3 const& vec) { return -1 * vec; });
537  }
538 
539  if (xyz.size() != dircos.size())
540  mf::LogError("EMShowerAlg") << "Problem converting pma::Track3D to recob::Track";
541 
542  track = std::make_unique<recob::Track>(
545  recob::Track::Flags_t(xyz.size()),
546  false),
547  0,
548  -1.,
549  0,
552  -1);
553 
554  return track;
555 }
556 
557 std::unique_ptr<recob::Track> shower::EMShowerAlg::ConstructTrack(
558  detinfo::DetectorPropertiesData const& detProp,
559  std::vector<art::Ptr<recob::Hit>> const& track1,
560  std::vector<art::Ptr<recob::Hit>> const& track2) const
561 {
562  std::map<int, TVector2> showerCentreMap;
563  return ConstructTrack(detProp, track1, track2, showerCentreMap);
564 }
565 
567  detinfo::DetectorPropertiesData const& detProp,
568  std::vector<art::Ptr<recob::Hit>> const& trackHits,
569  std::unique_ptr<recob::Track> const& track) const
570 {
571  assert(not empty(trackHits));
572  if (!track) return -999;
573 
574  recob::Hit const& firstHit = *trackHits.front();
575 
576  // Get the pitch
577  double pitch = 0;
578  try {
579  pitch = lar::util::TrackPitchInView(*track, firstHit.View());
580  }
581  catch (...) {
582  }
583 
584  // Deal with large pitches
585  if (pitch > fdEdxTrackLength) {
586  return fCalorimetryAlg.dEdx_AREA(clockData, detProp, firstHit, pitch);
587  }
588 
589  double totalCharge = 0, totalDistance = 0, avHitTime = 0;
590  unsigned int nHits = 0;
591 
592  for (auto const& hit : trackHits | views::transform(to_element)) {
593  if (totalDistance + pitch < fdEdxTrackLength) {
594  totalDistance += pitch;
595  totalCharge += hit.Integral();
596  avHitTime += hit.PeakTime();
597  ++nHits;
598  }
599  }
600 
601  avHitTime /= (double)nHits;
602 
603  double const dQdx = totalCharge / totalDistance;
604  return fCalorimetryAlg.dEdx_AREA(clockData, detProp, dQdx, avHitTime, firstHit.WireID().Plane);
605 }
606 
608  detinfo::DetectorPropertiesData const& detProp,
609  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap,
610  std::unique_ptr<recob::Track>& initialTrack,
611  std::map<int, std::vector<art::Ptr<recob::Hit>>>& initialTrackHits) const
612 {
613 
617 
618  // Now find the hits belonging to the track
619  if (fDebug > 1)
620  std::cout << " ------------------ Finding initial track hits "
621  "-------------------- \n";
622  initialTrackHits = FindShowerStart_(showerHitsMap);
623  if (fDebug > 1) {
624  std::cout << "Here are the initial shower hits... \n";
625  for (auto const& [plane, hitPtrs] : initialTrackHits) {
626  std::cout << " Plane " << plane << '\n';
627  for (auto const& hit : hitPtrs | views::transform(to_element)) {
628  std::cout << " Hit is (" << HitCoordinates_(hit).X() << " (real hit " << hit.WireID()
629  << "), " << HitCoordinates_(hit).Y() << ")\n";
630  }
631  }
632  }
633  if (fDebug > 1)
634  std::cout << " ------------------ End finding initial track hits "
635  "-------------------- \n";
636 
637  // Now we have the track hits -- can make a track!
638  if (fDebug > 1) std::cout << " ------------------ Finding initial track -------------------- \n";
639  initialTrack = MakeInitialTrack_(detProp, initialTrackHits, showerHitsMap);
640  if (initialTrack and fDebug > 1) {
641  std::cout << "The track start is (" << initialTrack->Vertex().X() << ", "
642  << initialTrack->Vertex().Y() << ", " << initialTrack->Vertex().Z() << ")\n";
643  std::cout << "The track direction is (" << initialTrack->VertexDirection().X() << ", "
644  << initialTrack->VertexDirection().Y() << ", " << initialTrack->VertexDirection().Z()
645  << ")\n";
646  }
647  if (fDebug > 1)
648  std::cout << " ------------------ End finding initial track "
649  "-------------------- \n";
650 }
651 
652 std::vector<art::Ptr<recob::Hit>> shower::EMShowerAlg::FindOrderOfHits_(
653  detinfo::DetectorPropertiesData const& detProp,
655  bool perpendicular) const
656 {
657  // Find the charge-weighted centre (in [cm]) of this shower
658  TVector2 centre = ShowerCenter_(detProp, hits);
659 
660  // Find a rough shower 'direction'
661  TVector2 direction = ShowerDirection_(detProp, hits);
662 
663  if (perpendicular) direction = direction.Rotate(TMath::Pi() / 2);
664 
665  // Find how far each hit (projected onto this axis) is from the centre
666  TVector2 pos;
667  std::map<double, art::Ptr<recob::Hit>> hitProjection;
668  for (auto const& hitPtr : hits) {
669  pos = HitPosition_(detProp, *hitPtr) - centre;
670  hitProjection[direction * pos] = hitPtr;
671  }
672 
673  // Get a vector of hits in order of the shower
674  std::vector<art::Ptr<recob::Hit>> showerHits;
675  cet::transform_all(
676  hitProjection, std::back_inserter(showerHits), [](auto const& pr) { return pr.second; });
677 
678  // Make gradient plot
679  if (fMakeGradientPlot) {
680  std::map<int, TGraph*> graphs;
681  for (auto const& hit : showerHits | views::transform(to_element)) {
682  int tpc = hit.WireID().TPC;
683  if (graphs[tpc] == nullptr) graphs[tpc] = new TGraph();
684  graphs[tpc]->SetPoint(
685  graphs[tpc]->GetN(), HitPosition_(detProp, hit).X(), HitPosition_(detProp, hit).Y());
686  }
687  TMultiGraph* multigraph = new TMultiGraph();
688  for (auto const [color, graph] : graphs) {
689  graph->SetMarkerColor(color);
690  graph->SetMarkerStyle(8);
691  graph->SetMarkerSize(2);
692  multigraph->Add(graph);
693  }
694  TCanvas* canvas = new TCanvas();
695  multigraph->Draw("AP");
696  TLine line;
697  line.SetLineColor(2);
698  line.DrawLine(centre.X() - 1000 * direction.X(),
699  centre.Y() - 1000 * direction.Y(),
700  centre.X() + 1000 * direction.X(),
701  centre.Y() + 1000 * direction.Y());
702  canvas->SaveAs("Gradient.png");
703  delete canvas;
704  delete multigraph;
705  }
706 
707  return showerHits;
708 }
709 
710 std::vector<std::vector<int>> shower::EMShowerAlg::FindShowers(
711  std::map<int, std::vector<int>> const& trackToClusters) const
712 {
713  // Showers are vectors of clusters
714  std::vector<std::vector<int>> showers;
715 
716  // Loop over all tracks
717  for (auto const& clusters : trackToClusters | views::values) {
718 
719  // Find which showers already made are associated with this track
720  std::vector<int> matchingShowers;
721  for (unsigned int shower = 0; shower < showers.size(); ++shower)
722  for (int const cluster : clusters) {
723  if (cet::search_all(showers[shower], cluster) and
724  not cet::search_all(matchingShowers, shower)) {
725  matchingShowers.push_back(shower);
726  }
727  }
728 
729  // THINK THERE PROBABLY CAN BE MORE THAN ONE!
730  // IN FACT, THIS WOULD BE A SUCCESS OF THE SHOWERING METHOD!
731  // // Shouldn't be more than one
732  // if (matchingShowers.size() > 1)
733  // mf::LogInfo("EMShowerAlg") << "Number of showers this track matches is " << matchingShowers.size() << std::endl;
734 
735  // New shower
736  if (matchingShowers.size() < 1) showers.push_back(clusters);
737 
738  // Add to existing shower
739  else {
740  for (int const cluster : clusters) {
741  if (not cet::search_all(showers.at(matchingShowers[0]), cluster))
742  showers.at(matchingShowers.at(0)).push_back(cluster);
743  }
744  }
745  }
746 
747  return showers;
748 }
749 
750 std::map<int, std::vector<art::Ptr<recob::Hit>>> shower::EMShowerAlg::FindShowerStart_(
751  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& orderedShowerMap) const
752 {
753 
754  std::map<int, std::vector<art::Ptr<recob::Hit>>> initialHitsMap;
755 
756  for (auto const& [plane, orderedShower] : orderedShowerMap) {
757  std::vector<art::Ptr<recob::Hit>> initialHits;
758 
759  // Find if the shower is traveling along ticks or wires
760  bool wireDirection = true;
761  std::vector<int> wires;
762  for (auto const& hit : orderedShower | views::transform(to_element))
763  wires.push_back(std::round(HitCoordinates_(hit).X()));
764 
765  cet::sort_all(wires);
766  if (std::abs(*wires.begin() - std::round(HitCoordinates_(**orderedShower.begin()).X())) > 3 and
767  std::abs(*wires.rbegin() - std::round(HitCoordinates_(**orderedShower.begin()).X())) > 3)
768  wireDirection = false;
769 
770  // Deal with showers traveling along wires
771  if (wireDirection) {
772  bool increasing = HitCoordinates_(**orderedShower.rbegin()).X() >
773  HitCoordinates_(**orderedShower.begin()).X();
774  std::map<int, std::vector<art::Ptr<recob::Hit>>> wireHitMap;
775  int multipleWires = 0;
776  for (auto const& hitPtr : orderedShower)
777  wireHitMap[std::round(HitCoordinates_(*hitPtr).X())].push_back(hitPtr);
778 
779  for (auto const& hitPtr : orderedShower) {
780  int wire = std::round(HitCoordinates_(*hitPtr).X());
781  if (wireHitMap[wire].size() > 1) {
782  ++multipleWires;
783  if (multipleWires > 5) break;
784  continue;
785  }
786  else if ((increasing and wireHitMap[wire + 1].size() > 1 and
787  (wireHitMap[wire + 2].size() > 1 or wireHitMap[wire + 3].size() > 1)) or
788  (!increasing and wireHitMap[wire - 1].size() > 1 and
789  (wireHitMap[wire - 2].size() > 1 or wireHitMap[wire - 3].size() > 1))) {
790  if ((increasing and
791  (wireHitMap[wire + 4].size() < 2 and wireHitMap[wire + 5].size() < 2 and
792  wireHitMap[wire + 6].size() < 2 and wireHitMap[wire + 9].size() > 1)) or
793  (!increasing and
794  (wireHitMap[wire - 4].size() < 2 and wireHitMap[wire - 5].size() < 2 and
795  wireHitMap[wire - 6].size() < 2) and
796  wireHitMap[wire - 9].size() > 1))
797  initialHits.push_back(hitPtr);
798  else
799  break;
800  }
801  else
802  initialHits.push_back(hitPtr);
803  }
804  if (!initialHits.size()) initialHits.push_back(*orderedShower.begin());
805  }
806 
807  // Deal with showers travelling along ticks
808  else {
809  bool increasing = HitCoordinates_(**orderedShower.rbegin()).Y() >
810  HitCoordinates_(**orderedShower.begin()).Y();
811  std::map<int, std::vector<art::Ptr<recob::Hit>>> tickHitMap;
812  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hitIt = orderedShower.begin();
813  hitIt != orderedShower.end();
814  ++hitIt)
815  tickHitMap[std::round(HitCoordinates_(**hitIt).Y())].push_back(*hitIt);
816 
817  for (auto const& hitPtr : orderedShower) {
818  int const tick = std::round(HitCoordinates_(*hitPtr).Y());
819  if ((increasing and (tickHitMap[tick + 1].size() or tickHitMap[tick + 2].size() or
820  tickHitMap[tick + 3].size() or tickHitMap[tick + 4].size() or
821  tickHitMap[tick + 5].size())) or
822  (!increasing and (tickHitMap[tick - 1].size() or tickHitMap[tick - 2].size() or
823  tickHitMap[tick - 3].size() or tickHitMap[tick - 4].size() or
824  tickHitMap[tick - 5].size())))
825  break;
826  else
827  initialHits.push_back(hitPtr);
828  }
829  if (initialHits.empty()) initialHits.push_back(*orderedShower.begin());
830  }
831 
832  // Need at least two hits to make a track
833  if (initialHits.size() == 1 and orderedShower.size() > 2)
834  initialHits.push_back(orderedShower[1]);
835 
836  // Quality check -- make sure there isn't a single hit in a different TPC (artefact of tracking failure)
837  std::vector<art::Ptr<recob::Hit>> newInitialHits;
838  std::map<int, int> tpcHitMap;
839  std::vector<int> singleHitTPCs;
840  for (auto const& hit : initialHits | views::transform(to_element))
841  ++tpcHitMap[hit.WireID().TPC];
842 
843  for (auto const [tpc, count] : tpcHitMap)
844  if (count == 1) singleHitTPCs.push_back(tpc);
845 
846  if (singleHitTPCs.size()) {
847  if (fDebug > 2)
848  for (int const tpc : singleHitTPCs)
849  std::cout << "Removed hits in TPC " << tpc << '\n';
850 
851  for (auto const& hitPtr : initialHits)
852  if (not cet::search_all(singleHitTPCs, hitPtr->WireID().TPC))
853  newInitialHits.push_back(hitPtr);
854  if (!newInitialHits.size()) newInitialHits.push_back(*orderedShower.begin());
855  }
856  else
857  newInitialHits = initialHits;
858 
859  initialHitsMap[plane] = newInitialHits;
860  }
861 
862  return initialHitsMap;
863 }
864 
865 std::map<int, std::map<int, bool>> shower::EMShowerAlg::GetPlanePermutations_(
866  const detinfo::DetectorPropertiesData& detProp,
867  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
868 {
869 
870  // The map to return
871  std::map<int, std::map<int, bool>> permutations;
872 
873  // Get the properties of the shower hits across the planes which will be used to
874  // determine the likelihood of a particular reorientation permutation
875  // -- relative width in the wire direction (if showers travel along the wire
876  // direction in a particular plane)
877  // -- the RMS gradients (how likely it is the RMS of the hit positions from
878  // central axis increases along a particular orientation)
879 
880  // Find the RMS, RMS gradient and wire widths
881  std::map<int, double> planeRMSGradients, planeRMS;
882  for (auto const& [plane, hitPtrs] : showerHitsMap) {
883  planeRMS[plane] = ShowerHitRMS_(detProp, hitPtrs);
884  planeRMSGradients[plane] = ShowerHitRMSGradient_(detProp, hitPtrs);
885  }
886 
887  // Order these backwards so they can be used to discriminate between planes
888  std::map<double, int> gradientMap;
889  for (int const plane : showerHitsMap | views::keys)
890  gradientMap[std::abs(planeRMSGradients.at(plane))] = plane;
891 
892  std::map<double, int> wireWidthMap = RelativeWireWidth_(showerHitsMap);
893 
894  if (fDebug > 1)
895  for (auto const [gradient, plane] : wireWidthMap)
896  std::cout << "Plane " << plane << " has relative width in wire of " << gradient
897  << "; and an RMS gradient of " << planeRMSGradients[plane] << '\n';
898 
899  // Find the correct permutations
900  int perm = 0;
901  std::vector<std::map<int, bool>> usedPermutations;
902 
903  // Most likely is to not change anything
904  for (int const plane : showerHitsMap | views::keys)
905  permutations[perm][plane] = 0;
906  ++perm;
907 
908  // Use properties of the shower to determine the middle cases
909  for (int const plane : wireWidthMap | views::values) {
910  std::map<int, bool> permutation;
911  permutation[plane] = true;
912  for (int const plane2 : wireWidthMap | views::values)
913  if (plane != plane2) permutation[plane2] = false;
914 
915  if (not cet::search_all(usedPermutations, permutation)) {
916  permutations[perm] = permutation;
917  usedPermutations.push_back(permutation);
918  ++perm;
919  }
920  }
921  for (int const plane : wireWidthMap | views::reverse | views::values) {
922  std::map<int, bool> permutation;
923  permutation[plane] = false;
924  for (int const plane2 : wireWidthMap | views::values)
925  if (plane != plane2) permutation[plane2] = true;
926 
927  if (not cet::search_all(usedPermutations, permutation)) {
928  permutations[perm] = permutation;
929  usedPermutations.push_back(permutation);
930  ++perm;
931  }
932  }
933 
934  // Least likely is to change everything
935  for (int const plane : showerHitsMap | views::keys)
936  permutations[perm][plane] = 1;
937  ++perm;
938 
939  if (fDebug > 1) {
940  std::cout << "Here are the permutations!\n";
941  for (auto const& [index, permutation] : permutations) {
942  std::cout << "Permutation " << index << '\n';
943  for (auto const [plane, value] : permutation)
944  std::cout << " Plane " << plane << " has value " << value << '\n';
945  }
946  }
947 
948  return permutations;
949 }
950 
951 std::unique_ptr<recob::Track> shower::EMShowerAlg::MakeInitialTrack_(
952  detinfo::DetectorPropertiesData const& detProp,
953  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& initialHitsMap,
954  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& showerHitsMap) const
955 {
956  // Can't do much with just one view
957  if (initialHitsMap.size() < 2) {
958  mf::LogInfo("EMShowerAlg") << "Only one useful view for this shower; nothing can be done\n";
959  return std::unique_ptr<recob::Track>();
960  }
961 
962  std::vector<std::pair<int, int>> initialHitsSize;
963  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator initialHitIt =
964  initialHitsMap.begin();
965  initialHitIt != initialHitsMap.end();
966  ++initialHitIt)
967  initialHitsSize.push_back(std::make_pair(initialHitIt->first, initialHitIt->second.size()));
968 
969  // Sort the planes by number of hits
970  std::sort(initialHitsSize.begin(),
971  initialHitsSize.end(),
972  [](std::pair<int, int> const& size1, std::pair<int, int> const& size2) {
973  return size1.second > size2.second;
974  });
975 
976  // Pick the two planes to use to make the track
977  // -- if more than two planes, can choose the two views
978  // more accurately
979  // -- if not, just use the two views available
980 
981  std::vector<int> planes;
982 
983  if (initialHitsSize.size() > 2 and !CheckShowerHits_(detProp, showerHitsMap)) {
984  int planeToIgnore = WorstPlane_(showerHitsMap);
985  if (fDebug > 1)
986  std::cout << "Igoring plane " << planeToIgnore << " in creation of initial track\n";
987  for (std::vector<std::pair<int, int>>::iterator hitsSizeIt = initialHitsSize.begin();
988  hitsSizeIt != initialHitsSize.end() and planes.size() < 2;
989  ++hitsSizeIt) {
990  if (hitsSizeIt->first == planeToIgnore) continue;
991  planes.push_back(hitsSizeIt->first);
992  }
993  }
994  else
995  planes = {initialHitsSize.at(0).first, initialHitsSize.at(1).first};
996 
997  return ConstructTrack(detProp, initialHitsMap.at(planes.at(0)), initialHitsMap.at(planes.at(1)));
998 }
999 
1001  detinfo::DetectorClocksData const& clockData,
1002  detinfo::DetectorPropertiesData const& detProp,
1004  std::unique_ptr<recob::Track> const& initialTrack,
1005  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& initialHitsMap) const
1006 {
1007 
1008  // Find the shower hits on each plane
1009  std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
1010  for (auto const& hitPtr : hits)
1011  planeHitsMap[hitPtr->View()].push_back(hitPtr);
1012 
1013  int bestPlane = -1;
1014  unsigned int highestNumberOfHits = 0;
1015  std::vector<double> totalEnergy, totalEnergyError, dEdx, dEdxError;
1016 
1017  // Look at each of the planes
1018  for (unsigned int plane = 0; plane < fGeom->MaxPlanes(); ++plane) {
1019 
1020  // If there's hits on this plane...
1021  if (planeHitsMap.count(plane)) {
1022  dEdx.push_back(FinddEdx_(clockData, detProp, initialHitsMap.at(plane), initialTrack));
1023  totalEnergy.push_back(
1024  fShowerEnergyAlg.ShowerEnergy(clockData, detProp, planeHitsMap.at(plane), plane));
1025  if (planeHitsMap.at(plane).size() > highestNumberOfHits and initialHitsMap.count(plane)) {
1026  bestPlane = plane;
1027  highestNumberOfHits = planeHitsMap.at(plane).size();
1028  }
1029  }
1030 
1031  // If not...
1032  else {
1033  dEdx.push_back(0);
1034  totalEnergy.push_back(0);
1035  }
1036  }
1037 
1038  TVector3 direction, directionError, showerStart, showerStartError;
1039  if (initialTrack) {
1040  direction = initialTrack->VertexDirection<TVector3>();
1041  showerStart = initialTrack->Vertex<TVector3>();
1042  }
1043 
1044  if (fDebug > 0) {
1045  std::cout << "Best plane is " << bestPlane;
1046  std::cout << "\ndE/dx for each plane is: " << dEdx[0] << ", " << dEdx[1] << " and " << dEdx[2];
1047  std::cout << "\nTotal energy for each plane is: " << totalEnergy[0] << ", " << totalEnergy[1]
1048  << " and " << totalEnergy[2];
1049  std::cout << "\nThe shower start is (" << showerStart.X() << ", " << showerStart.Y() << ", "
1050  << showerStart.Z() << ")\n";
1051  std::cout << "The shower direction is (" << direction.X() << ", " << direction.Y() << ", "
1052  << direction.Z() << ")\n";
1053  }
1054 
1055  return recob::Shower(direction,
1056  directionError,
1057  showerStart,
1058  showerStartError,
1059  totalEnergy,
1060  totalEnergyError,
1061  dEdx,
1062  dEdxError,
1063  bestPlane);
1064 }
1065 
1067  detinfo::DetectorPropertiesData const& detProp,
1070  int& iok) const
1071 {
1072  iok = 1;
1073 
1074  // Find the shower hits on each plane
1075  std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
1076  for (auto const& hitPtr : hits)
1077  planeHitsMap[hitPtr->WireID().Plane].push_back(hitPtr);
1078 
1079  std::vector<std::vector<art::Ptr<recob::Hit>>> initialTrackHits(3);
1080 
1081  int pl0 = -1;
1082  int pl1 = -1;
1083  unsigned maxhits0 = 0;
1084  unsigned maxhits1 = 0;
1085 
1086  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator planeHits = planeHitsMap.begin();
1087  planeHits != planeHitsMap.end();
1088  ++planeHits) {
1089 
1090  std::vector<art::Ptr<recob::Hit>> showerHits;
1091  OrderShowerHits_(detProp, planeHits->second, showerHits, vertex);
1092  FindInitialTrackHits(showerHits, vertex, initialTrackHits[planeHits->first]);
1093  if ((planeHits->second).size() > maxhits0) {
1094  if (pl0 != -1) {
1095  maxhits1 = maxhits0;
1096  pl1 = pl0;
1097  }
1098  pl0 = planeHits->first;
1099  maxhits0 = (planeHits->second).size();
1100  }
1101  else if ((planeHits->second).size() > maxhits1) {
1102  pl1 = planeHits->first;
1103  maxhits1 = (planeHits->second).size();
1104  }
1105  }
1106  if (pl0 != -1 && pl1 != -1 && initialTrackHits[pl0].size() >= 2 &&
1107  initialTrackHits[pl1].size() >= 2 &&
1108  initialTrackHits[pl0][0]->WireID().TPC == initialTrackHits[pl1][0]->WireID().TPC) {
1109  double xyz[3];
1110  vertex->XYZ(xyz);
1111  TVector3 vtx(xyz);
1112  pma::Track3D* pmatrack =
1113  fProjectionMatchingAlg.buildSegment(detProp, initialTrackHits[pl0], initialTrackHits[pl1]);
1114  std::vector<TVector3> spts;
1115 
1116  for (size_t i = 0; i < pmatrack->size(); ++i) {
1117  if ((*pmatrack)[i]->IsEnabled()) {
1118  TVector3 p3d = (*pmatrack)[i]->Point3D();
1119  spts.push_back(p3d);
1120  }
1121  }
1122  if (spts.size() >= 2) { // at least two space points
1123  TVector3 shwxyz, shwxyzerr;
1124  TVector3 shwdir, shwdirerr;
1125  std::vector<double> totalEnergy, totalEnergyError, dEdx, dEdxError;
1126  int bestPlane = pl0;
1127  double minpitch = 1000;
1128  std::vector<TVector3> dirs;
1129  if ((spts[0] - vtx).Mag() < (spts.back() - vtx).Mag()) {
1130  shwxyz = spts[0];
1131  size_t i = 5;
1132  if (spts.size() - 1 < 5) i = spts.size() - 1;
1133  shwdir = spts[i] - spts[0];
1134  shwdir = shwdir.Unit();
1135  }
1136  else {
1137  shwxyz = spts.back();
1138  size_t i = 0;
1139  if (spts.size() > 6) i = spts.size() - 6;
1140  shwdir = spts[i] - spts[spts.size() - 1];
1141  shwdir = shwdir.Unit();
1142  }
1143  for (unsigned int plane = 0; plane < fGeom->MaxPlanes(); ++plane) {
1144  if (planeHitsMap.find(plane) != planeHitsMap.end()) {
1145  totalEnergy.push_back(
1146  fShowerEnergyAlg.ShowerEnergy(clockData, detProp, planeHitsMap[plane], plane));
1147  }
1148  else {
1149  totalEnergy.push_back(0);
1150  }
1151  if (initialTrackHits[plane].size()) {
1152  double fdEdx = 0;
1153  double totQ = 0;
1154  double avgT = 0;
1155  double pitch = 0;
1156  double wirepitch = fGeom->WirePitch(initialTrackHits[plane][0]->WireID().planeID());
1157  double angleToVert =
1159  initialTrackHits[plane][0]->WireID().planeID()) -
1160  0.5 * TMath::Pi();
1161  double cosgamma = std::abs(sin(angleToVert) * shwdir.Y() + cos(angleToVert) * shwdir.Z());
1162  if (cosgamma > 0) pitch = wirepitch / cosgamma;
1163  if (pitch) {
1164  if (pitch < minpitch) {
1165  minpitch = pitch;
1166  bestPlane = plane;
1167  }
1168  int nhits = 0;
1169  std::vector<float> vQ;
1170  for (auto const& hit : initialTrackHits[plane]) {
1171  int w1 = hit->WireID().Wire;
1172  int w0 = initialTrackHits[plane][0]->WireID().Wire;
1173  if (std::abs((w1 - w0) * pitch) < fdEdxTrackLength) {
1174  vQ.push_back(hit->Integral());
1175  totQ += hit->Integral();
1176  avgT += hit->PeakTime();
1177  ++nhits;
1178  }
1179  }
1180  if (totQ) {
1181  double dQdx = TMath::Median(vQ.size(), &vQ[0]) / pitch;
1182  fdEdx = fCalorimetryAlg.dEdx_AREA(
1183  clockData, detProp, dQdx, avgT / nhits, initialTrackHits[plane][0]->WireID().Plane);
1184  }
1185  }
1186  dEdx.push_back(fdEdx);
1187  }
1188  else {
1189  dEdx.push_back(0);
1190  }
1191  }
1192  iok = 0;
1193  if (fDebug > 1) {
1194  std::cout << "Best plane is " << bestPlane;
1195  std::cout << "\ndE/dx for each plane is: " << dEdx[0] << ", " << dEdx[1] << " and "
1196  << dEdx[2];
1197  std::cout << "\nTotal energy for each plane is: " << totalEnergy[0] << ", "
1198  << totalEnergy[1] << " and " << totalEnergy[2];
1199  std::cout << "\nThe shower start is (" << shwxyz.X() << ", " << shwxyz.Y() << ", "
1200  << shwxyz.Z() << ")\n";
1201  shwxyz.Print();
1202  }
1203 
1204  return recob::Shower(shwdir,
1205  shwdirerr,
1206  shwxyz,
1207  shwxyzerr,
1208  totalEnergy,
1209  totalEnergyError,
1210  dEdx,
1211  dEdxError,
1212  bestPlane);
1213  }
1214  }
1215  return recob::Shower();
1216 }
1217 
1218 std::vector<recob::SpacePoint> shower::EMShowerAlg::MakeSpacePoints(
1219  detinfo::DetectorPropertiesData const& detProp,
1220  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& showerHits,
1221  std::vector<std::vector<art::Ptr<recob::Hit>>>& hitAssns) const
1222 {
1223  // Space points to return
1224  std::vector<recob::SpacePoint> spacePoints;
1225 
1226  // Make space points
1227  // Use the following procedure:
1228  // -- Consider hits plane by plane
1229  // -- For each hit on the first plane, consider the 3D point made by combining with each hit from the second plane
1230  // -- Project this 3D point back into the two planes
1231  // -- Determine how close to a the original hits this point lies
1232  // -- If close enough, make a 3D space point from this point
1233  // -- Discard these used hits in future iterations, along with hits in the
1234  // third plane (if exists) close to the projection of the point into this
1235  // plane
1236 
1237  // Container to hold used hits
1238  std::vector<int> usedHits;
1239 
1240  // Look through plane by plane
1241  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitIt =
1242  showerHits.begin();
1243  showerHitIt != showerHits.end();
1244  ++showerHitIt) {
1245 
1246  // Find the other planes with hits
1247  std::vector<int> otherPlanes;
1248  for (unsigned int otherPlane = 0; otherPlane < fGeom->MaxPlanes(); ++otherPlane)
1249  if ((int)otherPlane != showerHitIt->first and showerHits.count(otherPlane))
1250  otherPlanes.push_back(otherPlane);
1251 
1252  // Can't make space points if we only have one view
1253  if (otherPlanes.size() == 0) return spacePoints;
1254 
1255  // Look at all hits on this plane
1256  for (std::vector<art::Ptr<recob::Hit>>::const_iterator planeHitIt = showerHitIt->second.begin();
1257  planeHitIt != showerHitIt->second.end();
1258  ++planeHitIt) {
1259 
1260  if (std::find(usedHits.begin(), usedHits.end(), planeHitIt->key()) != usedHits.end())
1261  continue;
1262 
1263  // Make a 3D point with every hit on the second plane
1264  const std::vector<art::Ptr<recob::Hit>> otherPlaneHits = showerHits.at(otherPlanes.at(0));
1265  for (std::vector<art::Ptr<recob::Hit>>::const_iterator otherPlaneHitIt =
1266  otherPlaneHits.begin();
1267  otherPlaneHitIt != otherPlaneHits.end() and
1268  std::find(usedHits.begin(), usedHits.end(), planeHitIt->key()) == usedHits.end();
1269  ++otherPlaneHitIt) {
1270 
1271  if ((*otherPlaneHitIt)->WireID().TPC != (*planeHitIt)->WireID().TPC or
1272  std::find(usedHits.begin(), usedHits.end(), otherPlaneHitIt->key()) != usedHits.end())
1273  continue;
1274 
1275  auto const point = Construct3DPoint_(detProp, *planeHitIt, *otherPlaneHitIt);
1276  std::vector<art::Ptr<recob::Hit>> pointHits;
1277  bool truePoint = false;
1278 
1279  if (otherPlanes.size() > 1) {
1280 
1281  TVector2 projThirdPlane = Project3DPointOntoPlane_(detProp, point, otherPlanes.at(1));
1282  const std::vector<art::Ptr<recob::Hit>> otherOtherPlaneHits =
1283  showerHits.at(otherPlanes.at(1));
1284 
1285  for (std::vector<art::Ptr<recob::Hit>>::const_iterator otherOtherPlaneHitIt =
1286  otherOtherPlaneHits.begin();
1287  otherOtherPlaneHitIt != otherOtherPlaneHits.end() and !truePoint;
1288  ++otherOtherPlaneHitIt) {
1289 
1290  if ((*otherOtherPlaneHitIt)->WireID().TPC == (*planeHitIt)->WireID().TPC and
1291  (projThirdPlane - HitPosition_(detProp, **otherOtherPlaneHitIt)).Mod() <
1292  fSpacePointSize) {
1293 
1294  truePoint = true;
1295 
1296  // Remove hits used to make the point
1297  usedHits.push_back(planeHitIt->key());
1298  usedHits.push_back(otherPlaneHitIt->key());
1299  usedHits.push_back(otherOtherPlaneHitIt->key());
1300 
1301  pointHits.push_back(*planeHitIt);
1302  pointHits.push_back(*otherPlaneHitIt);
1303  pointHits.push_back(*otherOtherPlaneHitIt);
1304  }
1305  }
1306  }
1307 
1308  else if ((Project3DPointOntoPlane_(detProp, point, (*planeHitIt)->WireID().Plane) -
1309  HitPosition_(detProp, **planeHitIt))
1310  .Mod() < fSpacePointSize and
1311  (Project3DPointOntoPlane_(detProp, point, (*otherPlaneHitIt)->WireID().Plane) -
1312  HitPosition_(detProp, **otherPlaneHitIt))
1313  .Mod() < fSpacePointSize) {
1314 
1315  truePoint = true;
1316 
1317  usedHits.push_back(planeHitIt->key());
1318  usedHits.push_back(otherPlaneHitIt->key());
1319 
1320  pointHits.push_back(*planeHitIt);
1321  pointHits.push_back(*otherPlaneHitIt);
1322  }
1323 
1324  // Make space point
1325  if (truePoint) {
1326  double xyz[3] = {point.X(), point.Y(), point.Z()};
1327  double xyzerr[6] = {fSpacePointSize,
1332  fSpacePointSize};
1333  double chisq = 0.;
1334  spacePoints.emplace_back(xyz, xyzerr, chisq);
1335  hitAssns.push_back(pointHits);
1336  }
1337 
1338  } // end loop over second plane hits
1339 
1340  } // end loop over first plane hits
1341 
1342  } // end loop over planes
1343 
1344  if (fDebug > 0) {
1345  std::cout << "-------------------- Space points -------------------\n";
1346  std::cout << "There are " << spacePoints.size() << " space points:\n";
1347  if (fDebug > 1)
1348  for (std::vector<recob::SpacePoint>::const_iterator spacePointIt = spacePoints.begin();
1349  spacePointIt != spacePoints.end();
1350  ++spacePointIt) {
1351  const double* xyz = spacePointIt->XYZ();
1352  std::cout << " Space point (" << xyz[0] << ", " << xyz[1] << ", " << xyz[2] << ")\n";
1353  }
1354  }
1355 
1356  return spacePoints;
1357 }
1358 
1359 std::map<int, std::vector<art::Ptr<recob::Hit>>> shower::EMShowerAlg::OrderShowerHits(
1360  detinfo::DetectorPropertiesData const& detProp,
1362  int desired_plane) const
1363 {
1368 
1369  // ------------- Put hits in order ------------
1370 
1371  // Find the shower hits on each plane
1372  std::map<int, std::vector<art::Ptr<recob::Hit>>> showerHitsMap;
1373  for (auto const& hitPtr : shower) {
1374  showerHitsMap[hitPtr->WireID().Plane].push_back(hitPtr);
1375  }
1376 
1377  // Order the hits, get the RMS and the RMS gradient for the hits in this plane
1378  std::map<int, double> planeRMSGradients, planeRMS;
1379  for (auto const& [plane, hits] : showerHitsMap) {
1380  if (desired_plane != plane and desired_plane != -1) continue;
1381  std::vector<art::Ptr<recob::Hit>> orderedHits = FindOrderOfHits_(detProp, hits);
1382  planeRMS[plane] = ShowerHitRMS_(detProp, orderedHits);
1383  planeRMSGradients[plane] = ShowerHitRMSGradient_(detProp, orderedHits);
1384  showerHitsMap[plane] = orderedHits;
1385  }
1386 
1387  if (fDebug > 1) {
1388  for (auto const [plane, shower_hit_rms] : planeRMS) {
1389  std::cout << "Plane " << plane << " has RMS " << shower_hit_rms << " and RMS gradient "
1390  << planeRMSGradients.at(plane) << '\n';
1391  }
1392  }
1393 
1394  if (fDebug > 2) {
1395  std::cout << "\nHits in order; after ordering, before reversing...\n";
1396  for (auto const& [plane, hitPtrs] : showerHitsMap) {
1397  std::cout << " Plane " << plane << '\n';
1398  for (auto const& hit : hitPtrs | views::transform(to_element)) {
1399  std::cout << " Hit at (" << HitCoordinates_(hit).X() << ", " << HitCoordinates_(hit).Y()
1400  << ") -- real wire " << hit.WireID() << ", hit position ("
1401  << HitPosition_(detProp, hit).X() << ", " << HitPosition_(detProp, hit).Y()
1402  << ")\n";
1403  }
1404  }
1405  }
1406 
1407  // ------------- Check between the views to ensure consistency of ordering -------------
1408 
1409  // Check between the views to make sure there isn't a poorly formed shower in just one view
1410  // First, determine the average RMS and RMS gradient across the other planes
1411  std::map<int, double> planeOtherRMS, planeOtherRMSGradients;
1412  for (std::map<int, double>::iterator planeRMSIt = planeRMS.begin(); planeRMSIt != planeRMS.end();
1413  ++planeRMSIt) {
1414  planeOtherRMS[planeRMSIt->first] = 0;
1415  planeOtherRMSGradients[planeRMSIt->first] = 0;
1416  int nOtherPlanes = 0;
1417  for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane) {
1418  if (plane != planeRMSIt->first and planeRMS.count(plane)) {
1419  planeOtherRMS[planeRMSIt->first] += planeRMS.at(plane);
1420  planeOtherRMSGradients[planeRMSIt->first] += planeRMSGradients.at(plane);
1421  ++nOtherPlanes;
1422  }
1423  }
1424  planeOtherRMS[planeRMSIt->first] /= (double)nOtherPlanes;
1425  planeOtherRMSGradients[planeRMSIt->first] /= (double)nOtherPlanes;
1426  }
1427 
1428  // Look to see if one plane has a particularly high RMS (compared to the
1429  // others) whilst having a similar gradient
1430  for (auto const& [plane, hitPtrs] : showerHitsMap) {
1431  if (planeRMS.at(plane) > planeOtherRMS.at(plane) * 2) {
1432  if (fDebug > 1) std::cout << "Plane " << plane << " was perpendicular... recalculating\n";
1433  std::vector<art::Ptr<recob::Hit>> orderedHits =
1434  this->FindOrderOfHits_(detProp, hitPtrs, true);
1435  showerHitsMap[plane] = orderedHits;
1436  planeRMSGradients[plane] = this->ShowerHitRMSGradient_(detProp, orderedHits);
1437  }
1438  }
1439 
1440  // ------------- Orient the shower correctly ---------------
1441 
1442  if (fDebug > 1) {
1443  std::cout << "Before reversing, here are the start and end points...\n";
1444  for (auto const& [plane, hitPtrs] : showerHitsMap) {
1445  std::cout << " Plane " << plane << " has start (" << HitCoordinates_(*hitPtrs.front()).X()
1446  << ", " << HitCoordinates_(*hitPtrs.front()).Y() << ") (real wire "
1447  << hitPtrs.front()->WireID() << ") and end ("
1448  << HitCoordinates_(*hitPtrs.back()).X() << ", "
1449  << HitCoordinates_(*hitPtrs.back()).Y() << ") (real wire "
1450  << hitPtrs.back()->WireID() << ")\n";
1451  }
1452  }
1453 
1454  // Use the RMS gradient information to get an initial ordering
1455  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
1456  showerHitsMap.begin();
1457  showerHitsIt != showerHitsMap.end();
1458  ++showerHitsIt) {
1459  double gradient = planeRMSGradients.at(showerHitsIt->first);
1460  if (gradient < 0) std::reverse(showerHitsIt->second.begin(), showerHitsIt->second.end());
1461  }
1462 
1463  if (fDebug > 2) {
1464  std::cout << "\nHits in order; after reversing, before checking isolated hits...\n";
1465  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
1466  showerHitsMap.begin();
1467  showerHitsIt != showerHitsMap.end();
1468  ++showerHitsIt) {
1469  std::cout << " Plane " << showerHitsIt->first << '\n';
1470  for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsIt->second.begin();
1471  hitIt != showerHitsIt->second.end();
1472  ++hitIt)
1473  std::cout << " Hit at (" << HitCoordinates_(**hitIt).X() << ", "
1474  << HitCoordinates_(**hitIt).Y() << ") -- real wire " << (*hitIt)->WireID()
1475  << ", hit position (" << HitPosition_(detProp, **hitIt).X() << ", "
1476  << HitPosition_(detProp, **hitIt).Y() << ")\n";
1477  }
1478  }
1479 
1480  CheckIsolatedHits_(showerHitsMap);
1481 
1482  if (fDebug > 2) {
1483  std::cout << "\nHits in order; after checking isolated hits...\n";
1484  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
1485  showerHitsMap.begin();
1486  showerHitsIt != showerHitsMap.end();
1487  ++showerHitsIt) {
1488  std::cout << " Plane " << showerHitsIt->first << '\n';
1489  for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsIt->second.begin();
1490  hitIt != showerHitsIt->second.end();
1491  ++hitIt)
1492  std::cout << " Hit at (" << HitCoordinates_(**hitIt).X() << ", "
1493  << HitCoordinates_(**hitIt).Y() << ") -- real wire " << (*hitIt)->WireID()
1494  << ", hit position (" << HitPosition_(detProp, **hitIt).X() << ", "
1495  << HitPosition_(detProp, **hitIt).Y() << ")\n";
1496  }
1497  }
1498 
1499  if (fDebug > 1) {
1500  std::cout << "After reversing and checking isolated hits, here are the "
1501  "start and end points...\n";
1502  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
1503  showerHitsMap.begin();
1504  showerHitsIt != showerHitsMap.end();
1505  ++showerHitsIt)
1506  std::cout << " Plane " << showerHitsIt->first << " has start ("
1507  << HitCoordinates_(*showerHitsIt->second.front()).X() << ", "
1508  << HitCoordinates_(*showerHitsIt->second.front()).Y() << ") (real wire "
1509  << showerHitsIt->second.front()->WireID() << ") and end ("
1510  << HitCoordinates_(*showerHitsIt->second.back()).X() << ", "
1511  << HitCoordinates_(*showerHitsIt->second.back()).Y() << ")\n";
1512  }
1513 
1514  // Check for views in which the shower travels almost along the wire planes
1515  // (shown by a small relative wire width)
1516  std::map<double, int> wireWidths = RelativeWireWidth_(showerHitsMap);
1517  std::vector<int> badPlanes;
1518  if (fDebug > 1) std::cout << "Here are the relative wire widths... \n";
1519  for (auto const [relative_wire_width, plane] : wireWidths) {
1520  if (fDebug > 1)
1521  std::cout << "Plane " << plane << " has relative wire width " << relative_wire_width << '\n';
1522  if (relative_wire_width < 1 / (double)wireWidths.size()) badPlanes.push_back(plane);
1523  }
1524 
1525  std::map<int, std::vector<art::Ptr<recob::Hit>>> ignoredPlanes;
1526  if (badPlanes.size() == 1) {
1527  int const badPlane = badPlanes[0];
1528  if (fDebug > 1) std::cout << "Ignoring plane " << badPlane << " when orientating\n";
1529  ignoredPlanes[badPlane] = showerHitsMap.at(badPlane);
1530  showerHitsMap.erase(badPlane);
1531  }
1532 
1533  // Consider all possible permutations of planes (0/1, oriented
1534  // correctly/incorrectly)
1535  std::map<int, std::map<int, bool>> permutations = GetPlanePermutations_(detProp, showerHitsMap);
1536 
1537  // Go through all permutations until we have a satisfactory orientation
1538  auto const originalShowerHitsMap = showerHitsMap;
1539 
1540  int n = 0;
1541  while (!CheckShowerHits_(detProp, showerHitsMap) and
1542  n < TMath::Power(2, (int)showerHitsMap.size())) {
1543  if (fDebug > 1) std::cout << "Permutation " << n << '\n';
1544  for (int const plane : showerHitsMap | views::keys) {
1545  auto hits = originalShowerHitsMap.at(plane);
1546  if (permutations[n][plane] == 1) { std::reverse(hits.begin(), hits.end()); }
1547  showerHitsMap[plane] = hits;
1548  }
1549  ++n;
1550  }
1551 
1552  // Go back to original if still no consistency
1553  if (!CheckShowerHits_(detProp, showerHitsMap)) showerHitsMap = originalShowerHitsMap;
1554 
1555  if (fDebug > 2) {
1556  std::cout << "End of OrderShowerHits: here are the order of hits:\n";
1557  for (auto const& [plane, hitPtrs] : showerHitsMap) {
1558  std::cout << " Plane " << plane << '\n';
1559  for (auto const& hit : hitPtrs | views::transform(to_element)) {
1560  std::cout << " Hit (" << HitCoordinates_(hit).X() << " (real wire " << hit.WireID()
1561  << "), " << HitCoordinates_(hit).Y() << ") -- pos ("
1562  << HitPosition_(detProp, hit).X() << ", " << HitPosition_(detProp, hit).Y()
1563  << ")\n";
1564  }
1565  }
1566  }
1567 
1568  if (ignoredPlanes.size())
1569  showerHitsMap[ignoredPlanes.begin()->first] = ignoredPlanes.begin()->second;
1570 
1571  return showerHitsMap;
1572 }
1573 
1576  std::vector<art::Ptr<recob::Hit>>& showerHits,
1577  art::Ptr<recob::Vertex> const& vertex) const
1578 {
1579  showerHits = FindOrderOfHits_(detProp, shower);
1580 
1581  // Find TPC for the vertex
1582  auto const& xyz = vertex->position();
1583  geo::TPCID tpc = fGeom->FindTPCAtPosition(xyz);
1584  if (!tpc.isValid && showerHits.size()) tpc = geo::TPCID(showerHits[0]->WireID());
1585 
1586  // Find hits in the same TPC
1587  art::Ptr<recob::Hit> hit0, hit1;
1588  for (auto& hit : showerHits) {
1589  if (hit->WireID().TPC == tpc.TPC) {
1590  if (hit0.isNull()) { hit0 = hit; }
1591  hit1 = hit;
1592  }
1593  }
1594  if (hit0.isNull() || hit1.isNull()) return;
1595  TVector2 coord0 = TVector2(hit0->WireID().Wire, hit0->PeakTime());
1596  TVector2 coord1 = TVector2(hit1->WireID().Wire, hit1->PeakTime());
1597  TVector2 coordvtx = TVector2(fGeom->WireCoordinate(xyz, hit0->WireID().planeID()),
1598  detProp.ConvertXToTicks(xyz.X(), hit0->WireID().planeID()));
1599  if ((coord1 - coordvtx).Mod() < (coord0 - coordvtx).Mod()) {
1600  std::reverse(showerHits.begin(), showerHits.end());
1601  }
1602 }
1603 
1606  std::vector<art::Ptr<recob::Hit>>& trackHits) const
1607 {
1608  // Find TPC for the vertex
1609  auto const& xyz = vertex->position();
1610  geo::TPCID tpc = fGeom->FindTPCAtPosition(xyz);
1611 
1612  // vertex cannot be projected into a TPC, find the TPC that has the most hits
1613  if (!tpc.isValid) {
1614  std::map<geo::TPCID, unsigned int> tpcmap;
1615  unsigned maxhits = 0;
1616  for (auto const& hit : showerHits) {
1617  ++tpcmap[geo::TPCID(hit->WireID())];
1618  }
1619  for (auto const& t : tpcmap) {
1620  if (t.second > maxhits) {
1621  maxhits = t.second;
1622  tpc = t.first;
1623  }
1624  }
1625  }
1626  if (!tpc.isValid) return;
1627 
1628  double parm[2];
1629  int fitok = 0;
1630  std::vector<double> wfit;
1631  std::vector<double> tfit;
1632  std::vector<double> cfit;
1633 
1634  for (size_t i = 0; i < fNfitpass; ++i) {
1635 
1636  // Fit a straight line through hits
1637  unsigned int nhits = 0;
1638  for (auto& hit : showerHits) {
1639  if (hit->WireID().TPC == tpc.TPC) {
1640  TVector2 coord = HitCoordinates_(*hit);
1641  if (i == 0 ||
1642  (std::abs((coord.Y() - (parm[0] + coord.X() * parm[1])) * cos(atan(parm[1]))) <
1643  fToler[i - 1]) ||
1644  fitok == 1) {
1645  ++nhits;
1646  if (nhits == fNfithits[i] + 1) break;
1647  wfit.push_back(coord.X());
1648  tfit.push_back(coord.Y());
1649  cfit.push_back(1.);
1650  if (i == fNfitpass - 1) { trackHits.push_back(hit); }
1651  }
1652  }
1653  }
1654 
1655  if (i < fNfitpass - 1 && wfit.size()) {
1656  fitok = WeightedFit(wfit.size(), &wfit[0], &tfit[0], &cfit[0], &parm[0]);
1657  }
1658  wfit.clear();
1659  tfit.clear();
1660  cfit.clear();
1661  }
1662 }
1663 
1665 {
1666  return TVector2(GlobalWire_(hit.WireID()), hit.PeakTime());
1667 }
1668 
1670  recob::Hit const& hit) const
1671 {
1672  geo::PlaneID planeID = hit.WireID().planeID();
1673  return HitPosition_(detProp, HitCoordinates_(hit), planeID);
1674 }
1675 
1677  TVector2 const& pos,
1678  geo::PlaneID planeID) const
1679 {
1680  return TVector2(pos.X() * fGeom->WirePitch(planeID), detProp.ConvertTicksToX(pos.Y(), planeID));
1681 }
1682 
1684 {
1685  double globalWire = -999;
1686 
1687  // Induction
1688  if (fGeom->SignalType(wireID) == geo::kInduction) {
1689  auto const wireCenter = fGeom->WireIDToWireGeo(wireID).GetCenter();
1690  globalWire = fGeom->WireCoordinate(wireCenter,
1691  geo::PlaneID{wireID.Cryostat,
1692  wireID.TPC % 2, // 0 or 1
1693  wireID.Plane});
1694  }
1695 
1696  // Collection
1697  else {
1698  // FOR COLLECTION WIRES, HARD CODE THE GEOMETRY FOR GIVEN DETECTORS
1699  // THIS _SHOULD_ BE TEMPORARY. GLOBAL WIRE SUPPORT IS BEING ADDED TO THE LARSOFT GEOMETRY AND SHOULD BE AVAILABLE SOON
1700  if (fDetector == "dune35t") {
1701  unsigned int nwires = fGeom->Nwires(geo::PlaneID{wireID.Cryostat, 0, wireID.Plane});
1702  if (wireID.TPC == 0 or wireID.TPC == 1)
1703  globalWire = wireID.Wire;
1704  else if (wireID.TPC == 2 or wireID.TPC == 3 or wireID.TPC == 4 or wireID.TPC == 5)
1705  globalWire = nwires + wireID.Wire;
1706  else if (wireID.TPC == 6 or wireID.TPC == 7)
1707  globalWire = (2 * nwires) + wireID.Wire;
1708  else
1709  mf::LogError("BlurredClusterAlg")
1710  << "Error when trying to find a global induction plane coordinate for TPC " << wireID.TPC
1711  << " (geometry" << fDetector << ")";
1712  }
1713  else if (fDetector == "dune10kt") {
1714  unsigned int nwires = fGeom->Nwires(geo::PlaneID{wireID.Cryostat, 0, wireID.Plane});
1715  // Detector geometry has four TPCs, two on top of each other, repeated along z...
1716  int block = wireID.TPC / 4;
1717  globalWire = (nwires * block) + wireID.Wire;
1718  }
1719  else {
1720  auto const wireCenter = fGeom->WireIDToWireGeo(wireID).GetCenter();
1721  globalWire = fGeom->WireCoordinate(wireCenter,
1722  geo::PlaneID{wireID.Cryostat,
1723  wireID.TPC % 2, // 0 or 1
1724  wireID.Plane});
1725  }
1726  }
1727 
1728  return globalWire;
1729 }
1730 
1732  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
1733 {
1734 
1735  // Get the wire widths
1736  std::map<int, int> planeWireLength;
1737  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitsIt =
1738  showerHitsMap.begin();
1739  showerHitsIt != showerHitsMap.end();
1740  ++showerHitsIt)
1741  planeWireLength[showerHitsIt->first] =
1742  std::abs(HitCoordinates_(*showerHitsIt->second.front()).X() -
1743  HitCoordinates_(*showerHitsIt->second.back()).X());
1744 
1745  // Find the relative wire width for each plane with respect to the others
1746  std::map<int, double> planeOtherWireLengths;
1747  for (std::map<int, int>::iterator planeWireLengthIt = planeWireLength.begin();
1748  planeWireLengthIt != planeWireLength.end();
1749  ++planeWireLengthIt) {
1750  double quad = 0.;
1751  for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane)
1752  if (plane != planeWireLengthIt->first and planeWireLength.count(plane))
1753  quad += cet::square(planeWireLength[plane]);
1754  quad = std::sqrt(quad);
1755  planeOtherWireLengths[planeWireLengthIt->first] =
1756  planeWireLength[planeWireLengthIt->first] / (double)quad;
1757  }
1758 
1759  // Order these backwards
1760  std::map<double, int> wireWidthMap;
1761  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitsIt =
1762  showerHitsMap.begin();
1763  showerHitsIt != showerHitsMap.end();
1764  ++showerHitsIt)
1765  wireWidthMap[planeOtherWireLengths.at(showerHitsIt->first)] = showerHitsIt->first;
1766 
1767  return wireWidthMap;
1768 }
1769 
1771  detinfo::DetectorPropertiesData const& detProp,
1772  const std::vector<art::Ptr<recob::Hit>>& showerHits) const
1773 {
1774 
1775  TVector2 pos;
1776  double weight = 1;
1777  double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1778  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hit = showerHits.begin();
1779  hit != showerHits.end();
1780  ++hit) {
1781  //++nhits;
1782  pos = HitPosition_(detProp, **hit);
1783  weight = cet::square((*hit)->Integral());
1784  sumweight += weight;
1785  sumx += weight * pos.X();
1786  sumy += weight * pos.Y();
1787  sumx2 += weight * pos.X() * pos.X();
1788  sumxy += weight * pos.X() * pos.Y();
1789  }
1790  double gradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1791  TVector2 direction = TVector2(1, gradient).Unit();
1792 
1793  return direction;
1794 }
1795 
1797  detinfo::DetectorPropertiesData const& detProp,
1798  std::vector<art::Ptr<recob::Hit>> const& showerHits) const
1799 {
1800 
1801  TVector2 pos, chargePoint = TVector2(0, 0);
1802  double totalCharge = 0;
1803  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hit = showerHits.begin();
1804  hit != showerHits.end();
1805  ++hit) {
1806  pos = HitPosition_(detProp, **hit);
1807  chargePoint += (*hit)->Integral() * pos;
1808  totalCharge += (*hit)->Integral();
1809  }
1810  TVector2 centre = chargePoint / totalCharge;
1811 
1812  return centre;
1813 }
1814 
1816  const std::vector<art::Ptr<recob::Hit>>& showerHits) const
1817 {
1818 
1819  TVector2 direction = ShowerDirection_(detProp, showerHits);
1820  TVector2 centre = ShowerCenter_(detProp, showerHits);
1821 
1822  std::vector<double> distanceToAxis;
1823  for (std::vector<art::Ptr<recob::Hit>>::const_iterator showerHitsIt = showerHits.begin();
1824  showerHitsIt != showerHits.end();
1825  ++showerHitsIt) {
1826  TVector2 proj = (HitPosition_(detProp, **showerHitsIt) - centre).Proj(direction) + centre;
1827  distanceToAxis.push_back((HitPosition_(detProp, **showerHitsIt) - proj).Mod());
1828  }
1829  double RMS = TMath::RMS(distanceToAxis.begin(), distanceToAxis.end());
1830 
1831  return RMS;
1832 }
1833 
1835  detinfo::DetectorPropertiesData const& detProp,
1836  const std::vector<art::Ptr<recob::Hit>>& showerHits) const
1837 {
1838  // Find a rough shower 'direction' and center
1839  TVector2 direction = ShowerDirection_(detProp, showerHits);
1840 
1841  // Bin the hits into discreet chunks
1842  int nShowerSegments = fNumShowerSegments;
1843  double lengthOfShower =
1844  (HitPosition_(detProp, *showerHits.back()) - HitPosition_(detProp, *showerHits.front())).Mod();
1845  double lengthOfSegment = lengthOfShower / (double)nShowerSegments;
1846  std::map<int, std::vector<art::Ptr<recob::Hit>>> showerSegments;
1847  std::map<int, double> segmentCharge;
1848  for (auto const& hitPtr : showerHits) {
1849  auto const& hit = *hitPtr;
1850  int const segment =
1851  (HitPosition_(detProp, hit) - HitPosition_(detProp, *showerHits.front())).Mod() /
1852  lengthOfSegment;
1853  showerSegments[segment].push_back(hitPtr);
1854  segmentCharge[segment] += hit.Integral();
1855  }
1856 
1857  TGraph* graph = new TGraph();
1858  std::vector<std::pair<int, double>> binVsRMS;
1859 
1860  // Loop over the bins to find the distribution of hits as the shower
1861  // progresses
1862  for (auto const& [segment, hitPtrs] : showerSegments) {
1863 
1864  // Get the mean position of the hits in this bin
1865  TVector2 meanPosition(0, 0);
1866  for (auto const& hit : hitPtrs | views::transform(to_element))
1867  meanPosition += HitPosition_(detProp, hit);
1868  meanPosition /= (double)hitPtrs.size();
1869 
1870  // Get the RMS of this bin
1871  std::vector<double> distanceToAxisBin;
1872  for (auto const& hit : hitPtrs | views::transform(to_element)) {
1873  TVector2 proj = (HitPosition_(detProp, hit) - meanPosition).Proj(direction) + meanPosition;
1874  distanceToAxisBin.push_back((HitPosition_(detProp, hit) - proj).Mod());
1875  }
1876 
1877  double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
1878  binVsRMS.emplace_back(segment, RMSBin);
1879  if (fMakeRMSGradientPlot) graph->SetPoint(graph->GetN(), segment, RMSBin);
1880  }
1881 
1882  // Get the gradient of the RMS-bin plot
1883  double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1884  for (auto const& [bin, RMSBin] : binVsRMS) {
1885  double weight = segmentCharge.at(bin);
1886  sumweight += weight;
1887  sumx += weight * bin;
1888  sumy += weight * RMSBin;
1889  sumx2 += weight * bin * bin;
1890  sumxy += weight * bin * RMSBin;
1891  }
1892  double RMSgradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1893 
1894  if (fMakeRMSGradientPlot) {
1895  TVector2 direction = TVector2(1, RMSgradient).Unit();
1896  TCanvas* canv = new TCanvas();
1897  graph->Draw();
1898  graph->GetXaxis()->SetTitle("Shower segment");
1899  graph->GetYaxis()->SetTitle("RMS of hit distribution");
1900  TVector2 centre = TVector2(graph->GetMean(1), graph->GetMean(2));
1901  TLine line;
1902  line.SetLineColor(2);
1903  line.DrawLine(centre.X() - 1000 * direction.X(),
1904  centre.Y() - 1000 * direction.Y(),
1905  centre.X() + 1000 * direction.X(),
1906  centre.Y() + 1000 * direction.Y());
1907  canv->SaveAs("RMSGradient.png");
1908  delete canv;
1909  }
1910  delete graph;
1911 
1912  return RMSgradient;
1913 }
1914 
1916  detinfo::DetectorPropertiesData const& detProp,
1917  geo::Point_t const& point,
1918  int plane,
1919  int cryostat) const
1920 {
1921 
1922  TVector2 wireTickPos = TVector2(-999., -999.);
1923 
1924  geo::TPCID tpcID = fGeom->FindTPCAtPosition(point);
1925  int tpc = 0;
1926  if (tpcID.isValid)
1927  tpc = tpcID.TPC;
1928  else
1929  return wireTickPos;
1930 
1931  // Construct wire ID for this point projected onto the plane
1932  geo::PlaneID planeID = geo::PlaneID(cryostat, tpc, plane);
1933  geo::WireID wireID;
1934  try {
1935  wireID = fGeom->NearestWireID(point, planeID);
1936  }
1937  catch (geo::InvalidWireError const& e) {
1938  wireID = e.suggestedWireID(); // pick the closest valid wire
1939  }
1940 
1941  wireTickPos = TVector2(GlobalWire_(wireID), detProp.ConvertXToTicks(point.X(), planeID));
1942 
1943  return HitPosition_(detProp, wireTickPos, planeID);
1944 }
1945 
1947  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
1948 {
1949  // Get the width of the shower in wire coordinate
1950  std::map<int, int> planeWireLength;
1951  std::map<int, double> planeOtherWireLengths;
1952  for (auto const& [plane, hits] : showerHitsMap)
1953  planeWireLength[plane] =
1954  std::abs(HitCoordinates_(*hits.front()).X() - HitCoordinates_(*hits.back()).X());
1955  for (std::map<int, int>::iterator planeWireLengthIt = planeWireLength.begin();
1956  planeWireLengthIt != planeWireLength.end();
1957  ++planeWireLengthIt) {
1958  double quad = 0.;
1959  for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane)
1960  if (plane != planeWireLengthIt->first and planeWireLength.count(plane))
1961  quad += cet::square(planeWireLength[plane]);
1962  quad = std::sqrt(quad);
1963  planeOtherWireLengths[planeWireLengthIt->first] =
1964  planeWireLength[planeWireLengthIt->first] / (double)quad;
1965  }
1966 
1967  if (fDebug > 1)
1968  for (auto const [plane, relative_width] : planeOtherWireLengths) {
1969  std::cout << "Plane " << plane << " has " << planeWireLength[plane]
1970  << " wire width and therefore has relative width in wire of " << relative_width
1971  << '\n';
1972  }
1973 
1974  std::map<double, int> wireWidthMap;
1975  for (int const plane : showerHitsMap | views::keys) {
1976  double wireWidth = planeWireLength.at(plane);
1977  wireWidthMap[wireWidth] = plane;
1978  }
1979 
1980  return wireWidthMap.begin()->second;
1981 }
1982 
1984  const Double_t* x,
1985  const Double_t* y,
1986  const Double_t* w,
1987  Double_t* parm) const
1988 {
1989 
1990  Double_t sumx = 0.;
1991  Double_t sumx2 = 0.;
1992  Double_t sumy = 0.;
1993  Double_t sumxy = 0.;
1994  Double_t sumw = 0.;
1995  Double_t eparm[2];
1996 
1997  parm[0] = 0.;
1998  parm[1] = 0.;
1999  eparm[0] = 0.;
2000  eparm[1] = 0.;
2001 
2002  for (Int_t i = 0; i < n; i++) {
2003  sumx += x[i] * w[i];
2004  sumx2 += x[i] * x[i] * w[i];
2005  sumy += y[i] * w[i];
2006  sumxy += x[i] * y[i] * w[i];
2007  sumw += w[i];
2008  }
2009 
2010  if (sumx2 * sumw - sumx * sumx == 0.) return 1;
2011  if (sumx2 - sumx * sumx / sumw == 0.) return 1;
2012 
2013  parm[0] = (sumy * sumx2 - sumx * sumxy) / (sumx2 * sumw - sumx * sumx);
2014  parm[1] = (sumxy - sumx * sumy / sumw) / (sumx2 - sumx * sumx / sumw);
2015 
2016  eparm[0] = sumx2 * (sumx2 * sumw - sumx * sumx);
2017  eparm[1] = (sumx2 - sumx * sumx / sumw);
2018 
2019  if (eparm[0] < 0. || eparm[1] < 0.) return 1;
2020 
2021  eparm[0] = sqrt(eparm[0]) / (sumx2 * sumw - sumx * sumx);
2022  eparm[1] = sqrt(eparm[1]) / (sumx2 - sumx * sumx / sumw);
2023 
2024  return 0;
2025 }
2026 
2028 {
2029  if (hits.empty()) return false;
2030  if (hits.size() > 2000) return true;
2031  if (hits.size() < 20) return true;
2032 
2033  std::map<int, int> hitmap;
2034  unsigned nhits = 0;
2035  for (auto const& hit : hits | views::transform(to_element)) {
2036  ++nhits;
2037  if (nhits > 2) ++hitmap[hit.WireID().Wire];
2038  if (nhits == 20) break;
2039  }
2040  if (float(nhits - 2) / hitmap.size() > 1.4)
2041  return false;
2042  else
2043  return true;
2044 }
Float_t x
Definition: compare.C:6
void XYZ(double *xyz) const
Legacy method to access vertex position, preserved to avoid breaking code. Please try to use Vertex::...
Definition: Vertex.cxx:34
intermediate_table::iterator iterator
TVector2 HitPosition_(detinfo::DetectorPropertiesData const &detProp, recob::Hit const &hit) const
Return the coordinates of this hit in units of cm.
geo::Point_t Construct3DPoint_(detinfo::DetectorPropertiesData const &detProp, art::Ptr< recob::Hit > const &hit1, art::Ptr< recob::Hit > const &hit2) const
Constructs a 3D point (in [cm]) to represent the hits given in two views.
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
std::vector< recob::SpacePoint > MakeSpacePoints(detinfo::DetectorPropertiesData const &detProp, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &hits, std::vector< std::vector< art::Ptr< recob::Hit >>> &hitAssns) const
Makes space points from the shower hits in each plane.
void FindInitialTrackHits(std::vector< art::Ptr< recob::Hit >> const &showerHits, art::Ptr< recob::Vertex > const &vertex, std::vector< art::Ptr< recob::Hit >> &trackHits) const
<Tingjun to="" document>="">
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:221
int WorstPlane_(const std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const
Returns the plane which is determined to be the least likely to be correct.
double z
z position of intersection
Definition: geo_types.h:792
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
shower::ShowerEnergyAlg const fShowerEnergyAlg
Definition: EMShowerAlg.h:278
constexpr to_element_t to_element
Definition: ToElement.h:25
pma::ProjectionMatchingAlg const fProjectionMatchingAlg
Definition: EMShowerAlg.h:280
double ShowerEnergy(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &hits, geo::PlaneID::PlaneID_t plane) const
recob::Shower MakeShower(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &hits, std::unique_ptr< recob::Track > const &initialTrack, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &initialTrackHits) const
Makes a recob::Shower object given the hits in the shower and the initial track-like part...
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:78
auto coord(Vector &v, unsigned int n) noexcept
Returns an object to manage the coordinate n of a vector.
Float_t y
Definition: compare.C:6
Float_t Y
Definition: plot.C:37
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:210
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:67
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
constexpr auto abs(T v)
Returns the absolute value of the argument.
Vector_t VertexDirection() const
Access to track direction at different points.
Definition: Track.h:166
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
void AssociateClustersAndTracks(std::vector< art::Ptr< recob::Cluster >> const &clusters, art::FindManyP< recob::Hit > const &fmh, art::FindManyP< recob::Track > const &fmt, std::map< int, std::vector< int >> &clusterToTracks, std::map< int, std::vector< int >> &trackToClusters) const
Map associated tracks and clusters together given their associated hits.
Definition: EMShowerAlg.cxx:59
std::map< int, std::map< int, bool > > GetPlanePermutations_(const detinfo::DetectorPropertiesData &detProp, const std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
double GlobalWire_(const geo::WireID &wireID) const
Find the global wire position.
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:276
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
intermediate_table::const_iterator const_iterator
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Definition: Cluster.h:717
Cluster finding and building.
double const fdEdxTrackLength
Definition: EMShowerAlg.h:266
unsigned int const fNfitpass
Definition: EMShowerAlg.h:270
TVector2 ShowerCenter_(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &showerHits) const
Returns the charge-weighted shower center.
std::vector< std::vector< int > > FindShowers(std::map< int, std::vector< int >> const &trackToClusters) const
Makes showers given a map between tracks and all clusters associated with them.
std::vector< double > const fToler
Definition: EMShowerAlg.h:272
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
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
bool const fMakeRMSGradientPlot
Definition: EMShowerAlg.h:286
void CheckIsolatedHits_(std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:166
std::vector< unsigned int > const fNfithits
Definition: EMShowerAlg.h:271
Int_t WeightedFit(const Int_t n, const Double_t *x, const Double_t *y, const Double_t *w, Double_t *parm) const
<Tingjun to="" document>="">
Collection of exceptions for Geometry system.
void hits()
Definition: readHits.C:15
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
if(nlines<=0)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Signal from induction planes.
Definition: geo_types.h:151
void FindInitialTrack(detinfo::DetectorPropertiesData const &detProp, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &hits, std::unique_ptr< recob::Track > &initialTrack, std::map< int, std::vector< art::Ptr< recob::Hit >>> &initialTrackHits) const
DebugStuff debug
Definition: DebugStruct.cxx:4
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
A trajectory in space reconstructed from hits.
decltype(auto) values(Coll &&coll)
Range-for loop helper iterating across the values of the specified collection.
std::vector< int > CheckShowerPlanes(std::vector< std::vector< int >> const &initialShowers, std::vector< art::Ptr< recob::Cluster >> const &clusters, art::FindManyP< recob::Hit > const &fmh) const
Takes the initial showers found and tries to resolve issues where one bad view ruins the event...
key_type key() const noexcept
Definition: Ptr.h:166
TVector2 Project3DPointOntoPlane_(detinfo::DetectorPropertiesData const &detProp, geo::Point_t const &point, int plane, int cryostat=0) const
Point_t const & Vertex() const
Access to track position at different points.
Definition: Track.h:158
std::unique_ptr< recob::Track > ConstructTrack(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &track1, std::vector< art::Ptr< recob::Hit >> const &track2, std::map< int, TVector2 > const &showerCentreMap) const
bool isNull() const noexcept
Definition: Ptr.h:211
std::map< int, std::vector< art::Ptr< recob::Hit > > > OrderShowerHits(detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &shower, int plane) const
Takes the hits associated with a shower and orders them so they follow the direction of the shower...
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:73
double ConvertXToTicks(double X, int p, int t, int c) const
double ShowerHitRMS_(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &showerHits) const
double const fMinTrackLength
Definition: EMShowerAlg.h:265
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
float bin[41]
Definition: plottest35.C:14
bool isCleanShower(std::vector< art::Ptr< recob::Hit >> const &hits) const
<Tingjun to="" document>="">
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
bool const fMakeGradientPlot
Definition: EMShowerAlg.h:285
double value
Definition: spectrum.C:18
EMShowerAlg(fhicl::ParameterSet const &pset, int debug=0)
Definition: EMShowerAlg.cxx:37
TVector2 ShowerDirection_(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &showerHits) const
Detector simulation of raw signals on wires.
bool CheckShowerHits_(detinfo::DetectorPropertiesData const &detProp, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &showerHitsMap) const
double ConvertTicksToX(double ticks, int p, int t, int c) const
void OrderShowerHits_(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &shower, std::vector< art::Ptr< recob::Hit >> &orderedShower, art::Ptr< recob::Vertex > const &vertex) const
std::unique_ptr< recob::Track > MakeInitialTrack_(detinfo::DetectorPropertiesData const &detProp, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &initialHitsMap, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &showerHitsMap) const
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::size_t color(std::string const &procname)
double FinddEdx_(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &trackHits, std::unique_ptr< recob::Track > const &track) const
Finds dE/dx for the track given a set of hits.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
double weight
Definition: plottest35.C:25
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
WireID NearestWireID(Point_t const &point, PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
Contains all timing reference information for the detector.
double y
y position of intersection
Definition: geo_types.h:791
TVector2 HitCoordinates_(recob::Hit const &hit) const
Return the coordinates of this hit in global wire/tick space.
Float_t proj
Definition: plot.C:35
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double TrackPitchInView(recob::Track const &track, geo::View_t view, size_t trajectory_point=0U)
Returns the projected length of track on a wire pitch step [cm].
Definition: TrackUtils.cxx:75
int const fNumShowerSegments
Definition: EMShowerAlg.h:287
calo::CalorimetryAlg const fCalorimetryAlg
Definition: EMShowerAlg.h:279
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
std::vector< art::Ptr< recob::Hit > > FindOrderOfHits_(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &hits, bool perpendicular=false) const
Exception thrown on invalid wire number.
Definition: Exceptions.h:39
double const fSpacePointSize
Definition: EMShowerAlg.h:267
Char_t n[5]
IDparameter< geo::TPCID > TPCID
Member type of validated geo::TPCID parameter.
std::map< double, int > RelativeWireWidth_(const std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
constexpr PlaneID const & planeID() const
Definition: geo_types.h:620
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
Double_t sum
Definition: plot.C:31
art::ServiceHandle< geo::Geometry const > fGeom
Definition: EMShowerAlg.h:275
Float_t e
Definition: plot.C:35
size_t size() const
Definition: PmaTrack3D.h:89
geo::WireID suggestedWireID() const
Returns a better wire ID.
Definition: Exceptions.h:87
recob::tracking::Plane Plane
Definition: TrackState.h:17
Float_t X
Definition: plot.C:37
Float_t track
Definition: plot.C:35
Float_t w
Definition: plot.C:20
Utility functions to extract information from recob::Track
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
std::string const fDetector
Definition: EMShowerAlg.h:282
const Point_t & position() const
Return vertex 3D position.
Definition: Vertex.h:64
std::map< int, std::vector< art::Ptr< recob::Hit > > > FindShowerStart_(std::map< int, std::vector< art::Ptr< recob::Hit >>> const &orderedShowerMap) const
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
double ShowerHitRMSGradient_(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &showerHits) const
Returns the gradient of the RMS vs shower segment graph.
vertex reconstruction