#include "EMShowerAlg.h"

Public Member Functions
	EMShowerAlg (fhicl::ParameterSet const &pset, int debug=0)

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. More...

void	AssociateClustersAndTracks (std::vector< art::Ptr< recob::Cluster >> const &clusters, art::FindManyP< recob::Hit > const &fmh, art::FindManyP< recob::Track > const &fmt, std::vector< int > const &clustersToIgnore, 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, whilst ignoring certain clusters. More...

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. More...

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

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) const

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

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. More...

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. More...

recob::Shower	MakeShower (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &hits, art::Ptr< recob::Vertex > const &vertex, int &iok) const
	<Tingjun to="" document>=""> More...

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. More...

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. More...

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>=""> More...

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>=""> More...

bool	isCleanShower (std::vector< art::Ptr< recob::Hit >> const &hits) const
	<Tingjun to="" document>=""> More...

Public Attributes
int	fDebug

Private Member Functions
void	CheckIsolatedHits_ (std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const

bool	CheckShowerHits_ (detinfo::DetectorPropertiesData const &detProp, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &showerHitsMap) const

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. More...

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. More...

std::vector< art::Ptr< recob::Hit > >	FindOrderOfHits_ (detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &hits, bool perpendicular=false) const

std::map< int, std::vector< art::Ptr< recob::Hit > > >	FindShowerStart_ (std::map< int, std::vector< art::Ptr< recob::Hit >>> const &orderedShowerMap) const

std::map< int, std::map< int, bool > >	GetPlanePermutations_ (const detinfo::DetectorPropertiesData &detProp, const std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const

double	GlobalWire_ (const geo::WireID &wireID) const
	Find the global wire position. More...

TVector2	HitCoordinates_ (recob::Hit const &hit) const
	Return the coordinates of this hit in global wire/tick space. More...

TVector2	HitPosition_ (detinfo::DetectorPropertiesData const &detProp, recob::Hit const &hit) const
	Return the coordinates of this hit in units of cm. More...

TVector2	HitPosition_ (detinfo::DetectorPropertiesData const &detProp, TVector2 const &pos, geo::PlaneID planeID) const
	Return the coordinates of this hit in units of cm. More...

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

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

TVector2	Project3DPointOntoPlane_ (detinfo::DetectorPropertiesData const &detProp, geo::Point_t const &point, int plane, int cryostat=0) const

std::map< double, int >	RelativeWireWidth_ (const std::map< int, std::vector< art::Ptr< recob::Hit >>> &showerHitsMap) const

TVector2	ShowerCenter_ (detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &showerHits) const
	Returns the charge-weighted shower center. More...

TVector2	ShowerDirection_ (detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &showerHits) const

double	ShowerHitRMS_ (detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &showerHits) const

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. More...

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. More...

Private Attributes
double const	fMinTrackLength

double const	fdEdxTrackLength

double const	fSpacePointSize

unsigned int const	fNfitpass

std::vector< unsigned int > const	fNfithits

std::vector< double > const	fToler

art::ServiceHandle< geo::Geometry const >	fGeom

shower::ShowerEnergyAlg const	fShowerEnergyAlg

calo::CalorimetryAlg const	fCalorimetryAlg

pma::ProjectionMatchingAlg const	fProjectionMatchingAlg

std::string const	fDetector

bool const	fMakeGradientPlot

bool const	fMakeRMSGradientPlot

int const	fNumShowerSegments

Detailed Description

Definition at line 60 of file EMShowerAlg.h.

Constructor & Destructor Documentation

shower::EMShowerAlg::EMShowerAlg	(	fhicl::ParameterSet const &	pset,
		int	debug = `0`
	)

Definition at line 37 of file EMShowerAlg.cxx.

References art::errors::Configuration, fCalorimetryAlg, fdEdxTrackLength, fDetector, fMakeGradientPlot, fMakeRMSGradientPlot, fMinTrackLength, fNfithits, fNfitpass, fNumShowerSegments, fProjectionMatchingAlg, fShowerEnergyAlg, fSpacePointSize, and fToler.

   : fDebug{debug}
   , fMinTrackLength{pset.get<double>("MinTrackLength")}
   , fdEdxTrackLength{pset.get<double>("dEdxTrackLength")}
   , fSpacePointSize{pset.get<double>("SpacePointSize")}
   , fNfitpass{pset.get<unsigned int>("Nfitpass")}
   , fNfithits{pset.get<std::vector<unsigned int>>("Nfithits")}
   , fToler{pset.get<std::vector<double>>("Toler")}
   , fShowerEnergyAlg(pset.get<fhicl::ParameterSet>("ShowerEnergyAlg"))
   , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
   , fProjectionMatchingAlg(pset.get<fhicl::ParameterSet>("ProjectionMatchingAlg"))
   , fDetector{pset.get<std::string>("Detector", "dune35t")} // tmp
   , fMakeGradientPlot{pset.get<bool>("MakeGradientPlot", false)}
   , fMakeRMSGradientPlot{pset.get<bool>("MakeRMSGradientPlot", false)}
   , fNumShowerSegments{pset.get<int>("NumShowerSegments", 4)}
 {
   if (fNfitpass != fNfithits.size() || fNfitpass != fToler.size()) {
     throw art::Exception(art::errors::Configuration)
       << "EMShowerAlg: fNfithits and fToler need to have size fNfitpass";
   }
 }

Member Function Documentation

void shower::EMShowerAlg::AssociateClustersAndTracks	(	std::vector< art::Ptr< recob::Cluster >> const &	clusters,
		art::FindManyP< recob::Hit > const &	fmh,
		art::FindManyP< recob::Track > const &	fmt,
		std::map< int, std::vector< int >> &	clusterToTracks,
		std::map< int, std::vector< int >> &	trackToClusters
	)		const

Map associated tracks and clusters together given their associated hits.

Definition at line 59 of file EMShowerAlg.cxx.

Referenced by shower::EMShower::produce().

 {
   std::vector<int> clustersToIgnore = {-999};
   AssociateClustersAndTracks(
     clusters, fmh, fmt, clustersToIgnore, clusterToTracks, trackToClusters);
 }

void shower::EMShowerAlg::AssociateClustersAndTracks	(	std::vector< art::Ptr< recob::Cluster >> const &	clusters,
		art::FindManyP< recob::Hit > const &	fmh,
		art::FindManyP< recob::Track > const &	fmt,
		std::vector< int > const &	clustersToIgnore,
		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, whilst ignoring certain clusters.

Definition at line 71 of file EMShowerAlg.cxx.

References fDebug, fMinTrackLength, and track.

 {
   // Look through all the clusters
   for (auto const& clusterPtr : clusters) {
 
     // Get the hits in this cluster
     auto const& clusterHits = fmh.at(clusterPtr.key());
 
     // Look at all these hits and find the associated tracks
     for (auto const& clusterHitPtr : clusterHits) {
       // Get the tracks associated with this hit
       auto const& clusterHitTracks = fmt.at(clusterHitPtr.key());
       if (clusterHitTracks.size() > 1) {
         std::cout << "More than one track associated with this hit!\n";
         continue;
       }
 
       if (clusterHitTracks.size() < 1) continue;
 
       auto const& clusterHitTrackPtr = clusterHitTracks[0];
       if (clusterHitTrackPtr->Length() < fMinTrackLength) {
         if (fDebug > 2)
           std::cout << "Track " << clusterHitTrackPtr->ID() << " is too short! ("
                     << clusterHitTrackPtr->Length() << ")\n";
         continue;
       }
 
       // Add this cluster to the track map
       int track = clusterHitTrackPtr.key();
       int cluster = clusterPtr.key();
       if (cet::search_all(clustersToIgnore, cluster)) continue;
       if (not cet::search_all(trackToClusters[track], cluster))
         trackToClusters[track].push_back(cluster);
       if (not cet::search_all(clusterToTracks[cluster], track))
         clusterToTracks[cluster].push_back(track);
     }
   }
 }

void shower::EMShowerAlg::CheckIsolatedHits_ ( std::map< int, std::vector< art::Ptr< recob::Hit >>> & showerHitsMap ) const

private

Checks the hits across the views in a given shower to determine if there is one in the incorrect TPC

Definition at line 116 of file EMShowerAlg.cxx.

References fGeom, hits(), geo::GeometryCore::MaxPlanes(), util::values(), and lar::dump::vector().

Referenced by OrderShowerHits().

 {
   std::map<int, std::vector<int>> firstTPC;
   for (auto const& [plane, hits] : showerHitsMap)
     firstTPC[hits.at(0)->WireID().TPC].push_back(plane);
 
   // If all in the same TPC then that's great!
   if (firstTPC.size() != 2) return;
 
   // If they are in more than two TPCs, not much we can do
   else if (firstTPC.size() > 2)
     return;
 
   // If we get to this point, there should be something we can do!
 
   // Find the problem plane
   int problemPlane = -1;
   for (auto const& planes : firstTPC | views::values)
     if (planes.size() == 1) problemPlane = planes[0];
 
   // Require three hits
   if (showerHitsMap.at(problemPlane).size() < 3) return;
 
   // and get the other planes with at least three hits
   std::vector<int> otherPlanes;
   for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane)
     if (plane != problemPlane and showerHitsMap.count(plane) and
         showerHitsMap.at(plane).size() >= 3)
       otherPlanes.push_back(plane);
 
   if (otherPlanes.size() == 0) return;
 
   // Look at the hits after the first one
   unsigned int wrongTPC = showerHitsMap.at(problemPlane).at(0)->WireID().TPC;
   unsigned int nHits = 0;
   for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsMap.at(problemPlane).begin();
        hitIt != showerHitsMap.at(problemPlane).end() and (*hitIt)->WireID().TPC == wrongTPC;
        ++hitIt)
     ++nHits;
 
   // If there are more than two hits in the 'wrong TPC', we can't be sure it is indeed wrong
   if (nHits > 2) return;
 
   // See if at least the next four times as many hits are in a different TPC
   std::map<int, int> otherTPCs;
   for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt =
          std::next(showerHitsMap.at(problemPlane).begin(), nHits);
        hitIt != showerHitsMap.at(problemPlane).end() and
        std::distance(std::next(showerHitsMap.at(problemPlane).begin(), nHits), hitIt) < 4 * nHits;
        ++hitIt)
     ++otherTPCs[(*hitIt)->WireID().TPC];
 
   if (otherTPCs.size() > 1) return;
 
   // If we get this far, we can move the problem hits from the front of the shower to the back
   std::map<int, int> tpcCount;
   for (int const otherPlane : otherPlanes)
     for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt =
            std::next(showerHitsMap.at(otherPlane).begin());
          hitIt != showerHitsMap.at(otherPlane).end() and
          hitIt != std::next(showerHitsMap.at(otherPlane).begin(), 2);
          ++hitIt)
       ++tpcCount[(*hitIt)->WireID().TPC];
 
   // Remove the first hit if it is in the wrong TPC
   if (tpcCount.size() == 1 and
       tpcCount.begin()->first ==
         (int)(*std::next(showerHitsMap.at(problemPlane).begin(), nHits))->WireID().TPC) {
     std::vector<art::Ptr<recob::Hit>> naughty_hits;
     for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsMap.at(problemPlane).begin();
          hitIt != std::next(showerHitsMap.at(problemPlane).begin(), nHits);
          ++hitIt) {
       naughty_hits.push_back(*hitIt);
       showerHitsMap.at(problemPlane).erase(hitIt);
     }
     for (auto const& naughty_hit : naughty_hits)
       showerHitsMap.at(problemPlane).push_back(naughty_hit);
   }
 }

bool shower::EMShowerAlg::CheckShowerHits_	(	detinfo::DetectorPropertiesData const &	detProp,
		std::map< int, std::vector< art::Ptr< recob::Hit >>> const &	showerHitsMap
	)		const

private

Takes the shower hits in all views and ensure the ordering is consistent Returns bool, indicating whether or not everything makes sense!

Definition at line 197 of file EMShowerAlg.cxx.

References util::abs(), Construct3DPoint_(), fDebug, HitPosition_(), hits(), Project3DPointOntoPlane_(), X, and Y.

