cvn::GCNFeatureUtils Class Reference

Class containing some utility functions for all things CVN. More...

#include <GCNFeatureUtils.h>

Public Member Functions
	GCNFeatureUtils ()
	~GCNFeatureUtils ()
int	GetTrackIDFromHit (detinfo::DetectorClocksData const &clockData, const recob::Hit &) const
	Get primary true G4 ID for hit. More...
std::pair< unsigned int, float >	GetClosestApproach (const recob::SpacePoint sp, const simb::MCParticle p) const
	Get closest trajectory point for true particle given reconstructed spacepoint. More...
unsigned int	GetSpacePointNeighbours (const recob::SpacePoint &sp, art::Event const &evt, const float rangeCut, const std::string &spLabel) const
	Get the number of neighbours within rangeCut cm of this space point. More...
unsigned int	GetSpacePointNeighbours (const recob::SpacePoint &sp, art::Event const &evt, const float rangeCut, const std::vector< art::Ptr< recob::SpacePoint >> &sps) const
std::map< int, unsigned int >	GetAllNeighbours (art::Event const &evt, const float rangeCut, const std::string &spLabel) const
std::map< int, unsigned int >	GetAllNeighbours (art::Event const &evt, const float rangeCut, const std::vector< art::Ptr< recob::SpacePoint >> &sps) const
std::vector< std::map< int, unsigned int > >	GetNeighboursForRadii (art::Event const &evt, const std::vector< float > &rangeCuts, const std::string &spLabel) const
	Gets the number of nearest neigbours for each space point for a vector of cut values. Much more efficient that using the above functions multiple times. More...
std::vector< std::map< int, unsigned int > >	GetNeighboursForRadii (art::Event const &evt, const std::vector< float > &rangeCuts, const std::vector< art::Ptr< recob::SpacePoint >> &sps) const
std::vector< std::map< int, unsigned int > >	GetNeighboursForRadii (art::Event const &evt, const std::vector< float > &rangeCuts, const std::map< unsigned int, art::Ptr< recob::SpacePoint >> &sps) const
std::map< int, int >	GetNearestNeighbours (art::Event const &evt, const std::string &spLabel) const
	Get the nearest neighbour map for the spacepoints. Returns a map of <nose_spacepoint_id,nearest_neighbour_spacepoint_id> More...
std::map< int, int >	GetNearestNeighbours (art::Event const &evt, const std::vector< art::Ptr< recob::SpacePoint >> &sps) const
std::map< int, std::pair< int, int > >	GetTwoNearestNeighbours (art::Event const &evt, const std::string &spLabel) const
	Get the two nearest neighbours to use for calcuation of angles between them and the node in question. More...
std::map< int, std::pair< int, int > >	GetTwoNearestNeighbours (art::Event const &evt, const std::vector< art::Ptr< recob::SpacePoint >> &sps) const
std::map< int, std::pair< int, int > >	GetTwoNearestNeighbours (art::Event const &evt, const std::map< unsigned int, art::Ptr< recob::SpacePoint >> &sps) const
void	GetAngleAndDotProduct (const recob::SpacePoint &baseNode, const recob::SpacePoint &n1, const recob::SpacePoint &n2, float &dotProduct, float &angle) const
	Get the angle and the dot product between the vector from the base node to its neighbours. More...
std::map< unsigned int, float >	GetSpacePointChargeMap (std::vector< art::Ptr< recob::SpacePoint >> const &spacePoints, std::vector< std::vector< art::Ptr< recob::Hit >>> const &sp2Hit) const
	Use the association between space points and hits to return a charge. More...
std::map< unsigned int, float >	GetSpacePointChargeMap (art::Event const &evt, const std::string &spLabel) const
std::map< unsigned int, float >	GetSpacePointMeanHitRMSMap (art::Event const &evt, const std::string &spLabel) const
std::map< unsigned int, int >	GetTrueG4ID (detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::SpacePoint >> const &spacePoints, std::vector< std::vector< art::Ptr< recob::Hit >>> const &sp2Hit) const
	Get the true G4 ID for each spacepoint using energy matching. More...
std::map< unsigned int, int >	GetTrueG4ID (detinfo::DetectorClocksData const &clockData, art::Event const &evt, const std::string &spLabel) const
std::map< unsigned int, int >	GetTrueG4IDFromHits (detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::SpacePoint >> const &spacePoints, std::vector< std::vector< art::Ptr< recob::Hit >>> const &sp2Hit) const
	Get the true G4 ID for each spacepoint using energy matching. More...
std::map< unsigned int, int >	GetTrueG4IDFromHits (detinfo::DetectorClocksData const &clockData, art::Event const &evt, const std::string &spLabel) const
std::map< unsigned int, int >	GetTruePDG (detinfo::DetectorClocksData const &clockData, art::Event const &evt, const std::string &spLabel, bool useAbsoluteTrackID, bool useHits) const
	Get the true pdg code for each spacepoint. More...
std::map< unsigned int, std::vector< float > >	Get2DFeatures (std::vector< art::Ptr< recob::SpacePoint >> const &spacePoints, std::vector< std::vector< art::Ptr< recob::Hit >>> const &sp2Hit) const
	Get 2D hit features for a given spacepoint. More...
std::vector< cvn::GCNGraph >	ExtractGraphsFromPixelMap (const cvn::PixelMap &pm, const float chargeThreshold) const
	Convert a pixel map into three 2D GCNGraph objects. More...
