shower::EMShowerAlg Class Reference

#include "EMShowerAlg.h"

Public Member Functions
	EMShowerAlg (fhicl::ParameterSet const &pset)
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)
	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)
	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)
	Takes the initial showers found and tries to resolve issues where one bad view ruins the event. More...
std::unique_ptr< recob::Track >	ConstructTrack (std::vector< art::Ptr< recob::Hit > > const &track1, std::vector< art::Ptr< recob::Hit > > const &track2, std::map< int, TVector2 > const &showerCentreMap)
std::unique_ptr< recob::Track >	ConstructTrack (std::vector< art::Ptr< recob::Hit > > const &track1, std::vector< art::Ptr< recob::Hit > > const &track2)
void	FindInitialTrack (const std::map< int, std::vector< art::Ptr< recob::Hit > > > &hits, std::unique_ptr< recob::Track > &initialTrack, std::map< int, std::vector< art::Ptr< recob::Hit > > > &initialTrackHits, int plane)
	Finds the initial track-like part of the shower and the hits in all views associated with it. More...
std::vector< std::vector< int > >	FindShowers (std::map< int, std::vector< int > > const &trackToClusters)
	Makes showers given a map between tracks and all clusters associated with them. More...
recob::Shower	MakeShower (art::PtrVector< recob::Hit > const &hits, std::unique_ptr< recob::Track > const &initialTrack, std::map< int, std::vector< art::Ptr< recob::Hit > > > const &initialTrackHits)
	Makes a recob::Shower object given the hits in the shower and the initial track-like part. More...
recob::Shower	MakeShower (art::PtrVector< recob::Hit > const &hits, art::Ptr< recob::Vertex > const &vertex, int &iok)
	<Tingjun to="" document>=""> More...
std::vector< recob::SpacePoint >	MakeSpacePoints (std::map< int, std::vector< art::Ptr< recob::Hit > > > hits, std::vector< std::vector< art::Ptr< recob::Hit > > > &hitAssns)
	Makes space points from the shower hits in each plane. More...
std::map< int, std::vector< art::Ptr< recob::Hit > > >	OrderShowerHits (const art::PtrVector< recob::Hit > &shower, int plane)
	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)
	<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)
	<Tingjun to="" document>=""> More...
bool	isCleanShower (std::vector< art::Ptr< recob::Hit > > const &hits)
	<Tingjun to="" document>=""> More...

Public Attributes
int	fDebug

Private Member Functions
void	CheckIsolatedHits (std::map< int, std::vector< art::Ptr< recob::Hit > > > &showerHitsMap)
	Checks the hits across the views in a given shower to determine if there is one in the incorrect TPC. More...
bool	CheckShowerHits (std::map< int, std::vector< art::Ptr< recob::Hit > > > const &showerHitsMap)
TVector3	Construct3DPoint (art::Ptr< recob::Hit > const &hit1, art::Ptr< recob::Hit > const &hit2)
	Constructs a 3D point (in [cm]) to represent the hits given in two views. More...
double	FinddEdx (std::vector< art::Ptr< recob::Hit > > const &trackHits, std::unique_ptr< recob::Track > const &track)
	Finds dE/dx for the track given a set of hits. More...
std::vector< art::Ptr< recob::Hit > >	FindOrderOfHits (std::vector< art::Ptr< recob::Hit > > const &hits, bool perpendicular=false)
std::map< int, std::vector< art::Ptr< recob::Hit > > >	FindShowerStart (std::map< int, std::vector< art::Ptr< recob::Hit > > > const &orderedShowerMap)
std::map< int, std::map< int, bool > >	GetPlanePermutations (const std::map< int, std::vector< art::Ptr< recob::Hit > > > &showerHitsMap)
double	GlobalWire (const geo::WireID &wireID)
	Find the global wire position. More...
TVector2	HitCoordinates (art::Ptr< recob::Hit > const &hit)
	Return the coordinates of this hit in global wire/tick space. More...
TVector2	HitPosition (art::Ptr< recob::Hit > const &hit)
	Return the coordinates of this hit in units of cm. More...
TVector2	HitPosition (TVector2 const &pos, geo::PlaneID planeID)
	Return the coordinates of this hit in units of cm. More...
std::unique_ptr< recob::Track >	MakeInitialTrack (std::map< int, std::vector< art::Ptr< recob::Hit > > > const &initialHitsMap, std::map< int, std::vector< art::Ptr< recob::Hit > > > const &showerHitsMap)
	Takes initial track hits from multiple views and forms a track object which best represents the start of the shower. More...
void	OrderShowerHits (std::vector< art::Ptr< recob::Hit > > const &shower, std::vector< art::Ptr< recob::Hit > > &orderedShower, art::Ptr< recob::Vertex > const &vertex)
	Takes the hits associated with a shower and orders then so they follow the direction of the shower. More...
TVector2	Project3DPointOntoPlane (TVector3 const &point, int plane, int cryostat=0)
std::map< double, int >	RelativeWireWidth (const std::map< int, std::vector< art::Ptr< recob::Hit > > > &showerHitsMap)
TVector2	ShowerCentre (const std::vector< art::Ptr< recob::Hit > > &showerHits)
	Returns the charge-weighted shower centre. More...
TVector2	ShowerDirection (const std::vector< art::Ptr< recob::Hit > > &showerHits)
	Returns a rough charge-weighted shower 'direction' given the hits in the shower. More...
double	ShowerHitRMS (const std::vector< art::Ptr< recob::Hit > > &showerHits)
	Returns the RMS of the hits from the central shower 'axis' along the length of the shower. More...
double	ShowerHitRMSGradient (const std::vector< art::Ptr< recob::Hit > > &showerHits, TVector2 trueStart=TVector2(0, 0))
	Returns the gradient of the RMS vs shower segment graph. More...
int	WorstPlane (const std::map< int, std::vector< art::Ptr< recob::Hit > > > &showerHitsMap)
	Returns the plane which is determined to be the least likely to be correct. More...
int	FindTrueParticle (const std::vector< art::Ptr< recob::Hit > > &showerHits)
int	FindParticleID (const art::Ptr< recob::Hit > &hit)
void	MakePicture ()

Private Attributes
double	fMinTrackLength
double	fdEdxTrackLength
double	fSpacePointSize
unsigned int	fNfitpass
std::vector< unsigned int >	fNfithits
std::vector< double >	fToler
art::ServiceHandle< geo::Geometry >	fGeom
detinfo::DetectorProperties const *	fDetProp
art::ServiceHandle< art::TFileService >	tfs
shower::ShowerEnergyAlg	fShowerEnergyAlg
calo::CalorimetryAlg	fCalorimetryAlg
pma::ProjectionMatchingAlg	fProjectionMatchingAlg
std::string	fDetector
art::ServiceHandle< cheat::BackTrackerService >	bt_serv
TH1I *	hTrueDirection
TProfile *	hNumHitsInSegment
TProfile *	hNumSegments
bool	fMakeGradientPlot
bool	fMakeRMSGradientPlot
int	fNumShowerSegments

Detailed Description

Definition at line 74 of file EMShowerAlg.h.

Constructor & Destructor Documentation

shower::EMShowerAlg::EMShowerAlg

(

fhicl::ParameterSet const &

pset

)

Definition at line 13 of file EMShowerAlg.cxx.

References art::errors::Configuration, fdEdxTrackLength, fDetector, fMakeGradientPlot, fMakeRMSGradientPlot, fMinTrackLength, fNfithits, fNfitpass, fNumShowerSegments, fSpacePointSize, fToler, fhicl::ParameterSet::get(), hNumHitsInSegment, hNumSegments, hTrueDirection, art::TFileDirectory::make(), and tfs.

                                                            : fDetProp(lar::providerFrom<detinfo::DetectorPropertiesService>()),
                                                                    fShowerEnergyAlg(pset.get<fhicl::ParameterSet>("ShowerEnergyAlg")),
                                                                    fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg")),
                                                                    fProjectionMatchingAlg(pset.get<fhicl::ParameterSet>("ProjectionMatchingAlg")) {

  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");
  if (fNfitpass!=fNfithits.size() ||
      fNfitpass!=fToler.size()) {
    throw art::Exception(art::errors::Configuration)
      << "EMShowerAlg: fNfithits and fToler need to have size fNfitpass";
  }
  fDetector = pset.get<std::string>("Detector","dune35t");

  hTrueDirection = tfs->make<TH1I>("trueDir","",2,0,2);

  //this->MakePicture();
  // tmp
  fMakeGradientPlot = pset.get<bool>("MakeGradientPlot",false);
  fMakeRMSGradientPlot = pset.get<bool>("MakeRMSGradientPlot",false);
  fNumShowerSegments = pset.get<int>("NumShowerSegments",4);
  hNumHitsInSegment = tfs->make<TProfile>("NumHitsInSegment","",100,0,100);
  hNumSegments = tfs->make<TProfile>("NumSegments","",100,0,100);

}

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
	)

Map associated tracks and clusters together given their associated hits.

Definition at line 43 of file EMShowerAlg.cxx.

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

                                                                                                   {

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

  return;

}

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
	)

Map associated tracks and clusters together given their associated hits, whilst ignoring certain clusters.

Definition at line 56 of file EMShowerAlg.cxx.

References evd::details::begin(), evd::details::end(), fDebug, fMinTrackLength, track, and lar::dump::vector().

                                                                                                   {

  // Look through all the clusters
  for (std::vector<art::Ptr<recob::Cluster> >::const_iterator clusterIt = clusters.begin(); clusterIt != clusters.end(); ++clusterIt) {

    // Get the hits in this cluster
    std::vector<art::Ptr<recob::Hit> > clusterHits = fmh.at(clusterIt->key());

    // Look at all these hits and find the associated tracks
    for (std::vector<art::Ptr<recob::Hit> >::iterator clusterHitIt = clusterHits.begin(); clusterHitIt != clusterHits.end(); ++clusterHitIt) {

      // Get the tracks associated with this hit
      std::vector<art::Ptr<recob::Track> > clusterHitTracks = fmt.at(clusterHitIt->key());
      if (clusterHitTracks.size() > 1) { std::cout << "More than one track associated with this hit!" << std::endl; continue; }
      if (clusterHitTracks.size() < 1) continue;
      if (clusterHitTracks.at(0)->Length() < fMinTrackLength) {
        if (fDebug > 2)
          std::cout << "Track " << clusterHitTracks.at(0)->ID() << " is too short! (" << clusterHitTracks.at(0)->Length() << ")" << std::endl;
        continue;
      }

      // Add this cluster to the track map
      int track = clusterHitTracks.at(0).key();
      //int trackID = clusterHitTracks.at(0)->ID();
      int cluster = (*clusterIt).key();
      if (std::find(clustersToIgnore.begin(), clustersToIgnore.end(), cluster) != clustersToIgnore.end())
        continue;
      if (std::find(trackToClusters[track].begin(), trackToClusters[track].end(), cluster) == trackToClusters[track].end())
        trackToClusters[track].push_back(cluster);
      if (std::find(clusterToTracks[cluster].begin(), clusterToTracks[cluster].end(), track) == clusterToTracks[cluster].end())
        clusterToTracks[cluster].push_back(track);

    }

  }

  return;

}

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

private

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

Definition at line 101 of file EMShowerAlg.cxx.

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