Referenced by MakeInitialTrack_(), and OrderShowerHits().

 {
   bool consistencyCheck = true;
 
   if (showerHitsMap.size() < 2) { consistencyCheck = true; }
   else if (showerHitsMap.size() == 2) {
 
     // With two views, we can check:
     //  -- timing between views is consistent
     //  -- the 3D start point makes sense when projected back onto the individual planes
 
     std::vector<art::Ptr<recob::Hit>> startHits;
     std::vector<int> planes;
     for (auto const& [plane, hits] : showerHitsMap) {
       startHits.push_back(hits.front());
       planes.push_back(plane);
     }
 
     auto const showerStartPos = Construct3DPoint_(detProp, startHits.at(0), startHits.at(1));
     TVector2 proj1 = Project3DPointOntoPlane_(detProp, showerStartPos, planes.at(0));
     TVector2 proj2 = Project3DPointOntoPlane_(detProp, showerStartPos, planes.at(1));
 
     double timingDifference = std::abs(startHits.at(0)->PeakTime() - startHits.at(1)->PeakTime());
     double projectionDifference = ((HitPosition_(detProp, *startHits.at(0)) - proj1).Mod() +
                                    (HitPosition_(detProp, *startHits.at(1)) - proj2).Mod()) /
                                   (double)2;
 
     if (timingDifference > 40 or projectionDifference > 5 or showerStartPos.X() == -9999 or
         showerStartPos.Y() == -9999 or showerStartPos.Z() == -9999)
       consistencyCheck = false;
 
     if (fDebug > 1)
       std::cout << "Timing difference is " << timingDifference << " and projection distance is "
                 << projectionDifference << " (start is (" << showerStartPos.X() << ", "
                 << showerStartPos.Y() << ", " << showerStartPos.Z() << ")\n";
   }
   else if (showerHitsMap.size() == 3) {
 
     // With three views, we can check:
     //  -- the timing between views is consistent
     //  -- the 3D start point formed by two views and projected back into the third is close to the start point in that view
 
     std::map<int, art::Ptr<recob::Hit>> start2DMap;
     for (auto const& [plane, hits] : showerHitsMap) {
       start2DMap[plane] = hits.front();
     }
 
     std::map<int, double> projDiff;
     std::map<int, double> timingDiff;
 
     for (int plane = 0; plane < 3; ++plane) {
 
       std::vector<int> otherPlanes;
       for (int otherPlane = 0; otherPlane < 3; ++otherPlane)
         if (otherPlane != plane) otherPlanes.push_back(otherPlane);
 
       auto const showerStartPos = Construct3DPoint_(
         detProp, start2DMap.at(otherPlanes.at(0)), start2DMap.at(otherPlanes.at(1)));
       TVector2 showerStartProj = Project3DPointOntoPlane_(detProp, showerStartPos, plane);
 
       if (fDebug > 2) {
         std::cout << "Plane... " << plane;
         std::cout << "\nStart position in this plane is "
                   << HitPosition_(detProp, *start2DMap.at(plane)).X() << ", "
                   << HitPosition_(detProp, *start2DMap.at(plane)).Y() << ")\n";
         std::cout << "Shower start from other two planes is (" << showerStartPos.X() << ", "
                   << showerStartPos.Y() << ", " << showerStartPos.Z() << ")\n";
         std::cout << "Projecting the other two planes gives position (" << showerStartProj.X()
                   << ", " << showerStartProj.Y() << ")\n";
       }
 
       double projDiff =
         std::abs((showerStartProj - HitPosition_(detProp, *start2DMap.at(plane))).Mod());
       double timeDiff = TMath::Max(
         std::abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(0))->PeakTime()),
         std::abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(1))->PeakTime()));
 
       if (fDebug > 1)
         std::cout << "Plane " << plane << " has projDiff " << projDiff << " and timeDiff "
                   << timeDiff << '\n';
 
       if (projDiff > 5 or timeDiff > 40) consistencyCheck = false;
     }
   }
 
   if (fDebug > 1) std::cout << "Consistency check is " << consistencyCheck << '\n';
 
   return consistencyCheck;
 }

std::vector< int > shower::EMShowerAlg::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.

Definition at line 289 of file EMShowerAlg.cxx.

References fDebug, hits(), art::Ptr< T >::key(), geo::PlaneID::Plane, recob::Cluster::Plane(), sum, lar::to_element, util::to_string(), and lar::dump::vector().

Referenced by shower::EMShower::produce().

 {
   std::vector<int> clustersToIgnore;
 
   // Look at each shower
   for (auto initialShowerIt = initialShowers.cbegin(); initialShowerIt != initialShowers.cend();
        ++initialShowerIt) {
 
     if (std::distance(initialShowers.cbegin(), initialShowerIt) > 0) continue;
 
     // Map the clusters and cluster hits to each view
     std::map<int, std::vector<art::Ptr<recob::Cluster>>> planeClusters;
     std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHits;
     for (int const clusterIndex : *initialShowerIt) {
       art::Ptr<recob::Cluster> const& cluster = clusters.at(clusterIndex);
       planeClusters[cluster->Plane().Plane].push_back(cluster);
       for (auto const& hit : fmh.at(cluster.key()))
         planeHits[hit->WireID().Plane].push_back(hit);
     }
 
     TFile* outFile = new TFile("chargeDistributions.root", "RECREATE");
     std::map<int, TH1D*> chargeDist;
     for (auto const& [plane, clusterPtrs] : planeClusters) {
       for (auto const& clusterPtr : clusterPtrs) {
         chargeDist[plane] = new TH1D(std::string("chargeDist_Plane" + std::to_string(plane) +
                                                  "_Cluster" + std::to_string(clusterPtr.key()))
                                        .c_str(),
                                      "",
                                      150,
                                      0,
                                      1000);
         auto const& hits = fmh.at(clusterPtr.key());
         for (auto const& hit : hits | views::transform(to_element)) {
           chargeDist[plane]->Fill(hit.Integral());
         }
         outFile->cd();
         chargeDist[plane]->Write();
       }
     }
     outFile->Close();
     delete outFile;
 
     // Can't do much with fewer than three views
     if (planeClusters.size() < 3) continue;
 
     // Look at how many clusters each plane has, and the proportion of hits each one uses
     std::map<int, std::vector<double>> planeClusterSizes;
     for (std::map<int, std::vector<art::Ptr<recob::Cluster>>>::iterator planeClustersIt =
            planeClusters.begin();
          planeClustersIt != planeClusters.end();
          ++planeClustersIt) {
       for (std::vector<art::Ptr<recob::Cluster>>::iterator planeClusterIt =
              planeClustersIt->second.begin();
            planeClusterIt != planeClustersIt->second.end();
            ++planeClusterIt) {
         std::vector<art::Ptr<recob::Hit>> hits = fmh.at(planeClusterIt->key());
         planeClusterSizes[planeClustersIt->first].push_back(
           (double)hits.size() / (double)planeHits.at(planeClustersIt->first).size());
       }
     }
 
     // Find the average hit fraction across all clusters in the plane
     std::map<int, double> planeClustersAvSizes;
     for (auto const& [plane, cluster_sizes] : planeClusterSizes) {
       double const average = accumulate(cluster_sizes, 0.) / cluster_sizes.size();
       planeClustersAvSizes[plane] = average;
     }
 
     // Now decide if there is one plane which is ruining the reconstruction
     // If two planes have a low average cluster fraction and one high, this plane likely merges two particle deposits together
     int badPlane = -1;
     for (auto const [plane, avg] : planeClustersAvSizes) {
 
       // Get averages from other planes and add in quadrature
       std::vector<double> otherAverages;
       for (auto const [other_plane, other_avg] : planeClustersAvSizes)
         if (plane != other_plane) otherAverages.push_back(other_avg);
 
       double const sumSquareAvgsOtherPlanes = accumulate(
         otherAverages, 0., [](double sum, double avg) { return sum + cet::square(avg); });
       double const quadOtherPlanes = std::sqrt(sumSquareAvgsOtherPlanes);
 
       // If this plane has an average higher than the quadratic sum of the
       // others, it may be bad
       if (avg >= quadOtherPlanes) badPlane = plane;
     }
 
     if (badPlane != -1) {
       if (fDebug > 1) std::cout << "Bad plane is " << badPlane << '\n';
       for (auto const& cluster : planeClusters.at(badPlane))
         clustersToIgnore.push_back(cluster.key());
     }
   }
 
   return clustersToIgnore;
 }

geo::Point_t shower::EMShowerAlg::Construct3DPoint_	(	detinfo::DetectorPropertiesData const &	detProp,
		art::Ptr< recob::Hit > const &	hit1,
		art::Ptr< recob::Hit > const &	hit2
	)		const

private

Constructs a 3D point (in [cm]) to represent the hits given in two views.

Definition at line 389 of file EMShowerAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), fGeom, recob::Hit::PeakTime(), geo::WireID::planeID(), recob::Hit::WireID(), geo::GeometryCore::WireIDsIntersect(), x, geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by CheckShowerHits_(), ConstructTrack(), and MakeSpacePoints().

 {
 
   // x is average of the two x's
   double x = (detProp.ConvertTicksToX(hit1->PeakTime(), hit1->WireID().planeID()) +
               detProp.ConvertTicksToX(hit2->PeakTime(), hit2->WireID().planeID())) /
              (double)2;
 
   // y and z got from the wire interections
   geo::WireIDIntersection intersection;
   fGeom->WireIDsIntersect(hit1->WireID(), hit2->WireID(), intersection);
 
   return {x, intersection.y, intersection.z};
 }

std::unique_ptr< recob::Track > shower::EMShowerAlg::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

Constructs a recob::Track from sets of hits in two views. Intended to be used to construct the initial first part of a shower. All PMA methods taken from the pma tracking algorithm (R. Sulej and D. Stefan). This implementation also orients the track in the correct direction if a map of shower centres (units [cm]) in each view is provided.

Definition at line 406 of file EMShowerAlg.cxx.

References pma::ProjectionMatchingAlg::buildSegment(), Construct3DPoint_(), recob::tracking::convertCollToPoint(), recob::tracking::convertCollToVector(), util::end(), fDebug, fProjectionMatchingAlg, HitCoordinates_(), if(), Project3DPointOntoPlane_(), util::size(), lar::to_element, geo::vect::toPoint(), and track.