std::map< unsigned int, unsigned int >	Get2DGraphNeighbourMap (const cvn::GCNGraph &g, const unsigned int npixel) const
	Get the neighbours map <graph node, neighbours> for the three 2D graph in 2 box (npixel+1) around the pixel. More...
std::vector< float >	GetNodeGroundTruth (detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::SpacePoint >> const &spacePoints, std::vector< std::vector< art::Ptr< recob::Hit >>> const &spToHit, float distCut, std::vector< std::vector< float >> *dirTruth=nullptr) const
	Get ground truth for spacepoint deghosting graph network. More...
std::map< unsigned int, unsigned int >	GetParticleFlowMap (const std::set< unsigned int > &particles) const
	Get hierarchy map from set of particles. More...

Static Public Member Functions
static std::vector< ptruth >	GetParticleTree (const cvn::GCNGraph *g)

Private Attributes
geoalgo::GeoAlgo	fGeoAlgo

Detailed Description

Class containing some utility functions for all things CVN.

Definition at line 36 of file GCNFeatureUtils.h.

Constructor & Destructor Documentation

cvn::GCNFeatureUtils::GCNFeatureUtils

(

)

Definition at line 41 of file GCNFeatureUtils.cxx.

{}

cvn::GCNFeatureUtils::~GCNFeatureUtils

(

)

Definition at line 43 of file GCNFeatureUtils.cxx.

{}

Member Function Documentation

vector< GCNGraph > cvn::GCNFeatureUtils::ExtractGraphsFromPixelMap	(	const cvn::PixelMap &	pm,
		const float	chargeThreshold
	)		const

Convert a pixel map into three 2D GCNGraph objects.

Definition at line 542 of file GCNFeatureUtils.cxx.

                                                                                                                  {

    // Each pixel map has three vectors of length (nWires*nTDCs)
    // Each value is the hit charge, and we will make GCNGraph for each view
    vector<vector<float>> allViews;
    allViews.push_back(pm.fPEX);
    allViews.push_back(pm.fPEY);
    allViews.push_back(pm.fPEZ);

    const unsigned int nWires = pm.fNWire;
    const unsigned int nTDCs = pm.fNTdc;

    vector<GCNGraph> outputGraphs;

    for(unsigned int v = 0; v < allViews.size(); ++v){

      GCNGraph newGraph;

      for(unsigned int w = 0; w < nWires; ++w){

        for(unsigned int t = 0; t < nTDCs; ++t){

          const unsigned int index = w*nTDCs + t;
          const float charge = allViews[v][index];

          // If the charge is very small then ignore this pixel
          if(charge < chargeThreshold) continue;

          // Put the positons into a vector
          vector<float> pos = {static_cast<float>(w),static_cast<float>(t)};
          // and the charge
          vector<float> features = {charge};

          GCNGraphNode newNode(pos,features);
          newGraph.AddNode(newNode);
        } // loop over TDCs
      } // loop over wires
      outputGraphs.push_back(newGraph);
    } // loop over views

    return outputGraphs;
  }

std::map< unsigned int, std::vector< float > > cvn::GCNFeatureUtils::Get2DFeatures	(	std::vector< art::Ptr< recob::SpacePoint >> const &	spacePoints,
		std::vector< std::vector< art::Ptr< recob::Hit >>> const &	sp2Hit
	)		const

Get 2D hit features for a given spacepoint.

Definition at line 508 of file GCNFeatureUtils.cxx.

                                                                {

    // For each spacepoint, we want a size-nine float vector containing wire,
    // time and charge for each of the associated hits. These features assume
    // each spacepoint has a maximum of one hit associated from each plane,
    // and will throw an exception if it finds a spacepoint derived from
    // more than one hit on the same plane. Be warned!

    map<unsigned int, vector<float>> ret;
    for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
      vector<float> feat(9, 0.);
      // Loop over hits
      for (Ptr<Hit> hit : sp2Hit[spIdx]) {
        int offset = 3 * hit->WireID().Plane;
        // Throw an error if this plane's info has already been filled
        if (feat[offset+2] != 0) {
          std::ostringstream err;
          err << "2D features for plane " << hit->WireID().Plane << " have already been "
          << "filled.";
          throw std::runtime_error(err.str());
        }
        feat[offset] = hit->WireID().Wire;
        feat[offset+1] = hit->PeakTime();
        feat[offset+2] = hit->Integral();
      } // for hit
      ret[spacePoints[spIdx]->ID()] = feat;
    } // for spacepoint
    return ret;

  } // function GCNFeatureUtils::Get2DFeatures

std::map< unsigned int, unsigned int > cvn::GCNFeatureUtils::Get2DGraphNeighbourMap	(	const cvn::GCNGraph &	g,
		const unsigned int	npixel
	)		const

Get the neighbours map <graph node, neighbours> for the three 2D graph in 2 box (npixel+1) around the pixel.