Referenced by OrderShowerHits().

                                                                                                  {

  std::map<int,std::vector<int> > firstTPC;
  for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt)
    firstTPC[showerHitsIt->second.at(0)->WireID().TPC].push_back(showerHitsIt->first);

  // If all in the same TPC then that's great!
  if (firstTPC.size() == 1)
    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 (std::map<int,std::vector<int> >::iterator firstTPCIt = firstTPC.begin(); firstTPCIt != firstTPC.end(); ++firstTPCIt)
    if (firstTPCIt->second.size() == 1)
      problemPlane = firstTPCIt->second.at(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 (std::vector<int>::iterator otherPlaneIt = otherPlanes.begin(); otherPlaneIt != otherPlanes.end(); ++otherPlaneIt)
    for (std::vector<art::Ptr<recob::Hit> >::iterator hitIt = std::next(showerHitsMap.at(*otherPlaneIt).begin());
         hitIt != showerHitsMap.at(*otherPlaneIt).end() and hitIt != std::next(showerHitsMap.at(*otherPlaneIt).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 (std::vector<art::Ptr<recob::Hit> >::iterator naughty_hitIt = naughty_hits.begin(); naughty_hitIt != naughty_hits.end(); ++naughty_hitIt)
      showerHitsMap.at(problemPlane).push_back(*naughty_hitIt);
  }

  return;

}

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

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 183 of file EMShowerAlg.cxx.

References Construct3DPoint(), fDebug, HitPosition(), Project3DPointOntoPlane(), lar::dump::vector(), 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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
      startHits.push_back(showerHitsIt->second.front());
      planes.push_back(showerHitsIt->first);
    }

    TVector3 showerStartPos = Construct3DPoint(startHits.at(0), startHits.at(1));
    TVector2 proj1 = Project3DPointOntoPlane(showerStartPos, planes.at(0));
    TVector2 proj2 = Project3DPointOntoPlane(showerStartPos, planes.at(1));

    double timingDifference = TMath::Abs( startHits.at(0)->PeakTime() - startHits.at(1)->PeakTime() );
    double projectionDifference = ( (HitPosition(startHits.at(0)) - proj1).Mod() + (HitPosition(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() << ")" << std::endl;

  }

  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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitIt = showerHitsMap.begin(); showerHitIt != showerHitsMap.end(); ++showerHitIt)
      start2DMap[showerHitIt->first] = showerHitIt->second.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);

      TVector3 showerStartPos = Construct3DPoint(start2DMap.at(otherPlanes.at(0)), start2DMap.at(otherPlanes.at(1)));
      TVector2 showerStartProj = Project3DPointOntoPlane(showerStartPos, plane);

      if (fDebug > 2) {
        std::cout << "Plane... " << plane << std::endl;
        std::cout << "Start position in this plane is " << HitPosition(start2DMap.at(plane)).X() << ", " << HitPosition(start2DMap.at(plane)).Y() << ")" << std::endl;
        std::cout << "Shower start from other two planes is (" << showerStartPos.X() << ", " << showerStartPos.Y() << ", " << showerStartPos.Z() << ")" << std::endl;
        std::cout << "Projecting the other two planes gives position (" << showerStartProj.X() << ", " << showerStartProj.Y() << ")" << std::endl;
      }

      double projDiff = TMath::Abs((showerStartProj-HitPosition(start2DMap.at(plane))).Mod());
      double timeDiff = TMath::Max(TMath::Abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(0))->PeakTime()),
                                   TMath::Abs(start2DMap.at(plane)->PeakTime() - start2DMap.at(otherPlanes.at(1))->PeakTime()));

      if (fDebug > 1)
        std::cout << "Plane " << plane << " has projDiff " << projDiff << " and timeDiff " << timeDiff << std::endl;

      if (projDiff > 5 or timeDiff > 40)
        consistencyCheck = false;

    }

  }

  if (fDebug > 1)
    std::cout << "Consistency check is " << consistencyCheck << std::endl;

  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
	)

Takes the initial showers found and tries to resolve issues where one bad view ruins the event.

Definition at line 271 of file EMShowerAlg.cxx.

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

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

                                                                                           {

  std::vector<int> clustersToIgnore;

  // // Look at each shower
  // for (std::vector<std::vector<int> >::const_iterator initialShowerIt = initialShowers.begin(); initialShowerIt != initialShowers.end(); ++initialShowerIt) {

  //   // 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> > > planeClusterHits;
  //   for (std::vector<int>::const_iterator clusterIt = initialShowerIt->begin(); clusterIt != initialShowerIt->end(); ++clusterIt) {
  //     art::Ptr<recob::Cluster> cluster = clusters.at(*clusterIt);
  //     std::vector<art::Ptr<recob::Hit> > hits = fmh.at(cluster.key());
  //     planeClusters[cluster->Plane().Plane].push_back(cluster);
  //     for (std::vector<art::Ptr<recob::Hit> >::iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt)
  //    planeClusterHits[cluster.key()].push_back(*hitIt);
  //   }

  //   // Can't do much with fewer than three views
  //   if (planeClusters.size() < 3)
  //     continue;

  //   // Look at the average RMS of clusters in each view
  //   std::map<int,double> avRMS;
  //   for (std::map<int,std::vector<art::Ptr<recob::Cluster> > >::iterator planeClusterIt = planeClusters.begin(); planeClusterIt != planeClusters.end(); ++planeClusterIt) {
  //     std::cout << "Plane " << planeClusterIt->first << std::endl;
  //     for (std::vector<art::Ptr<recob::Cluster> >::iterator clusterIt = planeClusterIt->second.begin(); clusterIt != planeClusterIt->second.end(); ++clusterIt) {
  //    double rms = ShowerHitRMS(planeClusterHits.at(clusterIt->key()));
  //    std::cout << "Cluster " << clusterIt->key() << " has RMS " << rms << std::endl;
  //     }
  //   }
  // }


  // Look at each shower
  for (std::vector<std::vector<int> >::const_iterator initialShowerIt = initialShowers.begin(); initialShowerIt != initialShowers.end(); ++initialShowerIt) {

    if (std::distance(initialShowers.begin(),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 (std::vector<int>::const_iterator clusterIt = initialShowerIt->begin(); clusterIt != initialShowerIt->end(); ++clusterIt) {
      art::Ptr<recob::Cluster> cluster = clusters.at(*clusterIt);
      std::vector<art::Ptr<recob::Hit> > hits = fmh.at(cluster.key());
      planeClusters[cluster->Plane().Plane].push_back(cluster);
      for (std::vector<art::Ptr<recob::Hit> >::iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt)
        planeHits[(*hitIt)->WireID().Plane].push_back(*hitIt);
    }

    TFile* outFile = new TFile("chargeDistributions.root","RECREATE");
    std::map<int,TH1D*> chargeDist;
    for (std::map<int,std::vector<art::Ptr<recob::Cluster> > >::iterator planeIt = planeClusters.begin(); planeIt != planeClusters.end(); ++planeIt) {
      for (std::vector<art::Ptr<recob::Cluster> >::iterator clusterIt = planeIt->second.begin(); clusterIt != planeIt->second.end(); ++clusterIt) {
        chargeDist[planeIt->first] = new TH1D(std::string("chargeDist_Plane"+std::to_string(planeIt->first)+"_Cluster"+std::to_string(clusterIt->key())).c_str(),"",150,0,1000);
        std::vector<art::Ptr<recob::Hit> > hits = fmh.at(clusterIt->key());
        for (std::vector<art::Ptr<recob::Hit> >::iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt)
          chargeDist[planeIt->first]->Fill((*hitIt)->Integral());
        outFile->cd();
        chargeDist[planeIt->first]->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 (std::map<int,std::vector<double> >::iterator planeClusterSizesIt = planeClusterSizes.begin(); planeClusterSizesIt != planeClusterSizes.end(); ++planeClusterSizesIt) {
      double average = 0;
      for (std::vector<double>::iterator planeClusterSizeIt = planeClusterSizesIt->second.begin(); planeClusterSizeIt != planeClusterSizesIt->second.end(); ++planeClusterSizeIt)
        average += *planeClusterSizeIt;
      average /= planeClusterSizesIt->second.size();
      planeClustersAvSizes[planeClusterSizesIt->first] = 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 (std::map<int,double>::iterator clusterAvSizeIt = planeClustersAvSizes.begin(); clusterAvSizeIt != planeClustersAvSizes.end(); ++clusterAvSizeIt) {

      // Get averages from other planes and add in quadrature
      std::vector<double> otherAverages;
      for (std::map<int,double>::iterator otherClustersAvSizeIt = planeClustersAvSizes.begin(); otherClustersAvSizeIt != planeClustersAvSizes.end(); ++otherClustersAvSizeIt)
        if (clusterAvSizeIt->first != otherClustersAvSizeIt->first)
          otherAverages.push_back(otherClustersAvSizeIt->second);
      double quadOtherPlanes = 0;
      for (std::vector<double>::iterator otherAvIt = otherAverages.begin(); otherAvIt != otherAverages.end(); ++otherAvIt)
        quadOtherPlanes += TMath::Power(*otherAvIt,2);
      quadOtherPlanes = TMath::Sqrt(quadOtherPlanes);

      // If this plane has an average higher than the quadratic sum of the others, it may be bad
      if (clusterAvSizeIt->second >= quadOtherPlanes)
        badPlane = clusterAvSizeIt->first;

    }

    if (badPlane != -1) {
      if (fDebug > 1)
        std::cout << "Bad plane is " << badPlane << std::endl;
      for (std::vector<art::Ptr<recob::Cluster> >::iterator clusterIt = planeClusters.at(badPlane).begin(); clusterIt != planeClusters.at(badPlane).end(); ++clusterIt)
        clustersToIgnore.push_back(clusterIt->key());
    }

  }

  return clustersToIgnore;

}

TVector3 shower::EMShowerAlg::Construct3DPoint ( art::Ptr< recob::Hit > const &  hit1, art::Ptr< recob::Hit > const &  hit2  )

private

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

Definition at line 396 of file EMShowerAlg.cxx.

References detinfo::DetectorProperties::ConvertTicksToX(), fDetProp, 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 = (fDetProp->ConvertTicksToX(hit1->PeakTime(), hit1->WireID().planeID()) + fDetProp->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 TVector3(x, intersection.y, intersection.z);

}

std::unique_ptr< recob::Track > shower::EMShowerAlg::ConstructTrack	(	std::vector< art::Ptr< recob::Hit > > const &	track1,
		std::vector< art::Ptr< recob::Hit > > const &	track2,
		std::map< int, TVector2 > const &	showerCentreMap
	)

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 409 of file EMShowerAlg.cxx.

References pma::ProjectionMatchingAlg::buildSegment(), Construct3DPoint(), evd::details::end(), fDebug, fProjectionMatchingAlg, HitCoordinates(), if(), Project3DPointOntoPlane(), track, and lar::dump::vector().

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 (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits1.begin(); hitIt != hits1.end(); ++hitIt)
    ++tpcMap[(*hitIt)->WireID().TPC];
  for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits2.begin(); hitIt != hits2.end(); ++hitIt)
    ++tpcMap[(*hitIt)->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 (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits1.begin(); hitIt != hits1.end(); ++hitIt)
      if ((*hitIt)->WireID().TPC == firstTPC1) track1.push_back(*hitIt);
    for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits2.begin(); hitIt != hits2.end(); ++hitIt)
      if ((*hitIt)->WireID().TPC == firstTPC2) track2.push_back(*hitIt);    
  }
  else {
    track1 = hits1;
    track2 = hits2;
  }

  if (fDebug > 1) {
    std::cout << "About to make a track from these hits:" << std::endl;
    for (std::vector<art::Ptr<recob::Hit> >::const_iterator hit1 = track1.begin(); hit1 != track1.end(); ++hit1)
      std::cout << "Hit (" << HitCoordinates(*hit1).X() << ", " << HitCoordinates(*hit1).Y() << ") (real wire " << (*hit1)->WireID().Wire << ") in TPC " << (*hit1)->WireID().TPC << std::endl;
    for (std::vector<art::Ptr<recob::Hit> >::const_iterator hit2 = track2.begin(); hit2 != track2.end(); ++hit2)
      std::cout << "Hit (" << HitCoordinates(*hit2).X() << ", " << HitCoordinates(*hit2).Y() << ") (real wire " << (*hit2)->WireID().Wire << ") in TPC " << (*hit2)->WireID().TPC << std::endl;
  }

  TVector3 trackStart = Construct3DPoint(track1.at(0), track2.at(0));
  pma::Track3D* pmatrack = fProjectionMatchingAlg.buildSegment(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;
  std::vector<std::vector<double> > dEdx; // Right now, not finding the dE/dx for these tracks.  Can extend if needed.

  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();
      // TVector3 dc((*pmatrack)[i+1]->Point3D());
      // dc -= (*pmatrack)[i]->Point3D();
      // dc *= 1.0 / dc.Mag();
      dircos.push_back(dc);
    }
    else dircos.push_back(dircos.back());

  }

  // Orient the track correctly
  std::map<int,double> distanceToVertex, distanceToEnd;
  TVector3 vertex = *xyz.begin(), end = *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(vertex, showerCentreIt->first);
    TVector2 endProj    = Project3DPointOntoPlane(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 << std::endl;

  // 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>(xyz, dircos, dEdx);

  return track;

}

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

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 532 of file EMShowerAlg.cxx.

References ConstructTrack().

                                                                                                              {

  std::map<int,TVector2> showerCentreMap;

  return this->ConstructTrack(track1, track2, showerCentreMap);

}

double shower::EMShowerAlg::FinddEdx ( std::vector< art::Ptr< recob::Hit > > const &  trackHits, std::unique_ptr< recob::Track > const &  track  )

private

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

Definition at line 541 of file EMShowerAlg.cxx.

References calo::CalorimetryAlg::dEdx_AREA(), fCalorimetryAlg, fdEdxTrackLength, lar::util::TrackPitchInView(), and lar::dump::vector().

Referenced by MakeShower().

                                                                                                                        {

  double totalCharge = 0, totalDistance = 0, avHitTime = 0;
  unsigned int nHits = 0;

  if (!track)
    return -999;

  // size_t trajectory_point = 0;
  // double wirePitch   = fGeom->WirePitch(trackHits.at(0)->View(), 1, 0);
  // double angleToVert = fGeom->WireAngleToVertical(trackHits.at(0)->View(), 1, 0) - 0.5*::util::pi<>();
  // const TVector3& dir = track->DirectionAtPoint(trajectory_point);
  // //(sin(angleToVert),cos(angleToVert)) is the direction perpendicular to wire
  // double cosgamma = std::abs(std::sin(angleToVert)*dir.Y() +
  //                         std::cos(angleToVert)*dir.Z());

  // if(cosgamma < 1.e-5)
  //   throw cet::exception("Track") << "cosgamma is basically 0, that can't be right\n";
  // double pitch =  wirePitch/cosgamma;

  // Get the pitch
  double pitch = 0;
  try { pitch = lar::util::TrackPitchInView(*track, trackHits.at(0)->View()); }
  catch(...) { pitch = 0; }

  // Deal with large pitches
  if (pitch > fdEdxTrackLength) {
    double dEdx = fCalorimetryAlg.dEdx_AREA(*trackHits.begin(), pitch);
    return dEdx;
  }

  for (std::vector<art::Ptr<recob::Hit> >::const_iterator trackHitIt = trackHits.begin(); trackHitIt != trackHits.end(); ++trackHitIt) {
    if (totalDistance + pitch < fdEdxTrackLength) {
      totalDistance += pitch;
      totalCharge += (*trackHitIt)->Integral();
      avHitTime += (*trackHitIt)->PeakTime();
      ++nHits;
    }
  }

  avHitTime /= (double)nHits;

  double dQdx = totalCharge / totalDistance;
  double dEdx = fCalorimetryAlg.dEdx_AREA(dQdx, avHitTime, trackHits.at(0)->WireID().Plane);

  return dEdx;

}

void shower::EMShowerAlg::FindInitialTrack	(	const std::map< int, std::vector< art::Ptr< recob::Hit > > > &	hits,
		std::unique_ptr< recob::Track > &	initialTrack,
		std::map< int, std::vector< art::Ptr< recob::Hit > > > &	initialTrackHits,
		int	plane
	)

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: – put the hits in the correct order in each view – find the initial track-like hits in each view – use these to construct a track

Definition at line 590 of file EMShowerAlg.cxx.

References fDebug, FindShowerStart(), HitCoordinates(), MakeInitialTrack(), lar::dump::vector(), recob::Track::Vertex(), and recob::Track::VertexDirection().

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

                                                                                                                   {


  // // First, order the hits into the correct shower order in each plane
  // if (fDebug > 1)
  //   std::cout << " ------------------ Ordering shower hits -------------------- " << std::endl;
  // std::map<int,std::vector<art::Ptr<recob::Hit> > > showerHitsMap = OrderShowerHits(hits, plane);
  // if (fDebug > 1)
  //   std::cout << " ------------------ End ordering shower hits -------------------- " << std::endl;

  // Now find the hits belonging to the track
  if (fDebug > 1)
    std::cout << " ------------------ Finding initial track hits -------------------- " << std::endl;
  initialTrackHits = FindShowerStart(showerHitsMap);
  if (fDebug > 1) {
    std::cout << "Here are the initial shower hits... " << std::endl;
    for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator initialTrackHitsIt = initialTrackHits.begin(); initialTrackHitsIt != initialTrackHits.end(); ++initialTrackHitsIt) {
      std::cout << "  Plane " << initialTrackHitsIt->first << std::endl;
      for (std::vector<art::Ptr<recob::Hit> >::iterator initialTrackHitIt = initialTrackHitsIt->second.begin(); initialTrackHitIt != initialTrackHitsIt->second.end(); ++initialTrackHitIt)
        std::cout << "    Hit is (" << HitCoordinates(*initialTrackHitIt).X() << " (real hit " << (*initialTrackHitIt)->WireID() << "), " << HitCoordinates(*initialTrackHitIt).Y() << ")" << std::endl;
    }
  }
  if (fDebug > 1)
    std::cout << " ------------------ End finding initial track hits -------------------- " << std::endl;

  // Now we have the track hits -- can make a track!
  if (fDebug > 1)
    std::cout << " ------------------ Finding initial track -------------------- " << std::endl;
  initialTrack = MakeInitialTrack(initialTrackHits, showerHitsMap);
  if (initialTrack and fDebug > 1) {
    std::cout << "The track start is (" << initialTrack->Vertex().X() << ", " << initialTrack->Vertex().Y() << ", " << initialTrack->Vertex().Z() << ")" << std::endl;
    std::cout << "The track direction is (" << initialTrack->VertexDirection().X() << ", " << initialTrack->VertexDirection().Y() << ", " << initialTrack->VertexDirection().Z() << ")" << std::endl;
  }
  if (fDebug > 1)
    std::cout << " ------------------ End finding initial track -------------------- " << std::endl;

  // // Fill correct or incorrect direction histogram
  // std::map<int,int> trackHits;
  // for (art::PtrVector<recob::Hit>::const_iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt)
  //   ++trackHits[FindTrackID(*hitIt)];
  // int trueTrack = -9999;
  // for (std::map<int,int>::iterator trackHitIt = trackHits.begin(); trackHitIt != trackHits.end(); ++trackHitIt)
  //   if (trackHitIt->second/(double)hits.size() > 0.8)
  //     trueTrack = trackHitIt->first;
  // if (trueTrack != -9999) {
  //   const simb::MCParticle* trueParticle = backtracker->TrackIDToParticle(trueTrack);
  //   TVector3 trueStart = trueParticle->Position().Vect();
  //   if (initialTrack) {
  //     TVector3 recoStart = initialTrack->Vertex();
  //     if ((recoStart-trueStart).Mag() < 5)
  //    hTrueDirection->Fill(1);
  //     else
  //    hTrueDirection->Fill(0);
  //   }
  // }
  // else
  //   hTrueDirection->Fill(0);

  return;

}

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
	)