Referenced by ConstructTrack(), and MakeInitialTrack_().

 {
 
   std::unique_ptr<recob::Track> track;
 
   std::vector<art::Ptr<recob::Hit>> track1, track2;
 
   // Check the TPCs
   if ((*hits1.begin())->WireID().TPC != (*hits2.begin())->WireID().TPC) {
     mf::LogWarning("EMShowerAlg") << "Warning: attempting to construct a track from two different "
                                      "TPCs.  Returning a null track.";
     return track;
   }
 
   // Check for tracks crossing TPC boundaries
   std::map<int, int> tpcMap;
   for (auto const& hit : hits1)
     ++tpcMap[hit->WireID().TPC];
   for (auto const& hit : hits2)
     ++tpcMap[hit->WireID().TPC];
   if (tpcMap.size() > 1) {
     mf::LogWarning("EMShowerAlg")
       << "Warning: attempting to construct a track which crosses more than one TPC -- PMTrack "
          "can't handle this right now.  Returning a track made just from hits in the first TPC it "
          "traverses.";
     unsigned int firstTPC1 = hits1.at(0)->WireID().TPC, firstTPC2 = hits2.at(0)->WireID().TPC;
     for (auto const& hit : hits1)
       if (hit->WireID().TPC == firstTPC1) track1.push_back(hit);
     for (auto const& hit : hits2)
       if (hit->WireID().TPC == firstTPC2) track2.push_back(hit);
   }
   else {
     track1 = hits1;
     track2 = hits2;
   }
 
   if (fDebug > 1) {
     std::cout << "About to make a track from these hits:\n";
     auto print_hits = [this](auto const& track) {
       for (auto const& hit : track | views::transform(to_element)) {
         std::cout << "Hit (" << HitCoordinates_(hit).X() << ", " << HitCoordinates_(hit).Y()
                   << ") (real wire " << hit.WireID().Wire << ") in TPC " << hit.WireID().TPC
                   << '\n';
       }
     };
     print_hits(track1);
     print_hits(track2);
   }
 
   auto const trackStart = Construct3DPoint_(detProp, track1.at(0), track2.at(0));
   pma::Track3D* pmatrack = fProjectionMatchingAlg.buildSegment(detProp, track1, track2, trackStart);
 
   if (!pmatrack) {
     mf::LogInfo("EMShowerAlg") << "Skipping this event because not enough hits in two views";
     return track;
   }
 
   std::vector<TVector3> xyz, dircos;
 
   for (unsigned int i = 0; i < pmatrack->size(); i++) {
 
     xyz.push_back((*pmatrack)[i]->Point3D());
 
     if (i < pmatrack->size() - 1) {
       size_t j = i + 1;
       double mag = 0.0;
       TVector3 dc(0., 0., 0.);
       while ((mag == 0.0) and (j < pmatrack->size())) {
         dc = (*pmatrack)[j]->Point3D();
         dc -= (*pmatrack)[i]->Point3D();
         mag = dc.Mag();
         ++j;
       }
       if (mag > 0.0)
         dc *= 1.0 / mag;
       else if (!dircos.empty())
         dc = dircos.back();
       dircos.push_back(dc);
     }
     else
       dircos.push_back(dircos.back());
   }
 
   // Orient the track correctly
   std::map<int, double> distanceToVertex, distanceToEnd;
   using geo::vect::toPoint;
   geo::Point_t const vertex = toPoint(*xyz.begin());
   geo::Point_t const end = toPoint(*xyz.rbegin());
 
   // Loop over all the planes and find the distance from the vertex and end
   // projections to the centre in each plane
   for (std::map<int, TVector2>::const_iterator showerCentreIt = showerCentreMap.begin();
        showerCentreIt != showerCentreMap.end();
        ++showerCentreIt) {
 
     // Project the vertex and the end point onto this plane
     TVector2 vertexProj = Project3DPointOntoPlane_(detProp, vertex, showerCentreIt->first);
     TVector2 endProj = Project3DPointOntoPlane_(detProp, end, showerCentreIt->first);
 
     // Find the distance of each to the centre of the cluster
     distanceToVertex[showerCentreIt->first] = (vertexProj - showerCentreIt->second).Mod();
     distanceToEnd[showerCentreIt->first] = (endProj - showerCentreIt->second).Mod();
   }
 
   // Find the average distance to the vertex and the end across the planes
   double avDistanceToVertex = 0, avDistanceToEnd = 0;
   for (std::map<int, double>::iterator distanceToVertexIt = distanceToVertex.begin();
        distanceToVertexIt != distanceToVertex.end();
        ++distanceToVertexIt)
     avDistanceToVertex += distanceToVertexIt->second;
   avDistanceToVertex /= distanceToVertex.size();
 
   for (std::map<int, double>::iterator distanceToEndIt = distanceToEnd.begin();
        distanceToEndIt != distanceToEnd.end();
        ++distanceToEndIt)
     avDistanceToEnd += distanceToEndIt->second;
   if (distanceToEnd.size() != 0) avDistanceToEnd /= distanceToEnd.size();
 
   if (fDebug > 2)
     std::cout << "Distance to vertex is " << avDistanceToVertex << " and distance to end is "
               << avDistanceToEnd << '\n';
 
   // Change order if necessary
   if (avDistanceToEnd > avDistanceToVertex) {
     std::reverse(xyz.begin(), xyz.end());
     std::transform(
       dircos.begin(), dircos.end(), dircos.begin(), [](TVector3 const& vec) { return -1 * vec; });
   }
 
   if (xyz.size() != dircos.size())
     mf::LogError("EMShowerAlg") << "Problem converting pma::Track3D to recob::Track";
 
   track = std::make_unique<recob::Track>(
     recob::TrackTrajectory(recob::tracking::convertCollToPoint(xyz),
                            recob::tracking::convertCollToVector(dircos),
                            recob::Track::Flags_t(xyz.size()),
                            false),
     0,
     -1.,
     0,
     recob::tracking::SMatrixSym55(),
     recob::tracking::SMatrixSym55(),
     -1);
 
   return track;
 }

std::unique_ptr< recob::Track > shower::EMShowerAlg::ConstructTrack	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::Hit >> const &	track1,
		std::vector< art::Ptr< recob::Hit >> const &	track2
	)		const

Constructs a recob::Track from sets of hits in two views. Intended to be used to construct the initial first part of a shower. All methods taken from the pma tracking algorithm (R. Sulej and D. Stefan).

Definition at line 557 of file EMShowerAlg.cxx.

References ConstructTrack().

 {
   std::map<int, TVector2> showerCentreMap;
   return ConstructTrack(detProp, track1, track2, showerCentreMap);
 }

double shower::EMShowerAlg::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

private

Finds dE/dx for the track given a set of hits.

Definition at line 566 of file EMShowerAlg.cxx.

References calo::CalorimetryAlg::dEdx_AREA(), util::empty(), fCalorimetryAlg, fdEdxTrackLength, geo::PlaneID::Plane, lar::to_element, lar::util::TrackPitchInView(), recob::Hit::View(), and recob::Hit::WireID().

Referenced by MakeShower().

 {
   assert(not empty(trackHits));
   if (!track) return -999;
 
   recob::Hit const& firstHit = *trackHits.front();
 
   // Get the pitch
   double pitch = 0;
   try {
     pitch = lar::util::TrackPitchInView(*track, firstHit.View());
   }
   catch (...) {
   }
 
   // Deal with large pitches
   if (pitch > fdEdxTrackLength) {
     return fCalorimetryAlg.dEdx_AREA(clockData, detProp, firstHit, pitch);
   }
 
   double totalCharge = 0, totalDistance = 0, avHitTime = 0;
   unsigned int nHits = 0;
 
   for (auto const& hit : trackHits | views::transform(to_element)) {
     if (totalDistance + pitch < fdEdxTrackLength) {
       totalDistance += pitch;
       totalCharge += hit.Integral();
       avHitTime += hit.PeakTime();
       ++nHits;
     }
   }
 
   avHitTime /= (double)nHits;
 
   double const dQdx = totalCharge / totalDistance;
   return fCalorimetryAlg.dEdx_AREA(clockData, detProp, dQdx, avHitTime, firstHit.WireID().Plane);
 }

void shower::EMShowerAlg::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

Finds the initial track-like part of the shower and the hits in all views associated with it

Finding the initial track requires three stages: – find the initial track-like hits in each view – use these to construct a track

Definition at line 607 of file EMShowerAlg.cxx.

References fDebug, FindShowerStart_(), HitCoordinates_(), MakeInitialTrack_(), lar::to_element, recob::Track::Vertex(), and recob::Track::VertexDirection().

Referenced by shower::EMShower::produce().

 {
 
 
   // Now find the hits belonging to the track
   if (fDebug > 1)
     std::cout << " ------------------ Finding initial track hits "
                  "-------------------- \n";
   initialTrackHits = FindShowerStart_(showerHitsMap);
   if (fDebug > 1) {
     std::cout << "Here are the initial shower hits... \n";
     for (auto const& [plane, hitPtrs] : initialTrackHits) {
       std::cout << "  Plane " << plane << '\n';
       for (auto const& hit : hitPtrs | views::transform(to_element)) {
         std::cout << "    Hit is (" << HitCoordinates_(hit).X() << " (real hit " << hit.WireID()
                   << "), " << HitCoordinates_(hit).Y() << ")\n";
       }
     }
   }
   if (fDebug > 1)
     std::cout << " ------------------ End finding initial track hits "
                  "-------------------- \n";
 
   // Now we have the track hits -- can make a track!
   if (fDebug > 1) std::cout << " ------------------ Finding initial track -------------------- \n";
   initialTrack = MakeInitialTrack_(detProp, initialTrackHits, showerHitsMap);
   if (initialTrack and fDebug > 1) {
     std::cout << "The track start is (" << initialTrack->Vertex().X() << ", "
               << initialTrack->Vertex().Y() << ", " << initialTrack->Vertex().Z() << ")\n";
     std::cout << "The track direction is (" << initialTrack->VertexDirection().X() << ", "
               << initialTrack->VertexDirection().Y() << ", " << initialTrack->VertexDirection().Z()
               << ")\n";
   }
   if (fDebug > 1)
     std::cout << " ------------------ End finding initial track "
                  "-------------------- \n";
 }

void shower::EMShowerAlg::FindInitialTrackHits	(	std::vector< art::Ptr< recob::Hit >> const &	showerHits,
		art::Ptr< recob::Vertex > const &	vertex,
		std::vector< art::Ptr< recob::Hit >> &	trackHits
	)		const

Definition at line 1604 of file EMShowerAlg.cxx.

References util::abs(), geo::vect::coord(), fGeom, geo::GeometryCore::FindTPCAtPosition(), fNfithits, fNfitpass, fToler, HitCoordinates_(), geo::CryostatID::isValid, recob::Vertex::position(), geo::TPCID::TPC, and WeightedFit().

Referenced by MakeShower().

 {
   // Find TPC for the vertex
   auto const& xyz = vertex->position();
   geo::TPCID tpc = fGeom->FindTPCAtPosition(xyz);
 
   // vertex cannot be projected into a TPC, find the TPC that has the most hits
   if (!tpc.isValid) {
     std::map<geo::TPCID, unsigned int> tpcmap;
     unsigned maxhits = 0;
     for (auto const& hit : showerHits) {
       ++tpcmap[geo::TPCID(hit->WireID())];
     }
     for (auto const& t : tpcmap) {
       if (t.second > maxhits) {
         maxhits = t.second;
         tpc = t.first;
       }
     }
   }
   if (!tpc.isValid) return;
 
   double parm[2];
   int fitok = 0;
   std::vector<double> wfit;
   std::vector<double> tfit;
   std::vector<double> cfit;
 
   for (size_t i = 0; i < fNfitpass; ++i) {
 
     // Fit a straight line through hits
     unsigned int nhits = 0;
     for (auto& hit : showerHits) {
       if (hit->WireID().TPC == tpc.TPC) {
         TVector2 coord = HitCoordinates_(*hit);
         if (i == 0 ||
             (std::abs((coord.Y() - (parm[0] + coord.X() * parm[1])) * cos(atan(parm[1]))) <
              fToler[i - 1]) ||
             fitok == 1) {
           ++nhits;
           if (nhits == fNfithits[i] + 1) break;
           wfit.push_back(coord.X());
           tfit.push_back(coord.Y());
           cfit.push_back(1.);
           if (i == fNfitpass - 1) { trackHits.push_back(hit); }
         }
       }
     }
 
     if (i < fNfitpass - 1 && wfit.size()) {
       fitok = WeightedFit(wfit.size(), &wfit[0], &tfit[0], &cfit[0], &parm[0]);
     }
     wfit.clear();
     tfit.clear();
     cfit.clear();
   }
 }

std::vector< art::Ptr< recob::Hit > > shower::EMShowerAlg::FindOrderOfHits_	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::Hit >> const &	hits,
		bool	perpendicular = `false`
	)		const

private

Orders hits along the best fit line through the charge-weighted centre of the hits. Orders along the line perpendicular to the least squares line if perpendicular is set to true.

Definition at line 652 of file EMShowerAlg.cxx.

References color(), fMakeGradientPlot, HitPosition_(), hits(), ShowerCenter_(), ShowerDirection_(), lar::to_element, X, and Y.

Referenced by OrderShowerHits(), and OrderShowerHits_().

 {
   // Find the charge-weighted centre (in [cm]) of this shower
   TVector2 centre = ShowerCenter_(detProp, hits);
 
   // Find a rough shower 'direction'
   TVector2 direction = ShowerDirection_(detProp, hits);
 
   if (perpendicular) direction = direction.Rotate(TMath::Pi() / 2);
 
   // Find how far each hit (projected onto this axis) is from the centre
   TVector2 pos;
   std::map<double, art::Ptr<recob::Hit>> hitProjection;
   for (auto const& hitPtr : hits) {
     pos = HitPosition_(detProp, *hitPtr) - centre;
     hitProjection[direction * pos] = hitPtr;
   }
 
   // Get a vector of hits in order of the shower
   std::vector<art::Ptr<recob::Hit>> showerHits;
   cet::transform_all(
     hitProjection, std::back_inserter(showerHits), [](auto const& pr) { return pr.second; });
 
   // Make gradient plot
   if (fMakeGradientPlot) {
     std::map<int, TGraph*> graphs;
     for (auto const& hit : showerHits | views::transform(to_element)) {
       int tpc = hit.WireID().TPC;
       if (graphs[tpc] == nullptr) graphs[tpc] = new TGraph();
       graphs[tpc]->SetPoint(
         graphs[tpc]->GetN(), HitPosition_(detProp, hit).X(), HitPosition_(detProp, hit).Y());
     }
     TMultiGraph* multigraph = new TMultiGraph();
     for (auto const [color, graph] : graphs) {
       graph->SetMarkerColor(color);
       graph->SetMarkerStyle(8);
       graph->SetMarkerSize(2);
       multigraph->Add(graph);
     }
     TCanvas* canvas = new TCanvas();
     multigraph->Draw("AP");
     TLine line;
     line.SetLineColor(2);
     line.DrawLine(centre.X() - 1000 * direction.X(),
                   centre.Y() - 1000 * direction.Y(),
                   centre.X() + 1000 * direction.X(),
                   centre.Y() + 1000 * direction.Y());
     canvas->SaveAs("Gradient.png");
     delete canvas;
     delete multigraph;
   }
 
   return showerHits;
 }