Definition at line 585 of file GCNFeatureUtils.cxx.

                                                                                                                             {

    map<unsigned int,unsigned int> neighbourMap;

    // Loop over the nodes and get the wire and tdc
    for(unsigned int n1 = 0; n1 < g.GetNumberOfNodes(); ++n1){

      neighbourMap[n1] = 0;
      const GCNGraphNode &node1 = g.GetNode(n1);
      const unsigned int w1 = static_cast<unsigned int>(node1.GetPosition()[0]);
      const unsigned int t1 = static_cast<unsigned int>(node1.GetPosition()[1]);

      // Loop over the nodes again to compare w,t coordinates
      for(unsigned int n2 = 0; n2 < g.GetNumberOfNodes(); ++n2){

        if(n1 == n2) continue;

        const GCNGraphNode &node2 = g.GetNode(n2);
        const unsigned int w2 = static_cast<unsigned int>(node2.GetPosition()[0]);
        const unsigned int t2 = static_cast<unsigned int>(node2.GetPosition()[1]);

        // Check if the box around the node for neighbours
        // In this example npixel = 2.
        //
        // |---|---|---|---|---|---|---|
        // |   |   |   |   |   |   |   |
        // |---|---|---|---|---|---|---|
        // |   | x | x | x | x | x |   |
        // |---|---|---|---|---|---|---|
        // |   | x | x | x | x | x |   |
        // |---|---|---|---|---|---|---|
        // |   | x | x |n1 | x | x |   | t
        // |---|---|---|---|---|---|---|
        // |   | x | x | x | x | x |   |
        // |---|---|---|---|---|---|---|
        // |   | x | x | x | x | x |   |
        // |---|---|---|---|---|---|---|
        // |   |   |   |   |   |   |   |
        // |---|---|---|---|---|---|---|
        //               w

        if(w2 <= w1 + npixel && w2 >= w1 - npixel){
          if(t2 <= t1 + npixel && t2 >= t1 - npixel){
            neighbourMap[n1]++;
          }
        }

      }
    }

    return neighbourMap;

  }

std::map< int, unsigned int > cvn::GCNFeatureUtils::GetAllNeighbours	(	art::Event const &	evt,
		const float	rangeCut,
		const std::string &	spLabel
	)		const

Get a map of the number of neighbours for each space point ID. Using this function is much less wasteful than repeated calls to the above function

Definition at line 107 of file GCNFeatureUtils.cxx.

                                                                                                                                     {

    // Get the space points from the event and make the map
    vector<Ptr<SpacePoint>> spacePoints;
    auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
    if (spacePointHandle) {
      art::fill_ptr_vector(spacePoints, spacePointHandle);
    }
    return GetAllNeighbours(evt, rangeCut, spacePoints);
  }

std::map< int, unsigned int > cvn::GCNFeatureUtils::GetAllNeighbours	(	art::Event const &	evt,
		const float	rangeCut,
		const std::vector< art::Ptr< recob::SpacePoint >> &	sps
	)		const

Definition at line 118 of file GCNFeatureUtils.cxx.

                                                                                                                                                          {

    map<int,unsigned int> neighbourMap;

    for(const Ptr<SpacePoint> sp0 : sps){
      // We want an entry even if it ends up being zero
      neighbourMap[sp0->ID()] = 0;

      for(const Ptr<SpacePoint> sp1 : sps){

        if(sp0->ID() == sp1->ID()) continue;

        // For some reason we have to use arrays
        const double *p0 = sp0->XYZ();
        const double *p1 = sp1->XYZ();

        // Get the distance between the points
        const float dx = p1[0] - p0[0];
        const float dy = p1[1] - p0[1];
        const float dz = p1[2] - p0[2];
        const float dist = sqrt(dx*dx + dy*dy + dz*dz);

        if(dist < rangeCut){
          ++neighbourMap[sp0->ID()];
        }
      }
    }
    return neighbourMap;
  } // function GetAllNeighbours

void cvn::GCNFeatureUtils::GetAngleAndDotProduct	(	const recob::SpacePoint &	baseNode,
		const recob::SpacePoint &	n1,
		const recob::SpacePoint &	n2,
		float &	dotProduct,
		float &	angle
	)		const

Get the angle and the dot product between the vector from the base node to its neighbours.

Definition at line 306 of file GCNFeatureUtils.cxx.

                                                                                      {
    TVector3 basePos(baseNode.XYZ());
    TVector3 neighbour1Pos(n1.XYZ());
    TVector3 neighbour2Pos(n2.XYZ());
    TVector3 baseToNeighbour1 = neighbour1Pos - basePos;
    TVector3 baseToNeighbour2 = neighbour2Pos - basePos;
    dotProduct = baseToNeighbour1.Dot(baseToNeighbour2);
    angle = baseToNeighbour1.Angle(baseToNeighbour2);
    return;
  }

pair< unsigned int, float > cvn::GCNFeatureUtils::GetClosestApproach	(	const recob::SpacePoint	sp,
		const simb::MCParticle	p
	)		const

Get closest trajectory point for true particle given reconstructed spacepoint.

Definition at line 65 of file GCNFeatureUtils.cxx.

                              {

    // Spacepoint 3D position
    geoalgo::Point_t spPos(sp.XYZ()[0], sp.XYZ()[1], sp.XYZ()[2]);

    // Loop over trajectory segments and find point of closest approach
    unsigned int id = std::numeric_limits<unsigned int>::max();
    float dist = std::numeric_limits<float>::max();
    simb::MCTrajectory traj = p.Trajectory();
    for (size_t it = 1; it < traj.size(); ++it) {
      geoalgo::Point_t p1(TVector3(traj.Position(it-1).Vect()));
      geoalgo::Point_t p2(TVector3(traj.Position(it).Vect()));
      geoalgo::LineSegment_t ls(p1, p2);
      float distTmp = fGeoAlgo.SqDist(spPos, ls);
      if (distTmp < dist) {
        dist = distTmp;
        id = it-1;
      }
    }
    return std::make_pair(id, dist);

  } // function GCNFeatureUtils::GetClosestApproach