<Tingjun to="" document>="">

Definition at line 1539 of file EMShowerAlg.cxx.

References geo::vect::coord(), fGeom, geo::GeometryCore::FindTPCAtPosition(), fNfithits, fNfitpass, fToler, HitCoordinates(), geo::CryostatID::isValid, geo::TPCID::TPC, WeightedFit(), and recob::Vertex::XYZ().

Referenced by MakeShower().

                                                                                       {

  // Find TPC for the vertex
  //std::cout<<"here"<<std::endl;
  double xyz[3];
  vertex->XYZ(xyz);
  //std::cout<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<std::endl;
  geo::TPCID tpc = fGeom->FindTPCAtPosition(xyz);
  //std::cout<<tpc<<std::endl;
  //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;
      }
    }
  }
  //std::cout<<tpc<<std::endl;
    //if (!tpc.isValid&&showerHits.size()) tpc = geo::TPCID(showerHits[0]->WireID());
  if (!tpc.isValid) return;
  //std::cout<<"here 1"<<std::endl;

  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){
      //std::cout<<i<<" "<<hit->WireID()<<" "<<tpc<<std::endl;
      if (hit->WireID().TPC==tpc.TPC){
        TVector2 coord = HitCoordinates(hit);
        //std::cout<<i<<" "<<hit->WireID()<<" "<<hit->PeakTime()<<std::endl;
        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(hit->Integral());
          cfit.push_back(1.);
          if (i==fNfitpass-1) {
            trackHits.push_back(hit);
          }
        //std::cout<<*hit<<std::endl;
//
//<<hit->PeakTime()<<" "<<std::abs((coord.Y()-(parm[0]+coord.X()*parm[1]))*cos(atan(parm[1])))<<std::endl;
        }
      }
    }
  
    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 ( std::vector< art::Ptr< recob::Hit > > const &  hits, bool  perpendicular = false  )

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 658 of file EMShowerAlg.cxx.

References fMakeGradientPlot, HitPosition(), hits(), ShowerCentre(), ShowerDirection(), and lar::dump::vector().

