LArSoft  v10_04_05
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 #include <limits>
35 
36 using lar::to_element;
37 using namespace ranges;
38 
40  : fDebug{debug}
41  , fMinTrackLength{pset.get<double>("MinTrackLength")}
42  , fdEdxTrackLength{pset.get<double>("dEdxTrackLength")}
43  , fSpacePointSize{pset.get<double>("SpacePointSize")}
44  , fNfitpass{pset.get<unsigned int>("Nfitpass")}
45  , fNfithits{pset.get<std::vector<unsigned int>>("Nfithits")}
46  , fToler{pset.get<std::vector<double>>("Toler")}
47  , fShowerEnergyAlg(pset.get<fhicl::ParameterSet>("ShowerEnergyAlg"))
48  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
49  , fProjectionMatchingAlg(pset.get<fhicl::ParameterSet>("ProjectionMatchingAlg"))
50  , fDetector{pset.get<std::string>("Detector", "dune35t")} // tmp
51  , fMakeGradientPlot{pset.get<bool>("MakeGradientPlot", false)}
52  , fMakeRMSGradientPlot{pset.get<bool>("MakeRMSGradientPlot", false)}
53  , fNumShowerSegments{pset.get<int>("NumShowerSegments", 4)}
54 {
55  if (fNfitpass != fNfithits.size() || fNfitpass != fToler.size()) {
57  << "EMShowerAlg: fNfithits and fToler need to have size fNfitpass";
58  }
59 }
60 
62  std::vector<art::Ptr<recob::Cluster>> const& clusters,
63  art::FindManyP<recob::Hit> const& fmh,
65  std::map<int, std::vector<int>>& clusterToTracks,
66  std::map<int, std::vector<int>>& trackToClusters) const
67 {
68  std::vector<int> clustersToIgnore = {-999};
70  clusters, fmh, fmt, clustersToIgnore, clusterToTracks, trackToClusters);
71 }
72 
74  std::vector<art::Ptr<recob::Cluster>> const& clusters,
75  art::FindManyP<recob::Hit> const& fmh,
77  std::vector<int> const& clustersToIgnore,
78  std::map<int, std::vector<int>>& clusterToTracks,
79  std::map<int, std::vector<int>>& trackToClusters) const
80 {
81  // Look through all the clusters
82  for (auto const& clusterPtr : clusters) {
83 
84  // Get the hits in this cluster
85  auto const& clusterHits = fmh.at(clusterPtr.key());
86 
87  // Look at all these hits and find the associated tracks
88  for (auto const& clusterHitPtr : clusterHits) {
89  // Get the tracks associated with this hit
90  auto const& clusterHitTracks = fmt.at(clusterHitPtr.key());
91  if (clusterHitTracks.size() > 1) {
92  std::cout << "More than one track associated with this hit!\n";
93  continue;
94  }
95 
96  if (clusterHitTracks.size() < 1) continue;
97 
98  auto const& clusterHitTrackPtr = clusterHitTracks[0];
99  if (clusterHitTrackPtr->Length() < fMinTrackLength) {
100  if (fDebug > 2)
101  std::cout << "Track " << clusterHitTrackPtr->ID() << " is too short! ("
102  << clusterHitTrackPtr->Length() << ")\n";
103  continue;
104  }
105 
106  // Add this cluster to the track map
107  int track = clusterHitTrackPtr.key();
108  int cluster = clusterPtr.key();
109  if (cet::search_all(clustersToIgnore, cluster)) continue;
110  if (not cet::search_all(trackToClusters[track], cluster))
111  trackToClusters[track].push_back(cluster);
112  if (not cet::search_all(clusterToTracks[cluster], track))
113  clusterToTracks[cluster].push_back(track);
114  }
115  }
116 }
117 
119  std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
120 {
121  std::map<int, std::vector<int>> firstTPC;
122  for (auto const& [plane, hits] : showerHitsMap)
123  firstTPC[hits.at(0)->WireID().TPC].push_back(plane);
124 
125  // If all in the same TPC then that's great!
126  if (firstTPC.size() != 2) return;
127 
128  // If they are in more than two TPCs, not much we can do
129  else if (firstTPC.size() > 2)
130  return;
131 
132  // If we get to this point, there should be something we can do!
133 
134  // Find the problem plane
135  int problemPlane = -1;
136  for (auto const& planes : firstTPC | views::values)
137  if (planes.size() == 1) problemPlane = planes[0];
138 
139  // Require three hits
140  if (showerHitsMap.at(problemPlane).size() < 3) return;
141 
142  // and get the other planes with at least three hits
143  std::vector<int> otherPlanes;
144  for (int plane = 0; plane < (int)fWireReadoutGeom->MaxPlanes(); ++plane)
145  if (plane != problemPlane and showerHitsMap.count(plane) and
146  showerHitsMap.at(plane).size() >= 3)
147  otherPlanes.push_back(plane);
148 
149  if (otherPlanes.size() == 0) return;
150 
151  // Look at the hits after the first one
152  unsigned int wrongTPC = showerHitsMap.at(problemPlane).at(0)->WireID().TPC;
153  unsigned int nHits = 0;
154  for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsMap.at(problemPlane).begin();
155  hitIt != showerHitsMap.at(problemPlane).end() and (*hitIt)->WireID().TPC == wrongTPC;
156  ++hitIt)
157  ++nHits;
158 
159  // If there are more than two hits in the 'wrong TPC', we can't be sure it is indeed wrong
160  if (nHits > 2) return;
161 
162  // See if at least the next four times as many hits are in a different TPC
163  std::map<int, int> otherTPCs;
165  std::next(showerHitsMap.at(problemPlane).begin(), nHits);
166  hitIt != showerHitsMap.at(problemPlane).end() and
167  std::distance(std::next(showerHitsMap.at(problemPlane).begin(), nHits), hitIt) < 4 * nHits;
168  ++hitIt)
169  ++otherTPCs[(*hitIt)->WireID().TPC];
170 
171  if (otherTPCs.size() > 1) return;
172 
173  // If we get this far, we can move the problem hits from the front of the shower to the back
174  std::map<int, int> tpcCount;
175  for (int const otherPlane : otherPlanes)
177  std::next(showerHitsMap.at(otherPlane).begin());
178  hitIt != showerHitsMap.at(otherPlane).end() and
179  hitIt != std::next(showerHitsMap.at(otherPlane).begin(), 2);
180  ++hitIt)
181  ++tpcCount[(*hitIt)->WireID().TPC];
182 
183  // Remove the first hit if it is in the wrong TPC
184  if (tpcCount.size() == 1 and
185  tpcCount.begin()->first ==
186  (int)(*std::next(showerHitsMap.at(problemPlane).begin(), nHits))->WireID().TPC) {
187  std::vector<art::Ptr<recob::Hit>> naughty_hits;
188  for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsMap.at(problemPlane).begin();
189  hitIt != std::next(showerHitsMap.at(problemPlane).begin(), nHits);
190  ++hitIt) {
191  naughty_hits.push_back(*hitIt);
192  showerHitsMap.at(problemPlane).erase(hitIt);
193  }
194  for (auto const& naughty_hit : naughty_hits)
195  showerHitsMap.at(problemPlane).push_back(naughty_hit);
196  }
197 }
198 
200  detinfo::DetectorPropertiesData const& detProp,
201  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& showerHitsMap) const
202 {
203  bool consistencyCheck = true;
204 
205  if (showerHitsMap.size() < 2) { consistencyCheck = true; }
206  else if (showerHitsMap.size() == 2) {
207 
208  // With two views, we can check:
209  // -- timing between views is consistent
210  // -- the 3D start point makes sense when projected back onto the individual planes
211 
212  std::vector<art::Ptr<recob::Hit>> startHits;
213  std::vector<int> planes;
214  for (auto const& [plane, hits] : showerHitsMap) {
215  startHits.push_back(hits.front());
216  planes.push_back(plane);
217  }
218 
219  auto const showerStartPos = Construct3DPoint_(detProp, startHits.at(0), startHits.at(1));
220  TVector2 proj1 = Project3DPointOntoPlane_(detProp, showerStartPos, planes.at(0));
221  TVector2 proj2 = Project3DPointOntoPlane_(detProp, showerStartPos, planes.at(1));
222 
223  double timingDifference = std::abs(startHits.at(0)->PeakTime() - startHits.at(1)->PeakTime());
224  double projectionDifference = ((HitPosition_(detProp, *startHits.at(0)) - proj1).Mod() +
225  (HitPosition_(detProp, *startHits.at(1)) - proj2).Mod()) /
226  (double)2;
227 
228  if (timingDifference > 40 or projectionDifference > 5 or showerStartPos.X() == -9999 or
229  showerStartPos.Y() == -9999 or showerStartPos.Z() == -9999)
230  consistencyCheck = false;
231 
232  if (fDebug > 1)
233  std::cout << "Timing difference is " << timingDifference << " and projection distance is "
234  << projectionDifference << " (start is (" << showerStartPos.X() << ", "
235  << showerStartPos.Y() << ", " << showerStartPos.Z() << ")\n";
236  }
237  else if (showerHitsMap.size() == 3) {
238 
239  // With three views, we can check:
240  // -- the timing between views is consistent
241  // -- the 3D start point formed by two views and projected back into the third is close to the start point in that view
242 
243  std::map<int, art::Ptr<recob::Hit>> start2DMap;
244  for (auto const& [plane, hits] : showerHitsMap) {
245  start2DMap[plane] = hits.front();
246  }
247 
248  std::map<int, double> projDiff;
249  std::map<int, double> timingDiff;
250 
251  for (int plane = 0; plane < 3; ++plane) {
252 
253  std::vector<int> otherPlanes;
254  for (int otherPlane = 0; otherPlane < 3; ++otherPlane)
255  if (otherPlane != plane) otherPlanes.push_back(otherPlane);
256 
257  auto const showerStartPos = Construct3DPoint_(
258  detProp, start2DMap.at(otherPlanes.at(0)), start2DMap.at(otherPlanes.at(1)));
259  TVector2 showerStartProj = Project3DPointOntoPlane_(detProp, showerStartPos, plane);
260 
261  if (fDebug > 2) {
262  std::cout << "Plane... " << plane;
263  std::cout << "\nStart position in this plane is "
264  << HitPosition_(detProp, *start2DMap.at(plane)).X() << ", "
265  << HitPosition_(detProp, *start2DMap.at(plane)).Y() << ")\n";
266  std::cout << "Shower start from other two planes is (" << showerStartPos.X() << ", "
267  << showerStartPos.Y() << ", " << showerStartPos.Z() << ")\n";
268  std::cout << "Projecting the other two planes gives position (" << showerStartProj.X()
269  << ", " << showerStartProj.Y() << ")\n";
270  }
271 
272  double projDiff =
273  std::abs((showerStartProj - HitPosition_(detProp, *start2DMap.at(plane))).Mod());
274  double timeDiff = TMath::Max(
275  std::abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(0))->PeakTime()),
276  std::abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(1))->PeakTime()));
277 
278  if (fDebug > 1)
279  std::cout << "Plane " << plane << " has projDiff " << projDiff << " and timeDiff "
280  << timeDiff << '\n';
281 
282  if (projDiff > 5 or timeDiff > 40) consistencyCheck = false;
283  }
284  }
285 
286  if (fDebug > 1) std::cout << "Consistency check is " << consistencyCheck << '\n';
287 
288  return consistencyCheck;
289 }
290 
292  std::vector<std::vector<int>> const& initialShowers,
293  std::vector<art::Ptr<recob::Cluster>> const& clusters,
294  art::FindManyP<recob::Hit> const& fmh) const
295 {
296  std::vector<int> clustersToIgnore;
297 
298  // Look at each shower
299  for (auto initialShowerIt = initialShowers.cbegin(); initialShowerIt != initialShowers.cend();
300  ++initialShowerIt) {
301 
302  if (std::distance(initialShowers.cbegin(), initialShowerIt) > 0) continue;
303 
304  // Map the clusters and cluster hits to each view
305  std::map<int, std::vector<art::Ptr<recob::Cluster>>> planeClusters;
306  std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHits;
307  for (int const clusterIndex : *initialShowerIt) {
308  art::Ptr<recob::Cluster> const& cluster = clusters.at(clusterIndex);
309  planeClusters[cluster->Plane().Plane].push_back(cluster);
310  for (auto const& hit : fmh.at(cluster.key()))
311  planeHits[hit->WireID().Plane].push_back(hit);
312  }
313 
314  TFile* outFile = new TFile("chargeDistributions.root", "RECREATE");
315  std::map<int, TH1D*> chargeDist;
316  for (auto const& [plane, clusterPtrs] : planeClusters) {
317  for (auto const& clusterPtr : clusterPtrs) {
318  chargeDist[plane] = new TH1D(std::string("chargeDist_Plane" + std::to_string(plane) +
319  "_Cluster" + std::to_string(clusterPtr.key()))
320  .c_str(),
321  "",
322  150,
323  0,
324  1000);
325  auto const& hits = fmh.at(clusterPtr.key());
326  for (auto const& hit : hits | views::transform(to_element)) {
327  chargeDist[plane]->Fill(hit.Integral());
328  }
329  outFile->cd();
330  chargeDist[plane]->Write();
331  }
332  }
333  outFile->Close();
334  delete outFile;
335 
336  // Can't do much with fewer than three views
337  if (planeClusters.size() < 3) continue;
338 
339  // Look at how many clusters each plane has, and the proportion of hits each one uses
340  std::map<int, std::vector<double>> planeClusterSizes;
341  for (std::map<int, std::vector<art::Ptr<recob::Cluster>>>::iterator planeClustersIt =
342  planeClusters.begin();
343  planeClustersIt != planeClusters.end();
344  ++planeClustersIt) {
345  for (std::vector<art::Ptr<recob::Cluster>>::iterator planeClusterIt =
346  planeClustersIt->second.begin();
347  planeClusterIt != planeClustersIt->second.end();
348  ++planeClusterIt) {
349  std::vector<art::Ptr<recob::Hit>> hits = fmh.at(planeClusterIt->key());
350  planeClusterSizes[planeClustersIt->first].push_back(
351  (double)hits.size() / (double)planeHits.at(planeClustersIt->first).size());
352  }
353  }
354 
355  // Find the average hit fraction across all clusters in the plane
356  std::map<int, double> planeClustersAvSizes;
357  for (auto const& [plane, cluster_sizes] : planeClusterSizes) {
358  double const average = accumulate(cluster_sizes, 0.) / cluster_sizes.size();
359  planeClustersAvSizes[plane] = average;
360  }
361 
362  // Now decide if there is one plane which is ruining the reconstruction
363  // If two planes have a low average cluster fraction and one high, this plane likely merges two particle deposits together
364  int badPlane = -1;
365  for (auto const [plane, avg] : planeClustersAvSizes) {
366 
367  // Get averages from other planes and add in quadrature
368  std::vector<double> otherAverages;
369  for (auto const [other_plane, other_avg] : planeClustersAvSizes)
370  if (plane != other_plane) otherAverages.push_back(other_avg);
371 
372  double const sumSquareAvgsOtherPlanes = accumulate(
373  otherAverages, 0., [](double sum, double avg) { return sum + cet::square(avg); });
374  double const quadOtherPlanes = std::sqrt(sumSquareAvgsOtherPlanes);
375 
376  // If this plane has an average higher than the quadratic sum of the
377  // others, it may be bad
378  if (avg >= quadOtherPlanes) badPlane = plane;
379  }
380 
381  if (badPlane != -1) {
382  if (fDebug > 1) std::cout << "Bad plane is " << badPlane << '\n';
383  for (auto const& cluster : planeClusters.at(badPlane))
384  clustersToIgnore.push_back(cluster.key());
385  }
386  }
387 
388  return clustersToIgnore;
389 }
390 
392  art::Ptr<recob::Hit> const& hit1,
393  art::Ptr<recob::Hit> const& hit2) const
394 {
395  // x is average of the two x's
396  double x = (detProp.ConvertTicksToX(hit1->PeakTime(), hit1->WireID().planeID()) +
397  detProp.ConvertTicksToX(hit2->PeakTime(), hit2->WireID().planeID())) /
398  (double)2;
399 
400  // y and z got from the wire interections
401  auto intersection = fWireReadoutGeom->WireIDsIntersect(hit1->WireID(), hit2->WireID())
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  std::unique_ptr<recob::Track> track;
413 
414  std::vector<art::Ptr<recob::Hit>> track1, track2;
415 
416  // Check the TPCs
417  if ((*hits1.begin())->WireID().TPC != (*hits2.begin())->WireID().TPC) {
418  mf::LogWarning("EMShowerAlg") << "Warning: attempting to construct a track from two different "
419  "TPCs. Returning a null track.";
420  return track;
421  }
422 
423  // Check for tracks crossing TPC boundaries
424  std::map<int, int> tpcMap;
425  for (auto const& hit : hits1)
426  ++tpcMap[hit->WireID().TPC];
427  for (auto const& hit : hits2)
428  ++tpcMap[hit->WireID().TPC];
429  if (tpcMap.size() > 1) {
430  mf::LogWarning("EMShowerAlg")
431  << "Warning: attempting to construct a track which crosses more than one TPC -- PMTrack "
432  "can't handle this right now. Returning a track made just from hits in the first TPC it "
433  "traverses.";
434  unsigned int firstTPC1 = hits1.at(0)->WireID().TPC, firstTPC2 = hits2.at(0)->WireID().TPC;
435  for (auto const& hit : hits1)
436  if (hit->WireID().TPC == firstTPC1) track1.push_back(hit);
437  for (auto const& hit : hits2)
438  if (hit->WireID().TPC == firstTPC2) track2.push_back(hit);
439  }
440  else {
441  track1 = hits1;
442  track2 = hits2;
443  }
444 
445  if (fDebug > 1) {
446  std::cout << "About to make a track from these hits:\n";
447  auto print_hits = [this](auto const& track) {
448  for (auto const& hit : track | views::transform(to_element)) {
449  std::cout << "Hit (" << HitCoordinates_(hit).X() << ", " << HitCoordinates_(hit).Y()
450  << ") (real wire " << hit.WireID().Wire << ") in TPC " << hit.WireID().TPC
451  << '\n';
452  }
453  };
454  print_hits(track1);
455  print_hits(track2);
456  }
457 
458  auto const trackStart = Construct3DPoint_(detProp, track1.at(0), track2.at(0));
459  pma::Track3D* pmatrack = fProjectionMatchingAlg.buildSegment(detProp, track1, track2, trackStart);
460 
461  if (!pmatrack) {
462  mf::LogInfo("EMShowerAlg") << "Skipping this event because not enough hits in two views";
463  return track;
464  }
465 
466  std::vector<TVector3> xyz, dircos;
467 
468  for (unsigned int i = 0; i < pmatrack->size(); i++) {
469 
470  xyz.push_back((*pmatrack)[i]->Point3D());
471 
472  if (i < pmatrack->size() - 1) {
473  size_t j = i + 1;
474  double mag = 0.0;
475  TVector3 dc(0., 0., 0.);
476  while ((mag == 0.0) and (j < pmatrack->size())) {
477  dc = (*pmatrack)[j]->Point3D();
478  dc -= (*pmatrack)[i]->Point3D();
479  mag = dc.Mag();
480  ++j;
481  }
482  if (mag > 0.0)
483  dc *= 1.0 / mag;
484  else if (!dircos.empty())
485  dc = dircos.back();
486  dircos.push_back(dc);
487  }
488  else
489  dircos.push_back(dircos.back());
490  }
491 
492  // Orient the track correctly
493  std::map<int, double> distanceToVertex, distanceToEnd;
494  using geo::vect::toPoint;
495  geo::Point_t const vertex = toPoint(*xyz.begin());
496  geo::Point_t const end = toPoint(*xyz.rbegin());
497 
498  // Loop over all the planes and find the distance from the vertex and end
499  // projections to the centre in each plane
500  for (std::map<int, TVector2>::const_iterator showerCentreIt = showerCentreMap.begin();
501  showerCentreIt != showerCentreMap.end();
502  ++showerCentreIt) {
503 
504  // Project the vertex and the end point onto this plane
505  TVector2 vertexProj = Project3DPointOntoPlane_(detProp, vertex, showerCentreIt->first);
506  TVector2 endProj = Project3DPointOntoPlane_(detProp, end, showerCentreIt->first);
507 
508  // Find the distance of each to the centre of the cluster
509  distanceToVertex[showerCentreIt->first] = (vertexProj - showerCentreIt->second).Mod();
510  distanceToEnd[showerCentreIt->first] = (endProj - showerCentreIt->second).Mod();
511  }
512 
513  // Find the average distance to the vertex and the end across the planes
514  double avDistanceToVertex = 0, avDistanceToEnd = 0;
515  for (std::map<int, double>::iterator distanceToVertexIt = distanceToVertex.begin();
516  distanceToVertexIt != distanceToVertex.end();
517  ++distanceToVertexIt)
518  avDistanceToVertex += distanceToVertexIt->second;
519  avDistanceToVertex /= distanceToVertex.size();
520 
521  for (std::map<int, double>::iterator distanceToEndIt = distanceToEnd.begin();
522  distanceToEndIt != distanceToEnd.end();
523  ++distanceToEndIt)
524  avDistanceToEnd += distanceToEndIt->second;
525  if (distanceToEnd.size() != 0) avDistanceToEnd /= distanceToEnd.size();
526 
527  if (fDebug > 2)
528  std::cout << "Distance to vertex is " << avDistanceToVertex << " and distance to end is "
529  << avDistanceToEnd << '\n';
530 
531  // Change order if necessary
532  if (avDistanceToEnd > avDistanceToVertex) {
533  std::reverse(xyz.begin(), xyz.end());
534  std::transform(
535  dircos.begin(), dircos.end(), dircos.begin(), [](TVector3 const& vec) { return -1 * vec; });
536  }
537 
538  if (xyz.size() != dircos.size())
539  mf::LogError("EMShowerAlg") << "Problem converting pma::Track3D to recob::Track";
540 
541  track = std::make_unique<recob::Track>(
544  recob::Track::Flags_t(xyz.size()),
545  false),
546  0,
547  -1.,
548  0,
551  -1);
552 
553  return track;
554 }
555 
556 std::unique_ptr<recob::Track> shower::EMShowerAlg::ConstructTrack(
557  detinfo::DetectorPropertiesData const& detProp,
558  std::vector<art::Ptr<recob::Hit>> const& track1,
559  std::vector<art::Ptr<recob::Hit>> const& track2) const
560 {
561  std::map<int, TVector2> showerCentreMap;
562  return ConstructTrack(detProp, track1, track2, showerCentreMap);
563 }
564 
566  detinfo::DetectorPropertiesData const& detProp,
567  std::vector<art::Ptr<recob::Hit>> const& trackHits,
568  std::unique_ptr<recob::Track> const& track) const
569 {
570  assert(not empty(trackHits));
571  if (!track) return -999;
572 
573  recob::Hit const& firstHit = *trackHits.front();
574 
575  // Get the pitch
576  double pitch = 0;
577  try {
578  pitch = lar::util::TrackPitchInView(*track, firstHit.View());
579  }
580  catch (...) {
581  }
582 
583  // Deal with large pitches
584  if (pitch > fdEdxTrackLength) {
585  return fCalorimetryAlg.dEdx_AREA(clockData, detProp, firstHit, pitch);
586  }
587 
588  double totalCharge = 0, totalDistance = 0, avHitTime = 0;
589  unsigned int nHits = 0;
590 
591  for (auto const& hit : trackHits | views::transform(to_element)) {
592  if (totalDistance + pitch < fdEdxTrackLength) {
593  totalDistance += pitch;
594  totalCharge += hit.Integral();
595  avHitTime += hit.PeakTime();
596  ++nHits;
597  }
598  }
599 
600  avHitTime /= (double)nHits;
601 
602  double const dQdx = totalCharge / totalDistance;
603  return fCalorimetryAlg.dEdx_AREA(clockData, detProp, dQdx, avHitTime, firstHit.WireID().Plane);
604 }
605 
607  detinfo::DetectorPropertiesData const& detProp,
608  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap,
609  std::unique_ptr<recob::Track>& initialTrack,
610  std::map<int, std::vector<art::Ptr<recob::Hit>>>& initialTrackHits) const
611 {
612 
616 
617  // Now find the hits belonging to the track
618  if (fDebug > 1)
619  std::cout << " ------------------ Finding initial track hits "
620  "-------------------- \n";
621  initialTrackHits = FindShowerStart_(showerHitsMap);
622  if (fDebug > 1) {
623  std::cout << "Here are the initial shower hits... \n";
624  for (auto const& [plane, hitPtrs] : initialTrackHits) {
625  std::cout << " Plane " << plane << '\n';
626  for (auto const& hit : hitPtrs | views::transform(to_element)) {
627  std::cout << " Hit is (" << HitCoordinates_(hit).X() << " (real hit " << hit.WireID()
628  << "), " << HitCoordinates_(hit).Y() << ")\n";
629  }
630  }
631  }
632  if (fDebug > 1)
633  std::cout << " ------------------ End finding initial track hits "
634  "-------------------- \n";
635 
636  // Now we have the track hits -- can make a track!
637  if (fDebug > 1) std::cout << " ------------------ Finding initial track -------------------- \n";
638  initialTrack = MakeInitialTrack_(detProp, initialTrackHits, showerHitsMap);
639  if (initialTrack and fDebug > 1) {
640  std::cout << "The track start is (" << initialTrack->Vertex().X() << ", "
641  << initialTrack->Vertex().Y() << ", " << initialTrack->Vertex().Z() << ")\n";
642  std::cout << "The track direction is (" << initialTrack->VertexDirection().X() << ", "
643  << initialTrack->VertexDirection().Y() << ", " << initialTrack->VertexDirection().Z()
644  << ")\n";
645  }
646  if (fDebug > 1)
647  std::cout << " ------------------ End finding initial track "
648  "-------------------- \n";
649 }
650 
651 std::vector<art::Ptr<recob::Hit>> shower::EMShowerAlg::FindOrderOfHits_(
652  detinfo::DetectorPropertiesData const& detProp,
654  bool perpendicular) const
655 {
656  // Find the charge-weighted centre (in [cm]) of this shower
657  TVector2 centre = ShowerCenter_(detProp, hits);
658 
659  // Find a rough shower 'direction'
660  TVector2 direction = ShowerDirection_(detProp, hits);
661 
662  if (perpendicular) direction = direction.Rotate(TMath::Pi() / 2);
663 
664  // Find how far each hit (projected onto this axis) is from the centre
665  TVector2 pos;
666  std::map<double, art::Ptr<recob::Hit>> hitProjection;
667  for (auto const& hitPtr : hits) {
668  pos = HitPosition_(detProp, *hitPtr) - centre;
669  hitProjection[direction * pos] = hitPtr;
670  }
671 
672  // Get a vector of hits in order of the shower
673  std::vector<art::Ptr<recob::Hit>> showerHits;
674  cet::transform_all(
675  hitProjection, std::back_inserter(showerHits), [](auto const& pr) { return pr.second; });
676 
677  // Make gradient plot
678  if (fMakeGradientPlot) {
679  std::map<int, TGraph*> graphs;
680  for (auto const& hit : showerHits | views::transform(to_element)) {
681  int tpc = hit.WireID().TPC;
682  if (graphs[tpc] == nullptr) graphs[tpc] = new TGraph();
683  graphs[tpc]->SetPoint(
684  graphs[tpc]->GetN(), HitPosition_(detProp, hit).X(), HitPosition_(detProp, hit).Y());
685  }
686  TMultiGraph* multigraph = new TMultiGraph();
687  for (auto const [color, graph] : graphs) {
688  graph->SetMarkerColor(color);
689  graph->SetMarkerStyle(8);
690  graph->SetMarkerSize(2);
691  multigraph->Add(graph);
692  }
693  TCanvas* canvas = new TCanvas();
694  multigraph->Draw("AP");
695  TLine line;
696  line.SetLineColor(2);
697  line.DrawLine(centre.X() - 1000 * direction.X(),
698  centre.Y() - 1000 * direction.Y(),
699  centre.X() + 1000 * direction.X(),
700  centre.Y() + 1000 * direction.Y());
701  canvas->SaveAs("Gradient.png");
702  delete canvas;
703  delete multigraph;
704  }
705 
706  return showerHits;
707 }
708 
709 std::vector<std::vector<int>> shower::EMShowerAlg::FindShowers(
710  std::map<int, std::vector<int>> const& trackToClusters) const
711 {
712  // Showers are vectors of clusters
713  std::vector<std::vector<int>> showers;
714 
715  // Loop over all tracks
716  for (auto const& clusters : trackToClusters | views::values) {
717 
718  // Find which showers already made are associated with this track
719  std::vector<int> matchingShowers;
720  for (unsigned int shower = 0; shower < showers.size(); ++shower)
721  for (int const cluster : clusters) {
722  if (cet::search_all(showers[shower], cluster) and
723  not cet::search_all(matchingShowers, shower)) {
724  matchingShowers.push_back(shower);
725  }
726  }
727 
728  // THINK THERE PROBABLY CAN BE MORE THAN ONE!
729  // IN FACT, THIS WOULD BE A SUCCESS OF THE SHOWERING METHOD!
730  // // Shouldn't be more than one
731  // if (matchingShowers.size() > 1)
732  // mf::LogInfo("EMShowerAlg") << "Number of showers this track matches is " << matchingShowers.size() << std::endl;
733 
734  // New shower
735  if (matchingShowers.size() < 1) showers.push_back(clusters);
736 
737  // Add to existing shower
738  else {
739  for (int const cluster : clusters) {
740  if (not cet::search_all(showers.at(matchingShowers[0]), cluster))
741  showers.at(matchingShowers.at(0)).push_back(cluster);
742  }
743  }
744  }
745 
746  return showers;
747 }
748 
749 std::map<int, std::vector<art::Ptr<recob::Hit>>> shower::EMShowerAlg::FindShowerStart_(
750  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& orderedShowerMap) const
751 {
752 
753  std::map<int, std::vector<art::Ptr<recob::Hit>>> initialHitsMap;
754 
755  for (auto const& [plane, orderedShower] : orderedShowerMap) {
756  std::vector<art::Ptr<recob::Hit>> initialHits;
757 
758  // Find if the shower is traveling along ticks or wires
759  bool wireDirection = true;
760  std::vector<int> wires;
761  for (auto const& hit : orderedShower | views::transform(to_element))
762  wires.push_back(std::round(HitCoordinates_(hit).X()));
763 
764  cet::sort_all(wires);
765  if (std::abs(*wires.begin() - std::round(HitCoordinates_(**orderedShower.begin()).X())) > 3 and
766  std::abs(*wires.rbegin() - std::round(HitCoordinates_(**orderedShower.begin()).X())) > 3)
767  wireDirection = false;
768 
769  // Deal with showers traveling along wires
770  if (wireDirection) {
771  bool increasing = HitCoordinates_(**orderedShower.rbegin()).X() >
772  HitCoordinates_(**orderedShower.begin()).X();
773  std::map<int, std::vector<art::Ptr<recob::Hit>>> wireHitMap;
774  int multipleWires = 0;
775  for (auto const& hitPtr : orderedShower)
776  wireHitMap[std::round(HitCoordinates_(*hitPtr).X())].push_back(hitPtr);
777 
778  for (auto const& hitPtr : orderedShower) {
779  int wire = std::round(HitCoordinates_(*hitPtr).X());
780  if (wireHitMap[wire].size() > 1) {
781  ++multipleWires;
782  if (multipleWires > 5) break;
783  continue;
784  }
785  else if ((increasing and wireHitMap[wire + 1].size() > 1 and
786  (wireHitMap[wire + 2].size() > 1 or wireHitMap[wire + 3].size() > 1)) or
787  (!increasing and wireHitMap[wire - 1].size() > 1 and
788  (wireHitMap[wire - 2].size() > 1 or wireHitMap[wire - 3].size() > 1))) {
789  if ((increasing and
790  (wireHitMap[wire + 4].size() < 2 and wireHitMap[wire + 5].size() < 2 and
791  wireHitMap[wire + 6].size() < 2 and wireHitMap[wire + 9].size() > 1)) or
792  (!increasing and
793  (wireHitMap[wire - 4].size() < 2 and wireHitMap[wire - 5].size() < 2 and
794  wireHitMap[wire - 6].size() < 2) and
795  wireHitMap[wire - 9].size() > 1))
796  initialHits.push_back(hitPtr);
797  else
798  break;
799  }
800  else
801  initialHits.push_back(hitPtr);
802  }
803  if (!initialHits.size()) initialHits.push_back(*orderedShower.begin());
804  }
805 
806  // Deal with showers travelling along ticks
807  else {
808  bool increasing = HitCoordinates_(**orderedShower.rbegin()).Y() >
809  HitCoordinates_(**orderedShower.begin()).Y();
810  std::map<int, std::vector<art::Ptr<recob::Hit>>> tickHitMap;
811  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hitIt = orderedShower.begin();
812  hitIt != orderedShower.end();
813  ++hitIt)
814  tickHitMap[std::round(HitCoordinates_(**hitIt).Y())].push_back(*hitIt);
815 
816  for (auto const& hitPtr : orderedShower) {
817  int const tick = std::round(HitCoordinates_(*hitPtr).Y());
818  if ((increasing and (tickHitMap[tick + 1].size() or tickHitMap[tick + 2].size() or
819  tickHitMap[tick + 3].size() or tickHitMap[tick + 4].size() or
820  tickHitMap[tick + 5].size())) or
821  (!increasing and (tickHitMap[tick - 1].size() or tickHitMap[tick - 2].size() or
822  tickHitMap[tick - 3].size() or tickHitMap[tick - 4].size() or
823  tickHitMap[tick - 5].size())))
824  break;
825  else
826  initialHits.push_back(hitPtr);
827  }
828  if (initialHits.empty()) initialHits.push_back(*orderedShower.begin());
829  }
830 
831  // Need at least two hits to make a track
832  if (initialHits.size() == 1 and orderedShower.size() > 2)
833  initialHits.push_back(orderedShower[1]);
834 
835  // Quality check -- make sure there isn't a single hit in a different TPC (artefact of tracking failure)
836  std::vector<art::Ptr<recob::Hit>> newInitialHits;
837  std::map<int, int> tpcHitMap;
838  std::vector<int> singleHitTPCs;
839  for (auto const& hit : initialHits | views::transform(to_element))
840  ++tpcHitMap[hit.WireID().TPC];
841 
842  for (auto const [tpc, count] : tpcHitMap)
843  if (count == 1) singleHitTPCs.push_back(tpc);
844 
845  if (singleHitTPCs.size()) {
846  if (fDebug > 2)
847  for (int const tpc : singleHitTPCs)
848  std::cout << "Removed hits in TPC " << tpc << '\n';
849 
850  for (auto const& hitPtr : initialHits)
851  if (not cet::search_all(singleHitTPCs, hitPtr->WireID().TPC))
852  newInitialHits.push_back(hitPtr);
853  if (!newInitialHits.size()) newInitialHits.push_back(*orderedShower.begin());
854  }
855  else
856  newInitialHits = initialHits;
857 
858  initialHitsMap[plane] = newInitialHits;
859  }
860 
861  return initialHitsMap;
862 }
863 
864 std::map<int, std::map<int, bool>> shower::EMShowerAlg::GetPlanePermutations_(
865  const detinfo::DetectorPropertiesData& detProp,
866  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
867 {
868 
869  // The map to return
870  std::map<int, std::map<int, bool>> permutations;
871 
872  // Get the properties of the shower hits across the planes which will be used to
873  // determine the likelihood of a particular reorientation permutation
874  // -- relative width in the wire direction (if showers travel along the wire
875  // direction in a particular plane)
876  // -- the RMS gradients (how likely it is the RMS of the hit positions from
877  // central axis increases along a particular orientation)
878 
879  // Find the RMS, RMS gradient and wire widths
880  std::map<int, double> planeRMSGradients, planeRMS;
881  for (auto const& [plane, hitPtrs] : showerHitsMap) {
882  planeRMS[plane] = ShowerHitRMS_(detProp, hitPtrs);
883  planeRMSGradients[plane] = ShowerHitRMSGradient_(detProp, hitPtrs);
884  }
885 
886  // Order these backwards so they can be used to discriminate between planes
887  std::map<double, int> gradientMap;
888  for (int const plane : showerHitsMap | views::keys)
889  gradientMap[std::abs(planeRMSGradients.at(plane))] = plane;
890 
891  std::map<double, int> wireWidthMap = RelativeWireWidth_(showerHitsMap);
892 
893  if (fDebug > 1)
894  for (auto const [gradient, plane] : wireWidthMap)
895  std::cout << "Plane " << plane << " has relative width in wire of " << gradient
896  << "; and an RMS gradient of " << planeRMSGradients[plane] << '\n';
897 
898  // Find the correct permutations
899  int perm = 0;
900  std::vector<std::map<int, bool>> usedPermutations;
901 
902  // Most likely is to not change anything
903  for (int const plane : showerHitsMap | views::keys)
904  permutations[perm][plane] = 0;
905  ++perm;
906 
907  // Use properties of the shower to determine the middle cases
908  for (int const plane : wireWidthMap | views::values) {
909  std::map<int, bool> permutation;
910  permutation[plane] = true;
911  for (int const plane2 : wireWidthMap | views::values)
912  if (plane != plane2) permutation[plane2] = false;
913 
914  if (not cet::search_all(usedPermutations, permutation)) {
915  permutations[perm] = permutation;
916  usedPermutations.push_back(permutation);
917  ++perm;
918  }
919  }
920  for (int const plane : wireWidthMap | views::reverse | views::values) {
921  std::map<int, bool> permutation;
922  permutation[plane] = false;
923  for (int const plane2 : wireWidthMap | views::values)
924  if (plane != plane2) permutation[plane2] = true;
925 
926  if (not cet::search_all(usedPermutations, permutation)) {
927  permutations[perm] = permutation;
928  usedPermutations.push_back(permutation);
929  ++perm;
930  }
931  }
932 
933  // Least likely is to change everything
934  for (int const plane : showerHitsMap | views::keys)
935  permutations[perm][plane] = 1;
936  ++perm;
937 
938  if (fDebug > 1) {
939  std::cout << "Here are the permutations!\n";
940  for (auto const& [index, permutation] : permutations) {
941  std::cout << "Permutation " << index << '\n';
942  for (auto const [plane, value] : permutation)
943  std::cout << " Plane " << plane << " has value " << value << '\n';
944  }
945  }
946 
947  return permutations;
948 }
949 
950 std::unique_ptr<recob::Track> shower::EMShowerAlg::MakeInitialTrack_(
951  detinfo::DetectorPropertiesData const& detProp,
952  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& initialHitsMap,
953  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& showerHitsMap) const
954 {
955  // Can't do much with just one view
956  if (initialHitsMap.size() < 2) {
957  mf::LogInfo("EMShowerAlg") << "Only one useful view for this shower; nothing can be done\n";
958  return std::unique_ptr<recob::Track>();
959  }
960 
961  std::vector<std::pair<int, int>> initialHitsSize;
962  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator initialHitIt =
963  initialHitsMap.begin();
964  initialHitIt != initialHitsMap.end();
965  ++initialHitIt)
966  initialHitsSize.push_back(std::make_pair(initialHitIt->first, initialHitIt->second.size()));
967 
968  // Sort the planes by number of hits
969  std::sort(initialHitsSize.begin(),
970  initialHitsSize.end(),
971  [](std::pair<int, int> const& size1, std::pair<int, int> const& size2) {
972  return size1.second > size2.second;
973  });
974 
975  // Pick the two planes to use to make the track
976  // -- if more than two planes, can choose the two views
977  // more accurately
978  // -- if not, just use the two views available
979 
980  std::vector<int> planes;
981 
982  if (initialHitsSize.size() > 2 and !CheckShowerHits_(detProp, showerHitsMap)) {
983  int planeToIgnore = WorstPlane_(showerHitsMap);
984  if (fDebug > 1)
985  std::cout << "Igoring plane " << planeToIgnore << " in creation of initial track\n";
986  for (std::vector<std::pair<int, int>>::iterator hitsSizeIt = initialHitsSize.begin();
987  hitsSizeIt != initialHitsSize.end() and planes.size() < 2;
988  ++hitsSizeIt) {
989  if (hitsSizeIt->first == planeToIgnore) continue;
990  planes.push_back(hitsSizeIt->first);
991  }
992  }
993  else
994  planes = {initialHitsSize.at(0).first, initialHitsSize.at(1).first};
995 
996  return ConstructTrack(detProp, initialHitsMap.at(planes.at(0)), initialHitsMap.at(planes.at(1)));
997 }
998 
1000  detinfo::DetectorClocksData const& clockData,
1001  detinfo::DetectorPropertiesData const& detProp,
1003  std::unique_ptr<recob::Track> const& initialTrack,
1004  std::map<int, std::vector<art::Ptr<recob::Hit>>> const& initialHitsMap) const
1005 {
1006 
1007  // Find the shower hits on each plane
1008  std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
1009  for (auto const& hitPtr : hits)
1010  planeHitsMap[hitPtr->View()].push_back(hitPtr);
1011 
1012  int bestPlane = -1;
1013  unsigned int highestNumberOfHits = 0;
1014  std::vector<double> totalEnergy, totalEnergyError, dEdx, dEdxError;
1015 
1016  // Look at each of the planes
1017  for (unsigned int plane = 0; plane < fWireReadoutGeom->MaxPlanes(); ++plane) {
1018 
1019  // If there's hits on this plane...
1020  if (planeHitsMap.count(plane)) {
1021  dEdx.push_back(FinddEdx_(clockData, detProp, initialHitsMap.at(plane), initialTrack));
1022  totalEnergy.push_back(
1023  fShowerEnergyAlg.ShowerEnergy(clockData, detProp, planeHitsMap.at(plane), plane));
1024  if (planeHitsMap.at(plane).size() > highestNumberOfHits and initialHitsMap.count(plane)) {
1025  bestPlane = plane;
1026  highestNumberOfHits = planeHitsMap.at(plane).size();
1027  }
1028  }
1029 
1030  // If not...
1031  else {
1032  dEdx.push_back(0);
1033  totalEnergy.push_back(0);
1034  }
1035  }
1036 
1037  TVector3 direction, directionError, showerStart, showerStartError;
1038  if (initialTrack) {
1039  direction = initialTrack->VertexDirection<TVector3>();
1040  showerStart = initialTrack->Vertex<TVector3>();
1041  }
1042 
1043  if (fDebug > 0) {
1044  std::cout << "Best plane is " << bestPlane;
1045  std::cout << "\ndE/dx for each plane is: " << dEdx[0] << ", " << dEdx[1] << " and " << dEdx[2];
1046  std::cout << "\nTotal energy for each plane is: " << totalEnergy[0] << ", " << totalEnergy[1]
1047  << " and " << totalEnergy[2];
1048  std::cout << "\nThe shower start is (" << showerStart.X() << ", " << showerStart.Y() << ", "
1049  << showerStart.Z() << ")\n";
1050  std::cout << "The shower direction is (" << direction.X() << ", " << direction.Y() << ", "
1051  << direction.Z() << ")\n";
1052  }
1053 
1054  return recob::Shower(direction,
1055  directionError,
1056  showerStart,
1057  showerStartError,
1058  totalEnergy,
1059  totalEnergyError,
1060  dEdx,
1061  dEdxError,
1062  bestPlane);
1063 }
1064 
1066  detinfo::DetectorPropertiesData const& detProp,
1069  int& iok) const
1070 {
1071  iok = 1;
1072 
1073  // Find the shower hits on each plane
1074  std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
1075  for (auto const& hitPtr : hits)
1076  planeHitsMap[hitPtr->WireID().Plane].push_back(hitPtr);
1077 
1078  std::vector<std::vector<art::Ptr<recob::Hit>>> initialTrackHits(3);
1079 
1080  int pl0 = -1;
1081  int pl1 = -1;
1082  unsigned maxhits0 = 0;
1083  unsigned maxhits1 = 0;
1084 
1085  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator planeHits = planeHitsMap.begin();
1086  planeHits != planeHitsMap.end();
1087  ++planeHits) {
1088 
1089  std::vector<art::Ptr<recob::Hit>> showerHits;
1090  OrderShowerHits_(detProp, planeHits->second, showerHits, vertex);
1091  FindInitialTrackHits(showerHits, vertex, initialTrackHits[planeHits->first]);
1092  if ((planeHits->second).size() > maxhits0) {
1093  if (pl0 != -1) {
1094  maxhits1 = maxhits0;
1095  pl1 = pl0;
1096  }
1097  pl0 = planeHits->first;
1098  maxhits0 = (planeHits->second).size();
1099  }
1100  else if ((planeHits->second).size() > maxhits1) {
1101  pl1 = planeHits->first;
1102  maxhits1 = (planeHits->second).size();
1103  }
1104  }
1105  if (pl0 != -1 && pl1 != -1 && initialTrackHits[pl0].size() >= 2 &&
1106  initialTrackHits[pl1].size() >= 2 &&
1107  initialTrackHits[pl0][0]->WireID().TPC == initialTrackHits[pl1][0]->WireID().TPC) {
1108  double xyz[3];
1109  vertex->XYZ(xyz);
1110  TVector3 vtx(xyz);
1111  pma::Track3D* pmatrack =
1112  fProjectionMatchingAlg.buildSegment(detProp, initialTrackHits[pl0], initialTrackHits[pl1]);
1113  std::vector<TVector3> spts;
1114 
1115  for (size_t i = 0; i < pmatrack->size(); ++i) {
1116  if ((*pmatrack)[i]->IsEnabled()) {
1117  TVector3 p3d = (*pmatrack)[i]->Point3D();
1118  spts.push_back(p3d);
1119  }
1120  }
1121  if (spts.size() >= 2) { // at least two space points
1122  TVector3 shwxyz, shwxyzerr;
1123  TVector3 shwdir, shwdirerr;
1124  std::vector<double> totalEnergy, totalEnergyError, dEdx, dEdxError;
1125  int bestPlane = pl0;
1126  double minpitch = 1000;
1127  std::vector<TVector3> dirs;
1128  if ((spts[0] - vtx).Mag() < (spts.back() - vtx).Mag()) {
1129  shwxyz = spts[0];
1130  size_t i = 5;
1131  if (spts.size() - 1 < 5) i = spts.size() - 1;
1132  shwdir = spts[i] - spts[0];
1133  shwdir = shwdir.Unit();
1134  }
1135  else {
1136  shwxyz = spts.back();
1137  size_t i = 0;
1138  if (spts.size() > 6) i = spts.size() - 6;
1139  shwdir = spts[i] - spts[spts.size() - 1];
1140  shwdir = shwdir.Unit();
1141  }
1142  for (unsigned int plane = 0; plane < fWireReadoutGeom->MaxPlanes(); ++plane) {
1143  if (planeHitsMap.find(plane) != planeHitsMap.end()) {
1144  totalEnergy.push_back(
1145  fShowerEnergyAlg.ShowerEnergy(clockData, detProp, planeHitsMap[plane], plane));
1146  }
1147  else {
1148  totalEnergy.push_back(0);
1149  }
1150  if (initialTrackHits[plane].size()) {
1151  double fdEdx = 0;
1152  double totQ = 0;
1153  double avgT = 0;
1154  double pitch = 0;
1155  double wirepitch =
1156  fWireReadoutGeom->Plane(initialTrackHits[plane][0]->WireID().planeID()).WirePitch();
1157  double angleToVert = fWireReadoutGeom->WireAngleToVertical(
1158  fWireReadoutGeom->Plane(geo::PlaneID{0, 0, plane}).View(),
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 < fWireReadoutGeom->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)fWireReadoutGeom->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  auto const& planeID = hit0->WireID().parentID();
1598  TVector2 coordvtx = TVector2(fWireReadoutGeom->Plane(planeID).WireCoordinate(xyz),
1599  detProp.ConvertXToTicks(xyz.X(), planeID));
1600  if ((coord1 - coordvtx).Mod() < (coord0 - coordvtx).Mod()) {
1601  std::reverse(showerHits.begin(), showerHits.end());
1602  }
1603 }
1604 
1607  std::vector<art::Ptr<recob::Hit>>& trackHits) const
1608 {
1609  // Find TPC for the vertex
1610  auto const& xyz = vertex->position();
1611  geo::TPCID tpc = fGeom->FindTPCAtPosition(xyz);
1612 
1613  // vertex cannot be projected into a TPC, find the TPC that has the most hits
1614  if (!tpc.isValid) {
1615  std::map<geo::TPCID, unsigned int> tpcmap;
1616  unsigned maxhits = 0;
1617  for (auto const& hit : showerHits) {
1618  ++tpcmap[geo::TPCID(hit->WireID())];
1619  }
1620  for (auto const& t : tpcmap) {
1621  if (t.second > maxhits) {
1622  maxhits = t.second;
1623  tpc = t.first;
1624  }
1625  }
1626  }
1627  if (!tpc.isValid) return;
1628 
1629  double parm[2];
1630  int fitok = 0;
1631  std::vector<double> wfit;
1632  std::vector<double> tfit;
1633  std::vector<double> cfit;
1634 
1635  for (size_t i = 0; i < fNfitpass; ++i) {
1636 
1637  // Fit a straight line through hits
1638  unsigned int nhits = 0;
1639  for (auto& hit : showerHits) {
1640  if (hit->WireID().TPC == tpc.TPC) {
1641  TVector2 coord = HitCoordinates_(*hit);
1642  if (i == 0 ||
1643  (std::abs((coord.Y() - (parm[0] + coord.X() * parm[1])) * cos(atan(parm[1]))) <
1644  fToler[i - 1]) ||
1645  fitok == 1) {
1646  ++nhits;
1647  if (nhits == fNfithits[i] + 1) break;
1648  wfit.push_back(coord.X());
1649  tfit.push_back(coord.Y());
1650  cfit.push_back(1.);
1651  if (i == fNfitpass - 1) { trackHits.push_back(hit); }
1652  }
1653  }
1654  }
1655 
1656  if (i < fNfitpass - 1 && wfit.size()) {
1657  fitok = WeightedFit(wfit.size(), &wfit[0], &tfit[0], &cfit[0], &parm[0]);
1658  }
1659  wfit.clear();
1660  tfit.clear();
1661  cfit.clear();
1662  }
1663 }
1664 
1666 {
1667  return TVector2(GlobalWire_(hit.WireID()), hit.PeakTime());
1668 }
1669 
1671  recob::Hit const& hit) const
1672 {
1673  geo::PlaneID planeID = hit.WireID().planeID();
1674  return HitPosition_(detProp, HitCoordinates_(hit), planeID);
1675 }
1676 
1678  TVector2 const& pos,
1679  geo::PlaneID planeID) const
1680 {
1681  return {pos.X() * fWireReadoutGeom->Plane(planeID).WirePitch(),
1682  detProp.ConvertTicksToX(pos.Y(), planeID)};
1683 }
1684 
1686 {
1687  double globalWire = -999;
1688 
1689  // Induction
1690  if (fWireReadoutGeom->SignalType(wireID) == geo::kInduction) {
1691  auto const wireCenter = fWireReadoutGeom->Wire(wireID).GetCenter();
1692  globalWire = fWireReadoutGeom
1693  ->Plane({wireID.Cryostat,
1694  wireID.TPC % 2, // 0 or 1
1695  wireID.Plane})
1696  .WireCoordinate(wireCenter);
1697  }
1698 
1699  // Collection
1700  else {
1701  // FOR COLLECTION WIRES, HARD CODE THE GEOMETRY FOR GIVEN DETECTORS
1702  // THIS _SHOULD_ BE TEMPORARY. GLOBAL WIRE SUPPORT IS BEING ADDED TO THE LARSOFT GEOMETRY AND SHOULD BE AVAILABLE SOON
1703  if (fDetector == "dune35t") {
1704  unsigned int nwires =
1705  fWireReadoutGeom->Nwires(geo::PlaneID{wireID.Cryostat, 0, wireID.Plane});
1706  if (wireID.TPC == 0 or wireID.TPC == 1)
1707  globalWire = wireID.Wire;
1708  else if (wireID.TPC == 2 or wireID.TPC == 3 or wireID.TPC == 4 or wireID.TPC == 5)
1709  globalWire = nwires + wireID.Wire;
1710  else if (wireID.TPC == 6 or wireID.TPC == 7)
1711  globalWire = (2 * nwires) + wireID.Wire;
1712  else
1713  mf::LogError("BlurredClusterAlg")
1714  << "Error when trying to find a global induction plane coordinate for TPC " << wireID.TPC
1715  << " (geometry" << fDetector << ")";
1716  }
1717  else if (fDetector == "dune10kt") {
1718  unsigned int nwires =
1719  fWireReadoutGeom->Nwires(geo::PlaneID{wireID.Cryostat, 0, wireID.Plane});
1720  // Detector geometry has four TPCs, two on top of each other, repeated along z...
1721  int block = wireID.TPC / 4;
1722  globalWire = (nwires * block) + wireID.Wire;
1723  }
1724  else {
1725  auto const wireCenter = fWireReadoutGeom->Wire(wireID).GetCenter();
1726  globalWire = fWireReadoutGeom
1727  ->Plane({wireID.Cryostat,
1728  wireID.TPC % 2, // 0 or 1
1729  wireID.Plane})
1730  .WireCoordinate(wireCenter);
1731  }
1732  }
1733 
1734  return globalWire;
1735 }
1736 
1738  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
1739 {
1740 
1741  // Get the wire widths
1742  std::map<int, int> planeWireLength;
1743  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitsIt =
1744  showerHitsMap.begin();
1745  showerHitsIt != showerHitsMap.end();
1746  ++showerHitsIt)
1747  planeWireLength[showerHitsIt->first] =
1748  std::abs(HitCoordinates_(*showerHitsIt->second.front()).X() -
1749  HitCoordinates_(*showerHitsIt->second.back()).X());
1750 
1751  // Find the relative wire width for each plane with respect to the others
1752  std::map<int, double> planeOtherWireLengths;
1753  for (std::map<int, int>::iterator planeWireLengthIt = planeWireLength.begin();
1754  planeWireLengthIt != planeWireLength.end();
1755  ++planeWireLengthIt) {
1756  double quad = 0.;
1757  for (int plane = 0; plane < (int)fWireReadoutGeom->MaxPlanes(); ++plane)
1758  if (plane != planeWireLengthIt->first and planeWireLength.count(plane))
1759  quad += cet::square(planeWireLength[plane]);
1760  quad = std::sqrt(quad);
1761  planeOtherWireLengths[planeWireLengthIt->first] =
1762  planeWireLength[planeWireLengthIt->first] / (double)quad;
1763  }
1764 
1765  // Order these backwards
1766  std::map<double, int> wireWidthMap;
1767  for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitsIt =
1768  showerHitsMap.begin();
1769  showerHitsIt != showerHitsMap.end();
1770  ++showerHitsIt)
1771  wireWidthMap[planeOtherWireLengths.at(showerHitsIt->first)] = showerHitsIt->first;
1772 
1773  return wireWidthMap;
1774 }
1775 
1777  detinfo::DetectorPropertiesData const& detProp,
1778  const std::vector<art::Ptr<recob::Hit>>& showerHits) const
1779 {
1780 
1781  TVector2 pos;
1782  double weight = 1;
1783  double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1784  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hit = showerHits.begin();
1785  hit != showerHits.end();
1786  ++hit) {
1787  //++nhits;
1788  pos = HitPosition_(detProp, **hit);
1789  weight = cet::square((*hit)->Integral());
1790  sumweight += weight;
1791  sumx += weight * pos.X();
1792  sumy += weight * pos.Y();
1793  sumx2 += weight * pos.X() * pos.X();
1794  sumxy += weight * pos.X() * pos.Y();
1795  }
1796  double gradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1797  TVector2 direction = TVector2(1, gradient).Unit();
1798 
1799  return direction;
1800 }
1801 
1803  detinfo::DetectorPropertiesData const& detProp,
1804  std::vector<art::Ptr<recob::Hit>> const& showerHits) const
1805 {
1806 
1807  TVector2 pos, chargePoint = TVector2(0, 0);
1808  double totalCharge = 0;
1809  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hit = showerHits.begin();
1810  hit != showerHits.end();
1811  ++hit) {
1812  pos = HitPosition_(detProp, **hit);
1813  chargePoint += (*hit)->Integral() * pos;
1814  totalCharge += (*hit)->Integral();
1815  }
1816  TVector2 centre = chargePoint / totalCharge;
1817 
1818  return centre;
1819 }
1820 
1822  const std::vector<art::Ptr<recob::Hit>>& showerHits) const
1823 {
1824 
1825  TVector2 direction = ShowerDirection_(detProp, showerHits);
1826  TVector2 centre = ShowerCenter_(detProp, showerHits);
1827 
1828  std::vector<double> distanceToAxis;
1829  for (std::vector<art::Ptr<recob::Hit>>::const_iterator showerHitsIt = showerHits.begin();
1830  showerHitsIt != showerHits.end();
1831  ++showerHitsIt) {
1832  TVector2 proj = (HitPosition_(detProp, **showerHitsIt) - centre).Proj(direction) + centre;
1833  distanceToAxis.push_back((HitPosition_(detProp, **showerHitsIt) - proj).Mod());
1834  }
1835  double RMS = TMath::RMS(distanceToAxis.begin(), distanceToAxis.end());
1836 
1837  return RMS;
1838 }
1839 
1841  detinfo::DetectorPropertiesData const& detProp,
1842  const std::vector<art::Ptr<recob::Hit>>& showerHits) const
1843 {
1844  // Find a rough shower 'direction' and center
1845  TVector2 direction = ShowerDirection_(detProp, showerHits);
1846 
1847  // Bin the hits into discreet chunks
1848  int nShowerSegments = fNumShowerSegments;
1849  double lengthOfShower =
1850  (HitPosition_(detProp, *showerHits.back()) - HitPosition_(detProp, *showerHits.front())).Mod();
1851  double lengthOfSegment = lengthOfShower / (double)nShowerSegments;
1852  std::map<int, std::vector<art::Ptr<recob::Hit>>> showerSegments;
1853  std::map<int, double> segmentCharge;
1854  for (auto const& hitPtr : showerHits) {
1855  auto const& hit = *hitPtr;
1856  int const segment =
1857  (HitPosition_(detProp, hit) - HitPosition_(detProp, *showerHits.front())).Mod() /
1858  lengthOfSegment;
1859  showerSegments[segment].push_back(hitPtr);
1860  segmentCharge[segment] += hit.Integral();
1861  }
1862 
1863  TGraph* graph = new TGraph();
1864  std::vector<std::pair<int, double>> binVsRMS;
1865 
1866  // Loop over the bins to find the distribution of hits as the shower
1867  // progresses
1868  for (auto const& [segment, hitPtrs] : showerSegments) {
1869 
1870  // Get the mean position of the hits in this bin
1871  TVector2 meanPosition(0, 0);
1872  for (auto const& hit : hitPtrs | views::transform(to_element))
1873  meanPosition += HitPosition_(detProp, hit);
1874  meanPosition /= (double)hitPtrs.size();
1875 
1876  // Get the RMS of this bin
1877  std::vector<double> distanceToAxisBin;
1878  for (auto const& hit : hitPtrs | views::transform(to_element)) {
1879  TVector2 proj = (HitPosition_(detProp, hit) - meanPosition).Proj(direction) + meanPosition;
1880  distanceToAxisBin.push_back((HitPosition_(detProp, hit) - proj).Mod());
1881  }
1882 
1883  double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
1884  binVsRMS.emplace_back(segment, RMSBin);
1885  if (fMakeRMSGradientPlot) graph->SetPoint(graph->GetN(), segment, RMSBin);
1886  }
1887 
1888  // Get the gradient of the RMS-bin plot
1889  double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1890  for (auto const& [bin, RMSBin] : binVsRMS) {
1891  double weight = segmentCharge.at(bin);
1892  sumweight += weight;
1893  sumx += weight * bin;
1894  sumy += weight * RMSBin;
1895  sumx2 += weight * bin * bin;
1896  sumxy += weight * bin * RMSBin;
1897  }
1898  double RMSgradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1899 
1900  if (fMakeRMSGradientPlot) {
1901  TVector2 direction = TVector2(1, RMSgradient).Unit();
1902  TCanvas* canv = new TCanvas();
1903  graph->Draw();
1904  graph->GetXaxis()->SetTitle("Shower segment");
1905  graph->GetYaxis()->SetTitle("RMS of hit distribution");
1906  TVector2 centre = TVector2(graph->GetMean(1), graph->GetMean(2));
1907  TLine line;
1908  line.SetLineColor(2);
1909  line.DrawLine(centre.X() - 1000 * direction.X(),
1910  centre.Y() - 1000 * direction.Y(),
1911  centre.X() + 1000 * direction.X(),
1912  centre.Y() + 1000 * direction.Y());
1913  canv->SaveAs("RMSGradient.png");
1914  delete canv;
1915  }
1916  delete graph;
1917 
1918  return RMSgradient;
1919 }
1920 
1922  detinfo::DetectorPropertiesData const& detProp,
1923  geo::Point_t const& point,
1924  int plane,
1925  int cryostat) const
1926 {
1927  TVector2 wireTickPos{-999., -999.};
1928 
1929  geo::TPCID tpcID = fGeom->FindTPCAtPosition(point);
1930  int tpc = 0;
1931  if (tpcID.isValid)
1932  tpc = tpcID.TPC;
1933  else
1934  return wireTickPos;
1935 
1936  // Construct wire ID for this point projected onto the plane
1937  geo::PlaneID const planeID(cryostat, tpc, plane);
1938  geo::WireID wireID;
1939  try {
1940  wireID = fWireReadoutGeom->Plane(planeID).NearestWireID(point);
1941  }
1942  catch (geo::InvalidWireError const& e) {
1943  wireID = e.suggestedWireID(); // pick the closest valid wire
1944  }
1945 
1946  wireTickPos = TVector2(GlobalWire_(wireID), detProp.ConvertXToTicks(point.X(), planeID));
1947 
1948  return HitPosition_(detProp, wireTickPos, planeID);
1949 }
1950 
1952  const std::map<int, std::vector<art::Ptr<recob::Hit>>>& showerHitsMap) const
1953 {
1954  // Get the width of the shower in wire coordinate
1955  std::map<int, int> planeWireLength;
1956  std::map<int, double> planeOtherWireLengths;
1957  for (auto const& [plane, hits] : showerHitsMap)
1958  planeWireLength[plane] =
1959  std::abs(HitCoordinates_(*hits.front()).X() - HitCoordinates_(*hits.back()).X());
1960  for (std::map<int, int>::iterator planeWireLengthIt = planeWireLength.begin();
1961  planeWireLengthIt != planeWireLength.end();
1962  ++planeWireLengthIt) {
1963  double quad = 0.;
1964  for (int plane = 0; plane < (int)fWireReadoutGeom->MaxPlanes(); ++plane)
1965  if (plane != planeWireLengthIt->first and planeWireLength.count(plane))
1966  quad += cet::square(planeWireLength[plane]);
1967  quad = std::sqrt(quad);
1968  planeOtherWireLengths[planeWireLengthIt->first] =
1969  planeWireLength[planeWireLengthIt->first] / (double)quad;
1970  }
1971 
1972  if (fDebug > 1)
1973  for (auto const [plane, relative_width] : planeOtherWireLengths) {
1974  std::cout << "Plane " << plane << " has " << planeWireLength[plane]
1975  << " wire width and therefore has relative width in wire of " << relative_width
1976  << '\n';
1977  }
1978 
1979  std::map<double, int> wireWidthMap;
1980  for (int const plane : showerHitsMap | views::keys) {
1981  double wireWidth = planeWireLength.at(plane);
1982  wireWidthMap[wireWidth] = plane;
1983  }
1984 
1985  return wireWidthMap.begin()->second;
1986 }
1987 
1989  const Double_t* x,
1990  const Double_t* y,
1991  const Double_t* w,
1992  Double_t* parm) const
1993 {
1994 
1995  Double_t sumx = 0.;
1996  Double_t sumx2 = 0.;
1997  Double_t sumy = 0.;
1998  Double_t sumxy = 0.;
1999  Double_t sumw = 0.;
2000  Double_t eparm[2];
2001 
2002  parm[0] = 0.;
2003  parm[1] = 0.;
2004  eparm[0] = 0.;
2005  eparm[1] = 0.;
2006 
2007  for (Int_t i = 0; i < n; i++) {
2008  sumx += x[i] * w[i];
2009  sumx2 += x[i] * x[i] * w[i];
2010  sumy += y[i] * w[i];
2011  sumxy += x[i] * y[i] * w[i];
2012  sumw += w[i];
2013  }
2014 
2015  if (sumx2 * sumw - sumx * sumx == 0.) return 1;
2016  if (sumx2 - sumx * sumx / sumw == 0.) return 1;
2017 
2018  parm[0] = (sumy * sumx2 - sumx * sumxy) / (sumx2 * sumw - sumx * sumx);
2019  parm[1] = (sumxy - sumx * sumy / sumw) / (sumx2 - sumx * sumx / sumw);
2020 
2021  eparm[0] = sumx2 * (sumx2 * sumw - sumx * sumx);
2022  eparm[1] = (sumx2 - sumx * sumx / sumw);
2023 
2024  if (eparm[0] < 0. || eparm[1] < 0.) return 1;
2025 
2026  eparm[0] = sqrt(eparm[0]) / (sumx2 * sumw - sumx * sumx);
2027  eparm[1] = sqrt(eparm[1]) / (sumx2 - sumx * sumx / sumw);
2028 
2029  return 0;
2030 }
2031 
2033 {
2034  if (hits.empty()) return false;
2035  if (hits.size() > 2000) return true;
2036  if (hits.size() < 20) return true;
2037 
2038  std::map<int, int> hitmap;
2039  unsigned nhits = 0;
2040  for (auto const& hit : hits | views::transform(to_element)) {
2041  ++nhits;
2042  if (nhits > 2) ++hitmap[hit.WireID().Wire];
2043  if (nhits == 20) break;
2044  }
2045  if (float(nhits - 2) / hitmap.size() > 1.4)
2046  return false;
2047  else
2048  return true;
2049 }
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::WireReadoutGeom const * fWireReadoutGeom
Definition: EMShowerAlg.h:275
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.
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:219
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.
shower::ShowerEnergyAlg const fShowerEnergyAlg
Definition: EMShowerAlg.h:279
constexpr to_element_t to_element
Definition: ToElement.h:25
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
pma::ProjectionMatchingAlg const fProjectionMatchingAlg
Definition: EMShowerAlg.h:281
WireID suggestedWireID() const
Returns a better wire ID.
Definition: Exceptions.h:85
double ShowerEnergy(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &hits, geo::PlaneID::PlaneID_t plane) const
double WireCoordinate(Point_t const &point) const
Returns the coordinate of the point on the plane, in wire units.
Definition: PlaneGeo.h:704
WireID NearestWireID(Point_t const &pos) const
Returns the ID of wire closest to the specified position.
Definition: PlaneGeo.cxx:485
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:364
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:194
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).
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:195
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:61
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
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:286
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
intermediate_table::const_iterator const_iterator
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:290
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
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:265
unsigned int const fNfitpass
Definition: EMShowerAlg.h:269
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.
WireGeo const & Wire(WireID const &wireid) const
Returns the specified wire.
std::vector< double > const fToler
Definition: EMShowerAlg.h:271
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:287
void CheckIsolatedHits_(std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:155
std::vector< unsigned int > const fNfithits
Definition: EMShowerAlg.h:270
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>="">
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
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.
if(nlines<=0)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Signal from induction planes.
Definition: geo_types.h:147
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
static constexpr auto first()
Definition: geo_types.h:328
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
Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
double const fMinTrackLength
Definition: EMShowerAlg.h:264
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:306
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
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:2671
bool const fMakeGradientPlot
Definition: EMShowerAlg.h:286
double value
Definition: spectrum.C:18
EMShowerAlg(fhicl::ParameterSet const &pset, int debug=0)
Definition: EMShowerAlg.cxx:39
TVector2 ShowerDirection_(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &showerHits) const
Detector simulation of raw signals on wires.
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two 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:226
Contains all timing reference information for the detector.
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
constexpr ParentID_t const & parentID() const
Return the parent ID of this one (a plane ID).
Definition: geo_types.h:459
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:76
int const fNumShowerSegments
Definition: EMShowerAlg.h:288
calo::CalorimetryAlg const fCalorimetryAlg
Definition: EMShowerAlg.h:280
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:37
static constexpr WireIDIntersection invalid()
Definition: geo_types.h:594
double const fSpacePointSize
Definition: EMShowerAlg.h:266
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:471
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
Double_t sum
Definition: plot.C:31
art::ServiceHandle< geo::Geometry const > fGeom
Definition: EMShowerAlg.h:274
Float_t e
Definition: plot.C:35
size_t size() const
Definition: PmaTrack3D.h:89
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
std::string const fDetector
Definition: EMShowerAlg.h:283
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 WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
Definition: PlaneGeo.h:312
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