std::map< int, int > cvn::GCNFeatureUtils::GetNearestNeighbours	(	art::Event const &	evt,
		const std::string &	spLabel
	)		const

Get the nearest neighbour map for the spacepoints. Returns a map of <nose_spacepoint_id,nearest_neighbour_spacepoint_id>

Definition at line 201 of file GCNFeatureUtils.cxx.

                                   {
    std::vector<art::Ptr<recob::SpacePoint>> allSpacePoints;
    auto spacePointHandle = evt.getHandle<std::vector<recob::SpacePoint>>(spLabel);
    if (spacePointHandle) {
      art::fill_ptr_vector(allSpacePoints, spacePointHandle);
    }
    return GetNearestNeighbours(evt,allSpacePoints);
  }

std::map< int, int > cvn::GCNFeatureUtils::GetNearestNeighbours	(	art::Event const &	evt,
		const std::vector< art::Ptr< recob::SpacePoint >> &	sps
	)		const

Definition at line 211 of file GCNFeatureUtils.cxx.

                                                        {
    std::map<int,int> closestID;
    std::map<int,float> closestDistance;
    // Now loop over all of the space points and simultaneously fill our maps
    // Get the space points from the event and make the map
    for(const Ptr<SpacePoint> sp0 : sps){
      // We want an entry even if it ends up being zero
      closestID[sp0->ID()] = 0;
      closestDistance[sp0->ID()] = 99999;
      for(const Ptr<SpacePoint> sp1 : sps){
        if(sp0->ID() == sp1->ID()) continue;
        // For some reason we have to use arrays
        const double *p0 = sp0->XYZ();
        const double *p1 = sp1->XYZ();
        // Get the distance between the points
        const float dx = p1[0] - p0[0];
        const float dy = p1[1] - p0[1];
        const float dz = p1[2] - p0[2];
        const float dist = sqrt(dx*dx + dy*dy + dz*dz);
        if(dist < closestDistance[sp0->ID()]){
          closestDistance[sp0->ID()] = dist;
          closestID[sp0->ID()] = sp1->ID();
        }
      }
    }
    return closestID;
  }

std::vector< std::map< int, unsigned int > > cvn::GCNFeatureUtils::GetNeighboursForRadii	(	art::Event const &	evt,
		const std::vector< float > &	rangeCuts,
		const std::string &	spLabel
	)		const

Gets the number of nearest neigbours for each space point for a vector of cut values. Much more efficient that using the above functions multiple times.

Definition at line 149 of file GCNFeatureUtils.cxx.

                                                                                                                                                                  {
    // Get the space points from the event and make the map
    vector<Ptr<SpacePoint>> allSpacePoints;
    auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
    if (spacePointHandle) {
      art::fill_ptr_vector(allSpacePoints, spacePointHandle);
    }
    return GetNeighboursForRadii(evt,rangeCuts,allSpacePoints);
  }

std::vector< std::map< int, unsigned int > > cvn::GCNFeatureUtils::GetNeighboursForRadii	(	art::Event const &	evt,
		const std::vector< float > &	rangeCuts,
		const std::vector< art::Ptr< recob::SpacePoint >> &	sps
	)		const

Definition at line 168 of file GCNFeatureUtils.cxx.

                                                                                           {
    std::vector<std::map<int,unsigned int>> result;
    // Initialise a map for each range cut value
    for(unsigned int m = 0; m < rangeCuts.size(); ++m){
      result.push_back(map<int,unsigned int>());
    }
    for(const Ptr<SpacePoint> sp0 : sps){
      // We want an entry even if it ends up being zero
      for(auto &m : result){
        m[sp0->ID()] = 0;
      }
      for(const Ptr<SpacePoint> sp1 : sps){
        if(sp0->ID() == sp1->ID()) continue;
        // For some reason we have to use arrays
        const double *p0 = sp0->XYZ();
        const double *p1 = sp1->XYZ();
        // Get the distance between the points
        const float dx = p1[0] - p0[0];
        const float dy = p1[1] - p0[1];
        const float dz = p1[2] - p0[2];
        const float dist = sqrt(dx*dx + dy*dy + dz*dz);
        // Fill the maps if we satify the criteria
        for(unsigned int r = 0; r < rangeCuts.size(); ++r){
          if(dist < rangeCuts[r]){
            result[r][sp0->ID()] = result[r][sp0->ID()] + 1;
          }
        }
      }
    }
    return result;
  }

std::vector< std::map< int, unsigned int > > cvn::GCNFeatureUtils::GetNeighboursForRadii	(	art::Event const &	evt,
		const std::vector< float > &	rangeCuts,
		const std::map< unsigned int, art::Ptr< recob::SpacePoint >> &	sps
	)		const

Definition at line 159 of file GCNFeatureUtils.cxx.

                                                                                                     {
    std::vector<art::Ptr<recob::SpacePoint>> vec;
    for(auto m : sps){
      vec.push_back(m.second);
    }
    return GetNeighboursForRadii(evt, rangeCuts, vec);
  }

std::vector< float > cvn::GCNFeatureUtils::GetNodeGroundTruth	(	detinfo::DetectorClocksData const &	clockData,
		std::vector< art::Ptr< recob::SpacePoint >> const &	spacePoints,
		std::vector< std::vector< art::Ptr< recob::Hit >>> const &	spToHit,
		float	distCut,
		std::vector< std::vector< float >> *	dirTruth = nullptr
	)		const