std::vector< std::vector< int > > shower::EMShowerAlg::FindShowers ( std::map< int, std::vector< int >> const & trackToClusters ) const

Makes showers given a map between tracks and all clusters associated with them.

Definition at line 710 of file EMShowerAlg.cxx.

References util::values().

Referenced by shower::EMShower::produce().

 {
   // Showers are vectors of clusters
   std::vector<std::vector<int>> showers;
 
   // Loop over all tracks
   for (auto const& clusters : trackToClusters | views::values) {
 
     // Find which showers already made are associated with this track
     std::vector<int> matchingShowers;
     for (unsigned int shower = 0; shower < showers.size(); ++shower)
       for (int const cluster : clusters) {
         if (cet::search_all(showers[shower], cluster) and
             not cet::search_all(matchingShowers, shower)) {
           matchingShowers.push_back(shower);
         }
       }
 
     // THINK THERE PROBABLY CAN BE MORE THAN ONE!
     // IN FACT, THIS WOULD BE A SUCCESS OF THE SHOWERING METHOD!
     // // Shouldn't be more than one
     // if (matchingShowers.size() > 1)
     //   mf::LogInfo("EMShowerAlg") << "Number of showers this track matches is " << matchingShowers.size() << std::endl;
 
     // New shower
     if (matchingShowers.size() < 1) showers.push_back(clusters);
 
     // Add to existing shower
     else {
       for (int const cluster : clusters) {
         if (not cet::search_all(showers.at(matchingShowers[0]), cluster))
           showers.at(matchingShowers.at(0)).push_back(cluster);
       }
     }
   }
 
   return showers;
 }

std::map< int, std::vector< art::Ptr< recob::Hit > > > shower::EMShowerAlg::FindShowerStart_ ( std::map< int, std::vector< art::Ptr< recob::Hit >>> const & orderedShowerMap ) const

private

Takes a map of the shower hits on each plane (ordered from what has been decided to be the start) Returns a map of the initial track-like part of the shower on each plane

Definition at line 750 of file EMShowerAlg.cxx.

References util::abs(), fDebug, HitCoordinates_(), util::size(), lar::to_element, lar::dump::vector(), X, and Y.

Referenced by FindInitialTrack().

 {
 
   std::map<int, std::vector<art::Ptr<recob::Hit>>> initialHitsMap;
 
   for (auto const& [plane, orderedShower] : orderedShowerMap) {
     std::vector<art::Ptr<recob::Hit>> initialHits;
 
     // Find if the shower is traveling along ticks or wires
     bool wireDirection = true;
     std::vector<int> wires;
     for (auto const& hit : orderedShower | views::transform(to_element))
       wires.push_back(std::round(HitCoordinates_(hit).X()));
 
     cet::sort_all(wires);
     if (std::abs(*wires.begin() - std::round(HitCoordinates_(**orderedShower.begin()).X())) > 3 and
         std::abs(*wires.rbegin() - std::round(HitCoordinates_(**orderedShower.begin()).X())) > 3)
       wireDirection = false;
 
     // Deal with showers traveling along wires
     if (wireDirection) {
       bool increasing = HitCoordinates_(**orderedShower.rbegin()).X() >
                         HitCoordinates_(**orderedShower.begin()).X();
       std::map<int, std::vector<art::Ptr<recob::Hit>>> wireHitMap;
       int multipleWires = 0;
       for (auto const& hitPtr : orderedShower)
         wireHitMap[std::round(HitCoordinates_(*hitPtr).X())].push_back(hitPtr);
 
       for (auto const& hitPtr : orderedShower) {
         int wire = std::round(HitCoordinates_(*hitPtr).X());
         if (wireHitMap[wire].size() > 1) {
           ++multipleWires;
           if (multipleWires > 5) break;
           continue;
         }
         else if ((increasing and wireHitMap[wire + 1].size() > 1 and
                   (wireHitMap[wire + 2].size() > 1 or wireHitMap[wire + 3].size() > 1)) or
                  (!increasing and wireHitMap[wire - 1].size() > 1 and
                   (wireHitMap[wire - 2].size() > 1 or wireHitMap[wire - 3].size() > 1))) {
           if ((increasing and
                (wireHitMap[wire + 4].size() < 2 and wireHitMap[wire + 5].size() < 2 and
                 wireHitMap[wire + 6].size() < 2 and wireHitMap[wire + 9].size() > 1)) or
               (!increasing and
                (wireHitMap[wire - 4].size() < 2 and wireHitMap[wire - 5].size() < 2 and
                 wireHitMap[wire - 6].size() < 2) and
                wireHitMap[wire - 9].size() > 1))
             initialHits.push_back(hitPtr);
           else
             break;
         }
         else
           initialHits.push_back(hitPtr);
       }
       if (!initialHits.size()) initialHits.push_back(*orderedShower.begin());
     }
 
     // Deal with showers travelling along ticks
     else {
       bool increasing = HitCoordinates_(**orderedShower.rbegin()).Y() >
                         HitCoordinates_(**orderedShower.begin()).Y();
       std::map<int, std::vector<art::Ptr<recob::Hit>>> tickHitMap;
       for (std::vector<art::Ptr<recob::Hit>>::const_iterator hitIt = orderedShower.begin();
            hitIt != orderedShower.end();
            ++hitIt)
         tickHitMap[std::round(HitCoordinates_(**hitIt).Y())].push_back(*hitIt);
 
       for (auto const& hitPtr : orderedShower) {
         int const tick = std::round(HitCoordinates_(*hitPtr).Y());
         if ((increasing and (tickHitMap[tick + 1].size() or tickHitMap[tick + 2].size() or
                              tickHitMap[tick + 3].size() or tickHitMap[tick + 4].size() or
                              tickHitMap[tick + 5].size())) or
             (!increasing and (tickHitMap[tick - 1].size() or tickHitMap[tick - 2].size() or
                               tickHitMap[tick - 3].size() or tickHitMap[tick - 4].size() or
                               tickHitMap[tick - 5].size())))
           break;
         else
           initialHits.push_back(hitPtr);
       }
       if (initialHits.empty()) initialHits.push_back(*orderedShower.begin());
     }
 
     // Need at least two hits to make a track
     if (initialHits.size() == 1 and orderedShower.size() > 2)
       initialHits.push_back(orderedShower[1]);
 
     // Quality check -- make sure there isn't a single hit in a different TPC (artefact of tracking failure)
     std::vector<art::Ptr<recob::Hit>> newInitialHits;
     std::map<int, int> tpcHitMap;
     std::vector<int> singleHitTPCs;
     for (auto const& hit : initialHits | views::transform(to_element))
       ++tpcHitMap[hit.WireID().TPC];
 
     for (auto const [tpc, count] : tpcHitMap)
       if (count == 1) singleHitTPCs.push_back(tpc);
 
     if (singleHitTPCs.size()) {
       if (fDebug > 2)
         for (int const tpc : singleHitTPCs)
           std::cout << "Removed hits in TPC " << tpc << '\n';
 
       for (auto const& hitPtr : initialHits)
         if (not cet::search_all(singleHitTPCs, hitPtr->WireID().TPC))
           newInitialHits.push_back(hitPtr);
       if (!newInitialHits.size()) newInitialHits.push_back(*orderedShower.begin());
     }
     else
       newInitialHits = initialHits;
 
     initialHitsMap[plane] = newInitialHits;
   }
 
   return initialHitsMap;
 }

std::map< int, std::map< int, bool > > shower::EMShowerAlg::GetPlanePermutations_	(	const detinfo::DetectorPropertiesData &	detProp,
		const std::map< int, std::vector< art::Ptr< recob::Hit >>> &	showerHitsMap
	)		const

private

Takes all the shower hits, ready ordered, and returns information to help with the orientation of the shower in each view Returns map of most likely permutations of reorientation Starts at 0,0,0 (i.e. don't need to reorient any plane) and ends with 1,1,1 (i.e. every plane needs reorienting) Every permutation inbetween represent increasing less likely changes to satisfy the correct orientation criteria

Definition at line 865 of file EMShowerAlg.cxx.

References util::abs(), fDebug, RelativeWireWidth_(), ShowerHitRMS_(), ShowerHitRMSGradient_(), value, and util::values().

Referenced by OrderShowerHits().

 {
 
   // The map to return
   std::map<int, std::map<int, bool>> permutations;
 
   // Get the properties of the shower hits across the planes which will be used to
   // determine the likelihood of a particular reorientation permutation
   //   -- relative width in the wire direction (if showers travel along the wire
   //      direction in a particular plane)
   //   -- the RMS gradients (how likely it is the RMS of the hit positions from
   //      central axis increases along a particular orientation)
 
   // Find the RMS, RMS gradient and wire widths
   std::map<int, double> planeRMSGradients, planeRMS;
   for (auto const& [plane, hitPtrs] : showerHitsMap) {
     planeRMS[plane] = ShowerHitRMS_(detProp, hitPtrs);
     planeRMSGradients[plane] = ShowerHitRMSGradient_(detProp, hitPtrs);
   }
 
   // Order these backwards so they can be used to discriminate between planes
   std::map<double, int> gradientMap;
   for (int const plane : showerHitsMap | views::keys)
     gradientMap[std::abs(planeRMSGradients.at(plane))] = plane;
 
   std::map<double, int> wireWidthMap = RelativeWireWidth_(showerHitsMap);
 
   if (fDebug > 1)
     for (auto const [gradient, plane] : wireWidthMap)
       std::cout << "Plane " << plane << " has relative width in wire of " << gradient
                 << "; and an RMS gradient of " << planeRMSGradients[plane] << '\n';
 
   // Find the correct permutations
   int perm = 0;
   std::vector<std::map<int, bool>> usedPermutations;
 
   // Most likely is to not change anything
   for (int const plane : showerHitsMap | views::keys)
     permutations[perm][plane] = 0;
   ++perm;
 
   // Use properties of the shower to determine the middle cases
   for (int const plane : wireWidthMap | views::values) {
     std::map<int, bool> permutation;
     permutation[plane] = true;
     for (int const plane2 : wireWidthMap | views::values)
       if (plane != plane2) permutation[plane2] = false;
 
     if (not cet::search_all(usedPermutations, permutation)) {
       permutations[perm] = permutation;
       usedPermutations.push_back(permutation);
       ++perm;
     }
   }
   for (int const plane : wireWidthMap | views::reverse | views::values) {
     std::map<int, bool> permutation;
     permutation[plane] = false;
     for (int const plane2 : wireWidthMap | views::values)
       if (plane != plane2) permutation[plane2] = true;
 
     if (not cet::search_all(usedPermutations, permutation)) {
       permutations[perm] = permutation;
       usedPermutations.push_back(permutation);
       ++perm;
     }
   }
 
   // Least likely is to change everything
   for (int const plane : showerHitsMap | views::keys)
     permutations[perm][plane] = 1;
   ++perm;
 
   if (fDebug > 1) {
     std::cout << "Here are the permutations!\n";
     for (auto const& [index, permutation] : permutations) {
       std::cout << "Permutation " << index << '\n';
       for (auto const [plane, value] : permutation)
         std::cout << "  Plane " << plane << " has value " << value << '\n';
     }
   }
 
   return permutations;
 }

double shower::EMShowerAlg::GlobalWire_ ( const geo::WireID & wireID ) const

private

Find the global wire position.

Definition at line 1683 of file EMShowerAlg.cxx.

References geo::CryostatID::Cryostat, fDetector, fGeom, geo::WireGeo::GetCenter(), geo::kInduction, geo::GeometryCore::Nwires(), geo::PlaneID::Plane, geo::GeometryCore::SignalType(), geo::TPCID::TPC, geo::WireID::Wire, geo::GeometryCore::WireCoordinate(), and geo::GeometryCore::WireIDToWireGeo().

Referenced by HitCoordinates_(), and Project3DPointOntoPlane_().

 {
   double globalWire = -999;
 
   // Induction
   if (fGeom->SignalType(wireID) == geo::kInduction) {
     auto const wireCenter = fGeom->WireIDToWireGeo(wireID).GetCenter();
     globalWire = fGeom->WireCoordinate(wireCenter,
                                        geo::PlaneID{wireID.Cryostat,
                                                     wireID.TPC % 2, // 0 or 1
                                                     wireID.Plane});
   }
 
   // Collection
   else {
     // FOR COLLECTION WIRES, HARD CODE THE GEOMETRY FOR GIVEN DETECTORS
     // THIS _SHOULD_ BE TEMPORARY. GLOBAL WIRE SUPPORT IS BEING ADDED TO THE LARSOFT GEOMETRY AND SHOULD BE AVAILABLE SOON
     if (fDetector == "dune35t") {
       unsigned int nwires = fGeom->Nwires(geo::PlaneID{wireID.Cryostat, 0, wireID.Plane});
       if (wireID.TPC == 0 or wireID.TPC == 1)
         globalWire = wireID.Wire;
       else if (wireID.TPC == 2 or wireID.TPC == 3 or wireID.TPC == 4 or wireID.TPC == 5)
         globalWire = nwires + wireID.Wire;
       else if (wireID.TPC == 6 or wireID.TPC == 7)
         globalWire = (2 * nwires) + wireID.Wire;
       else
         mf::LogError("BlurredClusterAlg")
           << "Error when trying to find a global induction plane coordinate for TPC " << wireID.TPC
           << " (geometry" << fDetector << ")";
     }
     else if (fDetector == "dune10kt") {
       unsigned int nwires = fGeom->Nwires(geo::PlaneID{wireID.Cryostat, 0, wireID.Plane});
       // Detector geometry has four TPCs, two on top of each other, repeated along z...
       int block = wireID.TPC / 4;
       globalWire = (nwires * block) + wireID.Wire;
     }
     else {
       auto const wireCenter = fGeom->WireIDToWireGeo(wireID).GetCenter();
       globalWire = fGeom->WireCoordinate(wireCenter,
                                          geo::PlaneID{wireID.Cryostat,
                                                       wireID.TPC % 2, // 0 or 1
                                                       wireID.Plane});
     }
   }
 
   return globalWire;
 }

TVector2 shower::EMShowerAlg::HitCoordinates_ ( recob::Hit const & hit ) const

private

Return the coordinates of this hit in global wire/tick space.

Definition at line 1664 of file EMShowerAlg.cxx.

References GlobalWire_(), recob::Hit::PeakTime(), and recob::Hit::WireID().

Referenced by ConstructTrack(), FindInitialTrack(), FindInitialTrackHits(), FindShowerStart_(), HitPosition_(), OrderShowerHits(), RelativeWireWidth_(), and WorstPlane_().

 {
   return TVector2(GlobalWire_(hit.WireID()), hit.PeakTime());
 }

TVector2 shower::EMShowerAlg::HitPosition_	(	detinfo::DetectorPropertiesData const &	detProp,
		recob::Hit const &	hit
	)		const

private

Return the coordinates of this hit in units of cm.

Definition at line 1669 of file EMShowerAlg.cxx.

References HitCoordinates_(), geo::WireID::planeID(), and recob::Hit::WireID().

Referenced by CheckShowerHits_(), FindOrderOfHits_(), MakeSpacePoints(), OrderShowerHits(), Project3DPointOntoPlane_(), ShowerCenter_(), ShowerDirection_(), ShowerHitRMS_(), and ShowerHitRMSGradient_().

 {
   geo::PlaneID planeID = hit.WireID().planeID();
   return HitPosition_(detProp, HitCoordinates_(hit), planeID);
 }

TVector2 shower::EMShowerAlg::HitPosition_	(	detinfo::DetectorPropertiesData const &	detProp,
		TVector2 const &	pos,
		geo::PlaneID	planeID
	)		const

private

Return the coordinates of this hit in units of cm.

Definition at line 1676 of file EMShowerAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), fGeom, and geo::GeometryCore::WirePitch().

 {
   return TVector2(pos.X() * fGeom->WirePitch(planeID), detProp.ConvertTicksToX(pos.Y(), planeID));
 }