Referenced by OrderShowerHits().

                                                                                                                            {

  // Find the charge-weighted centre (in [cm]) of this shower
  TVector2 centre = ShowerCentre(hits);

  // Find a rough shower 'direction'
  TVector2 direction = ShowerDirection(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 (std::vector<art::Ptr<recob::Hit> >::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
    pos = HitPosition(*hit) - centre;
    hitProjection[direction*pos] = *hit;
  }

  // Get a vector of hits in order of the shower
  std::vector<art::Ptr<recob::Hit> > showerHits;
  std::transform(hitProjection.begin(), hitProjection.end(), std::back_inserter(showerHits), [](std::pair<double,art::Ptr<recob::Hit> > const& hit) { return hit.second; });

  // Make gradient plot
  if (fMakeGradientPlot) {
    std::map<int,TGraph*> graphs;
    for (std::vector<art::Ptr<recob::Hit> >::iterator hitIt = showerHits.begin(); hitIt != showerHits.end(); ++hitIt) {
      int tpc = (*hitIt)->WireID().TPC;
      if (graphs[tpc] == nullptr)
        graphs[tpc] = new TGraph();
      graphs[tpc]->SetPoint(graphs[tpc]->GetN(), HitPosition(*hitIt).X(), HitPosition(*hitIt).Y());
      //graphs[tpc]->SetPoint(graphs[tpc]->GetN(), HitCoordinates(*hitIt).X(), HitCoordinates(*hitIt).Y());
    }
    TMultiGraph* multigraph = new TMultiGraph();
    for (std::map<int,TGraph*>::iterator graphIt = graphs.begin(); graphIt != graphs.end(); ++graphIt) {
      graphIt->second->SetMarkerColor(graphIt->first);
      graphIt->second->SetMarkerStyle(8);
      graphIt->second->SetMarkerSize(2);
      multigraph->Add(graphIt->second);
    }
    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;

}

int shower::EMShowerAlg::FindParticleID ( const art::Ptr< recob::Hit > &  hit )

private

Returns the true track ID associated with this hit (if more than one, returns the one with highest energy)

Definition at line 2327 of file EMShowerAlg.cxx.

                                                                 {


  double particleEnergy = 0;
  int likelyTrackID = 0;
  std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(hit);
  for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
    if (trackIDs.at(idIt).energy > particleEnergy) {
      particleEnergy = trackIDs.at(idIt).energy;
      likelyTrackID = TMath::Abs(trackIDs.at(idIt).trackID);
    }
  }

  return likelyTrackID;

}

std::vector< std::vector< int > > shower::EMShowerAlg::FindShowers

(

std::map< int, std::vector< int > > const &

trackToClusters

)

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

Definition at line 711 of file EMShowerAlg.cxx.

References evd::details::begin(), and evd::details::end().

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

                                                                                                            {

  // Showers are vectors of clusters
  std::vector<std::vector<int> > showers;

  // Loop over all tracks 
  for (std::map<int,std::vector<int> >::const_iterator trackToClusterIt = trackToClusters.begin(); trackToClusterIt != trackToClusters.end(); ++ trackToClusterIt) {

    // Find which showers already made are associated with this track
    std::vector<int> matchingShowers;
    for (unsigned int shower = 0; shower < showers.size(); ++shower)
      for (std::vector<int>::const_iterator cluster = trackToClusterIt->second.begin(); cluster != trackToClusterIt->second.end(); ++cluster)
        if ( (std::find(showers.at(shower).begin(), showers.at(shower).end(), *cluster) != showers.at(shower).end()) and
             (std::find(matchingShowers.begin(), matchingShowers.end(), shower)) == matchingShowers.end() )
          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(trackToClusterIt->second);

    // Add to existing shower
    else {
      for (std::vector<int>::const_iterator cluster = trackToClusterIt->second.begin(); cluster != trackToClusterIt->second.end(); ++cluster)
        if (std::find(showers.at(matchingShowers.at(0)).begin(), showers.at(matchingShowers.at(0)).end(), *cluster) == showers.at(matchingShowers.at(0)).end())
          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 )

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 749 of file EMShowerAlg.cxx.

References fDebug, HitCoordinates(), lar::dump::vector(), X, and Y.

Referenced by FindInitialTrack().

                                                                                                                                              {

  std::map<int,std::vector<art::Ptr<recob::Hit> > > initialHitsMap;

  for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator orderedShowerIt = orderedShowerMap.begin(); orderedShowerIt != orderedShowerMap.end(); ++orderedShowerIt) {

    std::vector<art::Ptr<recob::Hit> > initialHits;
    const std::vector<art::Ptr<recob::Hit> > orderedShower = orderedShowerIt->second;

    // Find if the shower is travelling along ticks or wires
    bool wireDirection = true;
    std::vector<int> wires;
    for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = orderedShower.begin(); hitIt != orderedShower.end(); ++hitIt)
      wires.push_back(std::round(HitCoordinates(*hitIt).X()));
    std::sort(wires.begin(), wires.end());
    if (TMath::Abs(*wires.begin()-std::round(HitCoordinates(*orderedShower.begin()).X())) > 3 and
        TMath::Abs(*wires.rbegin()-std::round(HitCoordinates(*orderedShower.begin()).X())) > 3)
      wireDirection = false;

    // Deal with showers travelling 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 (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = orderedShower.begin(); hitIt != orderedShower.end(); ++hitIt)
        wireHitMap[std::round(HitCoordinates(*hitIt).X())].push_back(*hitIt);
      for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = orderedShower.begin(); hitIt != orderedShower.end(); ++hitIt) {
        int wire = std::round(HitCoordinates(*hitIt).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(*hitIt);
          else
            break;
        }
        else
          initialHits.push_back(*hitIt);
      }
      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 (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = orderedShower.begin(); hitIt != orderedShower.end(); ++hitIt) {
        int tick = std::round(HitCoordinates(*hitIt).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(*hitIt);
      }
      if (!initialHits.size()) 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.at(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 (std::vector<art::Ptr<recob::Hit> >::iterator initialHitIt = initialHits.begin(); initialHitIt != initialHits.end(); ++initialHitIt)
      ++tpcHitMap[(*initialHitIt)->WireID().TPC];
    for (std::map<int,int>::iterator tpcIt = tpcHitMap.begin(); tpcIt != tpcHitMap.end(); ++tpcIt)
      if (tpcIt->second == 1) singleHitTPCs.push_back(tpcIt->first);
    if (singleHitTPCs.size()) {
      if (fDebug > 2)
        for (std::vector<int>::iterator tpcIt = singleHitTPCs.begin(); tpcIt != singleHitTPCs.end(); ++tpcIt)
          std::cout << "Removed hits in TPC " << *tpcIt << std::endl;
      for (std::vector<art::Ptr<recob::Hit> >::iterator initialHitIt = initialHits.begin(); initialHitIt != initialHits.end(); ++initialHitIt)
        if (std::find(singleHitTPCs.begin(), singleHitTPCs.end(), (*initialHitIt)->WireID().TPC) == singleHitTPCs.end())
          newInitialHits.push_back(*initialHitIt);
      if (!newInitialHits.size()) newInitialHits.push_back(*orderedShower.begin());
    }
    else
      newInitialHits = initialHits;

    initialHitsMap[orderedShowerIt->first] = newInitialHits;

  }

  return initialHitsMap;

}

int shower::EMShowerAlg::FindTrueParticle ( const std::vector< art::Ptr< recob::Hit > > &  showerHits )

private

Returns the true particle most likely associated with this shower

Definition at line 2301 of file EMShowerAlg.cxx.

References recob::Hit::Integral(), and lar::dump::vector().

                                                                                      {


  // Make a map of the tracks which are associated with this shower and the charge each contributes
  std::map<int,double> trackMap;
  for (std::vector<art::Ptr<recob::Hit> >::const_iterator showerHitIt = showerHits.begin(); showerHitIt != showerHits.end(); ++showerHitIt) {
    art::Ptr<recob::Hit> hit = *showerHitIt;
    int trackID = FindParticleID(hit);
    trackMap[trackID] += hit->Integral();
  }

  // Pick the track with the highest charge as the 'true track'
  double highestCharge = 0;
  int showerTrack = 0;
  for (std::map<int,double>::iterator trackIt = trackMap.begin(); trackIt != trackMap.end(); ++trackIt) {
    if (trackIt->second > highestCharge) {
      highestCharge = trackIt->second;
      showerTrack  = trackIt->first;
    }
  }

  return showerTrack;

}

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

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 845 of file EMShowerAlg.cxx.

References fDebug, for(), if(), RelativeWireWidth(), ShowerHitRMS(), ShowerHitRMSGradient(), and lar::dump::vector().

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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
    planeRMS[showerHitsIt->first] = ShowerHitRMS(showerHitsIt->second);
    planeRMSGradients[showerHitsIt->first] = ShowerHitRMSGradient(showerHitsIt->second);
  }

  // Order these backwards so they can be used to discriminate between planes
  std::map<double,int> gradientMap;
  for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt)
    gradientMap[TMath::Abs(planeRMSGradients.at(showerHitsIt->first))] = showerHitsIt->first;

  std::map<double,int> wireWidthMap = RelativeWireWidth(showerHitsMap);

  if (fDebug > 1)
    for (std::map<double,int>::const_iterator wireWidthIt = wireWidthMap.begin(); wireWidthIt != wireWidthMap.end(); ++wireWidthIt)
      std::cout << "Plane " << wireWidthIt->second << " has relative width in wire of " << wireWidthIt->first << "; and an RMS gradient of " << planeRMSGradients[wireWidthIt->second] << std::endl;

  // Find the correct permutations
  int perm = 0;
  std::vector<std::map<int,bool> > usedPermutations;

  // Most likely is to not change anything
  for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt)
    permutations[perm][showerHitsIt->first] = 0;
  ++perm;

  // Use properties of the shower to determine the middle cases
  for (std::map<double,int>::iterator wireWidthIt = wireWidthMap.begin(); wireWidthIt != wireWidthMap.end(); ++wireWidthIt) {
    std::map<int,bool> permutation;
    permutation[wireWidthIt->second] = 1;
    for (std::map<double,int>::iterator wireWidth2It = wireWidthMap.begin(); wireWidth2It != wireWidthMap.end(); ++wireWidth2It)
      if (wireWidthIt->second != wireWidth2It->second)
        permutation[wireWidth2It->second] = 0;
    if (std::find(usedPermutations.begin(), usedPermutations.end(), permutation) == usedPermutations.end()) {
      permutations[perm] = permutation;
      usedPermutations.push_back(permutation);
      ++perm;
    }
  }
  for (std::map<double,int>::reverse_iterator wireWidthIt = wireWidthMap.rbegin(); wireWidthIt != wireWidthMap.rend(); ++wireWidthIt) {
    std::map<int,bool> permutation;
    permutation[wireWidthIt->second] = 0;
    for (std::map<double,int>::iterator wireWidth2It = wireWidthMap.begin(); wireWidth2It != wireWidthMap.end(); ++wireWidth2It)
      if (wireWidthIt->second != wireWidth2It->second)
        permutation[wireWidth2It->second] = 1;
    if (std::find(usedPermutations.begin(), usedPermutations.end(), permutation) == usedPermutations.end()) {
      permutations[perm] = permutation;
      usedPermutations.push_back(permutation);
      ++perm;
    }
  }

  // Least likely is to change everything
  for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt)
    permutations[perm][showerHitsIt->first] = 1;
  ++perm;

  if (fDebug > 1) {
    std::cout << "Here are the permutations!" << std::endl;
    for (std::map<int,std::map<int,bool> >::iterator permIt = permutations.begin(); permIt != permutations.end(); ++permIt) {
      std::cout << "Permutation " << permIt->first << std::endl;
      for (std::map<int,bool>::iterator planePermIt = permIt->second.begin(); planePermIt != permIt->second.end(); ++planePermIt)
        std::cout << "  Plane " << planePermIt->first << " has value " << planePermIt->second << std::endl;
    }
  }

  return permutations;

}

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

private

Find the global wire position.

Definition at line 1633 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) {
    double wireCentre[3];
    fGeom->WireIDToWireGeo(wireID).GetCenter(wireCentre);
    if (wireID.TPC % 2 == 0) globalWire = fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 0, wireID.Cryostat);
    else globalWire = fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 1, wireID.Cryostat);
  }

  // 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(wireID.Plane, 0, wireID.Cryostat);
      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(wireID.Plane, 0, wireID.Cryostat);
      // 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 {
      double wireCentre[3];
      fGeom->WireIDToWireGeo(wireID).GetCenter(wireCentre);
      if (wireID.TPC % 2 == 0) globalWire = fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 0, wireID.Cryostat);
      else globalWire = fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 1, wireID.Cryostat);
    }
  }

  return globalWire;

}

TVector2 shower::EMShowerAlg::HitCoordinates ( art::Ptr< recob::Hit > const &  hit )

private

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

Definition at line 1612 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 ( art::Ptr< recob::Hit > const &  hit )

private

Return the coordinates of this hit in units of cm.

Definition at line 1618 of file EMShowerAlg.cxx.

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

Referenced by CheckShowerHits(), FindOrderOfHits(), MakeSpacePoints(), OrderShowerHits(), Project3DPointOntoPlane(), ShowerCentre(), ShowerDirection(), ShowerHitRMS(), and ShowerHitRMSGradient().

                                                                   {

  geo::PlaneID planeID = hit->WireID().planeID();

  return HitPosition(HitCoordinates(hit), planeID);

}

TVector2 shower::EMShowerAlg::HitPosition ( TVector2 const &  pos, geo::PlaneID  planeID  )

private

Return the coordinates of this hit in units of cm.

Definition at line 1626 of file EMShowerAlg.cxx.

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

                                                                               {

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

}