Get ground truth for spacepoint deghosting graph network.

Definition at line 639 of file GCNFeatureUtils.cxx.

                                                {

    // Fetch cheat services
    ServiceHandle<BackTrackerService> bt;
    ServiceHandle<ParticleInventoryService> pi;

    vector<float> ret(spacePoints.size(), -1);
    if (dirTruth) {
      dirTruth->clear();
      dirTruth->resize(spacePoints.size());
    }
    vector<pair<size_t, float>> dist;
    vector<int> trueIDs(spacePoints.size(), -1);
    vector<int> closestTrajPoint(spacePoints.size(), -1);
    // vector<float> dist(spacePoints.size(), -1);

    // Debug counters!
    size_t nMismatched(0), nNoHit(0), nGood(0);//, nOutOfRange(0), nHitUsed(0), nGood(0);
    // Loop over each spacepoint
    for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
      // Check index and ID agree
      if (spIdx != (size_t)spacePoints[spIdx]->ID()) {
        throw art::Exception(art::errors::LogicError) << "Spacepoint index "
          << spIdx << " mismatched with spacepoint ID "
          << spacePoints[spIdx]->ID();
      }
      bool done = false;
      Ptr<SpacePoint> sp = spacePoints[spIdx];

      // Loop over each hit associated with this spacepoint. At this point we
      // ask three questions to figure out if this should be considered a "true"
      // spacepoint:
      //   1) Do all the hits associated with this spacepoint come from a true
      //      particle (ie. are not noise)?
      //   2) If yes to 1), do all the hits associated with this spacepoint
      //      derive the majority of their energy from the same true particle?
      //   3) Does the reconstructed position of the spacepoint fall within
      //      some FHICL-configurable distance to the true particle's
      //      trajectory?
      //   4) Is the spacepoint's hit already used by a different spacepoint?
      // If yes to all of these, the spacepoint is true. Otherwise, it's false.

      int trueParticleID = std::numeric_limits<int>::max();
      bool firstHit = true;
      for (art::Ptr<recob::Hit> hit : sp2Hit[spIdx]) {
        int trueID = GetTrackIDFromHit(clockData, *hit);
        if (trueID == std::numeric_limits<int>::max()) {
          ++nNoHit;
          done = true;
          break;
        }
        // Check whether this hit's true particle matches other hits from this spacepoint
        if (firstHit) {
          trueParticleID = trueID;
          firstHit = false;
        }
        else if (trueParticleID != trueID) {
          // If true IDs don't match, quit
          ++nMismatched;
          done = true;
          break;
        }
      } // for hit

      if (done) continue; // if this isn't a valid spacepoint then stop here
      if (trueParticleID == std::numeric_limits<int>::max()) {
        ++nNoHit;
        continue;
      }

      // Get the true particle, and find the closest MC trajectory point to spacepoint
      MCParticle p = pi->TrackIdToParticle(abs(trueParticleID));
      pair<unsigned int, float> closest = GetClosestApproach(*sp, p);
      dist.push_back(std::make_pair(spIdx, closest.second));
      trueIDs[spIdx] = abs(trueParticleID);
      closestTrajPoint[spIdx] = closest.first;
      ret[spIdx] = closest.second;
      ++nGood;
      //if (closest.second < distCut) {
        // ++nGood;
        //ret[spIdx] = -1;
      //}
      //else {
      //  ++nOutOfRange;
      //}

    } // for spIdx

    // Sort spacepoints in ascending order of distance to true trajectory
    //std::sort(std::begin(dist), std::end(dist),
    //  [](const auto& a, const auto&b) { return a.second < b.second; });
    //set<size_t> usedHits;
    // Loop over all truth-matched spacepoints
    //for (auto d : dist) {
    //  // Check whether any of the hits associated with this spacepoint have already been used
    //  for (Ptr<Hit> hit : sp2Hit[d.first]) {
    //    if (usedHits.count(hit.key())) {
    //      ret[d.first] = false; // Remove this spacepoint from the list of good ones
    //      ++nHitUsed;
    //      break;
    //    }
    //  }
    //  if (!ret[d.first]) continue; // If this spacepoint is bad, move on to the next one
    //  // None of the hits have been used, so we'll let this one through but
    //  // add its hits to the list so no other spacepoints can use them
    //  for (Ptr<Hit> hit : sp2Hit[d.first]) {
    //    usedHits.insert(hit.key());
    //  }
    //  ++nGood;

    //} // for spacepoint

    // Loop over spacepoints to set directionality
    for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
      // If this was a true node & we're storing direction, store it!
      if (ret[spIdx] && dirTruth) {
        MCParticle p = pi->TrackIdToParticle(trueIDs[spIdx]);
        TVector3 dir = p.Trajectory().Momentum(closestTrajPoint[spIdx]).Vect().Unit();
        for (size_t it = 0; it < 3; ++it) {
          dirTruth->at(spIdx).push_back(dir[it]);
        }
      }

    } // for spIdx

    std::cout << "There were " << nGood << " valid spacepoints, " << nNoHit//<< nHitUsed
//      << " spacepoints with a hit already used, " << nNoHit
      << " spacepoints with no associated hits or noise hits, and " << nMismatched
      << " spacepoints produced from mismatched hits." << std::endl;// and " << nOutOfRange