bool shower::EMShowerAlg::isCleanShower ( std::vector< art::Ptr< recob::Hit >> const & hits ) const

Definition at line 2027 of file EMShowerAlg.cxx.

References hits(), and lar::to_element.

 {
   if (hits.empty()) return false;
   if (hits.size() > 2000) return true;
   if (hits.size() < 20) return true;
 
   std::map<int, int> hitmap;
   unsigned nhits = 0;
   for (auto const& hit : hits | views::transform(to_element)) {
     ++nhits;
     if (nhits > 2) ++hitmap[hit.WireID().Wire];
     if (nhits == 20) break;
   }
   if (float(nhits - 2) / hitmap.size() > 1.4)
     return false;
   else
     return true;
 }

std::unique_ptr< recob::Track > shower::EMShowerAlg::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

private

Takes initial track hits from multiple views and forms a track object which best represents the start of the shower

Definition at line 951 of file EMShowerAlg.cxx.

References CheckShowerHits_(), ConstructTrack(), fDebug, lar::dump::vector(), and WorstPlane_().

Referenced by FindInitialTrack().

 {
   // Can't do much with just one view
   if (initialHitsMap.size() < 2) {
     mf::LogInfo("EMShowerAlg") << "Only one useful view for this shower; nothing can be done\n";
     return std::unique_ptr<recob::Track>();
   }
 
   std::vector<std::pair<int, int>> initialHitsSize;
   for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator initialHitIt =
          initialHitsMap.begin();
        initialHitIt != initialHitsMap.end();
        ++initialHitIt)
     initialHitsSize.push_back(std::make_pair(initialHitIt->first, initialHitIt->second.size()));
 
   // Sort the planes by number of hits
   std::sort(initialHitsSize.begin(),
             initialHitsSize.end(),
             [](std::pair<int, int> const& size1, std::pair<int, int> const& size2) {
               return size1.second > size2.second;
             });
 
   // Pick the two planes to use to make the track
   //   -- if more than two planes, can choose the two views
   //      more accurately
   //   -- if not, just use the two views available
 
   std::vector<int> planes;
 
   if (initialHitsSize.size() > 2 and !CheckShowerHits_(detProp, showerHitsMap)) {
     int planeToIgnore = WorstPlane_(showerHitsMap);
     if (fDebug > 1)
       std::cout << "Igoring plane " << planeToIgnore << " in creation of initial track\n";
     for (std::vector<std::pair<int, int>>::iterator hitsSizeIt = initialHitsSize.begin();
          hitsSizeIt != initialHitsSize.end() and planes.size() < 2;
          ++hitsSizeIt) {
       if (hitsSizeIt->first == planeToIgnore) continue;
       planes.push_back(hitsSizeIt->first);
     }
   }
   else
     planes = {initialHitsSize.at(0).first, initialHitsSize.at(1).first};
 
   return ConstructTrack(detProp, initialHitsMap.at(planes.at(0)), initialHitsMap.at(planes.at(1)));
 }

recob::Shower shower::EMShowerAlg::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.

Definition at line 1000 of file EMShowerAlg.cxx.

References tca::dEdx(), fDebug, fGeom, FinddEdx_(), fShowerEnergyAlg, geo::GeometryCore::MaxPlanes(), shower::ShowerEnergyAlg::ShowerEnergy(), recob::Track::Vertex(), and recob::Track::VertexDirection().

Referenced by shower::EMShower::produce().

 {
 
   // Find the shower hits on each plane
   std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
   for (auto const& hitPtr : hits)
     planeHitsMap[hitPtr->View()].push_back(hitPtr);
 
   int bestPlane = -1;
   unsigned int highestNumberOfHits = 0;
   std::vector<double> totalEnergy, totalEnergyError, dEdx, dEdxError;
 
   // Look at each of the planes
   for (unsigned int plane = 0; plane < fGeom->MaxPlanes(); ++plane) {
 
     // If there's hits on this plane...
     if (planeHitsMap.count(plane)) {
       dEdx.push_back(FinddEdx_(clockData, detProp, initialHitsMap.at(plane), initialTrack));
       totalEnergy.push_back(
         fShowerEnergyAlg.ShowerEnergy(clockData, detProp, planeHitsMap.at(plane), plane));
       if (planeHitsMap.at(plane).size() > highestNumberOfHits and initialHitsMap.count(plane)) {
         bestPlane = plane;
         highestNumberOfHits = planeHitsMap.at(plane).size();
       }
     }
 
     // If not...
     else {
       dEdx.push_back(0);
       totalEnergy.push_back(0);
     }
   }
 
   TVector3 direction, directionError, showerStart, showerStartError;
   if (initialTrack) {
     direction = initialTrack->VertexDirection<TVector3>();
     showerStart = initialTrack->Vertex<TVector3>();
   }
 
   if (fDebug > 0) {
     std::cout << "Best plane is " << bestPlane;
     std::cout << "\ndE/dx for each plane is: " << dEdx[0] << ", " << dEdx[1] << " and " << dEdx[2];
     std::cout << "\nTotal energy for each plane is: " << totalEnergy[0] << ", " << totalEnergy[1]
               << " and " << totalEnergy[2];
     std::cout << "\nThe shower start is (" << showerStart.X() << ", " << showerStart.Y() << ", "
               << showerStart.Z() << ")\n";
     std::cout << "The shower direction is (" << direction.X() << ", " << direction.Y() << ", "
               << direction.Z() << ")\n";
   }
 
   return recob::Shower(direction,
                        directionError,
                        showerStart,
                        showerStartError,
                        totalEnergy,
                        totalEnergyError,
                        dEdx,
                        dEdxError,
                        bestPlane);
 }

recob::Shower shower::EMShowerAlg::MakeShower	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		art::PtrVector< recob::Hit > const &	hits,
		art::Ptr< recob::Vertex > const &	vertex,
		int &	iok
	)		const

Definition at line 1066 of file EMShowerAlg.cxx.

