11 #include "cetlib/container_algorithms.h" 12 #include "cetlib/pow.h" 27 #include "TMathBase.h" 28 #include "TMultiGraph.h" 31 #include "range/v3/numeric.hpp" 32 #include "range/v3/view.hpp" 42 ,
fNfitpass{pset.get<
unsigned int>(
"Nfitpass")}
43 ,
fNfithits{pset.get<std::vector<unsigned int>>(
"Nfithits")}
44 ,
fToler{pset.get<std::vector<double>>(
"Toler")}
48 ,
fDetector{pset.get<std::string>(
"Detector",
"dune35t")}
55 <<
"EMShowerAlg: fNfithits and fToler need to have size fNfitpass";
63 std::map<
int, std::vector<int>>& clusterToTracks,
64 std::map<
int, std::vector<int>>& trackToClusters)
const 66 std::vector<int> clustersToIgnore = {-999};
68 clusters, fmh, fmt, clustersToIgnore, clusterToTracks, trackToClusters);
75 std::vector<int>
const& clustersToIgnore,
76 std::map<
int, std::vector<int>>& clusterToTracks,
77 std::map<
int, std::vector<int>>& trackToClusters)
const 80 for (
auto const& clusterPtr : clusters) {
83 auto const& clusterHits = fmh.at(clusterPtr.key());
86 for (
auto const& clusterHitPtr : clusterHits) {
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";
94 if (clusterHitTracks.size() < 1)
continue;
96 auto const& clusterHitTrackPtr = clusterHitTracks[0];
99 std::cout <<
"Track " << clusterHitTrackPtr->ID() <<
" is too short! (" 100 << clusterHitTrackPtr->Length() <<
")\n";
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);
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);
124 if (firstTPC.size() != 2)
return;
127 else if (firstTPC.size() > 2)
133 int problemPlane = -1;
135 if (planes.size() == 1) problemPlane = planes[0];
138 if (showerHitsMap.at(problemPlane).size() < 3)
return;
141 std::vector<int> otherPlanes;
143 if (plane != problemPlane and showerHitsMap.count(plane) and
144 showerHitsMap.at(plane).size() >= 3)
145 otherPlanes.push_back(plane);
147 if (otherPlanes.size() == 0)
return;
150 unsigned int wrongTPC = showerHitsMap.at(problemPlane).at(0)->WireID().TPC;
151 unsigned int nHits = 0;
153 hitIt != showerHitsMap.at(problemPlane).end() and (*hitIt)->WireID().TPC == wrongTPC;
158 if (nHits > 2)
return;
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;
167 ++otherTPCs[(*hitIt)->WireID().TPC];
169 if (otherTPCs.size() > 1)
return;
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);
179 ++tpcCount[(*hitIt)->WireID().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;
187 hitIt != std::next(showerHitsMap.at(problemPlane).begin(), nHits);
189 naughty_hits.push_back(*hitIt);
190 showerHitsMap.at(problemPlane).erase(hitIt);
192 for (
auto const& naughty_hit : naughty_hits)
193 showerHitsMap.at(problemPlane).push_back(naughty_hit);
201 bool consistencyCheck =
true;
203 if (showerHitsMap.size() < 2) { consistencyCheck =
true; }
204 else if (showerHitsMap.size() == 2) {
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);
217 auto const showerStartPos =
Construct3DPoint_(detProp, startHits.at(0), startHits.at(1));
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()) /
226 if (timingDifference > 40 or projectionDifference > 5 or showerStartPos.X() == -9999 or
227 showerStartPos.Y() == -9999 or showerStartPos.Z() == -9999)
228 consistencyCheck =
false;
231 std::cout <<
"Timing difference is " << timingDifference <<
" and projection distance is " 232 << projectionDifference <<
" (start is (" << showerStartPos.X() <<
", " 233 << showerStartPos.Y() <<
", " << showerStartPos.Z() <<
")\n";
235 else if (showerHitsMap.size() == 3) {
241 std::map<int, art::Ptr<recob::Hit>> start2DMap;
242 for (
auto const& [plane,
hits] : showerHitsMap) {
243 start2DMap[plane] =
hits.front();
246 std::map<int, double> projDiff;
247 std::map<int, double> timingDiff;
249 for (
int plane = 0; plane < 3; ++plane) {
251 std::vector<int> otherPlanes;
252 for (
int otherPlane = 0; otherPlane < 3; ++otherPlane)
253 if (otherPlane != plane) otherPlanes.push_back(otherPlane);
256 detProp, start2DMap.at(otherPlanes.at(0)), start2DMap.at(otherPlanes.at(1)));
260 std::cout <<
"Plane... " << plane;
261 std::cout <<
"\nStart position in this plane is " 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";
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()));
277 std::cout <<
"Plane " << plane <<
" has projDiff " << projDiff <<
" and timeDiff " 280 if (projDiff > 5 or timeDiff > 40) consistencyCheck =
false;
284 if (
fDebug > 1) std::cout <<
"Consistency check is " << consistencyCheck <<
'\n';
286 return consistencyCheck;
290 std::vector<std::vector<int>>
const& initialShowers,
294 std::vector<int> clustersToIgnore;
297 for (
auto initialShowerIt = initialShowers.cbegin(); initialShowerIt != initialShowers.cend();
300 if (std::distance(initialShowers.cbegin(), initialShowerIt) > 0)
continue;
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) {
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);
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) +
323 auto const&
hits = fmh.at(clusterPtr.key());
324 for (
auto const&
hit : hits | views::transform(
to_element)) {
325 chargeDist[plane]->Fill(
hit.Integral());
328 chargeDist[plane]->Write();
335 if (planeClusters.size() < 3)
continue;
338 std::map<int, std::vector<double>> planeClusterSizes;
340 planeClusters.begin();
341 planeClustersIt != planeClusters.end();
344 planeClustersIt->second.begin();
345 planeClusterIt != planeClustersIt->second.end();
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());
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;
363 for (
auto const [plane, avg] : planeClustersAvSizes) {
366 std::vector<double> otherAverages;
367 for (
auto const [other_plane, other_avg] : planeClustersAvSizes)
368 if (plane != other_plane) otherAverages.push_back(other_avg);
370 double const sumSquareAvgsOtherPlanes = accumulate(
371 otherAverages, 0., [](
double sum,
double avg) {
return sum + cet::square(avg); });
372 double const quadOtherPlanes = std::sqrt(sumSquareAvgsOtherPlanes);
376 if (avg >= quadOtherPlanes) badPlane = plane;
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());
386 return clustersToIgnore;
403 return {
x, intersection.
y, intersection.
z};
410 std::map<int, TVector2>
const& showerCentreMap)
const 413 std::unique_ptr<recob::Track>
track;
415 std::vector<art::Ptr<recob::Hit>> track1, track2;
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.";
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) {
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 " 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);
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)) {
451 <<
") (real wire " <<
hit.WireID().Wire <<
") in TPC " <<
hit.WireID().TPC
463 mf::LogInfo(
"EMShowerAlg") <<
"Skipping this event because not enough hits in two views";
467 std::vector<TVector3> xyz, dircos;
469 for (
unsigned int i = 0; i < pmatrack->size(); i++) {
471 xyz.push_back((*pmatrack)[i]->
Point3D());
473 if (i < pmatrack->
size() - 1) {
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();
485 else if (!dircos.empty())
487 dircos.push_back(dc);
490 dircos.push_back(dircos.back());
494 std::map<int, double> distanceToVertex, distanceToEnd;
502 showerCentreIt != showerCentreMap.end();
510 distanceToVertex[showerCentreIt->first] = (vertexProj - showerCentreIt->second).Mod();
511 distanceToEnd[showerCentreIt->first] = (endProj - showerCentreIt->second).Mod();
515 double avDistanceToVertex = 0, avDistanceToEnd = 0;
517 distanceToVertexIt != distanceToVertex.end();
518 ++distanceToVertexIt)
519 avDistanceToVertex += distanceToVertexIt->second;
520 avDistanceToVertex /= distanceToVertex.size();
523 distanceToEndIt != distanceToEnd.end();
525 avDistanceToEnd += distanceToEndIt->second;
526 if (distanceToEnd.size() != 0) avDistanceToEnd /= distanceToEnd.size();
529 std::cout <<
"Distance to vertex is " << avDistanceToVertex <<
" and distance to end is " 530 << avDistanceToEnd <<
'\n';
533 if (avDistanceToEnd > avDistanceToVertex) {
534 std::reverse(xyz.begin(), xyz.end());
536 dircos.begin(), dircos.end(), dircos.begin(), [](TVector3
const& vec) {
return -1 * vec; });
539 if (xyz.size() != dircos.size())
540 mf::LogError(
"EMShowerAlg") <<
"Problem converting pma::Track3D to recob::Track";
542 track = std::make_unique<recob::Track>(
562 std::map<int, TVector2> showerCentreMap;
569 std::unique_ptr<recob::Track>
const&
track)
const 571 assert(not
empty(trackHits));
572 if (!track)
return -999;
574 recob::Hit const& firstHit = *trackHits.front();
589 double totalCharge = 0, totalDistance = 0, avHitTime = 0;
590 unsigned int nHits = 0;
592 for (
auto const&
hit : trackHits | views::transform(
to_element)) {
594 totalDistance += pitch;
595 totalCharge +=
hit.Integral();
596 avHitTime +=
hit.PeakTime();
601 avHitTime /= (double)nHits;
603 double const dQdx = totalCharge / totalDistance;
610 std::unique_ptr<recob::Track>& initialTrack,
620 std::cout <<
" ------------------ Finding initial track hits " 621 "-------------------- \n";
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)) {
634 std::cout <<
" ------------------ End finding initial track hits " 635 "-------------------- \n";
638 if (
fDebug > 1) std::cout <<
" ------------------ Finding initial track -------------------- \n";
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() <<
", " 648 std::cout <<
" ------------------ End finding initial track " 649 "-------------------- \n";
655 bool perpendicular)
const 663 if (perpendicular) direction = direction.Rotate(TMath::Pi() / 2);
667 std::map<double, art::Ptr<recob::Hit>> hitProjection;
668 for (
auto const& hitPtr :
hits) {
670 hitProjection[direction * pos] = hitPtr;
674 std::vector<art::Ptr<recob::Hit>> showerHits;
676 hitProjection, std::back_inserter(showerHits), [](
auto const& pr) {
return pr.second; });
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(
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);
694 TCanvas* canvas =
new TCanvas();
695 multigraph->Draw(
"AP");
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");
711 std::map<
int, std::vector<int>>
const& trackToClusters)
const 714 std::vector<std::vector<int>> showers;
717 for (
auto const& clusters : trackToClusters |
views::values) {
720 std::vector<int> matchingShowers;
722 for (
int const cluster : clusters) {
724 not cet::search_all(matchingShowers, shower)) {
725 matchingShowers.push_back(shower);
736 if (matchingShowers.size() < 1) showers.push_back(clusters);
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);
754 std::map<int, std::vector<art::Ptr<recob::Hit>>> initialHitsMap;
756 for (
auto const& [plane, orderedShower] : orderedShowerMap) {
757 std::vector<art::Ptr<recob::Hit>> initialHits;
760 bool wireDirection =
true;
761 std::vector<int> wires;
762 for (
auto const&
hit : orderedShower | views::transform(
to_element))
765 cet::sort_all(wires);
768 wireDirection =
false;
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);
779 for (
auto const& hitPtr : orderedShower) {
781 if (wireHitMap[wire].
size() > 1) {
783 if (multipleWires > 5)
break;
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))) {
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
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);
802 initialHits.push_back(hitPtr);
804 if (!initialHits.size()) initialHits.push_back(*orderedShower.begin());
811 std::map<int, std::vector<art::Ptr<recob::Hit>>> tickHitMap;
813 hitIt != orderedShower.end();
815 tickHitMap[std::round(
HitCoordinates_(**hitIt).Y())].push_back(*hitIt);
817 for (
auto const& hitPtr : orderedShower) {
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())))
827 initialHits.push_back(hitPtr);
829 if (initialHits.empty()) initialHits.push_back(*orderedShower.begin());
833 if (initialHits.size() == 1 and orderedShower.size() > 2)
834 initialHits.push_back(orderedShower[1]);
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];
843 for (
auto const [tpc, count] : tpcHitMap)
844 if (count == 1) singleHitTPCs.push_back(tpc);
846 if (singleHitTPCs.size()) {
848 for (
int const tpc : singleHitTPCs)
849 std::cout <<
"Removed hits in TPC " << tpc <<
'\n';
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());
857 newInitialHits = initialHits;
859 initialHitsMap[plane] = newInitialHits;
862 return initialHitsMap;
871 std::map<int, std::map<int, bool>> permutations;
881 std::map<int, double> planeRMSGradients, planeRMS;
882 for (
auto const& [plane, hitPtrs] : showerHitsMap) {
888 std::map<double, int> gradientMap;
889 for (
int const plane : showerHitsMap | views::keys)
890 gradientMap[
std::abs(planeRMSGradients.at(plane))] = plane;
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';
901 std::vector<std::map<int, bool>> usedPermutations;
904 for (
int const plane : showerHitsMap | views::keys)
905 permutations[perm][plane] = 0;
910 std::map<int, bool> permutation;
911 permutation[plane] =
true;
913 if (plane != plane2) permutation[plane2] =
false;
915 if (not cet::search_all(usedPermutations, permutation)) {
916 permutations[perm] = permutation;
917 usedPermutations.push_back(permutation);
921 for (
int const plane : wireWidthMap | views::reverse |
views::values) {
922 std::map<int, bool> permutation;
923 permutation[plane] =
false;
925 if (plane != plane2) permutation[plane2] =
true;
927 if (not cet::search_all(usedPermutations, permutation)) {
928 permutations[perm] = permutation;
929 usedPermutations.push_back(permutation);
935 for (
int const plane : showerHitsMap | views::keys)
936 permutations[perm][plane] = 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';
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>();
962 std::vector<std::pair<int, int>> initialHitsSize;
964 initialHitsMap.begin();
965 initialHitIt != initialHitsMap.end();
967 initialHitsSize.push_back(std::make_pair(initialHitIt->first, initialHitIt->second.size()));
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;
981 std::vector<int> planes;
983 if (initialHitsSize.size() > 2 and !
CheckShowerHits_(detProp, showerHitsMap)) {
986 std::cout <<
"Igoring plane " << planeToIgnore <<
" in creation of initial track\n";
988 hitsSizeIt != initialHitsSize.end() and planes.size() < 2;
990 if (hitsSizeIt->first == planeToIgnore)
continue;
991 planes.push_back(hitsSizeIt->first);
995 planes = {initialHitsSize.at(0).first, initialHitsSize.at(1).first};
997 return ConstructTrack(detProp, initialHitsMap.at(planes.at(0)), initialHitsMap.at(planes.at(1)));
1004 std::unique_ptr<recob::Track>
const& initialTrack,
1009 std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
1010 for (
auto const& hitPtr : hits)
1011 planeHitsMap[hitPtr->View()].push_back(hitPtr);
1014 unsigned int highestNumberOfHits = 0;
1015 std::vector<double> totalEnergy, totalEnergyError,
dEdx, dEdxError;
1018 for (
unsigned int plane = 0; plane <
fGeom->
MaxPlanes(); ++plane) {
1021 if (planeHitsMap.count(plane)) {
1022 dEdx.push_back(
FinddEdx_(clockData, detProp, initialHitsMap.at(plane), initialTrack));
1023 totalEnergy.push_back(
1025 if (planeHitsMap.at(plane).size() > highestNumberOfHits and initialHitsMap.count(plane)) {
1027 highestNumberOfHits = planeHitsMap.at(plane).size();
1034 totalEnergy.push_back(0);
1038 TVector3 direction, directionError, showerStart, showerStartError;
1041 showerStart = initialTrack->
Vertex<TVector3>();
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";
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);
1079 std::vector<std::vector<art::Ptr<recob::Hit>>> initialTrackHits(3);
1083 unsigned maxhits0 = 0;
1084 unsigned maxhits1 = 0;
1087 planeHits != planeHitsMap.end();
1090 std::vector<art::Ptr<recob::Hit>> showerHits;
1093 if ((planeHits->second).size() > maxhits0) {
1095 maxhits1 = maxhits0;
1098 pl0 = planeHits->first;
1099 maxhits0 = (planeHits->second).
size();
1101 else if ((planeHits->second).size() > maxhits1) {
1102 pl1 = planeHits->first;
1103 maxhits1 = (planeHits->second).
size();
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) {
1114 std::vector<TVector3> spts;
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);
1122 if (spts.size() >= 2) {
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()) {
1132 if (spts.size() - 1 < 5) i = spts.size() - 1;
1133 shwdir = spts[i] - spts[0];
1134 shwdir = shwdir.Unit();
1137 shwxyz = spts.back();
1139 if (spts.size() > 6) i = spts.size() - 6;
1140 shwdir = spts[i] - spts[spts.size() - 1];
1141 shwdir = shwdir.Unit();
1143 for (
unsigned int plane = 0; plane <
fGeom->
MaxPlanes(); ++plane) {
1144 if (planeHitsMap.find(plane) != planeHitsMap.end()) {
1145 totalEnergy.push_back(
1149 totalEnergy.push_back(0);
1151 if (initialTrackHits[plane].
size()) {
1157 double angleToVert =
1159 initialTrackHits[plane][0]->WireID().planeID()) -
1161 double cosgamma =
std::abs(sin(angleToVert) * shwdir.Y() + cos(angleToVert) * shwdir.Z());
1162 if (cosgamma > 0) pitch = wirepitch / cosgamma;
1164 if (pitch < minpitch) {
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;
1174 vQ.push_back(
hit->Integral());
1175 totQ +=
hit->Integral();
1176 avgT +=
hit->PeakTime();
1181 double dQdx = TMath::Median(vQ.size(), &vQ[0]) / pitch;
1183 clockData, detProp, dQdx, avgT / nhits, initialTrackHits[plane][0]->
WireID().
Plane);
1186 dEdx.push_back(fdEdx);
1194 std::cout <<
"Best plane is " << bestPlane;
1195 std::cout <<
"\ndE/dx for each plane is: " << dEdx[0] <<
", " << dEdx[1] <<
" and " 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";
1224 std::vector<recob::SpacePoint> spacePoints;
1238 std::vector<int> usedHits;
1243 showerHitIt != showerHits.end();
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);
1253 if (otherPlanes.size() == 0)
return spacePoints;
1257 planeHitIt != showerHitIt->second.end();
1260 if (std::find(usedHits.begin(), usedHits.end(), planeHitIt->key()) != usedHits.end())
1264 const std::vector<art::Ptr<recob::Hit>> otherPlaneHits = showerHits.at(otherPlanes.at(0));
1266 otherPlaneHits.begin();
1267 otherPlaneHitIt != otherPlaneHits.end() and
1268 std::find(usedHits.begin(), usedHits.end(), planeHitIt->key()) == usedHits.end();
1269 ++otherPlaneHitIt) {
1271 if ((*otherPlaneHitIt)->WireID().TPC != (*planeHitIt)->WireID().TPC or
1272 std::find(usedHits.begin(), usedHits.end(), otherPlaneHitIt->key()) != usedHits.end())
1276 std::vector<art::Ptr<recob::Hit>> pointHits;
1277 bool truePoint =
false;
1279 if (otherPlanes.size() > 1) {
1282 const std::vector<art::Ptr<recob::Hit>> otherOtherPlaneHits =
1283 showerHits.at(otherPlanes.at(1));
1286 otherOtherPlaneHits.begin();
1287 otherOtherPlaneHitIt != otherOtherPlaneHits.end() and !truePoint;
1288 ++otherOtherPlaneHitIt) {
1290 if ((*otherOtherPlaneHitIt)->WireID().TPC == (*planeHitIt)->WireID().TPC and
1291 (projThirdPlane -
HitPosition_(detProp, **otherOtherPlaneHitIt)).Mod() <
1297 usedHits.push_back(planeHitIt->key());
1298 usedHits.push_back(otherPlaneHitIt->key());
1299 usedHits.push_back(otherOtherPlaneHitIt->key());
1301 pointHits.push_back(*planeHitIt);
1302 pointHits.push_back(*otherPlaneHitIt);
1303 pointHits.push_back(*otherOtherPlaneHitIt);
1317 usedHits.push_back(planeHitIt->key());
1318 usedHits.push_back(otherPlaneHitIt->key());
1320 pointHits.push_back(*planeHitIt);
1321 pointHits.push_back(*otherPlaneHitIt);
1326 double xyz[3] = {point.X(), point.Y(), point.Z()};
1334 spacePoints.emplace_back(xyz, xyzerr, chisq);
1335 hitAssns.push_back(pointHits);
1345 std::cout <<
"-------------------- Space points -------------------\n";
1346 std::cout <<
"There are " << spacePoints.size() <<
" space points:\n";
1349 spacePointIt != spacePoints.end();
1351 const double* xyz = spacePointIt->XYZ();
1352 std::cout <<
" Space point (" << xyz[0] <<
", " << xyz[1] <<
", " << xyz[2] <<
")\n";
1362 int desired_plane)
const 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);
1378 std::map<int, double> planeRMSGradients, planeRMS;
1379 for (
auto const& [plane,
hits] : showerHitsMap) {
1380 if (desired_plane != plane and desired_plane != -1)
continue;
1384 showerHitsMap[plane] = orderedHits;
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';
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)) {
1400 <<
") -- real wire " <<
hit.WireID() <<
", hit position (" 1411 std::map<int, double> planeOtherRMS, planeOtherRMSGradients;
1414 planeOtherRMS[planeRMSIt->first] = 0;
1415 planeOtherRMSGradients[planeRMSIt->first] = 0;
1416 int nOtherPlanes = 0;
1418 if (plane != planeRMSIt->first and planeRMS.count(plane)) {
1419 planeOtherRMS[planeRMSIt->first] += planeRMS.at(plane);
1420 planeOtherRMSGradients[planeRMSIt->first] += planeRMSGradients.at(plane);
1424 planeOtherRMS[planeRMSIt->first] /= (double)nOtherPlanes;
1425 planeOtherRMSGradients[planeRMSIt->first] /= (double)nOtherPlanes;
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 =
1435 showerHitsMap[plane] = orderedHits;
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()
1447 << hitPtrs.front()->WireID() <<
") and end (" 1450 << hitPtrs.back()->WireID() <<
")\n";
1456 showerHitsMap.begin();
1457 showerHitsIt != showerHitsMap.end();
1459 double gradient = planeRMSGradients.at(showerHitsIt->first);
1460 if (gradient < 0) std::reverse(showerHitsIt->second.begin(), showerHitsIt->second.end());
1464 std::cout <<
"\nHits in order; after reversing, before checking isolated hits...\n";
1466 showerHitsMap.begin();
1467 showerHitsIt != showerHitsMap.end();
1469 std::cout <<
" Plane " << showerHitsIt->first <<
'\n';
1471 hitIt != showerHitsIt->second.end();
1474 <<
HitCoordinates_(**hitIt).Y() <<
") -- real wire " << (*hitIt)->WireID()
1475 <<
", hit position (" <<
HitPosition_(detProp, **hitIt).X() <<
", " 1483 std::cout <<
"\nHits in order; after checking isolated hits...\n";
1485 showerHitsMap.begin();
1486 showerHitsIt != showerHitsMap.end();
1488 std::cout <<
" Plane " << showerHitsIt->first <<
'\n';
1490 hitIt != showerHitsIt->second.end();
1493 <<
HitCoordinates_(**hitIt).Y() <<
") -- real wire " << (*hitIt)->WireID()
1494 <<
", hit position (" <<
HitPosition_(detProp, **hitIt).X() <<
", " 1500 std::cout <<
"After reversing and checking isolated hits, here are the " 1501 "start and end points...\n";
1503 showerHitsMap.begin();
1504 showerHitsIt != showerHitsMap.end();
1506 std::cout <<
" Plane " << showerHitsIt->first <<
" has start (" 1509 << showerHitsIt->second.front()->WireID() <<
") and end (" 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) {
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);
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);
1538 auto const originalShowerHitsMap = showerHitsMap;
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;
1553 if (!
CheckShowerHits_(detProp, showerHitsMap)) showerHitsMap = originalShowerHitsMap;
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)) {
1568 if (ignoredPlanes.size())
1569 showerHitsMap[ignoredPlanes.begin()->first] = ignoredPlanes.begin()->second;
1571 return showerHitsMap;
1582 auto const& xyz = vertex->
position();
1588 for (
auto&
hit : showerHits) {
1589 if (
hit->WireID().TPC == tpc.
TPC) {
1599 if ((coord1 - coordvtx).Mod() < (coord0 - coordvtx).Mod()) {
1600 std::reverse(showerHits.begin(), showerHits.end());
1609 auto const& xyz = vertex->
position();
1614 std::map<geo::TPCID, unsigned int> tpcmap;
1615 unsigned maxhits = 0;
1616 for (
auto const&
hit : showerHits) {
1619 for (
auto const& t : tpcmap) {
1620 if (t.second > maxhits) {
1630 std::vector<double> wfit;
1631 std::vector<double> tfit;
1632 std::vector<double> cfit;
1634 for (
size_t i = 0; i <
fNfitpass; ++i) {
1637 unsigned int nhits = 0;
1638 for (
auto&
hit : showerHits) {
1639 if (
hit->WireID().TPC == tpc.
TPC) {
1642 (
std::abs((coord.Y() - (parm[0] + coord.X() * parm[1])) * cos(atan(parm[1]))) <
1647 wfit.push_back(coord.X());
1648 tfit.push_back(coord.Y());
1650 if (i == fNfitpass - 1) { trackHits.push_back(
hit); }
1655 if (i < fNfitpass - 1 && wfit.size()) {
1656 fitok =
WeightedFit(wfit.size(), &wfit[0], &tfit[0], &cfit[0], &parm[0]);
1677 TVector2
const& pos,
1685 double globalWire = -999;
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;
1710 <<
"Error when trying to find a global induction plane coordinate for TPC " << wireID.
TPC 1716 int block = wireID.
TPC / 4;
1717 globalWire = (nwires * block) + wireID.
Wire;
1736 std::map<int, int> planeWireLength;
1738 showerHitsMap.begin();
1739 showerHitsIt != showerHitsMap.end();
1741 planeWireLength[showerHitsIt->first] =
1746 std::map<int, double> planeOtherWireLengths;
1748 planeWireLengthIt != planeWireLength.end();
1749 ++planeWireLengthIt) {
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;
1760 std::map<double, int> wireWidthMap;
1762 showerHitsMap.begin();
1763 showerHitsIt != showerHitsMap.end();
1765 wireWidthMap[planeOtherWireLengths.at(showerHitsIt->first)] = showerHitsIt->first;
1767 return wireWidthMap;
1777 double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1779 hit != showerHits.end();
1783 weight = cet::square((*hit)->Integral());
1785 sumx += weight * pos.X();
1786 sumy += weight * pos.Y();
1787 sumx2 += weight * pos.X() * pos.X();
1788 sumxy += weight * pos.X() * pos.Y();
1790 double gradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1791 TVector2 direction = TVector2(1, gradient).Unit();
1801 TVector2 pos, chargePoint = TVector2(0, 0);
1802 double totalCharge = 0;
1804 hit != showerHits.end();
1807 chargePoint += (*hit)->Integral() * pos;
1808 totalCharge += (*hit)->Integral();
1810 TVector2 centre = chargePoint / totalCharge;
1822 std::vector<double> distanceToAxis;
1824 showerHitsIt != showerHits.end();
1826 TVector2
proj = (
HitPosition_(detProp, **showerHitsIt) - centre).Proj(direction) + centre;
1827 distanceToAxis.push_back((
HitPosition_(detProp, **showerHitsIt) - proj).Mod());
1829 double RMS = TMath::RMS(distanceToAxis.begin(), distanceToAxis.end());
1843 double lengthOfShower =
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;
1853 showerSegments[segment].push_back(hitPtr);
1854 segmentCharge[segment] +=
hit.Integral();
1857 TGraph* graph =
new TGraph();
1858 std::vector<std::pair<int, double>> binVsRMS;
1862 for (
auto const& [segment, hitPtrs] : showerSegments) {
1865 TVector2 meanPosition(0, 0);
1866 for (
auto const&
hit : hitPtrs | views::transform(
to_element))
1868 meanPosition /= (double)hitPtrs.size();
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());
1877 double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
1878 binVsRMS.emplace_back(segment, RMSBin);
1883 double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1884 for (
auto const& [
bin, RMSBin] : binVsRMS) {
1887 sumx += weight *
bin;
1888 sumy += weight * RMSBin;
1889 sumx2 += weight * bin *
bin;
1890 sumxy += weight * bin * RMSBin;
1892 double RMSgradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1895 TVector2 direction = TVector2(1, RMSgradient).Unit();
1896 TCanvas* canv =
new TCanvas();
1898 graph->GetXaxis()->SetTitle(
"Shower segment");
1899 graph->GetYaxis()->SetTitle(
"RMS of hit distribution");
1900 TVector2 centre = TVector2(graph->GetMean(1), graph->GetMean(2));
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");
1922 TVector2 wireTickPos = TVector2(-999., -999.);
1950 std::map<int, int> planeWireLength;
1951 std::map<int, double> planeOtherWireLengths;
1952 for (
auto const& [plane,
hits] : showerHitsMap)
1953 planeWireLength[plane] =
1956 planeWireLengthIt != planeWireLength.end();
1957 ++planeWireLengthIt) {
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;
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
1974 std::map<double, int> wireWidthMap;
1975 for (
int const plane : showerHitsMap | views::keys) {
1976 double wireWidth = planeWireLength.at(plane);
1977 wireWidthMap[wireWidth] = plane;
1980 return wireWidthMap.begin()->second;
1987 Double_t* parm)
const 1991 Double_t sumx2 = 0.;
1993 Double_t sumxy = 0.;
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];
2010 if (sumx2 * sumw - sumx * sumx == 0.)
return 1;
2011 if (sumx2 - sumx * sumx / sumw == 0.)
return 1;
2013 parm[0] = (sumy * sumx2 - sumx * sumxy) / (sumx2 * sumw - sumx * sumx);
2014 parm[1] = (sumxy - sumx * sumy / sumw) / (sumx2 - sumx * sumx / sumw);
2016 eparm[0] = sumx2 * (sumx2 * sumw - sumx * sumx);
2017 eparm[1] = (sumx2 - sumx * sumx / sumw);
2019 if (eparm[0] < 0. || eparm[1] < 0.)
return 1;
2021 eparm[0] = sqrt(eparm[0]) / (sumx2 * sumw - sumx * sumx);
2022 eparm[1] = sqrt(eparm[1]) / (sumx2 - sumx * sumx / sumw);
2029 if (
hits.empty())
return false;
2030 if (
hits.size() > 2000)
return true;
2031 if (
hits.size() < 20)
return true;
2033 std::map<int, int> hitmap;
2037 if (nhits > 2) ++hitmap[
hit.WireID().Wire];
2038 if (nhits == 20)
break;
2040 if (
float(nhits - 2) / hitmap.size() > 1.4)
void XYZ(double *xyz) const
Legacy method to access vertex position, preserved to avoid breaking code. Please try to use Vertex::...
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].
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
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
shower::ShowerEnergyAlg const fShowerEnergyAlg
constexpr to_element_t to_element
pma::ProjectionMatchingAlg const fProjectionMatchingAlg
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)
auto coord(Vector &v, unsigned int n) noexcept
Returns an object to manage the coordinate n of a vector.
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
TrackTrajectory::Flags_t Flags_t
::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.
CryostatID_t Cryostat
Index of cryostat.
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.
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.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Cluster finding and building.
double const fdEdxTrackLength
unsigned int const fNfitpass
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
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
bool const fMakeRMSGradientPlot
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.
std::vector< unsigned int > const fNfithits
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.
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.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Signal from induction planes.
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
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
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.
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
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.
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
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
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.
PlaneID_t Plane
Index of the plane within its TPC.
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)
bool const fMakeGradientPlot
EMShowerAlg(fhicl::ParameterSet const &pset, int debug=0)
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
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.
float PeakTime() const
Time of the signal peak, in tick units.
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
TVector2 HitCoordinates_(recob::Hit const &hit) const
Return the coordinates of this hit in global wire/tick space.
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].
int const fNumShowerSegments
calo::CalorimetryAlg const fCalorimetryAlg
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.
double const fSpacePointSize
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
constexpr PlaneID const & planeID() const
TPCID_t TPC
Index of the TPC within its cryostat.
art::ServiceHandle< geo::Geometry const > fGeom
geo::WireID suggestedWireID() const
Returns a better wire ID.
recob::tracking::Plane Plane
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
const Point_t & position() const
Return vertex 3D position.
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.
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.