//      << " spacepoints reconstructed too far away from the true particle trajectory."
//      << std::endl;

    return ret;

  } // function GCNFeatureUtils::GetNodeGroundTruth

std::map< unsigned int, unsigned int > cvn::GCNFeatureUtils::GetParticleFlowMap

(

const std::set< unsigned int > &

particles

)

const

Get hierarchy map from set of particles.

Definition at line 779 of file GCNFeatureUtils.cxx.

                                                 {

    ServiceHandle<ParticleInventoryService> pi;
    map<unsigned int, unsigned int> ret;
    for (int p : particles) {
      if (p == 0) continue; // No parent for the event primary
      ret[p] = pi->TrackIdToParticle(abs(p)).Mother();
    }
    return ret;

  } // function GCNFeatureUtils::GetParticleFlowMap

vector< ptruth > cvn::GCNFeatureUtils::GetParticleTree ( const cvn::GCNGraph *  g )

static

Definition at line 792 of file GCNFeatureUtils.cxx.

                          {

    ServiceHandle<ParticleInventoryService> pi;

    set<int> trackIDs;
    for (unsigned int i = 0; i < g->GetNumberOfNodes(); ++i) {
      trackIDs.insert(g->GetNode(i).GetGroundTruth()[0]);
    }

    set<int> allIDs = trackIDs; // Copy original so we can safely modify it

    vector<ptruth> ret;

    // Add invisible particles to hierarchy
    for (int id : trackIDs) {
      const MCParticle* p = pi->TrackIdToParticle_P(abs(id));
      while (p->Mother() != 0) {
        allIDs.insert(abs(p->Mother()));
        p = pi->TrackIdToParticle_P(abs(p->Mother()));
      }
    }

    for (int id : allIDs) {
      const MCParticle* p = pi->TrackIdToParticle_P(abs(id));

      ret.push_back(std::make_tuple(abs(id), p->PdgCode(), p->Mother(),
        p->P(), p->Vx(), p->Vy(), p->Vz(), p->EndX(), p->EndY(), p->EndZ(),
        p->Process(), p->EndProcess()));
    }

    return ret;

  } // function GCNFeatureUtils::GetParticleTree

std::map< unsigned int, float > cvn::GCNFeatureUtils::GetSpacePointChargeMap	(	std::vector< art::Ptr< recob::SpacePoint >> const &	spacePoints,
		std::vector< std::vector< art::Ptr< recob::Hit >>> const &	sp2Hit
	)		const

Use the association between space points and hits to return a charge.

Definition at line 319 of file GCNFeatureUtils.cxx.

                                                                {

    map<unsigned int, float> ret;

    for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
      float charge = 0.0;
      for (Ptr<Hit> hit : sp2Hit[spIdx]) {
        charge += hit->Integral();
      }
      ret[spacePoints[spIdx]->ID()] = charge;
    }

    return ret;

  } // function GetSpacePointChargeMap

std::map< unsigned int, float > cvn::GCNFeatureUtils::GetSpacePointChargeMap	(	art::Event const &	evt,
		const std::string &	spLabel
	)		const

Definition at line 337 of file GCNFeatureUtils.cxx.

                                                         {

    vector<Ptr<SpacePoint>> spacePoints;
    auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
    if (!spacePointHandle) {
      throw art::Exception(art::errors::LogicError)
        << "Could not find spacepoints with module label "
        << spLabel << "!";
    }
    art::fill_ptr_vector(spacePoints, spacePointHandle);
    art::FindManyP<Hit> fmp(spacePointHandle, evt, spLabel);
    vector<vector<Ptr<Hit>>> sp2Hit(spacePoints.size());
    for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) {
      sp2Hit[spIdx] = fmp.at(spIdx);
    } // for spacepoint

    return GetSpacePointChargeMap(spacePoints, sp2Hit);

  } // function GetSpacePointChargeMap

std::map< unsigned int, float > cvn::GCNFeatureUtils::GetSpacePointMeanHitRMSMap	(	art::Event const &	evt,
		const std::string &	spLabel
	)		const

Definition at line 358 of file GCNFeatureUtils.cxx.

                                                         {

    map<unsigned int, float> ret;

    // Get the hits associated to the space points
    auto allSP = evt.getValidHandle<vector<SpacePoint>>(spLabel);
    const art::FindManyP<Hit> sp2Hit(allSP, evt, spLabel);

    for (size_t itSP = 0; itSP < allSP->size(); ++itSP) {
      const SpacePoint& sp = allSP->at(itSP);
      float chargeRMS = 0.0;
      const vector<Ptr<Hit>> spHits = sp2Hit.at(itSP);
      unsigned int nHits = spHits.size();
      for (auto const hit : spHits) {
        chargeRMS += hit->RMS();
      }
      ret[sp.ID()] = chargeRMS / static_cast<float>(nHits);
    }

    return ret;

  } // function GetSpacePointMeanHitRMSMap

unsigned int cvn::GCNFeatureUtils::GetSpacePointNeighbours	(	const recob::SpacePoint &	sp,
		art::Event const &	evt,
		const float	rangeCut,
		const std::string &	spLabel
	)		const

Get the number of neighbours within rangeCut cm of this space point.

Definition at line 90 of file GCNFeatureUtils.cxx.

                                                                                                                        {

    return GetAllNeighbours(evt,rangeCut,spLabel).at(sp.ID());

  }