bool shower::EMShowerAlg::isCleanShower

(

std::vector< art::Ptr< recob::Hit > > const &

hits

)

<Tingjun to="" document>="">

Definition at line 2199 of file EMShowerAlg.cxx.

References hits().

                                                                             {

  if (!hits.size()) 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){
    ++nhits;
    if (nhits>2)
      ++hitmap[hit->WireID().Wire];
    if (nhits==20) break;
  }
  //std::cout<<hits.size()<<" "<<float(nhits-2)/hitmap.size()<<std::endl;
  if (float(nhits-2)/hitmap.size()>1.4) return false;
  else return true;
}

std::unique_ptr< recob::Track > shower::EMShowerAlg::MakeInitialTrack ( std::map< int, std::vector< art::Ptr< recob::Hit > > > const &  initialHitsMap, std::map< int, std::vector< art::Ptr< recob::Hit > > > const &  showerHitsMap  )

private

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

Definition at line 928 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" << std::endl;
    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(showerHitsMap)) {
    int planeToIgnore = WorstPlane(showerHitsMap);
    if (fDebug > 1)
      std::cout << "Igoring plane " << planeToIgnore << " in creation of initial track" << std::endl;
    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(initialHitsMap.at(planes.at(0)), initialHitsMap.at(planes.at(1)));

}

void shower::EMShowerAlg::MakePicture ( )

private

recob::Shower shower::EMShowerAlg::MakeShower	(	art::PtrVector< recob::Hit > const &	hits,
		std::unique_ptr< recob::Track > const &	initialTrack,
		std::map< int, std::vector< art::Ptr< recob::Hit > > > const &	initialTrackHits
	)

Makes a recob::Shower object given the hits in the shower and the initial track-like part.

Definition at line 968 of file EMShowerAlg.cxx.

References art::PtrVector< T >::begin(), art::PtrVector< T >::end(), fDebug, fGeom, FinddEdx(), fShowerEnergyAlg, geo::GeometryCore::MaxPlanes(), shower::ShowerEnergyAlg::ShowerEnergy(), recob::Track::Vertex(), and recob::Track::VertexDirection().

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

                                                                                                               {

  //return recob::Shower();

  // Find the shower hits on each plane
  std::map<int,std::vector<art::Ptr<recob::Hit> > > planeHitsMap;
  for (art::PtrVector<recob::Hit>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit)
    planeHitsMap[(*hit)->View()].push_back(*hit);

  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(initialHitsMap.at(plane), initialTrack));
      totalEnergy.push_back(fShowerEnergyAlg.ShowerEnergy(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();
    showerStart = initialTrack->Vertex();
  }

  if (fDebug > 0) {
    std::cout << "Best plane is " << bestPlane << std::endl;
    std::cout << "dE/dx for each plane is: " << dEdx[0] << ", " << dEdx[1] << " and " << dEdx[2] << std::endl;
    std::cout << "Total energy for each plane is: " << totalEnergy[0] << ", " << totalEnergy[1] << " and " << totalEnergy[2] << std::endl;
    std::cout << "The shower start is (" << showerStart.X() << ", " << showerStart.Y() << ", " << showerStart.Z() << ")" << std::endl;
    std::cout << "The shower direction is (" << direction.X() << ", " << direction.Y() << ", " << direction.Z() << ")" << std::endl;
  }

  return recob::Shower(direction, directionError, showerStart, showerStartError, totalEnergy, totalEnergyError, dEdx, dEdxError, bestPlane);

}

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

<Tingjun to="" document>="">

Definition at line 1022 of file EMShowerAlg.cxx.

References art::PtrVector< T >::begin(), pma::ProjectionMatchingAlg::buildSegment(), calo::CalorimetryAlg::dEdx_AREA(), art::PtrVector< T >::end(), fCalorimetryAlg, fDebug, fdEdxTrackLength, fGeom, FindInitialTrackHits(), fProjectionMatchingAlg, fShowerEnergyAlg, geo::GeometryCore::MaxPlanes(), OrderShowerHits(), geo::GeometryCore::Plane(), shower::ShowerEnergyAlg::ShowerEnergy(), pma::Track3D::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 (art::PtrVector<recob::Hit>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit)
    planeHitsMap[(*hit)->WireID().Plane].push_back(*hit);

  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(planeHits->second, showerHits, vertex);
    //if (!isCleanShower(showerHits)) continue;
    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();
    }

  }
  //std::cout<<pl0<<" "<<pl1<<std::endl;
//  if (pl0!=-1&&pl1!=-1) {
//    pl0 = 1;
//    pl1 = 2;
//  }
  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);
//    std::vector<art::Ptr<recob::Hit>> alltrackhits;
//    for (size_t i = 0; i<3; ++i){
//      for (auto const&hit : initialTrackHits[i]){
//      alltrackhits.push_back(hit);
//      }
//    }
    //std::cout<<"vertex "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<std::endl;
    //for (auto const&hit : initialTrackHits[pl0]) std::cout<<*hit<<std::endl;
    //for (auto const&hit : initialTrackHits[pl1]) std::cout<<*hit<<std::endl;
    pma::Track3D* pmatrack = fProjectionMatchingAlg.buildSegment(initialTrackHits[pl0], initialTrackHits[pl1]);
    //std::cout<<pmatrack->size()<<std::endl;
    //pma::Track3D* pmatrack = fProjectionMatchingAlg.buildSegment(alltrackhits);
    std::vector<TVector3> spts;

//  points are shifted now by tracking code, commented out shifts here
//    double xshift = pmatrack->GetXShift();
//    bool has_shift = (xshift != 0.0);
    for (size_t i = 0; i<pmatrack->size(); ++i){
      if ((*pmatrack)[i]->IsEnabled()){
        TVector3 p3d = (*pmatrack)[i]->Point3D();
//      if (has_shift) p3d.SetX(p3d.X() + xshift);
        //std::cout<<p3d.X()<<" "<<p3d.Y()<<" "<<p3d.Z()<<std::endl;
        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();
      }
      //std::cout<<shwxyz.X()<<" "<<shwxyz.Y()<<" "<<shwxyz.Z()<<std::endl;
      //std::cout<<shwdir.X()<<" "<<shwdir.Y()<<" "<<shwdir.Z()<<std::endl;
      for (unsigned int plane = 0; plane < fGeom->MaxPlanes(); ++plane) {
        if (planeHitsMap.find(plane)!=planeHitsMap.end()){
          totalEnergy.push_back(fShowerEnergyAlg.ShowerEnergy(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(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::cout<<"pitch = "<<pitch<<std::endl;
            std::vector<float> vQ;
            for (auto const& hit: initialTrackHits[plane]){
              //std::cout<<hit->WireID()<<" "<<hit->PeakTime()<<" "<<std::abs((hit->WireID().Wire-initialTrackHits[plane][0]->WireID().Wire)*pitch)<<" "<<fdEdxTrackLength<<std::endl;
              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;
                //std::cout<<hit->WireID()<<" "<<hit->PeakTime()<<" "<<hit->Integral()<<" "<<totQ<<" "<<avgT<<std::endl;
              }
            }
            if (totQ) {
              //double dQdx = totQ/(nhits*pitch);
              //std::sort(vQ.begin(), vQ.end());
              //double dQdx = vQ[vQ.size()/2]/pitch;
              double dQdx = TMath::Median(vQ.size(), &vQ[0])/pitch;
              fdEdx = fCalorimetryAlg.dEdx_AREA(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::endl;
        std::cout << "dE/dx for each plane is: " << dEdx[0] << ", " << dEdx[1] << " and " << dEdx[2] << std::endl;
        std::cout << "Total energy for each plane is: " << totalEnergy[0] << ", " << totalEnergy[1] << " and " << totalEnergy[2] << std::endl;
        std::cout << "The shower start is (" << shwxyz.X() << ", " << shwxyz.Y() << ", " << shwxyz.Z() << ")" << std::endl;
        shwxyz.Print();
      }

      return recob::Shower(shwdir, shwdirerr, shwxyz, shwxyzerr, totalEnergy, totalEnergyError, dEdx, dEdxError, bestPlane);
    }
  }
  return recob::Shower();
}

std::vector< recob::SpacePoint > shower::EMShowerAlg::MakeSpacePoints	(	std::map< int, std::vector< art::Ptr< recob::Hit > > >	hits,
		std::vector< std::vector< art::Ptr< recob::Hit > > > &	hitAssns
	)

Makes space points from the shower hits in each plane.

Definition at line 1187 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;

  // // Order by charge
  // for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator planeIt = showerHits.begin(); planeIt != showerHits.end(); ++planeIt)
  //   std::sort(planeIt->second.begin(), planeIt->second.end(), [](const art::Ptr<recob::Hit>& h1, const art::Ptr<recob::Hit>& h2) { return h1->Integral() > h2->Integral(); });

  // 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;

        TVector3 point = Construct3DPoint(*planeHitIt, *otherPlaneHitIt);
        std::vector<art::Ptr<recob::Hit> > pointHits;
        bool truePoint = false;

        if (otherPlanes.size() > 1) {

          TVector2 projThirdPlane = Project3DPointOntoPlane(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(*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(point, (*planeHitIt)->WireID().Plane) - HitPosition(*planeHitIt)).Mod() < fSpacePointSize and
                 (Project3DPointOntoPlane(point, (*otherPlaneHitIt)->WireID().Plane) - HitPosition(*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 -------------------" << std::endl;
    std::cout << "There are " << spacePoints.size() << " space points:" << std::endl;
    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] << ")" << std::endl;
      }
  }

  return spacePoints;

}

std::map< int, std::vector< art::Ptr< recob::Hit > > > shower::EMShowerAlg::OrderShowerHits	(	const art::PtrVector< recob::Hit > &	shower,
		int	plane
	)

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 1312 of file EMShowerAlg.cxx.

References art::PtrVector< T >::begin(), CheckIsolatedHits(), CheckShowerHits(), art::PtrVector< T >::end(), fDebug, fGeom, FindOrderOfHits(), GetPlanePermutations(), HitCoordinates(), HitPosition(), hits(), geo::GeometryCore::MaxPlanes(), n, RelativeWireWidth(), ShowerHitRMS(), ShowerHitRMSGradient(), lar::dump::vector(), X, and Y.

Referenced by MakeShower(), and shower::EMShower::produce().

                                                                                                                            {


  // ------------- Put hits in order ------------

  // // Find the true start for reconstruction validation purposes
  // std::vector<art::Ptr<recob::Hit> > showerHitsVector;
  // std::map<int,int> trueParticleHits;
  // for (art::PtrVector<recob::Hit>::const_iterator showerHitIt = shower.begin(); showerHitIt != shower.end(); ++showerHitIt) {
  //   showerHitsVector.push_back(*showerHitIt);
  //   ++trueParticleHits[FindParticleID(*showerHitIt)];
  // }
  // int trueParticleID = FindTrueParticle(showerHitsVector);
  // const simb::MCParticle* trueParticle = bt_serv->TrackIDToParticle(trueParticleID);
  // TVector3 trueStart3D = trueParticle->Position().Vect();
  // double purity = trueParticleHits[trueParticleID]/(double)shower.size();

  // Find the shower hits on each plane
  std::map<int,std::vector<art::Ptr<recob::Hit> > > showerHitsMap;
  for (art::PtrVector<recob::Hit>::const_iterator hit = shower.begin(); hit != shower.end(); ++hit)
    showerHitsMap[(*hit)->WireID().Plane].push_back(*hit);

  // if (purity < 0.9)
  //   return showerHitsMap;

  // Order the hits, get the RMS and the RMS gradient for the hits in this plane
  std::map<int,double> planeRMSGradients, planeRMS;
  for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
    if (plane != showerHitsIt->first and plane != -1)
      continue;
    std::vector<art::Ptr<recob::Hit> > orderedHits = FindOrderOfHits(showerHitsIt->second);
    planeRMS[showerHitsIt->first] = ShowerHitRMS(orderedHits);
    //TVector2 trueStart2D = Project3DPointOntoPlane(trueStart3D, showerHitsIt->first);
    planeRMSGradients[showerHitsIt->first] = ShowerHitRMSGradient(orderedHits);
    showerHitsMap[showerHitsIt->first] = orderedHits;
  }

  if (fDebug > 1)
    for (std::map<int,double>::iterator planeRMSIt = planeRMS.begin(); planeRMSIt != planeRMS.end(); ++planeRMSIt)
      std::cout << "Plane " << planeRMSIt->first << " has RMS " << planeRMSIt->second << " and RMS gradient " << planeRMSGradients.at(planeRMSIt->first) << std::endl;

  if (fDebug > 2) {
    std::cout << std::endl << "Hits in order; after ordering, before reversing..." << std::endl;
    for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
      std::cout << "  Plane " << showerHitsIt->first << std::endl;
      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(*hitIt).X() << ", " << HitPosition(*hitIt).Y() << ")" << std::endl;
    }
  }

  // ------------- 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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
    if (planeRMS.at(showerHitsIt->first) > planeOtherRMS.at(showerHitsIt->first) * 2) {
      //and TMath::Abs(planeRMSGradients.at(showerHitsIt->first) / planeOtherRMSGradients.at(showerHitsIt->first)) < 0.1) {
      if (fDebug > 1)
        std::cout << "Plane " << showerHitsIt->first << " was perpendicular... recalculating" << std::endl;
      std::vector<art::Ptr<recob::Hit> > orderedHits = this->FindOrderOfHits(showerHitsIt->second, true);
      showerHitsMap[showerHitsIt->first] = orderedHits;
      planeRMSGradients[showerHitsIt->first] = this->ShowerHitRMSGradient(orderedHits);
    }
  }

  // ------------- Orient the shower correctly ---------------

  if (fDebug > 1) {
    std::cout << "Before reversing, here are the start and end points..." << std::endl;
    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() << ") (real wire " << showerHitsIt->second.back()->WireID() << ")" << std::endl;
  }

  // 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 << std::endl << "Hits in order; after reversing, before checking isolated hits..." << std::endl;
    for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
      std::cout << "  Plane " << showerHitsIt->first << std::endl;
      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(*hitIt).X() << ", " << HitPosition(*hitIt).Y() << ")" << std::endl;
    }
  }

  CheckIsolatedHits(showerHitsMap);

  if (fDebug > 2) {
    std::cout << std::endl << "Hits in order; after checking isolated hits..." << std::endl;
    for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
      std::cout << "  Plane " << showerHitsIt->first << std::endl;
      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(*hitIt).X() << ", " << HitPosition(*hitIt).Y() << ")" << std::endl;
    }
  }

  if (fDebug > 1) {
    std::cout << "After reversing and checking isolated hits, here are the start and end points..." << std::endl;
    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() << ")" << std::endl;
  }

  // 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... " << std::endl;
  for (std::map<double,int>::iterator wireWidthIt = wireWidths.begin(); wireWidthIt != wireWidths.end(); ++wireWidthIt) {
    if (fDebug > 1)
      std::cout << "Plane " << wireWidthIt->second << " has relative wire width " << wireWidthIt->first << std::endl;
    if (wireWidthIt->first < 1/(double)wireWidths.size())
      badPlanes.push_back(wireWidthIt->second);
  }

  std::map<int,std::vector<art::Ptr<recob::Hit> > > ignoredPlanes;
  if (badPlanes.size() == 1) {
    int badPlane = badPlanes.at(0);
    if (fDebug > 1)
      std::cout << "Ignoring plane " << badPlane << " when orientating" << std::endl;
    ignoredPlanes[badPlane] = showerHitsMap.at(badPlane);
    showerHitsMap.erase(showerHitsMap.find(badPlane));
  }

  // Consider all possible permutations of planes (0/1, oriented correctly/incorrectly)
  std::map<int,std::map<int,bool> > permutations = GetPlanePermutations(showerHitsMap);

  // Go through all permutations until we have a satisfactory orientation
  std::map<int,std::vector<art::Ptr<recob::Hit> > > originalShowerHitsMap = showerHitsMap;
  int n = 0;
  while (!CheckShowerHits(showerHitsMap) and n < TMath::Power(2,(int)showerHitsMap.size())) {
    if (fDebug > 1)
      std::cout << "Permutation " << n << std::endl;
    for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
      std::vector<art::Ptr<recob::Hit> > hits = originalShowerHitsMap.at(showerHitsIt->first);
      if (permutations[n][showerHitsIt->first] == 1) {
        std::reverse(hits.begin(), hits.end());
        showerHitsMap[showerHitsIt->first] = hits;
      }
      else
        showerHitsMap[showerHitsIt->first] = hits;
    }
    ++n;      
  }

  // Go back to original if still no consistency
  if (!CheckShowerHits(showerHitsMap))
    showerHitsMap = originalShowerHitsMap;

  if (fDebug > 2) {
    std::cout << "End of OrderShowerHits: here are the order of hits:" << std::endl;
    for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator planeHitsIt = showerHitsMap.begin(); planeHitsIt != showerHitsMap.end(); ++planeHitsIt) {
      std::cout << "  Plane " << planeHitsIt->first << std::endl;
      for (std::vector<art::Ptr<recob::Hit> >::iterator hitIt = planeHitsIt->second.begin(); hitIt != planeHitsIt->second.end(); ++hitIt)
        std::cout << "    Hit (" << HitCoordinates(*hitIt).X() << " (real wire " << (*hitIt)->WireID() << "), " << HitCoordinates(*hitIt).Y() << ") -- pos (" << HitPosition(*hitIt).X() << ", " << HitPosition(*hitIt).Y() << ")" << std::endl;
    }
  }

  if (ignoredPlanes.size())
    showerHitsMap[ignoredPlanes.begin()->first] = ignoredPlanes.begin()->second;

  return showerHitsMap;

}

