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" 44 ,
fNfitpass{pset.get<
unsigned int>(
"Nfitpass")}
45 ,
fNfithits{pset.get<std::vector<unsigned int>>(
"Nfithits")}
46 ,
fToler{pset.get<std::vector<double>>(
"Toler")}
50 ,
fDetector{pset.get<std::string>(
"Detector",
"dune35t")}
57 <<
"EMShowerAlg: fNfithits and fToler need to have size fNfitpass";
65 std::map<
int, std::vector<int>>& clusterToTracks,
66 std::map<
int, std::vector<int>>& trackToClusters)
const 68 std::vector<int> clustersToIgnore = {-999};
70 clusters, fmh, fmt, clustersToIgnore, clusterToTracks, trackToClusters);
77 std::vector<int>
const& clustersToIgnore,
78 std::map<
int, std::vector<int>>& clusterToTracks,
79 std::map<
int, std::vector<int>>& trackToClusters)
const 82 for (
auto const& clusterPtr : clusters) {
85 auto const& clusterHits = fmh.at(clusterPtr.key());
88 for (
auto const& clusterHitPtr : clusterHits) {
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";
96 if (clusterHitTracks.size() < 1)
continue;
98 auto const& clusterHitTrackPtr = clusterHitTracks[0];
101 std::cout <<
"Track " << clusterHitTrackPtr->ID() <<
" is too short! (" 102 << clusterHitTrackPtr->Length() <<
")\n";
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);
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);
126 if (firstTPC.size() != 2)
return;
129 else if (firstTPC.size() > 2)
135 int problemPlane = -1;
137 if (planes.size() == 1) problemPlane = planes[0];
140 if (showerHitsMap.at(problemPlane).size() < 3)
return;
143 std::vector<int> otherPlanes;
145 if (plane != problemPlane and showerHitsMap.count(plane) and
146 showerHitsMap.at(plane).size() >= 3)
147 otherPlanes.push_back(plane);
149 if (otherPlanes.size() == 0)
return;
152 unsigned int wrongTPC = showerHitsMap.at(problemPlane).at(0)->WireID().TPC;
153 unsigned int nHits = 0;
155 hitIt != showerHitsMap.at(problemPlane).end() and (*hitIt)->WireID().TPC == wrongTPC;
160 if (nHits > 2)
return;
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;
169 ++otherTPCs[(*hitIt)->WireID().TPC];
171 if (otherTPCs.size() > 1)
return;
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);
181 ++tpcCount[(*hitIt)->WireID().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;
189 hitIt != std::next(showerHitsMap.at(problemPlane).begin(), nHits);
191 naughty_hits.push_back(*hitIt);
192 showerHitsMap.at(problemPlane).erase(hitIt);
194 for (
auto const& naughty_hit : naughty_hits)
195 showerHitsMap.at(problemPlane).push_back(naughty_hit);
203 bool consistencyCheck =
true;
205 if (showerHitsMap.size() < 2) { consistencyCheck =
true; }
206 else if (showerHitsMap.size() == 2) {
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);
219 auto const showerStartPos =
Construct3DPoint_(detProp, startHits.at(0), startHits.at(1));
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()) /
228 if (timingDifference > 40 or projectionDifference > 5 or showerStartPos.X() == -9999 or
229 showerStartPos.Y() == -9999 or showerStartPos.Z() == -9999)
230 consistencyCheck =
false;
233 std::cout <<
"Timing difference is " << timingDifference <<
" and projection distance is " 234 << projectionDifference <<
" (start is (" << showerStartPos.X() <<
", " 235 << showerStartPos.Y() <<
", " << showerStartPos.Z() <<
")\n";
237 else if (showerHitsMap.size() == 3) {
243 std::map<int, art::Ptr<recob::Hit>> start2DMap;
244 for (
auto const& [plane,
hits] : showerHitsMap) {
245 start2DMap[plane] =
hits.front();
248 std::map<int, double> projDiff;
249 std::map<int, double> timingDiff;
251 for (
int plane = 0; plane < 3; ++plane) {
253 std::vector<int> otherPlanes;
254 for (
int otherPlane = 0; otherPlane < 3; ++otherPlane)
255 if (otherPlane != plane) otherPlanes.push_back(otherPlane);
258 detProp, start2DMap.at(otherPlanes.at(0)), start2DMap.at(otherPlanes.at(1)));
262 std::cout <<
"Plane... " << plane;
263 std::cout <<
"\nStart position in this plane is " 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";
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()));
279 std::cout <<
"Plane " << plane <<
" has projDiff " << projDiff <<
" and timeDiff " 282 if (projDiff > 5 or timeDiff > 40) consistencyCheck =
false;
286 if (
fDebug > 1) std::cout <<
"Consistency check is " << consistencyCheck <<
'\n';
288 return consistencyCheck;
292 std::vector<std::vector<int>>
const& initialShowers,
296 std::vector<int> clustersToIgnore;
299 for (
auto initialShowerIt = initialShowers.cbegin(); initialShowerIt != initialShowers.cend();
302 if (std::distance(initialShowers.cbegin(), initialShowerIt) > 0)
continue;
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) {
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);
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) +
325 auto const&
hits = fmh.at(clusterPtr.key());
326 for (
auto const&
hit : hits | views::transform(
to_element)) {
327 chargeDist[plane]->Fill(
hit.Integral());
330 chargeDist[plane]->Write();
337 if (planeClusters.size() < 3)
continue;
340 std::map<int, std::vector<double>> planeClusterSizes;
342 planeClusters.begin();
343 planeClustersIt != planeClusters.end();
346 planeClustersIt->second.begin();
347 planeClusterIt != planeClustersIt->second.end();
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());
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;
365 for (
auto const [plane, avg] : planeClustersAvSizes) {
368 std::vector<double> otherAverages;
369 for (
auto const [other_plane, other_avg] : planeClustersAvSizes)
370 if (plane != other_plane) otherAverages.push_back(other_avg);
372 double const sumSquareAvgsOtherPlanes = accumulate(
373 otherAverages, 0., [](
double sum,
double avg) {
return sum + cet::square(avg); });
374 double const quadOtherPlanes = std::sqrt(sumSquareAvgsOtherPlanes);
378 if (avg >= quadOtherPlanes) badPlane = plane;
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());
388 return clustersToIgnore;
403 return {
x, intersection.y, intersection.z};
410 std::map<int, TVector2>
const& showerCentreMap)
const 412 std::unique_ptr<recob::Track>
track;
414 std::vector<art::Ptr<recob::Hit>> track1, track2;
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.";
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) {
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 " 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);
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)) {
450 <<
") (real wire " <<
hit.WireID().Wire <<
") in TPC " <<
hit.WireID().TPC
462 mf::LogInfo(
"EMShowerAlg") <<
"Skipping this event because not enough hits in two views";
466 std::vector<TVector3> xyz, dircos;
468 for (
unsigned int i = 0; i < pmatrack->size(); i++) {
470 xyz.push_back((*pmatrack)[i]->
Point3D());
472 if (i < pmatrack->
size() - 1) {
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();
484 else if (!dircos.empty())
486 dircos.push_back(dc);
489 dircos.push_back(dircos.back());
493 std::map<int, double> distanceToVertex, distanceToEnd;
501 showerCentreIt != showerCentreMap.end();
509 distanceToVertex[showerCentreIt->first] = (vertexProj - showerCentreIt->second).Mod();
510 distanceToEnd[showerCentreIt->first] = (endProj - showerCentreIt->second).Mod();
514 double avDistanceToVertex = 0, avDistanceToEnd = 0;
516 distanceToVertexIt != distanceToVertex.end();
517 ++distanceToVertexIt)
518 avDistanceToVertex += distanceToVertexIt->second;
519 avDistanceToVertex /= distanceToVertex.size();
522 distanceToEndIt != distanceToEnd.end();
524 avDistanceToEnd += distanceToEndIt->second;
525 if (distanceToEnd.size() != 0) avDistanceToEnd /= distanceToEnd.size();
528 std::cout <<
"Distance to vertex is " << avDistanceToVertex <<
" and distance to end is " 529 << avDistanceToEnd <<
'\n';
532 if (avDistanceToEnd > avDistanceToVertex) {
533 std::reverse(xyz.begin(), xyz.end());
535 dircos.begin(), dircos.end(), dircos.begin(), [](TVector3
const& vec) {
return -1 * vec; });
538 if (xyz.size() != dircos.size())
539 mf::LogError(
"EMShowerAlg") <<
"Problem converting pma::Track3D to recob::Track";
541 track = std::make_unique<recob::Track>(
561 std::map<int, TVector2> showerCentreMap;
568 std::unique_ptr<recob::Track>
const&
track)
const 570 assert(not
empty(trackHits));
571 if (!track)
return -999;
573 recob::Hit const& firstHit = *trackHits.front();
588 double totalCharge = 0, totalDistance = 0, avHitTime = 0;
589 unsigned int nHits = 0;
591 for (
auto const&
hit : trackHits | views::transform(
to_element)) {
593 totalDistance += pitch;
594 totalCharge +=
hit.Integral();
595 avHitTime +=
hit.PeakTime();
600 avHitTime /= (double)nHits;
602 double const dQdx = totalCharge / totalDistance;
609 std::unique_ptr<recob::Track>& initialTrack,
619 std::cout <<
" ------------------ Finding initial track hits " 620 "-------------------- \n";
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)) {
633 std::cout <<
" ------------------ End finding initial track hits " 634 "-------------------- \n";
637 if (
fDebug > 1) std::cout <<
" ------------------ Finding initial track -------------------- \n";
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() <<
", " 647 std::cout <<
" ------------------ End finding initial track " 648 "-------------------- \n";
654 bool perpendicular)
const 662 if (perpendicular) direction = direction.Rotate(TMath::Pi() / 2);
666 std::map<double, art::Ptr<recob::Hit>> hitProjection;
667 for (
auto const& hitPtr :
hits) {
669 hitProjection[direction * pos] = hitPtr;
673 std::vector<art::Ptr<recob::Hit>> showerHits;
675 hitProjection, std::back_inserter(showerHits), [](
auto const& pr) {
return pr.second; });
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(
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);
693 TCanvas* canvas =
new TCanvas();
694 multigraph->Draw(
"AP");
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");
710 std::map<
int, std::vector<int>>
const& trackToClusters)
const 713 std::vector<std::vector<int>> showers;
716 for (
auto const& clusters : trackToClusters |
views::values) {
719 std::vector<int> matchingShowers;
721 for (
int const cluster : clusters) {
723 not cet::search_all(matchingShowers, shower)) {
724 matchingShowers.push_back(shower);
735 if (matchingShowers.size() < 1) showers.push_back(clusters);
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);
753 std::map<int, std::vector<art::Ptr<recob::Hit>>> initialHitsMap;
755 for (
auto const& [plane, orderedShower] : orderedShowerMap) {
756 std::vector<art::Ptr<recob::Hit>> initialHits;
759 bool wireDirection =
true;
760 std::vector<int> wires;
761 for (
auto const&
hit : orderedShower | views::transform(
to_element))
764 cet::sort_all(wires);
767 wireDirection =
false;
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);
778 for (
auto const& hitPtr : orderedShower) {
780 if (wireHitMap[wire].
size() > 1) {
782 if (multipleWires > 5)
break;
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))) {
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
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);
801 initialHits.push_back(hitPtr);
803 if (!initialHits.size()) initialHits.push_back(*orderedShower.begin());
810 std::map<int, std::vector<art::Ptr<recob::Hit>>> tickHitMap;
812 hitIt != orderedShower.end();
814 tickHitMap[std::round(
HitCoordinates_(**hitIt).Y())].push_back(*hitIt);
816 for (
auto const& hitPtr : orderedShower) {
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())))
826 initialHits.push_back(hitPtr);
828 if (initialHits.empty()) initialHits.push_back(*orderedShower.begin());
832 if (initialHits.size() == 1 and orderedShower.size() > 2)
833 initialHits.push_back(orderedShower[1]);
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];
842 for (
auto const [tpc, count] : tpcHitMap)
843 if (count == 1) singleHitTPCs.push_back(tpc);
845 if (singleHitTPCs.size()) {
847 for (
int const tpc : singleHitTPCs)
848 std::cout <<
"Removed hits in TPC " << tpc <<
'\n';
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());
856 newInitialHits = initialHits;
858 initialHitsMap[plane] = newInitialHits;
861 return initialHitsMap;
870 std::map<int, std::map<int, bool>> permutations;
880 std::map<int, double> planeRMSGradients, planeRMS;
881 for (
auto const& [plane, hitPtrs] : showerHitsMap) {
887 std::map<double, int> gradientMap;
888 for (
int const plane : showerHitsMap | views::keys)
889 gradientMap[
std::abs(planeRMSGradients.at(plane))] = plane;
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';
900 std::vector<std::map<int, bool>> usedPermutations;
903 for (
int const plane : showerHitsMap | views::keys)
904 permutations[perm][plane] = 0;
909 std::map<int, bool> permutation;
910 permutation[plane] =
true;
912 if (plane != plane2) permutation[plane2] =
false;
914 if (not cet::search_all(usedPermutations, permutation)) {
915 permutations[perm] = permutation;
916 usedPermutations.push_back(permutation);
920 for (
int const plane : wireWidthMap | views::reverse |
views::values) {
921 std::map<int, bool> permutation;
922 permutation[plane] =
false;
924 if (plane != plane2) permutation[plane2] =
true;
926 if (not cet::search_all(usedPermutations, permutation)) {
927 permutations[perm] = permutation;
928 usedPermutations.push_back(permutation);
934 for (
int const plane : showerHitsMap | views::keys)
935 permutations[perm][plane] = 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';
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>();
961 std::vector<std::pair<int, int>> initialHitsSize;
963 initialHitsMap.begin();
964 initialHitIt != initialHitsMap.end();
966 initialHitsSize.push_back(std::make_pair(initialHitIt->first, initialHitIt->second.size()));
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;
980 std::vector<int> planes;
982 if (initialHitsSize.size() > 2 and !
CheckShowerHits_(detProp, showerHitsMap)) {
985 std::cout <<
"Igoring plane " << planeToIgnore <<
" in creation of initial track\n";
987 hitsSizeIt != initialHitsSize.end() and planes.size() < 2;
989 if (hitsSizeIt->first == planeToIgnore)
continue;
990 planes.push_back(hitsSizeIt->first);
994 planes = {initialHitsSize.at(0).first, initialHitsSize.at(1).first};
996 return ConstructTrack(detProp, initialHitsMap.at(planes.at(0)), initialHitsMap.at(planes.at(1)));
1003 std::unique_ptr<recob::Track>
const& initialTrack,
1008 std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
1009 for (
auto const& hitPtr : hits)
1010 planeHitsMap[hitPtr->View()].push_back(hitPtr);
1013 unsigned int highestNumberOfHits = 0;
1014 std::vector<double> totalEnergy, totalEnergyError,
dEdx, dEdxError;
1020 if (planeHitsMap.count(plane)) {
1021 dEdx.push_back(
FinddEdx_(clockData, detProp, initialHitsMap.at(plane), initialTrack));
1022 totalEnergy.push_back(
1024 if (planeHitsMap.at(plane).size() > highestNumberOfHits and initialHitsMap.count(plane)) {
1026 highestNumberOfHits = planeHitsMap.at(plane).size();
1033 totalEnergy.push_back(0);
1037 TVector3 direction, directionError, showerStart, showerStartError;
1040 showerStart = initialTrack->
Vertex<TVector3>();
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";
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);
1078 std::vector<std::vector<art::Ptr<recob::Hit>>> initialTrackHits(3);
1082 unsigned maxhits0 = 0;
1083 unsigned maxhits1 = 0;
1086 planeHits != planeHitsMap.end();
1089 std::vector<art::Ptr<recob::Hit>> showerHits;
1092 if ((planeHits->second).size() > maxhits0) {
1094 maxhits1 = maxhits0;
1097 pl0 = planeHits->first;
1098 maxhits0 = (planeHits->second).
size();
1100 else if ((planeHits->second).size() > maxhits1) {
1101 pl1 = planeHits->first;
1102 maxhits1 = (planeHits->second).
size();
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) {
1113 std::vector<TVector3> spts;
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);
1121 if (spts.size() >= 2) {
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()) {
1131 if (spts.size() - 1 < 5) i = spts.size() - 1;
1132 shwdir = spts[i] - spts[0];
1133 shwdir = shwdir.Unit();
1136 shwxyz = spts.back();
1138 if (spts.size() > 6) i = spts.size() - 6;
1139 shwdir = spts[i] - spts[spts.size() - 1];
1140 shwdir = shwdir.Unit();
1143 if (planeHitsMap.find(plane) != planeHitsMap.end()) {
1144 totalEnergy.push_back(
1148 totalEnergy.push_back(0);
1150 if (initialTrackHits[plane].
size()) {
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;
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) {
1600 if ((coord1 - coordvtx).Mod() < (coord0 - coordvtx).Mod()) {
1601 std::reverse(showerHits.begin(), showerHits.end());
1610 auto const& xyz = vertex->
position();
1615 std::map<geo::TPCID, unsigned int> tpcmap;
1616 unsigned maxhits = 0;
1617 for (
auto const&
hit : showerHits) {
1620 for (
auto const& t : tpcmap) {
1621 if (t.second > maxhits) {
1631 std::vector<double> wfit;
1632 std::vector<double> tfit;
1633 std::vector<double> cfit;
1635 for (
size_t i = 0; i <
fNfitpass; ++i) {
1638 unsigned int nhits = 0;
1639 for (
auto&
hit : showerHits) {
1640 if (
hit->WireID().TPC == tpc.
TPC) {
1643 (
std::abs((coord.Y() - (parm[0] + coord.X() * parm[1])) * cos(atan(parm[1]))) <
1648 wfit.push_back(coord.X());
1649 tfit.push_back(coord.Y());
1651 if (i == fNfitpass - 1) { trackHits.push_back(
hit); }
1656 if (i < fNfitpass - 1 && wfit.size()) {
1657 fitok =
WeightedFit(wfit.size(), &wfit[0], &tfit[0], &cfit[0], &parm[0]);
1678 TVector2
const& pos,
1687 double globalWire = -999;
1696 .WireCoordinate(wireCenter);
1704 unsigned int nwires =
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;
1714 <<
"Error when trying to find a global induction plane coordinate for TPC " << wireID.
TPC 1718 unsigned int nwires =
1721 int block = wireID.
TPC / 4;
1722 globalWire = (nwires * block) + wireID.
Wire;
1730 .WireCoordinate(wireCenter);
1742 std::map<int, int> planeWireLength;
1744 showerHitsMap.begin();
1745 showerHitsIt != showerHitsMap.end();
1747 planeWireLength[showerHitsIt->first] =
1752 std::map<int, double> planeOtherWireLengths;
1754 planeWireLengthIt != planeWireLength.end();
1755 ++planeWireLengthIt) {
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;
1766 std::map<double, int> wireWidthMap;
1768 showerHitsMap.begin();
1769 showerHitsIt != showerHitsMap.end();
1771 wireWidthMap[planeOtherWireLengths.at(showerHitsIt->first)] = showerHitsIt->first;
1773 return wireWidthMap;
1783 double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1785 hit != showerHits.end();
1789 weight = cet::square((*hit)->Integral());
1791 sumx += weight * pos.X();
1792 sumy += weight * pos.Y();
1793 sumx2 += weight * pos.X() * pos.X();
1794 sumxy += weight * pos.X() * pos.Y();
1796 double gradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1797 TVector2 direction = TVector2(1, gradient).Unit();
1807 TVector2 pos, chargePoint = TVector2(0, 0);
1808 double totalCharge = 0;
1810 hit != showerHits.end();
1813 chargePoint += (*hit)->Integral() * pos;
1814 totalCharge += (*hit)->Integral();
1816 TVector2 centre = chargePoint / totalCharge;
1828 std::vector<double> distanceToAxis;
1830 showerHitsIt != showerHits.end();
1832 TVector2
proj = (
HitPosition_(detProp, **showerHitsIt) - centre).Proj(direction) + centre;
1833 distanceToAxis.push_back((
HitPosition_(detProp, **showerHitsIt) - proj).Mod());
1835 double RMS = TMath::RMS(distanceToAxis.begin(), distanceToAxis.end());
1849 double lengthOfShower =
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;
1859 showerSegments[segment].push_back(hitPtr);
1860 segmentCharge[segment] +=
hit.Integral();
1863 TGraph* graph =
new TGraph();
1864 std::vector<std::pair<int, double>> binVsRMS;
1868 for (
auto const& [segment, hitPtrs] : showerSegments) {
1871 TVector2 meanPosition(0, 0);
1872 for (
auto const&
hit : hitPtrs | views::transform(
to_element))
1874 meanPosition /= (double)hitPtrs.size();
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());
1883 double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
1884 binVsRMS.emplace_back(segment, RMSBin);
1889 double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
1890 for (
auto const& [
bin, RMSBin] : binVsRMS) {
1893 sumx += weight *
bin;
1894 sumy += weight * RMSBin;
1895 sumx2 += weight * bin *
bin;
1896 sumxy += weight * bin * RMSBin;
1898 double RMSgradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
1901 TVector2 direction = TVector2(1, RMSgradient).Unit();
1902 TCanvas* canv =
new TCanvas();
1904 graph->GetXaxis()->SetTitle(
"Shower segment");
1905 graph->GetYaxis()->SetTitle(
"RMS of hit distribution");
1906 TVector2 centre = TVector2(graph->GetMean(1), graph->GetMean(2));
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");
1927 TVector2 wireTickPos{-999., -999.};
1955 std::map<int, int> planeWireLength;
1956 std::map<int, double> planeOtherWireLengths;
1957 for (
auto const& [plane,
hits] : showerHitsMap)
1958 planeWireLength[plane] =
1961 planeWireLengthIt != planeWireLength.end();
1962 ++planeWireLengthIt) {
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;
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
1979 std::map<double, int> wireWidthMap;
1980 for (
int const plane : showerHitsMap | views::keys) {
1981 double wireWidth = planeWireLength.at(plane);
1982 wireWidthMap[wireWidth] = plane;
1985 return wireWidthMap.begin()->second;
1992 Double_t* parm)
const 1996 Double_t sumx2 = 0.;
1998 Double_t sumxy = 0.;
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];
2015 if (sumx2 * sumw - sumx * sumx == 0.)
return 1;
2016 if (sumx2 - sumx * sumx / sumw == 0.)
return 1;
2018 parm[0] = (sumy * sumx2 - sumx * sumxy) / (sumx2 * sumw - sumx * sumx);
2019 parm[1] = (sumxy - sumx * sumy / sumw) / (sumx2 - sumx * sumx / sumw);
2021 eparm[0] = sumx2 * (sumx2 * sumw - sumx * sumx);
2022 eparm[1] = (sumx2 - sumx * sumx / sumw);
2024 if (eparm[0] < 0. || eparm[1] < 0.)
return 1;
2026 eparm[0] = sqrt(eparm[0]) / (sumx2 * sumw - sumx * sumx);
2027 eparm[1] = sqrt(eparm[1]) / (sumx2 - sumx * sumx / sumw);
2034 if (
hits.empty())
return false;
2035 if (
hits.size() > 2000)
return true;
2036 if (
hits.size() < 20)
return true;
2038 std::map<int, int> hitmap;
2042 if (nhits > 2) ++hitmap[
hit.WireID().Wire];
2043 if (nhits == 20)
break;
2045 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::WireReadoutGeom const * fWireReadoutGeom
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].
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
constexpr to_element_t to_element
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
pma::ProjectionMatchingAlg const fProjectionMatchingAlg
WireID suggestedWireID() const
Returns a better wire ID.
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.
WireID NearestWireID(Point_t const &pos) const
Returns the ID of wire closest to the specified position.
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.
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).
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.
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.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
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.
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.
WireGeo const & Wire(WireID const &wireid) const
Returns the specified wire.
std::vector< double > const fToler
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
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>="">
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
Collection of exceptions for Geometry system.
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
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
static constexpr auto first()
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
Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
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 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
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.
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.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
constexpr ParentID_t const & parentID() const
Return the parent ID of this one (a plane ID).
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
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.
static constexpr WireIDIntersection invalid()
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
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
TPCID_t TPC
Index of the TPC within its cryostat.
art::ServiceHandle< geo::Geometry const > fGeom
recob::tracking::Plane Plane
Utility functions to extract information from recob::Track
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 WirePitch() const
Return the wire pitch (in centimeters). It is assumed constant.
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.