unsigned int cvn::GCNFeatureUtils::GetSpacePointNeighbours	(	const recob::SpacePoint &	sp,
		art::Event const &	evt,
		const float	rangeCut,
		const std::vector< art::Ptr< recob::SpacePoint >> &	sps
	)		const

Definition at line 98 of file GCNFeatureUtils.cxx.

                                                                                                                                             {

    return GetAllNeighbours(evt,rangeCut,sps).at(sp.ID());

  }

int cvn::GCNFeatureUtils::GetTrackIDFromHit	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Hit &	hit
	)		const

Get primary true G4 ID for hit.

Definition at line 45 of file GCNFeatureUtils.cxx.

                                                                    {

    // Get the backtracker, then find the TrackIDE with the largest energy
    // deposit and return it
    art::ServiceHandle<cheat::BackTrackerService> bt;
    auto ides = bt->HitToTrackIDEs(clockData, hit);
    if (ides.empty()) return std::numeric_limits<int>::max(); // Return invalid number if no true tracks
    int id = 0;
    float energy = -1;
    for (auto ide : ides) {
      if (ide.energy > energy) {
        energy = ide.energy;
        id = ide.trackID;
      }
    }
    return id;

  } // function GCNFeatureUtils::GetTrackIDFromHit

std::map< unsigned int, int > cvn::GCNFeatureUtils::GetTrueG4ID	(	detinfo::DetectorClocksData const &	clockData,
		std::vector< art::Ptr< recob::SpacePoint >> const &	spacePoints,
		std::vector< std::vector< art::Ptr< recob::Hit >>> const &	sp2Hit
	)		const

Get the true G4 ID for each spacepoint using energy matching.

Definition at line 382 of file GCNFeatureUtils.cxx.

                                                                {

    ServiceHandle<BackTrackerService> bt;
    map<unsigned int, int> ret;

    for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
      // Use the backtracker to find the G4 IDs associated with these hits
      std::map<unsigned int, float> trueParticles;
      for (art::Ptr<recob::Hit> hit : sp2Hit[spIdx]) {
        auto ides = bt->HitToTrackIDEs(clockData, hit);
        for (auto ide : ides) {
          int id = ide.trackID;
          if (trueParticles.count(id)) trueParticles[id] += ide.energy;
          else trueParticles[id] = ide.energy;
        }
      }
      ret[spacePoints[spIdx]->ID()] = std::max_element(trueParticles.begin(), trueParticles.end(), [](const pair<unsigned int, float> &lhs,
        const pair<unsigned int, float> &rhs) { return lhs.second < rhs.second; })->first;
    }
    return ret;

  } // function GetTrueG4ID

std::map< unsigned int, int > cvn::GCNFeatureUtils::GetTrueG4ID	(	detinfo::DetectorClocksData const &	clockData,
		art::Event const &	evt,
		const std::string &	spLabel
	)		const

Definition at line 408 of file GCNFeatureUtils.cxx.

                                                         {

    vector<Ptr<SpacePoint>> spacePoints;
    auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
    if (!spacePointHandle) {

      throw art::Exception(art::errors::LogicError)
        << "Could not find spacepoints with module label "
        << spLabel << "!";
    }
    art::fill_ptr_vector(spacePoints, spacePointHandle);
    art::FindManyP<Hit> fmp(spacePointHandle, evt, spLabel);
    vector<vector<Ptr<Hit>>> sp2Hit(spacePoints.size());
    for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) {
      sp2Hit[spIdx] = fmp.at(spIdx);
    } // for spacepoint
    return GetTrueG4ID(clockData, spacePoints, sp2Hit);

  } // function GetTrueG4ID

std::map< unsigned int, int > cvn::GCNFeatureUtils::GetTrueG4IDFromHits	(	detinfo::DetectorClocksData const &	clockData,
		std::vector< art::Ptr< recob::SpacePoint >> const &	spacePoints,
		std::vector< std::vector< art::Ptr< recob::Hit >>> const &	sp2Hit
	)		const

Get the true G4 ID for each spacepoint using energy matching.

Definition at line 430 of file GCNFeatureUtils.cxx.

                                                                {

    ServiceHandle<BackTrackerService> bt;
    map<unsigned int, int> ret;

    for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
      // Use the backtracker to find the G4 IDs associated with these hits
      std::map<unsigned int, unsigned int> trueParticleHits;
      for (art::Ptr<recob::Hit> hit : sp2Hit[spIdx]) {
        auto ides = bt->HitToTrackIDEs(clockData, hit);
        for (auto ide : ides) {
          int id = ide.trackID;
          if (trueParticleHits.count(id)) trueParticleHits[id] += 1;
          else trueParticleHits[id] = 1;
        }
      }
      ret[spacePoints[spIdx]->ID()] = std::max_element(trueParticleHits.begin(), trueParticleHits.end(), [](const pair<unsigned int, unsigned> &lhs,
        const pair<unsigned int, unsigned int> &rhs) { return lhs.second < rhs.second; })->first;
    }
    return ret;

  } // function GetTrueG4IDFromHits

std::map< unsigned int, int > cvn::GCNFeatureUtils::GetTrueG4IDFromHits	(	detinfo::DetectorClocksData const &	clockData,
		art::Event const &	evt,
		const std::string &	spLabel
	)		const