void shower::EMShowerAlg::OrderShowerHits ( std::vector< art::Ptr< recob::Hit > > const &  shower, std::vector< art::Ptr< recob::Hit > > &  orderedShower, art::Ptr< recob::Vertex > const &  vertex  )

private

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

Definition at line 1502 of file EMShowerAlg.cxx.

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

                                                                              {

  showerHits = FindOrderOfHits(shower);

  // Find TPC for the vertex
  double xyz[3];
  vertex->XYZ(xyz);
  geo::TPCID tpc = fGeom->FindTPCAtPosition(xyz);
  if (!tpc.isValid&&showerHits.size()) tpc = geo::TPCID(showerHits[0]->WireID());
  //std::cout<<tpc<<std::endl;
  // 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[1], xyz[2], hit0->WireID().planeID()),
                               fDetProp->ConvertXToTicks(xyz[0],  hit0->WireID().planeID()));
//  std::cout<<coord0.X()<<" "<<coord0.Y()<<std::endl;
//  std::cout<<coord1.X()<<" "<<coord1.Y()<<std::endl;
//  std::cout<<coordvtx.X()<<" "<<coordvtx.Y()<<std::endl;
//  std::cout<<hit0->WireID()<<" "<<hit1->WireID()<<std::endl;
  if ((coord1-coordvtx).Mod()<(coord0-coordvtx).Mod()){
    std::reverse(showerHits.begin(), showerHits.end());
  }
  //std::cout<<showerHits[0]->WireID()<<" "<<showerHits.back()->WireID()<<std::endl;
}

TVector2 shower::EMShowerAlg::Project3DPointOntoPlane ( TVector3 const &  point, int  plane, int  cryostat = 0  )

private

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

Definition at line 2093 of file EMShowerAlg.cxx.

References detinfo::DetectorProperties::ConvertXToTicks(), e, fDetProp, 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.);

  double pointPosition[3] = {point.X(), point.Y(), point.Z()};

  geo::TPCID tpcID = fGeom->FindTPCAtPosition(pointPosition);
  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),
                         fDetProp->ConvertXToTicks(point.X(), planeID));

  // wireTickPos = TVector2(fGeom->WireCoordinate(point.Y(), point.Z(), planeID.Plane, tpc % 2, planeID.Cryostat),
  //                     fDetProp->ConvertXToTicks(point.X(), planeID.Plane, tpc % 2, planeID.Cryostat));

  //return wireTickPos;
  return HitPosition(wireTickPos, planeID);

}

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

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 1674 of file EMShowerAlg.cxx.

References 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] = TMath::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 += TMath::Power(planeWireLength[plane],2);
    quad = TMath::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::ShowerCentre ( const std::vector< art::Ptr< recob::Hit > > &  showerHits )

private

Returns the charge-weighted shower centre.

Definition at line 1726 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(*hit);
    chargePoint += (*hit)->Integral() * pos;
    totalCharge += (*hit)->Integral();
  }
  TVector2 centre = chargePoint / totalCharge;

  return centre;

}

TVector2 shower::EMShowerAlg::ShowerDirection ( const std::vector< art::Ptr< recob::Hit > > &  showerHits )

private

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

Definition at line 1701 of file EMShowerAlg.cxx.

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

Referenced by FindOrderOfHits(), ShowerHitRMS(), and ShowerHitRMSGradient().

                                                                                          {

  TVector2 pos;
  //double weight;
  double weight = 1;
  //int nhits = 0;
  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(*hit);
    weight = TMath::Power((*hit)->Integral(),2);
    sumweight += weight;
    sumx += weight * pos.X();
    sumy += weight * pos.Y();
    sumx2 += weight * pos.X() * pos.X();
    sumxy += weight * pos.X() * pos.Y();
  }
  //double gradient = (nhits * sumxy - sumx * sumy) / (nhits * sumx2 - sumx * sumx);
  double gradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);
  TVector2 direction = TVector2(1,gradient).Unit();

  return direction;

}

double shower::EMShowerAlg::ShowerHitRMS ( const std::vector< art::Ptr< recob::Hit > > &  showerHits )

private

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

Definition at line 1741 of file EMShowerAlg.cxx.

References HitPosition(), proj, ShowerCentre(), ShowerDirection(), and lar::dump::vector().

Referenced by GetPlanePermutations(), and OrderShowerHits().

                                                                                     {

  TVector2 direction = ShowerDirection(showerHits);
  TVector2 centre = ShowerCentre(showerHits);

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

  return RMS;

}

double shower::EMShowerAlg::ShowerHitRMSGradient ( const std::vector< art::Ptr< recob::Hit > > &  showerHits, TVector2  trueStart = TVector2(0,0)  )

private

Returns the gradient of the RMS vs shower segment graph.

Definition at line 1757 of file EMShowerAlg.cxx.

References fMakeRMSGradientPlot, fNumShowerSegments, HitPosition(), proj, ShowerDirection(), lar::dump::vector(), and weight.