References util::abs(), pma::ProjectionMatchingAlg::buildSegment(), tca::dEdx(), calo::CalorimetryAlg::dEdx_AREA(), fCalorimetryAlg, fDebug, fdEdxTrackLength, fGeom, FindInitialTrackHits(), fProjectionMatchingAlg, fShowerEnergyAlg, geo::GeometryCore::MaxPlanes(), OrderShowerHits_(), geo::GeometryCore::Plane(), shower::ShowerEnergyAlg::ShowerEnergy(), pma::Track3D::size(), util::size(), lar::dump::vector(), geo::PlaneGeo::View(), geo::GeometryCore::WireAngleToVertical(), geo::GeometryCore::WirePitch(), and recob::Vertex::XYZ().

 {
   iok = 1;
 
   // Find the shower hits on each plane
   std::map<int, std::vector<art::Ptr<recob::Hit>>> planeHitsMap;
   for (auto const& hitPtr : hits)
     planeHitsMap[hitPtr->WireID().Plane].push_back(hitPtr);
 
   std::vector<std::vector<art::Ptr<recob::Hit>>> initialTrackHits(3);
 
   int pl0 = -1;
   int pl1 = -1;
   unsigned maxhits0 = 0;
   unsigned maxhits1 = 0;
 
   for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator planeHits = planeHitsMap.begin();
        planeHits != planeHitsMap.end();
        ++planeHits) {
 
     std::vector<art::Ptr<recob::Hit>> showerHits;
     OrderShowerHits_(detProp, planeHits->second, showerHits, vertex);
     FindInitialTrackHits(showerHits, vertex, initialTrackHits[planeHits->first]);
     if ((planeHits->second).size() > maxhits0) {
       if (pl0 != -1) {
         maxhits1 = maxhits0;
         pl1 = pl0;
       }
       pl0 = planeHits->first;
       maxhits0 = (planeHits->second).size();
     }
     else if ((planeHits->second).size() > maxhits1) {
       pl1 = planeHits->first;
       maxhits1 = (planeHits->second).size();
     }
   }
   if (pl0 != -1 && pl1 != -1 && initialTrackHits[pl0].size() >= 2 &&
       initialTrackHits[pl1].size() >= 2 &&
       initialTrackHits[pl0][0]->WireID().TPC == initialTrackHits[pl1][0]->WireID().TPC) {
     double xyz[3];
     vertex->XYZ(xyz);
     TVector3 vtx(xyz);
     pma::Track3D* pmatrack =
       fProjectionMatchingAlg.buildSegment(detProp, initialTrackHits[pl0], initialTrackHits[pl1]);
     std::vector<TVector3> spts;
 
     for (size_t i = 0; i < pmatrack->size(); ++i) {
       if ((*pmatrack)[i]->IsEnabled()) {
         TVector3 p3d = (*pmatrack)[i]->Point3D();
         spts.push_back(p3d);
       }
     }
     if (spts.size() >= 2) { // at least two space points
       TVector3 shwxyz, shwxyzerr;
       TVector3 shwdir, shwdirerr;
       std::vector<double> totalEnergy, totalEnergyError, dEdx, dEdxError;
       int bestPlane = pl0;
       double minpitch = 1000;
       std::vector<TVector3> dirs;
       if ((spts[0] - vtx).Mag() < (spts.back() - vtx).Mag()) {
         shwxyz = spts[0];
         size_t i = 5;
         if (spts.size() - 1 < 5) i = spts.size() - 1;
         shwdir = spts[i] - spts[0];
         shwdir = shwdir.Unit();
       }
       else {
         shwxyz = spts.back();
         size_t i = 0;
         if (spts.size() > 6) i = spts.size() - 6;
         shwdir = spts[i] - spts[spts.size() - 1];
         shwdir = shwdir.Unit();
       }
       for (unsigned int plane = 0; plane < fGeom->MaxPlanes(); ++plane) {
         if (planeHitsMap.find(plane) != planeHitsMap.end()) {
           totalEnergy.push_back(
             fShowerEnergyAlg.ShowerEnergy(clockData, detProp, planeHitsMap[plane], plane));
         }
         else {
           totalEnergy.push_back(0);
         }
         if (initialTrackHits[plane].size()) {
           double fdEdx = 0;
           double totQ = 0;
           double avgT = 0;
           double pitch = 0;
           double wirepitch = fGeom->WirePitch(initialTrackHits[plane][0]->WireID().planeID());
           double angleToVert =
             fGeom->WireAngleToVertical(fGeom->Plane(geo::PlaneID{0, 0, plane}).View(),
                                        initialTrackHits[plane][0]->WireID().planeID()) -
             0.5 * TMath::Pi();
           double cosgamma = std::abs(sin(angleToVert) * shwdir.Y() + cos(angleToVert) * shwdir.Z());
           if (cosgamma > 0) pitch = wirepitch / cosgamma;
           if (pitch) {
             if (pitch < minpitch) {
               minpitch = pitch;
               bestPlane = plane;
             }
             int nhits = 0;
             std::vector<float> vQ;
             for (auto const& hit : initialTrackHits[plane]) {
               int w1 = hit->WireID().Wire;
               int w0 = initialTrackHits[plane][0]->WireID().Wire;
               if (std::abs((w1 - w0) * pitch) < fdEdxTrackLength) {
                 vQ.push_back(hit->Integral());
                 totQ += hit->Integral();
                 avgT += hit->PeakTime();
                 ++nhits;
               }
             }
             if (totQ) {
               double dQdx = TMath::Median(vQ.size(), &vQ[0]) / pitch;
               fdEdx = fCalorimetryAlg.dEdx_AREA(
                 clockData, detProp, dQdx, avgT / nhits, initialTrackHits[plane][0]->WireID().Plane);
             }
           }
           dEdx.push_back(fdEdx);
         }
         else {
           dEdx.push_back(0);
         }
       }
       iok = 0;
       if (fDebug > 1) {
         std::cout << "Best plane is " << bestPlane;
         std::cout << "\ndE/dx for each plane is: " << dEdx[0] << ", " << dEdx[1] << " and "
                   << dEdx[2];
         std::cout << "\nTotal energy for each plane is: " << totalEnergy[0] << ", "
                   << totalEnergy[1] << " and " << totalEnergy[2];
         std::cout << "\nThe shower start is (" << shwxyz.X() << ", " << shwxyz.Y() << ", "
                   << shwxyz.Z() << ")\n";
         shwxyz.Print();
       }
 
       return recob::Shower(shwdir,
                            shwdirerr,
                            shwxyz,
                            shwxyzerr,
                            totalEnergy,
                            totalEnergyError,
                            dEdx,
                            dEdxError,
                            bestPlane);
     }
   }
   return recob::Shower();
 }

std::vector< recob::SpacePoint > shower::EMShowerAlg::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.

Definition at line 1218 of file EMShowerAlg.cxx.

References Construct3DPoint_(), fDebug, fGeom, fSpacePointSize, HitPosition_(), geo::GeometryCore::MaxPlanes(), Project3DPointOntoPlane_(), and lar::dump::vector().

Referenced by shower::EMShower::produce().

 {
   // Space points to return
   std::vector<recob::SpacePoint> spacePoints;
 
   // Make space points
   // Use the following procedure:
   //  -- Consider hits plane by plane
   //  -- For each hit on the first plane, consider the 3D point made by combining with each hit from the second plane
   //  -- Project this 3D point back into the two planes
   //  -- Determine how close to a the original hits this point lies
   //  -- If close enough, make a 3D space point from this point
   //  -- Discard these used hits in future iterations, along with hits in the
   //     third plane (if exists) close to the projection of the point into this
   //     plane
 
   // Container to hold used hits
   std::vector<int> usedHits;
 
   // Look through plane by plane
   for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitIt =
          showerHits.begin();
        showerHitIt != showerHits.end();
        ++showerHitIt) {
 
     // Find the other planes with hits
     std::vector<int> otherPlanes;
     for (unsigned int otherPlane = 0; otherPlane < fGeom->MaxPlanes(); ++otherPlane)
       if ((int)otherPlane != showerHitIt->first and showerHits.count(otherPlane))
         otherPlanes.push_back(otherPlane);
 
     // Can't make space points if we only have one view
     if (otherPlanes.size() == 0) return spacePoints;
 
     // Look at all hits on this plane
     for (std::vector<art::Ptr<recob::Hit>>::const_iterator planeHitIt = showerHitIt->second.begin();
          planeHitIt != showerHitIt->second.end();
          ++planeHitIt) {
 
       if (std::find(usedHits.begin(), usedHits.end(), planeHitIt->key()) != usedHits.end())
         continue;
 
       // Make a 3D point with every hit on the second plane
       const std::vector<art::Ptr<recob::Hit>> otherPlaneHits = showerHits.at(otherPlanes.at(0));
       for (std::vector<art::Ptr<recob::Hit>>::const_iterator otherPlaneHitIt =
              otherPlaneHits.begin();
            otherPlaneHitIt != otherPlaneHits.end() and
            std::find(usedHits.begin(), usedHits.end(), planeHitIt->key()) == usedHits.end();
            ++otherPlaneHitIt) {
 
         if ((*otherPlaneHitIt)->WireID().TPC != (*planeHitIt)->WireID().TPC or
             std::find(usedHits.begin(), usedHits.end(), otherPlaneHitIt->key()) != usedHits.end())
           continue;
 
         auto const point = Construct3DPoint_(detProp, *planeHitIt, *otherPlaneHitIt);
         std::vector<art::Ptr<recob::Hit>> pointHits;
         bool truePoint = false;
 
         if (otherPlanes.size() > 1) {
 
           TVector2 projThirdPlane = Project3DPointOntoPlane_(detProp, point, otherPlanes.at(1));
           const std::vector<art::Ptr<recob::Hit>> otherOtherPlaneHits =
             showerHits.at(otherPlanes.at(1));
 
           for (std::vector<art::Ptr<recob::Hit>>::const_iterator otherOtherPlaneHitIt =
                  otherOtherPlaneHits.begin();
                otherOtherPlaneHitIt != otherOtherPlaneHits.end() and !truePoint;
                ++otherOtherPlaneHitIt) {
 
             if ((*otherOtherPlaneHitIt)->WireID().TPC == (*planeHitIt)->WireID().TPC and
                 (projThirdPlane - HitPosition_(detProp, **otherOtherPlaneHitIt)).Mod() <
                   fSpacePointSize) {
 
               truePoint = true;
 
               // Remove hits used to make the point
               usedHits.push_back(planeHitIt->key());
               usedHits.push_back(otherPlaneHitIt->key());
               usedHits.push_back(otherOtherPlaneHitIt->key());
 
               pointHits.push_back(*planeHitIt);
               pointHits.push_back(*otherPlaneHitIt);
               pointHits.push_back(*otherOtherPlaneHitIt);
             }
           }
         }
 
         else if ((Project3DPointOntoPlane_(detProp, point, (*planeHitIt)->WireID().Plane) -
                   HitPosition_(detProp, **planeHitIt))
                      .Mod() < fSpacePointSize and
                  (Project3DPointOntoPlane_(detProp, point, (*otherPlaneHitIt)->WireID().Plane) -
                   HitPosition_(detProp, **otherPlaneHitIt))
                      .Mod() < fSpacePointSize) {
 
           truePoint = true;
 
           usedHits.push_back(planeHitIt->key());
           usedHits.push_back(otherPlaneHitIt->key());
 
           pointHits.push_back(*planeHitIt);
           pointHits.push_back(*otherPlaneHitIt);
         }
 
         // Make space point
         if (truePoint) {
           double xyz[3] = {point.X(), point.Y(), point.Z()};
           double xyzerr[6] = {fSpacePointSize,
                               fSpacePointSize,
                               fSpacePointSize,
                               fSpacePointSize,
                               fSpacePointSize,
                               fSpacePointSize};
           double chisq = 0.;
           spacePoints.emplace_back(xyz, xyzerr, chisq);
           hitAssns.push_back(pointHits);
         }
 
       } // end loop over second plane hits
 
     } // end loop over first plane hits
 
   } // end loop over planes
 
   if (fDebug > 0) {
     std::cout << "-------------------- Space points -------------------\n";
     std::cout << "There are " << spacePoints.size() << " space points:\n";
     if (fDebug > 1)
       for (std::vector<recob::SpacePoint>::const_iterator spacePointIt = spacePoints.begin();
            spacePointIt != spacePoints.end();
            ++spacePointIt) {
         const double* xyz = spacePointIt->XYZ();
         std::cout << "  Space point (" << xyz[0] << ", " << xyz[1] << ", " << xyz[2] << ")\n";
       }
   }
 
   return spacePoints;
 }

std::map< int, std::vector< art::Ptr< recob::Hit > > > shower::EMShowerAlg::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.

Ordering the shower hits requires three stages: – putting all the hits in a given plane in some kind of order – use the properties of the hits in all three planes to check this order – orient the hits correctly using properties of the shower

Definition at line 1359 of file EMShowerAlg.cxx.

References CheckIsolatedHits_(), CheckShowerHits_(), fDebug, fGeom, FindOrderOfHits_(), GetPlanePermutations_(), HitCoordinates_(), HitPosition_(), hits(), geo::GeometryCore::MaxPlanes(), n, RelativeWireWidth_(), ShowerHitRMS_(), ShowerHitRMSGradient_(), lar::to_element, lar::dump::vector(), X, and Y.