Definition at line 456 of file GCNFeatureUtils.cxx.

                                                         {

    vector<Ptr<SpacePoint>> spacePoints;
    auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
    if (!spacePointHandle) {

      throw art::Exception(art::errors::LogicError)
        << "Could not find spacepoints with module label "
        << spLabel << "!";
    }
    art::fill_ptr_vector(spacePoints, spacePointHandle);
    art::FindManyP<Hit> fmp(spacePointHandle, evt, spLabel);
    vector<vector<Ptr<Hit>>> sp2Hit(spacePoints.size());
    for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) {
      sp2Hit[spIdx] = fmp.at(spIdx);
    } // for spacepoint
    return GetTrueG4IDFromHits(clockData, spacePoints, sp2Hit);

  } // function GetTrueG4IDFromHits

std::map< unsigned int, int > cvn::GCNFeatureUtils::GetTruePDG	(	detinfo::DetectorClocksData const &	clockData,
		art::Event const &	evt,
		const std::string &	spLabel,
		bool	useAbsoluteTrackID,
		bool	useHits
	)		const

Get the true pdg code for each spacepoint.

Definition at line 478 of file GCNFeatureUtils.cxx.

                                                                                                {

    std::map<unsigned int, int> idMap;
    if(useHits) idMap  = GetTrueG4IDFromHits(clockData, evt,spLabel);
    else idMap = GetTrueG4ID(clockData, evt,spLabel);

    map<unsigned int,int> pdgMap;

    ServiceHandle<ParticleInventoryService> pi;

    // Now we need to get the true pdg code for each GEANT track ID in the map
    for(const pair<unsigned int, int> m : idMap){
      int pdg = 0;
      if(m.second == 0) std::cout << "Getting particle with ID " << m.second << " for space point " << m.first << std::endl;
      else{
        int trackID = m.second;
        if(useAbsoluteTrackID || trackID >= 0){
          trackID = abs(trackID);
          pdg = pi->TrackIdToParticle_P(trackID)->PdgCode();
        }
        else pdg = 11; // Dummy value to flag EM activity
      }
      pdgMap.insert(std::make_pair(m.first,pdg));
    }

    return pdgMap;
  } // function GetTrueG4ID

std::map< int, std::pair< int, int > > cvn::GCNFeatureUtils::GetTwoNearestNeighbours	(	art::Event const &	evt,
		const std::string &	spLabel
	)		const

Get the two nearest neighbours to use for calcuation of angles between them and the node in question.

Definition at line 241 of file GCNFeatureUtils.cxx.

                                   {
    std::vector<art::Ptr<recob::SpacePoint>> allSpacePoints;
    auto spacePointHandle = evt.getHandle<std::vector<recob::SpacePoint>>(spLabel);
    if (spacePointHandle) {
      art::fill_ptr_vector(allSpacePoints, spacePointHandle);
    }
    return GetTwoNearestNeighbours(evt,allSpacePoints);
  }

std::map< int, std::pair< int, int > > cvn::GCNFeatureUtils::GetTwoNearestNeighbours	(	art::Event const &	evt,
		const std::vector< art::Ptr< recob::SpacePoint >> &	sps
	)		const

Definition at line 261 of file GCNFeatureUtils.cxx.

                                                        {
    std::map<int,int> closestID;
    std::map<int,int> secondID;
    // Now loop over all of the space points and simultaneously fill our maps
    // Get the space points from the event and make the map
    for(const Ptr<SpacePoint> sp0 : sps){
      // We want an entry even if it ends up being zero
      int thisSP = sp0->ID();
      int closest = -1;
      int second = -1;
      float closestDist = 99999;
      float secondDist = 99999;
      for(const Ptr<SpacePoint> sp1 : sps){
        if(thisSP == sp1->ID()) continue;
        // For some reason we have to use arrays
        const double *p0 = sp0->XYZ();
        const double *p1 = sp1->XYZ();
        // Get the distance between the points
        const float dx = p1[0] - p0[0];
        const float dy = p1[1] - p0[1];
        const float dz = p1[2] - p0[2];
        const float dist = sqrt(dx*dx + dy*dy + dz*dz);
        if(dist < closestDist){
          secondDist = closestDist;
          closestDist = dist;
          second = closest;
          closest = sp1->ID();
        }
        else if(dist < secondDist){
          secondDist = dist;
          second = sp1->ID();
        }
      }
      closestID.insert(std::make_pair(thisSP,closest));
      secondID.insert(std::make_pair(thisSP,second));
    }
    map<int,pair<int,int>> finalMap;
    for(unsigned int m = 0; m < closestID.size(); ++m){
      finalMap[m] = std::make_pair(closestID[m],secondID[m]);
    }
    return finalMap;
  }

std::map< int, std::pair< int, int > > cvn::GCNFeatureUtils::GetTwoNearestNeighbours	(	art::Event const &	evt,
		const std::map< unsigned int, art::Ptr< recob::SpacePoint >> &	sps
	)		const

Definition at line 251 of file GCNFeatureUtils.cxx.

                                                                   {

    vector<Ptr<SpacePoint>> vec;
    for(auto m : sps){
      vec.push_back(m.second);
    }
    return GetTwoNearestNeighbours(evt,vec);
  }

Member Data Documentation

geoalgo::GeoAlgo cvn::GCNFeatureUtils::fGeoAlgo

private

Definition at line 122 of file GCNFeatureUtils.h.

---

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

• dunereco/dunereco/CVN/func/GCNFeatureUtils.h
• dunereco/dunereco/CVN/func/GCNFeatureUtils.cxx