Referenced by GetPlanePermutations(), and OrderShowerHits().

                                                                                                                 {

  // Don't forget to clean up the header file!

  // Find a rough shower 'direction' and centre
  TVector2 direction = ShowerDirection(showerHits);

  //std::cout << std::endl << "Number of hits in this shower: " << showerHits.size() << std::endl;

  // for (int i = 0; i < 100; ++i) {

  //   // Bin the hits into discreet chunks
  //   //int nSegmentHits = fNumSegmentHits;
  //   int nSegmentHits = i;
  //   int nShowerSegments = showerHits.size() / (double)nSegmentHits;
  //   //int nShowerSegments = fNumShowerSegments;
  //   double lengthOfShower = (HitPosition(showerHits.back()) - HitPosition(showerHits.front())).Mod();
  //   double lengthOfSegment = lengthOfShower / (double)nShowerSegments;
  //   std::map<int,std::vector<art::Ptr<recob::Hit> > > showerSegments;
  //   std::map<int,double> segmentCharge;
  //   for (std::vector<art::Ptr<recob::Hit> >::const_iterator showerHitsIt = showerHits.begin(); showerHitsIt != showerHits.end(); ++showerHitsIt) {
  //     showerSegments[(int)(HitPosition(*showerHitsIt)-HitPosition(showerHits.front())).Mod() / lengthOfSegment].push_back(*showerHitsIt);
  //     segmentCharge[(int)(HitPosition(*showerHitsIt)-HitPosition(showerHits.front())).Mod() / lengthOfSegment] += (*showerHitsIt)->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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerSegmentIt = showerSegments.begin(); showerSegmentIt != showerSegments.end(); ++showerSegmentIt) {

  //     // Get the mean position of the hits in this bin
  //     TVector2 meanPosition(0,0);
  //     for (std::vector<art::Ptr<recob::Hit> >::iterator hitInSegmentIt = showerSegmentIt->second.begin(); hitInSegmentIt != showerSegmentIt->second.end(); ++hitInSegmentIt)
  //    meanPosition += HitPosition(*hitInSegmentIt);
  //     meanPosition /= (double)showerSegmentIt->second.size();

  //     // Get the RMS of this bin
  //     std::vector<double> distanceToAxisBin;
  //     for (std::vector<art::Ptr<recob::Hit> >::iterator hitInSegmentIt = showerSegmentIt->second.begin(); hitInSegmentIt != showerSegmentIt->second.end(); ++hitInSegmentIt) {
  //    TVector2 proj = (HitPosition(*hitInSegmentIt) - meanPosition).Proj(direction) + meanPosition;
  //    distanceToAxisBin.push_back((HitPosition(*hitInSegmentIt) - proj).Mod());
  //     }

  //     double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
  //     binVsRMS.push_back(std::make_pair(showerSegmentIt->first, RMSBin));
  //     if (fMakeRMSGradientPlot)
  //    graph->SetPoint(graph->GetN(), showerSegmentIt->first, RMSBin);

  //   }

  //   // Get the gradient of the RMS-bin plot

  //   // OLD
  //   int nhits1 = 0;
  //   double sumx1=0., sumy1=0., sumx21=0., sumxy1=0.;
  //   for (std::vector<std::pair<int,double> >::iterator binVsRMSIt = binVsRMS.begin(); binVsRMSIt != binVsRMS.end(); ++binVsRMSIt) {
  //     ++nhits1;
  //     sumx1 += binVsRMSIt->first;
  //     sumy1 += binVsRMSIt->second;
  //     sumx21 += binVsRMSIt->first * binVsRMSIt->first;
  //     sumxy1 += binVsRMSIt->first * binVsRMSIt->second;
  //   }
  //   double RMSgradient1 = (nhits1 * sumxy1 - sumx1 * sumy1) / (nhits1 * sumx21 - sumx1 * sumx1);

  //   // NEW
  //   double sumx=0., sumy=0., sumx2=0., sumxy=0., sumweight = 0.;
  //   for (std::vector<std::pair<int,double> >::iterator binVsRMSIt = binVsRMS.begin(); binVsRMSIt != binVsRMS.end(); ++binVsRMSIt) {
  //     //double weight = showerSegments.at(binVsRMSIt->first).size();
  //     double weight = TMath::Power(segmentCharge.at(binVsRMSIt->first),2);
  //     sumweight += weight;
  //     sumx += weight * binVsRMSIt->first;
  //     sumy += weight * binVsRMSIt->second;
  //     sumx2 += weight * binVsRMSIt->first * binVsRMSIt->first;
  //     sumxy += weight * binVsRMSIt->first * binVsRMSIt->second;
  //   }
  //   double RMSgradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);

  //   // PCA
  //   TPrincipal* pca = new TPrincipal(2,"");
  //   double point[2];
  //   for (std::vector<std::pair<int,double> >::iterator binVsRMSIt = binVsRMS.begin(); binVsRMSIt != binVsRMS.end(); ++binVsRMSIt) {
  //     if (fDebug > 2)
  //    std::cout << "Number of hits in segment " << binVsRMSIt->first << " is " << showerSegments.at(binVsRMSIt->first).size() << std::endl;
  //     for (int nhit = 0; nhit < (int)showerSegments.at(binVsRMSIt->first).size(); ++nhit) {
  //    point[0] = binVsRMSIt->first;
  //    point[1] = binVsRMSIt->second;
  //    pca->AddRow(point);
  //     }
  //   }
  //   pca->MakePrincipals();
  //   TVector2 RMSgradientPCA = TVector2( (*pca->GetEigenVectors())[0][0], (*pca->GetEigenVectors())[1][0] ).Unit();

  //   if (fDebug > 1)
  //     std::cout << "RMS gradient without weighting: " << RMSgradient1 << ", and with weighting: " << RMSgradient << "; PCA with weighting: " << RMSgradientPCA.Y()/RMSgradientPCA.X() << std::endl;

  //   if (fMakeRMSGradientPlot) {
  //     TVector2 direction = TVector2(1,RMSgradient).Unit();
  //     TVector2 direction1 = TVector2(1,RMSgradient1).Unit();
  //     TCanvas* canv = new TCanvas();
  //     graph->Draw();
  //     TVector2 centre = TVector2(graph->GetMean(1), graph->GetMean(2));
  //     TLine line, line1, line2;
  //     line.SetLineColor(3);
  //     line1.SetLineColor(4);
  //     line2.SetLineColor(5);
  //     line.DrawLine(centre.X()-1000*direction.X(),centre.Y()-1000*direction.Y(),centre.X()+1000*direction.X(),centre.Y()+1000*direction.Y());
  //     line1.DrawLine(centre.X()-1000*direction1.X(),centre.Y()-1000*direction1.Y(),centre.X()+1000*direction1.X(),centre.Y()+1000*direction1.Y());
  //     line2.DrawLine(centre.X()-1000*RMSgradientPCA.X(),centre.Y()-1000*RMSgradientPCA.Y(),centre.X()+1000*RMSgradientPCA.X(),centre.Y()+1000*RMSgradientPCA.Y());
  //     canv->SaveAs("RMSGradient.png");
  //     delete canv;
  //   }
  //   delete graph; delete pca;

  //   //std::cout << "Number of hits per bin " << i << " gives gradient of " << RMSgradient << std::endl;

  //   double distanceFromStart = 0, distanceFromEnd = 0;
  //   if (RMSgradient > 0) {
  //     distanceFromStart = (trueStart - HitPosition(showerHits.front())).Mod();
  //     distanceFromEnd = (trueStart - HitPosition(showerHits.back())).Mod();
  //   }
  //   if (RMSgradient < 0) {
  //    distanceFromStart = (trueStart - HitPosition(showerHits.back())).Mod();
  //    distanceFromEnd = (trueStart - HitPosition(showerHits.front())).Mod();
  //   }
  //   bool correct = distanceFromStart < distanceFromEnd;
  //   hNumHitsInSegment->Fill(i,correct);

  // }




  // for (int i = 0; i < 100; ++i) {

  //   // Bin the hits into discreet chunks
  //   int nShowerSegments = i;
  //   double lengthOfShower = (HitPosition(showerHits.back()) - HitPosition(showerHits.front())).Mod();
  //   double lengthOfSegment = lengthOfShower / (double)nShowerSegments;
  //   std::map<int,std::vector<art::Ptr<recob::Hit> > > showerSegments;
  //   std::map<int,double> segmentCharge;
  //   for (std::vector<art::Ptr<recob::Hit> >::const_iterator showerHitsIt = showerHits.begin(); showerHitsIt != showerHits.end(); ++showerHitsIt) {
  //     showerSegments[(int)(HitPosition(*showerHitsIt)-HitPosition(showerHits.front())).Mod() / lengthOfSegment].push_back(*showerHitsIt);
  //     segmentCharge[(int)(HitPosition(*showerHitsIt)-HitPosition(showerHits.front())).Mod() / lengthOfSegment] += (*showerHitsIt)->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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerSegmentIt = showerSegments.begin(); showerSegmentIt != showerSegments.end(); ++showerSegmentIt) {

  //     // Get the mean position of the hits in this bin
  //     TVector2 meanPosition(0,0);
  //     for (std::vector<art::Ptr<recob::Hit> >::iterator hitInSegmentIt = showerSegmentIt->second.begin(); hitInSegmentIt != showerSegmentIt->second.end(); ++hitInSegmentIt)
  //    meanPosition += HitPosition(*hitInSegmentIt);
  //     meanPosition /= (double)showerSegmentIt->second.size();

  //     // Get the RMS of this bin
  //     std::vector<double> distanceToAxisBin;
  //     for (std::vector<art::Ptr<recob::Hit> >::iterator hitInSegmentIt = showerSegmentIt->second.begin(); hitInSegmentIt != showerSegmentIt->second.end(); ++hitInSegmentIt) {
  //    TVector2 proj = (HitPosition(*hitInSegmentIt) - meanPosition).Proj(direction) + meanPosition;
  //    distanceToAxisBin.push_back((HitPosition(*hitInSegmentIt) - proj).Mod());
  //     }

  //     double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
  //     binVsRMS.push_back(std::make_pair(showerSegmentIt->first, RMSBin));
  //     if (fMakeRMSGradientPlot)
  //    graph->SetPoint(graph->GetN(), showerSegmentIt->first, RMSBin);

  //   }

  //   // Get the gradient of the RMS-bin plot

  //   // OLD
  //   int nhits1 = 0;
  //   double sumx1=0., sumy1=0., sumx21=0., sumxy1=0.;
  //   for (std::vector<std::pair<int,double> >::iterator binVsRMSIt = binVsRMS.begin(); binVsRMSIt != binVsRMS.end(); ++binVsRMSIt) {
  //     ++nhits1;
  //     sumx1 += binVsRMSIt->first;
  //     sumy1 += binVsRMSIt->second;
  //     sumx21 += binVsRMSIt->first * binVsRMSIt->first;
  //     sumxy1 += binVsRMSIt->first * binVsRMSIt->second;
  //   }
  //   double RMSgradient1 = (nhits1 * sumxy1 - sumx1 * sumy1) / (nhits1 * sumx21 - sumx1 * sumx1);

  //   // NEW
  //   double sumx=0., sumy=0., sumx2=0., sumxy=0., sumweight = 0.;
  //   for (std::vector<std::pair<int,double> >::iterator binVsRMSIt = binVsRMS.begin(); binVsRMSIt != binVsRMS.end(); ++binVsRMSIt) {
  //     //double weight = showerSegments.at(binVsRMSIt->first).size();
  //     double weight = TMath::Power(segmentCharge.at(binVsRMSIt->first),2);
  //     sumweight += weight;
  //     sumx += weight * binVsRMSIt->first;
  //     sumy += weight * binVsRMSIt->second;
  //     sumx2 += weight * binVsRMSIt->first * binVsRMSIt->first;
  //     sumxy += weight * binVsRMSIt->first * binVsRMSIt->second;
  //   }
  //   double RMSgradient = (sumweight * sumxy - sumx * sumy) / (sumweight * sumx2 - sumx * sumx);

  //   // PCA
  //   TPrincipal* pca = new TPrincipal(2,"");
  //   double point[2];
  //   for (std::vector<std::pair<int,double> >::iterator binVsRMSIt = binVsRMS.begin(); binVsRMSIt != binVsRMS.end(); ++binVsRMSIt) {
  //     if (fDebug > 2)
  //    std::cout << "Number of hits in segment " << binVsRMSIt->first << " is " << showerSegments.at(binVsRMSIt->first).size() << std::endl;
  //     for (int nhit = 0; nhit < (int)showerSegments.at(binVsRMSIt->first).size(); ++nhit) {
  //    point[0] = binVsRMSIt->first;
  //    point[1] = binVsRMSIt->second;
  //    pca->AddRow(point);
  //     }
  //   }
  //   pca->MakePrincipals();
  //   TVector2 RMSgradientPCA = TVector2( (*pca->GetEigenVectors())[0][0], (*pca->GetEigenVectors())[1][0] ).Unit();

  //   if (fDebug > 1)
  //     std::cout << "RMS gradient without weighting: " << RMSgradient1 << ", and with weighting: " << RMSgradient << "; PCA with weighting: " << RMSgradientPCA.Y()/RMSgradientPCA.X() << std::endl;

  //   if (fMakeRMSGradientPlot) {
  //     TVector2 direction = TVector2(1,RMSgradient).Unit();
  //     TVector2 direction1 = TVector2(1,RMSgradient1).Unit();
  //     TCanvas* canv = new TCanvas();
  //     graph->Draw();
  //     TVector2 centre = TVector2(graph->GetMean(1), graph->GetMean(2));
  //     TLine line, line1, line2;
  //     line.SetLineColor(3);
  //     line1.SetLineColor(4);
  //     line2.SetLineColor(5);
  //     line.DrawLine(centre.X()-1000*direction.X(),centre.Y()-1000*direction.Y(),centre.X()+1000*direction.X(),centre.Y()+1000*direction.Y());
  //     line1.DrawLine(centre.X()-1000*direction1.X(),centre.Y()-1000*direction1.Y(),centre.X()+1000*direction1.X(),centre.Y()+1000*direction1.Y());
  //     line2.DrawLine(centre.X()-1000*RMSgradientPCA.X(),centre.Y()-1000*RMSgradientPCA.Y(),centre.X()+1000*RMSgradientPCA.X(),centre.Y()+1000*RMSgradientPCA.Y());
  //     canv->SaveAs("RMSGradient.png");
  //     delete canv;
  //   }
  //   delete graph; delete pca;

  //   //std::cout << "Number of hits per bin " << i << " gives gradient of " << RMSgradient << std::endl;

  //   double distanceFromStart = 0, distanceFromEnd = 0;
  //   if (RMSgradient > 0) {
  //     distanceFromStart = (trueStart - HitPosition(showerHits.front())).Mod();
  //     distanceFromEnd = (trueStart - HitPosition(showerHits.back())).Mod();
  //   }
  //   if (RMSgradient < 0) {
  //    distanceFromStart = (trueStart - HitPosition(showerHits.back())).Mod();
  //    distanceFromEnd = (trueStart - HitPosition(showerHits.front())).Mod();
  //   }
  //   bool correct = distanceFromStart < distanceFromEnd;
  //   hNumSegments->Fill(i,correct);

  // }
















  // Bin the hits into discreet chunks
  int nShowerSegments = fNumShowerSegments;
  double lengthOfShower = (HitPosition(showerHits.back()) - HitPosition(showerHits.front())).Mod();
  double lengthOfSegment = lengthOfShower / (double)nShowerSegments;
  std::map<int,std::vector<art::Ptr<recob::Hit> > > showerSegments;
  std::map<int,double> segmentCharge;
  for (std::vector<art::Ptr<recob::Hit> >::const_iterator showerHitsIt = showerHits.begin(); showerHitsIt != showerHits.end(); ++showerHitsIt) {
    showerSegments[(int)(HitPosition(*showerHitsIt)-HitPosition(showerHits.front())).Mod() / lengthOfSegment].push_back(*showerHitsIt);
    segmentCharge[(int)(HitPosition(*showerHitsIt)-HitPosition(showerHits.front())).Mod() / lengthOfSegment] += (*showerHitsIt)->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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::iterator showerSegmentIt = showerSegments.begin(); showerSegmentIt != showerSegments.end(); ++showerSegmentIt) {

    // Get the mean position of the hits in this bin
    TVector2 meanPosition(0,0);
    for (std::vector<art::Ptr<recob::Hit> >::iterator hitInSegmentIt = showerSegmentIt->second.begin(); hitInSegmentIt != showerSegmentIt->second.end(); ++hitInSegmentIt)
      meanPosition += HitPosition(*hitInSegmentIt);
    meanPosition /= (double)showerSegmentIt->second.size();

    // Get the RMS of this bin
    std::vector<double> distanceToAxisBin;
    for (std::vector<art::Ptr<recob::Hit> >::iterator hitInSegmentIt = showerSegmentIt->second.begin(); hitInSegmentIt != showerSegmentIt->second.end(); ++hitInSegmentIt) {
      TVector2 proj = (HitPosition(*hitInSegmentIt) - meanPosition).Proj(direction) + meanPosition;
      distanceToAxisBin.push_back((HitPosition(*hitInSegmentIt) - proj).Mod());
    }

    double RMSBin = TMath::RMS(distanceToAxisBin.begin(), distanceToAxisBin.end());
    binVsRMS.push_back(std::make_pair(showerSegmentIt->first, RMSBin));
    if (fMakeRMSGradientPlot)
      graph->SetPoint(graph->GetN(), showerSegmentIt->first, RMSBin);

  }

  // Get the gradient of the RMS-bin plot
  double sumx=0., sumy=0., sumx2=0., sumxy=0., sumweight = 0.;
  for (std::vector<std::pair<int,double> >::iterator binVsRMSIt = binVsRMS.begin(); binVsRMSIt != binVsRMS.end(); ++binVsRMSIt) {
    //double weight = showerSegments.at(binVsRMSIt->first).size();
    //double weight = 1;
    double weight = segmentCharge.at(binVsRMSIt->first);
    sumweight += weight;
    sumx += weight * binVsRMSIt->first;
    sumy += weight * binVsRMSIt->second;
    sumx2 += weight * binVsRMSIt->first * binVsRMSIt->first;
    sumxy += weight * binVsRMSIt->first * binVsRMSIt->second;
  }
  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
	)