Referenced by shower::EMShower::produce().

 {
 
   // ------------- Put hits in order ------------
 
   // Find the shower hits on each plane
   std::map<int, std::vector<art::Ptr<recob::Hit>>> showerHitsMap;
   for (auto const& hitPtr : shower) {
     showerHitsMap[hitPtr->WireID().Plane].push_back(hitPtr);
   }
 
   // Order the hits, get the RMS and the RMS gradient for the hits in this plane
   std::map<int, double> planeRMSGradients, planeRMS;
   for (auto const& [plane, hits] : showerHitsMap) {
     if (desired_plane != plane and desired_plane != -1) continue;
     std::vector<art::Ptr<recob::Hit>> orderedHits = FindOrderOfHits_(detProp, hits);
     planeRMS[plane] = ShowerHitRMS_(detProp, orderedHits);
     planeRMSGradients[plane] = ShowerHitRMSGradient_(detProp, orderedHits);
     showerHitsMap[plane] = orderedHits;
   }
 
   if (fDebug > 1) {
     for (auto const [plane, shower_hit_rms] : planeRMS) {
       std::cout << "Plane " << plane << " has RMS " << shower_hit_rms << " and RMS gradient "
                 << planeRMSGradients.at(plane) << '\n';
     }
   }
 
   if (fDebug > 2) {
     std::cout << "\nHits in order; after ordering, before reversing...\n";
     for (auto const& [plane, hitPtrs] : showerHitsMap) {
       std::cout << "  Plane " << plane << '\n';
       for (auto const& hit : hitPtrs | views::transform(to_element)) {
         std::cout << "    Hit at (" << HitCoordinates_(hit).X() << ", " << HitCoordinates_(hit).Y()
                   << ") -- real wire " << hit.WireID() << ", hit position ("
                   << HitPosition_(detProp, hit).X() << ", " << HitPosition_(detProp, hit).Y()
                   << ")\n";
       }
     }
   }
 
   // ------------- Check between the views to ensure consistency of ordering -------------
 
   // Check between the views to make sure there isn't a poorly formed shower in just one view
   // First, determine the average RMS and RMS gradient across the other planes
   std::map<int, double> planeOtherRMS, planeOtherRMSGradients;
   for (std::map<int, double>::iterator planeRMSIt = planeRMS.begin(); planeRMSIt != planeRMS.end();
        ++planeRMSIt) {
     planeOtherRMS[planeRMSIt->first] = 0;
     planeOtherRMSGradients[planeRMSIt->first] = 0;
     int nOtherPlanes = 0;
     for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane) {
       if (plane != planeRMSIt->first and planeRMS.count(plane)) {
         planeOtherRMS[planeRMSIt->first] += planeRMS.at(plane);
         planeOtherRMSGradients[planeRMSIt->first] += planeRMSGradients.at(plane);
         ++nOtherPlanes;
       }
     }
     planeOtherRMS[planeRMSIt->first] /= (double)nOtherPlanes;
     planeOtherRMSGradients[planeRMSIt->first] /= (double)nOtherPlanes;
   }
 
   // Look to see if one plane has a particularly high RMS (compared to the
   // others) whilst having a similar gradient
   for (auto const& [plane, hitPtrs] : showerHitsMap) {
     if (planeRMS.at(plane) > planeOtherRMS.at(plane) * 2) {
       if (fDebug > 1) std::cout << "Plane " << plane << " was perpendicular... recalculating\n";
       std::vector<art::Ptr<recob::Hit>> orderedHits =
         this->FindOrderOfHits_(detProp, hitPtrs, true);
       showerHitsMap[plane] = orderedHits;
       planeRMSGradients[plane] = this->ShowerHitRMSGradient_(detProp, orderedHits);
     }
   }
 
   // ------------- Orient the shower correctly ---------------
 
   if (fDebug > 1) {
     std::cout << "Before reversing, here are the start and end points...\n";
     for (auto const& [plane, hitPtrs] : showerHitsMap) {
       std::cout << "  Plane " << plane << " has start (" << HitCoordinates_(*hitPtrs.front()).X()
                 << ", " << HitCoordinates_(*hitPtrs.front()).Y() << ") (real wire "
                 << hitPtrs.front()->WireID() << ") and end ("
                 << HitCoordinates_(*hitPtrs.back()).X() << ", "
                 << HitCoordinates_(*hitPtrs.back()).Y() << ") (real wire "
                 << hitPtrs.back()->WireID() << ")\n";
     }
   }
 
   // Use the RMS gradient information to get an initial ordering
   for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
          showerHitsMap.begin();
        showerHitsIt != showerHitsMap.end();
        ++showerHitsIt) {
     double gradient = planeRMSGradients.at(showerHitsIt->first);
     if (gradient < 0) std::reverse(showerHitsIt->second.begin(), showerHitsIt->second.end());
   }
 
   if (fDebug > 2) {
     std::cout << "\nHits in order; after reversing, before checking isolated hits...\n";
     for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
            showerHitsMap.begin();
          showerHitsIt != showerHitsMap.end();
          ++showerHitsIt) {
       std::cout << "  Plane " << showerHitsIt->first << '\n';
       for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsIt->second.begin();
            hitIt != showerHitsIt->second.end();
            ++hitIt)
         std::cout << "    Hit at (" << HitCoordinates_(**hitIt).X() << ", "
                   << HitCoordinates_(**hitIt).Y() << ") -- real wire " << (*hitIt)->WireID()
                   << ", hit position (" << HitPosition_(detProp, **hitIt).X() << ", "
                   << HitPosition_(detProp, **hitIt).Y() << ")\n";
     }
   }
 
   CheckIsolatedHits_(showerHitsMap);
 
   if (fDebug > 2) {
     std::cout << "\nHits in order; after checking isolated hits...\n";
     for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
            showerHitsMap.begin();
          showerHitsIt != showerHitsMap.end();
          ++showerHitsIt) {
       std::cout << "  Plane " << showerHitsIt->first << '\n';
       for (std::vector<art::Ptr<recob::Hit>>::iterator hitIt = showerHitsIt->second.begin();
            hitIt != showerHitsIt->second.end();
            ++hitIt)
         std::cout << "    Hit at (" << HitCoordinates_(**hitIt).X() << ", "
                   << HitCoordinates_(**hitIt).Y() << ") -- real wire " << (*hitIt)->WireID()
                   << ", hit position (" << HitPosition_(detProp, **hitIt).X() << ", "
                   << HitPosition_(detProp, **hitIt).Y() << ")\n";
     }
   }
 
   if (fDebug > 1) {
     std::cout << "After reversing and checking isolated hits, here are the "
                  "start and end points...\n";
     for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::iterator showerHitsIt =
            showerHitsMap.begin();
          showerHitsIt != showerHitsMap.end();
          ++showerHitsIt)
       std::cout << "  Plane " << showerHitsIt->first << " has start ("
                 << HitCoordinates_(*showerHitsIt->second.front()).X() << ", "
                 << HitCoordinates_(*showerHitsIt->second.front()).Y() << ") (real wire "
                 << showerHitsIt->second.front()->WireID() << ") and end ("
                 << HitCoordinates_(*showerHitsIt->second.back()).X() << ", "
                 << HitCoordinates_(*showerHitsIt->second.back()).Y() << ")\n";
   }
 
   // Check for views in which the shower travels almost along the wire planes
   // (shown by a small relative wire width)
   std::map<double, int> wireWidths = RelativeWireWidth_(showerHitsMap);
   std::vector<int> badPlanes;
   if (fDebug > 1) std::cout << "Here are the relative wire widths... \n";
   for (auto const [relative_wire_width, plane] : wireWidths) {
     if (fDebug > 1)
       std::cout << "Plane " << plane << " has relative wire width " << relative_wire_width << '\n';
     if (relative_wire_width < 1 / (double)wireWidths.size()) badPlanes.push_back(plane);
   }
 
   std::map<int, std::vector<art::Ptr<recob::Hit>>> ignoredPlanes;
   if (badPlanes.size() == 1) {
     int const badPlane = badPlanes[0];
     if (fDebug > 1) std::cout << "Ignoring plane " << badPlane << " when orientating\n";
     ignoredPlanes[badPlane] = showerHitsMap.at(badPlane);
     showerHitsMap.erase(badPlane);
   }
 
   // Consider all possible permutations of planes (0/1, oriented
   // correctly/incorrectly)
   std::map<int, std::map<int, bool>> permutations = GetPlanePermutations_(detProp, showerHitsMap);
 
   // Go through all permutations until we have a satisfactory orientation
   auto const originalShowerHitsMap = showerHitsMap;
 
   int n = 0;
   while (!CheckShowerHits_(detProp, showerHitsMap) and
          n < TMath::Power(2, (int)showerHitsMap.size())) {
     if (fDebug > 1) std::cout << "Permutation " << n << '\n';
     for (int const plane : showerHitsMap | views::keys) {
       auto hits = originalShowerHitsMap.at(plane);
       if (permutations[n][plane] == 1) { std::reverse(hits.begin(), hits.end()); }
       showerHitsMap[plane] = hits;
     }
     ++n;
   }
 
   // Go back to original if still no consistency
   if (!CheckShowerHits_(detProp, showerHitsMap)) showerHitsMap = originalShowerHitsMap;
 
   if (fDebug > 2) {
     std::cout << "End of OrderShowerHits: here are the order of hits:\n";
     for (auto const& [plane, hitPtrs] : showerHitsMap) {
       std::cout << "  Plane " << plane << '\n';
       for (auto const& hit : hitPtrs | views::transform(to_element)) {
         std::cout << "    Hit (" << HitCoordinates_(hit).X() << " (real wire " << hit.WireID()
                   << "), " << HitCoordinates_(hit).Y() << ") -- pos ("
                   << HitPosition_(detProp, hit).X() << ", " << HitPosition_(detProp, hit).Y()
                   << ")\n";
       }
     }
   }
 
   if (ignoredPlanes.size())
     showerHitsMap[ignoredPlanes.begin()->first] = ignoredPlanes.begin()->second;
 
   return showerHitsMap;
 }

void shower::EMShowerAlg::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

private

Takes the hits associated with a shower and orders then so they follow the direction of the shower

Definition at line 1574 of file EMShowerAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertXToTicks(), fGeom, FindOrderOfHits_(), geo::GeometryCore::FindTPCAtPosition(), art::Ptr< T >::isNull(), geo::CryostatID::isValid, recob::Hit::PeakTime(), geo::WireID::planeID(), recob::Vertex::position(), geo::TPCID::TPC, geo::WireID::Wire, geo::GeometryCore::WireCoordinate(), and recob::Hit::WireID().

Referenced by MakeShower().

 {
   showerHits = FindOrderOfHits_(detProp, shower);
 
   // Find TPC for the vertex
   auto const& xyz = vertex->position();
   geo::TPCID tpc = fGeom->FindTPCAtPosition(xyz);
   if (!tpc.isValid && showerHits.size()) tpc = geo::TPCID(showerHits[0]->WireID());
 
   // Find hits in the same TPC
   art::Ptr<recob::Hit> hit0, hit1;
   for (auto& hit : showerHits) {
     if (hit->WireID().TPC == tpc.TPC) {
       if (hit0.isNull()) { hit0 = hit; }
       hit1 = hit;
     }
   }
   if (hit0.isNull() || hit1.isNull()) return;
   TVector2 coord0 = TVector2(hit0->WireID().Wire, hit0->PeakTime());
   TVector2 coord1 = TVector2(hit1->WireID().Wire, hit1->PeakTime());
   TVector2 coordvtx = TVector2(fGeom->WireCoordinate(xyz, hit0->WireID().planeID()),
                                detProp.ConvertXToTicks(xyz.X(), hit0->WireID().planeID()));
   if ((coord1 - coordvtx).Mod() < (coord0 - coordvtx).Mod()) {
     std::reverse(showerHits.begin(), showerHits.end());
   }
 }

TVector2 shower::EMShowerAlg::Project3DPointOntoPlane_	(	detinfo::DetectorPropertiesData const &	detProp,
		geo::Point_t const &	point,
		int	plane,
		int	cryostat = `0`
	)		const

private

Projects a 3D point (units [cm]) onto a 2D plane Returns 2D point (units [cm])

Definition at line 1915 of file EMShowerAlg.cxx.

References detinfo::DetectorPropertiesData::ConvertXToTicks(), e, fGeom, geo::GeometryCore::FindTPCAtPosition(), GlobalWire_(), HitPosition_(), geo::CryostatID::isValid, geo::GeometryCore::NearestWireID(), geo::InvalidWireError::suggestedWireID(), and geo::TPCID::TPC.

Referenced by CheckShowerHits_(), ConstructTrack(), and MakeSpacePoints().

 {
 
   TVector2 wireTickPos = TVector2(-999., -999.);
 
   geo::TPCID tpcID = fGeom->FindTPCAtPosition(point);
   int tpc = 0;
   if (tpcID.isValid)
     tpc = tpcID.TPC;
   else
     return wireTickPos;
 
   // Construct wire ID for this point projected onto the plane
   geo::PlaneID planeID = geo::PlaneID(cryostat, tpc, plane);
   geo::WireID wireID;
   try {
     wireID = fGeom->NearestWireID(point, planeID);
   }
   catch (geo::InvalidWireError const& e) {
     wireID = e.suggestedWireID(); // pick the closest valid wire
   }
 
   wireTickPos = TVector2(GlobalWire_(wireID), detProp.ConvertXToTicks(point.X(), planeID));
 
   return HitPosition_(detProp, wireTickPos, planeID);
 }

std::map< double, int > shower::EMShowerAlg::RelativeWireWidth_ ( const std::map< int, std::vector< art::Ptr< recob::Hit >>> & showerHitsMap ) const

private

Determines the 'relative wire width', i.e. how spread a shower is across wires of each plane relative to the others If a shower travels along the wire directions in a particular view, it will have a smaller wire width in that view Returns a map relating these widths to each plane

