GCNGraphMakerProtoDUNE_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file GCNGraphMakerProtoDUNE_module.cc
3 // \brief Producer module for creating CVN Graph objects
4 // \author Leigh H. Whitehead - leigh.howard.whitehead@cern.ch
5 ////////////////////////////////////////////////////////////////////////
6 
7 // C/C++ includes
8 #include <iostream>
9 #include <sstream>
10 
11 // Framework includes
16 #include "fhiclcpp/ParameterSet.h"
20 
21 // LArSoft includes
26 
28 
29 namespace cvn {
30 
32  public:
33  explicit GCNGraphMakerProtoDUNE(fhicl::ParameterSet const& pset);
35 
36  void produce(art::Event& evt);
37  void beginJob();
38  void endJob();
39 
40 
41 
42  private:
43  /// Module label for input space points
45 
46  /// Minimum number of space points to produce a graph
47  unsigned short fMinClusterHits;
48 
49  /// Radii for calculating number of neighbours for any number of cut values
50  std::vector<float> fNeighbourRadii;
51 
52  /// Use only the beam slice as determined by Pandora
55 
56  /// Slicing module (typically Pandora)
58 
59  /// PFParticle module (typically Pandora)
61 
62  // Consider EM activity as different to the parent
63  bool fUseEM;
64 
65  // Use hits not energy for truth matching
67  };
68 
69 
70 
71  //.......................................................................
73  fSpacePointLabel (pset.get<std::string> ("SpacePointLabel")),
74  fMinClusterHits (pset.get<unsigned short> ("MinClusterHits")),
75  fNeighbourRadii (pset.get<std::vector<float>>("NeighbourRadii")),
76  fUseBeamSliceOnly(pset.get<bool>("UseBeamSliceOnly")),
77  fUseAllSlices(pset.get<bool>("UseAllSlices")),
78  fSliceLabel (pset.get<std::string>("SliceModuleLabel")),
79  fParticleLabel (pset.get<std::string>("ParticleModuleLabel")),
80  fUseEM (pset.get<bool>("UseEM",true)),
81  fUseHitsForTruthMatching (pset.get<bool>("UseHitsForTruthMatching",true))
82  {
83 
84  produces< std::vector<cvn::GCNGraph> >();
85 
86  }
87 
88  //......................................................................
90  {
91  //======================================================================
92  // Clean up any memory allocated by your module
93  //======================================================================
94  }
95 
96  //......................................................................
98  { }
99 
100  //......................................................................
102  {
103  }
104 
105  //......................................................................
107  {
108 
109 
110  // Create the Graph vector and fill it if we have enough hits
111  std::unique_ptr<std::vector<cvn::GCNGraph>> graphs(new std::vector<cvn::GCNGraph>);
112  // std::cout << "GCNGraphMakerProtoDUNE: checking if we have enough space points (" << pointList.size() << " / " << fMinClusterHits << ")" << std::endl;
113  // Vector for all of the space points
114  // Get the space points from the event
115  std::vector<art::Ptr<recob::SpacePoint>> allSpacePoints;
116  auto spacePointHandle = evt.getHandle<std::vector<recob::SpacePoint>>(fSpacePointLabel);
117  if (spacePointHandle){
118  art::fill_ptr_vector(allSpacePoints, spacePointHandle);
119  }
120 
121  // Graph space points for each slice we want to consider
122  std::map<unsigned int,std::map<unsigned int,art::Ptr<recob::SpacePoint>>> allGraphSpacePoints;
123 
124  // Get the utility to help us calculate features
125  cvn::GCNFeatureUtils graphUtil;
126 
127  // We can calculate the number of neighbours for each space point with some radius
128  // Store the number of neighbours for each spacepoint ID, for each slice. If we
129  // only want the beam slice, or all of the slices together, this will have one
130  // element
131  std::map<unsigned int,std::vector<std::map<int,unsigned int>>> neighbourMap;
132 
133  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
135  // We need to get a vector of space points from the beam slice
136  cvn::CVNProtoDUNEUtils pfpUtil;
137 
138  const unsigned short beamSliceIndex = pfpUtil.GetBeamSlice(evt,fParticleLabel);
139  if(beamSliceIndex > 500 && fUseBeamSliceOnly){
140  std::cerr << "Requested beam slice only yet one didn't exist, returning no graphs" << std::endl;
141  evt.put(std::move(graphs));
142  return;
143  }
144 
145  // Get a map of slice index to all of the PFParticles within it
146  const std::map<unsigned int,std::vector<const recob::PFParticle*>> particleSliceMap = pfpUtil.GetAllPFParticleSliceMap(evt, fParticleLabel);
147  for(const std::pair<unsigned int,std::vector<const recob::PFParticle*>> m : particleSliceMap){
148  unsigned int sliceID = m.first;
149  std::map<unsigned int,art::Ptr<recob::SpacePoint>> sliceSpacePoints;
150  // If we want the beam slice then make sure we have it
151  if(fUseBeamSliceOnly && (sliceID != beamSliceIndex)) continue;
152 
153  for(const recob::PFParticle* p : m.second){
154  // Get the SpacePoints associated to the PFParticle
155  const std::vector<const recob::SpacePoint*> particlePoints = pfpUtil.GetPFParticleSpacePoints(*p, evt, fParticleLabel);
156  //std::cout << "Got beam slice particle with " << particlePoints.size() << " space points" << std::endl;
157  for(const recob::SpacePoint* s : particlePoints){
158  sliceSpacePoints.insert(std::make_pair(s->ID(),allSpacePoints.at(s->ID())));
159  }
160  }
161  allGraphSpacePoints[sliceID] = sliceSpacePoints;
162  neighbourMap.insert(std::make_pair(sliceID,graphUtil.GetNeighboursForRadii(evt,fNeighbourRadii,sliceSpacePoints)));
163  }
164  }
165  else{
166  std::vector<art::Ptr<recob::SpacePoint>> eventSpacePoints;
167  spacePointHandle = evt.getHandle<std::vector<recob::SpacePoint>>(fSpacePointLabel);
168  if (spacePointHandle){
169  art::fill_ptr_vector(eventSpacePoints, spacePointHandle);
170  }
171  //Convert this vector to a map
172  std::map<unsigned int,art::Ptr<recob::SpacePoint>> mapVec;
173  for(art::Ptr<recob::SpacePoint> p : eventSpacePoints){
174  mapVec.insert(std::make_pair(p->ID(),p));
175  }
176  neighbourMap.insert(std::make_pair(0,graphUtil.GetNeighboursForRadii(evt,fNeighbourRadii,eventSpacePoints)));
177  allGraphSpacePoints.insert(std::make_pair(0,mapVec));
178  }
179 
180  // Get some maps we'll need to fill our features
181  std::cout << "Found all neighbours for " << neighbourMap.size() << " slices, building graphs..." << std::endl;
182 
183  // Function is linear in number of points so just do it once
184  std::map<unsigned int,float> chargeMap = graphUtil.GetSpacePointChargeMap(evt,fSpacePointLabel);
185 
186  // Mean hit RMS for each spacepoint
187  std::map<unsigned int, float> hitRMSMap = graphUtil.GetSpacePointMeanHitRMSMap(evt,fSpacePointLabel);
188 
189  // The true particle PDG code is needed for training node classifiers
190  std::map<unsigned int,int> trueIDMap = graphUtil.GetTruePDG(clockData, evt, fSpacePointLabel, !fUseEM, fUseHitsForTruthMatching);
191 
192  // Now we want to produce a graph for each one of the slices
193  for(const std::pair<unsigned int,std::map<unsigned int,art::Ptr<recob::SpacePoint>>> &sps : allGraphSpacePoints){
194  if(sps.second.size() >= fMinClusterHits){
195 
196  cvn::GCNGraph newGraph;
197 
198  std::map<int,std::pair<int,int>> twoNearest = graphUtil.GetTwoNearestNeighbours(evt,sps.second);
199 
200  std::cout << "Constructing graph for slice " << sps.first << " with " << sps.second.size() << " nodes." << std::endl;
201  for(const std::pair<unsigned int,art::Ptr<recob::SpacePoint>> &sp : sps.second){
202 
203  // Get the position
204  std::vector<float> position;
205  // Why does this use an array... we want a vector in any case
206  const double *pos = sp.second->XYZ();
207  for(unsigned int p = 0; p < 3; ++p) position.push_back(pos[p]);
208 
209  // Calculate some features
210  std::vector<float> features;
211 
212  // The neighbour map gives us our first feature(s)
213  for(unsigned int m = 0; m < fNeighbourRadii.size(); ++m){
214  features.push_back(neighbourMap.at(sps.first)[m].at(sp.second->ID()));
215  }
216 
217  // How about charge?
218  features.push_back(chargeMap.at(sp.second->ID()));
219 
220  // Now the hit width
221 // features.push_back(hitRMSMap.at(sp.second->ID()));
222 
223  // Angle and dot product between node and its two nearest neighbours
224  float angle = -999.;
225  float dotProduct = -999.;
226  int n1ID = twoNearest[sp.second->ID()].first;
227  int n2ID = twoNearest[sp.second->ID()].second;
228  const recob::SpacePoint n1 = *(sps.second.at(n1ID).get());
229  const recob::SpacePoint n2 = *(sps.second.at(n2ID).get());
230  graphUtil.GetAngleAndDotProduct(*(sp.second.get()),n1,n2,dotProduct,angle);
231  features.push_back(dotProduct);
232  features.push_back(angle);
233 
234  // We set the "ground truth" as the particle PDG code in this case
235  std::vector<float> truePDG;
236  truePDG.push_back(static_cast<float>(trueIDMap.at(sp.second->ID())));
237  newGraph.AddNode(position,features,truePDG);
238 // if(abs(trueIDMap.at(sp.second->ID())) != 13) std::cout << "Adding node " << sp.second->ID() << " with neighbours " << n1ID << " and " << n2ID << " and PDG = " << truePDG[0] << std::endl;
239  }
240 
241  std::cout << "GCNGraphMakerProtoDUNE: produced GCNGraph object with " << newGraph.GetNumberOfNodes() << " nodes" << std::endl;
242 
243  // Add out graph to the vector
244  graphs->push_back(newGraph);
245  }
246  }
247 
248  // Write our graph to the event
249  evt.put(std::move(graphs));
250  }
251 
253 } // end namespace cvn
254 ////////////////////////////////////////////////////////////////////////
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
GCNGraph, basic input for the GCN.
Definition: GCNGraph.h:18
STL namespace.
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.
void AddNode(std::vector< float > position, std::vector< float > features)
Add a new node.
Definition: GCNGraph.cxx:37
Utility class for truth labels.
const char features[]
Definition: feature_tests.c:2
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...
const std::map< unsigned int, std::vector< const recob::PFParticle * > > GetAllPFParticleSliceMap(art::Event const &evt, const std::string particleLabel) const
Get a map of slice index to all of the PFParticles within it.
bool fUseBeamSliceOnly
Use only the beam slice as determined by Pandora.
unsigned short fMinClusterHits
Minimum number of space points to produce a graph.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
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.
def move(depos, offset)
Definition: depos.py:107
const std::vector< const recob::SpacePoint * > GetPFParticleSpacePoints(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Get the SpacePoints associated to the PFParticle.
p
Definition: test.py:223
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
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.
GCNGraphMakerProtoDUNE(fhicl::ParameterSet const &pset)
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
std::string fParticleLabel
PFParticle module (typically Pandora)
unsigned short GetBeamSlice(art::Event const &evt, const std::string particleLabel) const
Try to get the slice tagged as beam. Returns 9999 if no beam slice was found.
GCNGraph for GCN.
std::string fSliceLabel
Slicing module (typically Pandora)
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.
Utilities for calculating feature values for the GCN.
const unsigned int GetNumberOfNodes() const
Get the number of nodes.
Definition: GCNGraph.cxx:54
TCEvent evt
Definition: DataStructs.cxx:7
std::string fSpacePointLabel
Module label for input space points.
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::map< unsigned int, float > GetSpacePointMeanHitRMSMap(art::Event const &evt, const std::string &spLabel) const
int bool
Definition: qglobal.h:345
static QCString * s
Definition: config.cpp:1042
std::vector< float > fNeighbourRadii
Radii for calculating number of neighbours for any number of cut values.
QTextStream & endl(QTextStream &s)
Class containing some utility functions for all things CVN.