<Tingjun to="" document>="">

Definition at line 2157 of file EMShowerAlg.cxx.

References n.

Referenced by FindInitialTrackHits().

                                                                                                                           {

  Double_t sumx=0.;
  Double_t sumx2=0.;
  Double_t sumy=0.;
  Double_t sumy2=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]; 
    sumy2 += y[i]*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 )

private

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

Definition at line 2127 of file EMShowerAlg.cxx.

References fDebug, fGeom, for(), HitCoordinates(), geo::GeometryCore::MaxPlanes(), lar::dump::vector(), 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 (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt)
    planeWireLength[showerHitsIt->first] = TMath::Abs(HitCoordinates(showerHitsIt->second.front()).X() - HitCoordinates(showerHitsIt->second.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 += TMath::Power(planeWireLength[plane],2);
    quad = TMath::Sqrt(quad);
    planeOtherWireLengths[planeWireLengthIt->first] = planeWireLength[planeWireLengthIt->first] / (double)quad;
  }

  if (fDebug > 1)
    for (std::map<int,double>::const_iterator planeOtherWireLengthIt = planeOtherWireLengths.begin(); planeOtherWireLengthIt != planeOtherWireLengths.end(); ++planeOtherWireLengthIt)
      std::cout << "Plane " << planeOtherWireLengthIt->first << " has " << planeWireLength[planeOtherWireLengthIt->first] << " wire width and therefore has relative width in wire of " << planeOtherWireLengthIt->second << std::endl;

  std::map<double,int> wireWidthMap;
  for (std::map<int,std::vector<art::Ptr<recob::Hit> > >::const_iterator showerHitsIt = showerHitsMap.begin(); showerHitsIt != showerHitsMap.end(); ++showerHitsIt) {
    double wireWidth = planeWireLength.at(showerHitsIt->first);
    wireWidthMap[wireWidth] = showerHitsIt->first;
  }

  return wireWidthMap.begin()->second;

}

Member Data Documentation

art::ServiceHandle<cheat::BackTrackerService> shower::EMShowerAlg::bt_serv

private

Definition at line 249 of file EMShowerAlg.h.

calo::CalorimetryAlg shower::EMShowerAlg::fCalorimetryAlg

private

Definition at line 240 of file EMShowerAlg.h.

Referenced by FinddEdx(), and MakeShower().

int shower::EMShowerAlg::fDebug

Definition at line 147 of file EMShowerAlg.h.

Referenced by AssociateClustersAndTracks(), CheckShowerHits(), CheckShowerPlanes(), ConstructTrack(), FindInitialTrack(), FindShowerStart(), GetPlanePermutations(), MakeInitialTrack(), MakeShower(), MakeSpacePoints(), OrderShowerHits(), shower::EMShower::reconfigure(), and WorstPlane().

double shower::EMShowerAlg::fdEdxTrackLength

private

Definition at line 226 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), FinddEdx(), and MakeShower().

std::string shower::EMShowerAlg::fDetector

private

Definition at line 243 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and GlobalWire().

detinfo::DetectorProperties const* shower::EMShowerAlg::fDetProp

private

Definition at line 235 of file EMShowerAlg.h.

Referenced by Construct3DPoint(), shower::HitPosition::ConvertWireTickToCm(), HitPosition(), OrderShowerHits(), and Project3DPointOntoPlane().

art::ServiceHandle<geo::Geometry> shower::EMShowerAlg::fGeom

private

Definition at line 234 of file EMShowerAlg.h.

Referenced by CheckIsolatedHits(), Construct3DPoint(), shower::HitPosition::ConvertWireTickToCm(), FindInitialTrackHits(), GlobalWire(), HitPosition(), MakeShower(), MakeSpacePoints(), OrderShowerHits(), Project3DPointOntoPlane(), RelativeWireWidth(), and WorstPlane().

bool shower::EMShowerAlg::fMakeGradientPlot

private

Definition at line 253 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and FindOrderOfHits().

bool shower::EMShowerAlg::fMakeRMSGradientPlot

private

Definition at line 253 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and ShowerHitRMSGradient().

double shower::EMShowerAlg::fMinTrackLength

private

Definition at line 225 of file EMShowerAlg.h.

Referenced by AssociateClustersAndTracks(), and EMShowerAlg().

std::vector<unsigned int> shower::EMShowerAlg::fNfithits

private

Definition at line 230 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and FindInitialTrackHits().

unsigned int shower::EMShowerAlg::fNfitpass

private

Definition at line 229 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and FindInitialTrackHits().

int shower::EMShowerAlg::fNumShowerSegments

private

Definition at line 254 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and ShowerHitRMSGradient().

pma::ProjectionMatchingAlg shower::EMShowerAlg::fProjectionMatchingAlg

private

Definition at line 241 of file EMShowerAlg.h.

Referenced by ConstructTrack(), and MakeShower().

shower::ShowerEnergyAlg shower::EMShowerAlg::fShowerEnergyAlg

private

Definition at line 239 of file EMShowerAlg.h.

Referenced by MakeShower().

double shower::EMShowerAlg::fSpacePointSize

private

Definition at line 227 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and MakeSpacePoints().

std::vector<double> shower::EMShowerAlg::fToler

private

Definition at line 231 of file EMShowerAlg.h.

Referenced by EMShowerAlg(), and FindInitialTrackHits().

TProfile* shower::EMShowerAlg::hNumHitsInSegment

private

Definition at line 251 of file EMShowerAlg.h.

Referenced by EMShowerAlg().

TProfile * shower::EMShowerAlg::hNumSegments

private

Definition at line 251 of file EMShowerAlg.h.

Referenced by EMShowerAlg().

TH1I* shower::EMShowerAlg::hTrueDirection

private

Definition at line 250 of file EMShowerAlg.h.

Referenced by EMShowerAlg().

art::ServiceHandle<art::TFileService> shower::EMShowerAlg::tfs

private

Definition at line 236 of file EMShowerAlg.h.

Referenced by EMShowerAlg().

---

The documentation for this class was generated from the following files:

• larreco/v06_64_02/source/larreco/RecoAlg/EMShowerAlg.h
• larreco/v06_64_02/source/larreco/RecoAlg/EMShowerAlg.cxx