Definition at line 1731 of file EMShowerAlg.cxx.

References util::abs(), fGeom, HitCoordinates_(), geo::GeometryCore::MaxPlanes(), lar::dump::vector(), and X.

Referenced by GetPlanePermutations_(), and OrderShowerHits().

 {
 
   // Get the wire widths
   std::map<int, int> planeWireLength;
   for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitsIt =
          showerHitsMap.begin();
        showerHitsIt != showerHitsMap.end();
        ++showerHitsIt)
     planeWireLength[showerHitsIt->first] =
       std::abs(HitCoordinates_(*showerHitsIt->second.front()).X() -
                HitCoordinates_(*showerHitsIt->second.back()).X());
 
   // Find the relative wire width for each plane with respect to the others
   std::map<int, double> planeOtherWireLengths;
   for (std::map<int, int>::iterator planeWireLengthIt = planeWireLength.begin();
        planeWireLengthIt != planeWireLength.end();
        ++planeWireLengthIt) {
     double quad = 0.;
     for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane)
       if (plane != planeWireLengthIt->first and planeWireLength.count(plane))
         quad += cet::square(planeWireLength[plane]);
     quad = std::sqrt(quad);
     planeOtherWireLengths[planeWireLengthIt->first] =
       planeWireLength[planeWireLengthIt->first] / (double)quad;
   }
 
   // Order these backwards
   std::map<double, int> wireWidthMap;
   for (std::map<int, std::vector<art::Ptr<recob::Hit>>>::const_iterator showerHitsIt =
          showerHitsMap.begin();
        showerHitsIt != showerHitsMap.end();
        ++showerHitsIt)
     wireWidthMap[planeOtherWireLengths.at(showerHitsIt->first)] = showerHitsIt->first;
 
   return wireWidthMap;
 }

TVector2 shower::EMShowerAlg::ShowerCenter_	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::Hit >> const &	showerHits
	)		const

private

Returns the charge-weighted shower center.

Definition at line 1796 of file EMShowerAlg.cxx.

References HitPosition_(), and lar::dump::vector().

Referenced by FindOrderOfHits_(), and ShowerHitRMS_().

 {
 
   TVector2 pos, chargePoint = TVector2(0, 0);
   double totalCharge = 0;
   for (std::vector<art::Ptr<recob::Hit>>::const_iterator hit = showerHits.begin();
        hit != showerHits.end();
        ++hit) {
     pos = HitPosition_(detProp, **hit);
     chargePoint += (*hit)->Integral() * pos;
     totalCharge += (*hit)->Integral();
   }
   TVector2 centre = chargePoint / totalCharge;
 
   return centre;
 }

TVector2 shower::EMShowerAlg::ShowerDirection_	(	detinfo::DetectorPropertiesData const &	detProp,
		const std::vector< art::Ptr< recob::Hit >> &	showerHits
	)		const

private

Returns a rough charge-weighted shower 'direction' given the hits in the shower

Definition at line 1770 of file EMShowerAlg.cxx.

References HitPosition_(), lar::dump::vector(), and weight.

Referenced by FindOrderOfHits_(), ShowerHitRMS_(), and ShowerHitRMSGradient_().

 {
 
   TVector2 pos;
   double weight = 1;
   double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
   for (std::vector<art::Ptr<recob::Hit>>::const_iterator hit = showerHits.begin();
        hit != showerHits.end();
        ++hit) {
     //++nhits;
     pos = HitPosition_(detProp, **hit);
     weight = cet::square((*hit)->Integral());
     sumweight += weight;
     sumx += weight * pos.X();
     sumy += weight * pos.Y();
     sumx2 += weight * pos.X() * pos.X();
     sumxy += weight * pos.X() * pos.Y();
   }
   double gradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
   TVector2 direction = TVector2(1, gradient).Unit();
 
   return direction;
 }

double shower::EMShowerAlg::ShowerHitRMS_	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::Hit >> const &	showerHits
	)		const

private

Returns the RMS of the hits from the central shower 'axis' along the length of the shower

Definition at line 1815 of file EMShowerAlg.cxx.

References HitPosition_(), proj, ShowerCenter_(), ShowerDirection_(), and lar::dump::vector().

Referenced by GetPlanePermutations_(), and OrderShowerHits().

 {
 
   TVector2 direction = ShowerDirection_(detProp, showerHits);
   TVector2 centre = ShowerCenter_(detProp, showerHits);
 
   std::vector<double> distanceToAxis;
   for (std::vector<art::Ptr<recob::Hit>>::const_iterator showerHitsIt = showerHits.begin();
        showerHitsIt != showerHits.end();
        ++showerHitsIt) {
     TVector2 proj = (HitPosition_(detProp, **showerHitsIt) - centre).Proj(direction) + centre;
     distanceToAxis.push_back((HitPosition_(detProp, **showerHitsIt) - proj).Mod());
   }
   double RMS = TMath::RMS(distanceToAxis.begin(), distanceToAxis.end());
 
   return RMS;
 }

double shower::EMShowerAlg::ShowerHitRMSGradient_	(	detinfo::DetectorPropertiesData const &	detProp,
		const std::vector< art::Ptr< recob::Hit >> &	showerHits
	)		const

private

Returns the gradient of the RMS vs shower segment graph.

Definition at line 1834 of file EMShowerAlg.cxx.

References bin, fMakeRMSGradientPlot, fNumShowerSegments, HitPosition_(), proj, ShowerDirection_(), lar::to_element, and weight.

Referenced by GetPlanePermutations_(), and OrderShowerHits().

 {
   // Find a rough shower 'direction' and center
   TVector2 direction = ShowerDirection_(detProp, showerHits);
 
   // Bin the hits into discreet chunks
   int nShowerSegments = fNumShowerSegments;
   double lengthOfShower =
     (HitPosition_(detProp, *showerHits.back()) - HitPosition_(detProp, *showerHits.front())).Mod();
   double lengthOfSegment = lengthOfShower / (double)nShowerSegments;
   std::map<int, std::vector<art::Ptr<recob::Hit>>> showerSegments;
   std::map<int, double> segmentCharge;
   for (auto const& hitPtr : showerHits) {
     auto const& hit = *hitPtr;
     int const segment =
       (HitPosition_(detProp, hit) - HitPosition_(detProp, *showerHits.front())).Mod() /
       lengthOfSegment;
     showerSegments[segment].push_back(hitPtr);
     segmentCharge[segment] += hit.Integral();
   }
 
   TGraph* graph = new TGraph();
   std::vector<std::pair<int, double>> binVsRMS;
 
   // Loop over the bins to find the distribution of hits as the shower
   // progresses
   for (auto const& [segment, hitPtrs] : showerSegments) {
 
     // Get the mean position of the hits in this bin
     TVector2 meanPosition(0, 0);
     for (auto const& hit : hitPtrs | views::transform(to_element))
       meanPosition += HitPosition_(detProp, hit);
     meanPosition /= (double)hitPtrs.size();
 
     // Get the RMS of this bin
     std::vector<double> distanceToAxisBin;
     for (auto const& hit : hitPtrs | views::transform(to_element)) {
       TVector2 proj = (HitPosition_(detProp, hit) - meanPosition).Proj(direction) + meanPosition;
       distanceToAxisBin.push_back((HitPosition_(detProp, hit) - proj).Mod());
     }
 
     double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
     binVsRMS.emplace_back(segment, RMSBin);
     if (fMakeRMSGradientPlot) graph->SetPoint(graph->GetN(), segment, RMSBin);
   }
 
   // Get the gradient of the RMS-bin plot
   double sumx = 0., sumy = 0., sumx2 = 0., sumxy = 0., sumweight = 0.;
   for (auto const& [bin, RMSBin] : binVsRMS) {
     double weight = segmentCharge.at(bin);
     sumweight += weight;
     sumx += weight * bin;
     sumy += weight * RMSBin;
     sumx2 += weight * bin * bin;
     sumxy += weight * bin * RMSBin;
   }
   double RMSgradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
 
   if (fMakeRMSGradientPlot) {
     TVector2 direction = TVector2(1, RMSgradient).Unit();
     TCanvas* canv = new TCanvas();
     graph->Draw();
     graph->GetXaxis()->SetTitle("Shower segment");
     graph->GetYaxis()->SetTitle("RMS of hit distribution");
     TVector2 centre = TVector2(graph->GetMean(1), graph->GetMean(2));
     TLine line;
     line.SetLineColor(2);
     line.DrawLine(centre.X() - 1000 * direction.X(),
                   centre.Y() - 1000 * direction.Y(),
                   centre.X() + 1000 * direction.X(),
                   centre.Y() + 1000 * direction.Y());
     canv->SaveAs("RMSGradient.png");
     delete canv;
   }
   delete graph;
 
   return RMSgradient;
 }

Int_t shower::EMShowerAlg::WeightedFit	(	const Int_t	n,
		const Double_t *	x,
		const Double_t *	y,
		const Double_t *	w,
		Double_t *	parm
	)		const

Definition at line 1983 of file EMShowerAlg.cxx.

References n.

Referenced by FindInitialTrackHits().

 {
 
   Double_t sumx = 0.;
   Double_t sumx2 = 0.;
   Double_t sumy = 0.;
   Double_t sumxy = 0.;
   Double_t sumw = 0.;
   Double_t eparm[2];
 
   parm[0] = 0.;
   parm[1] = 0.;
   eparm[0] = 0.;
   eparm[1] = 0.;
 
   for (Int_t i = 0; i < n; i++) {
     sumx += x[i] * w[i];
     sumx2 += x[i] * x[i] * w[i];
     sumy += y[i] * w[i];
     sumxy += x[i] * y[i] * w[i];
     sumw += w[i];
   }
 
   if (sumx2 * sumw - sumx * sumx == 0.) return 1;
   if (sumx2 - sumx * sumx / sumw == 0.) return 1;
 
   parm[0] = (sumy * sumx2 - sumx * sumxy) / (sumx2 * sumw - sumx * sumx);
   parm[1] = (sumxy - sumx * sumy / sumw) / (sumx2 - sumx * sumx / sumw);
 
   eparm[0] = sumx2 * (sumx2 * sumw - sumx * sumx);
   eparm[1] = (sumx2 - sumx * sumx / sumw);
 
   if (eparm[0] < 0. || eparm[1] < 0.) return 1;
 
   eparm[0] = sqrt(eparm[0]) / (sumx2 * sumw - sumx * sumx);
   eparm[1] = sqrt(eparm[1]) / (sumx2 - sumx * sumx / sumw);
 
   return 0;
 }

int shower::EMShowerAlg::WorstPlane_ ( const std::map< int, std::vector< art::Ptr< recob::Hit >>> & showerHitsMap ) const

private

Returns the plane which is determined to be the least likely to be correct.

Definition at line 1946 of file EMShowerAlg.cxx.

References util::abs(), fDebug, fGeom, HitCoordinates_(), hits(), geo::GeometryCore::MaxPlanes(), and X.

Referenced by MakeInitialTrack_().

 {
   // Get the width of the shower in wire coordinate
   std::map<int, int> planeWireLength;
   std::map<int, double> planeOtherWireLengths;
   for (auto const& [plane, hits] : showerHitsMap)
     planeWireLength[plane] =
       std::abs(HitCoordinates_(*hits.front()).X() - HitCoordinates_(*hits.back()).X());
   for (std::map<int, int>::iterator planeWireLengthIt = planeWireLength.begin();
        planeWireLengthIt != planeWireLength.end();
        ++planeWireLengthIt) {
     double quad = 0.;
     for (int plane = 0; plane < (int)fGeom->MaxPlanes(); ++plane)
       if (plane != planeWireLengthIt->first and planeWireLength.count(plane))
         quad += cet::square(planeWireLength[plane]);
     quad = std::sqrt(quad);
     planeOtherWireLengths[planeWireLengthIt->first] =
       planeWireLength[planeWireLengthIt->first] / (double)quad;
   }
 
   if (fDebug > 1)
     for (auto const [plane, relative_width] : planeOtherWireLengths) {
       std::cout << "Plane " << plane << " has " << planeWireLength[plane]
                 << " wire width and therefore has relative width in wire of " << relative_width
                 << '\n';
     }
 
   std::map<double, int> wireWidthMap;
   for (int const plane : showerHitsMap | views::keys) {
     double wireWidth = planeWireLength.at(plane);
     wireWidthMap[wireWidth] = plane;
   }
 
   return wireWidthMap.begin()->second;
 }