GCNFeatureUtils.cxx
Go to the documentation of this file.
1 #include <vector>
2 #include <iostream>
3 #include <ctime>
4 
13 
17 
18 #include "TVector3.h"
19 
20 using std::string;
21 using std::vector;
22 using std::pair;
23 using std::map;
24 using std::set;
25 
26 using art::Ptr;
27 using art::Event;
28 using art::ServiceHandle;
29 
32 
33 using simb::MCParticle;
34 
35 using recob::Hit;
36 using recob::SpacePoint;
37 
38 namespace cvn
39 {
40 
42 
44 
46  const recob::Hit& hit) const {
47 
48  // Get the backtracker, then find the TrackIDE with the largest energy
49  // deposit and return it
51  auto ides = bt->HitToTrackIDEs(clockData, hit);
52  if (ides.empty()) return std::numeric_limits<int>::max(); // Return invalid number if no true tracks
53  int id = 0;
54  float energy = -1;
55  for (auto ide : ides) {
56  if (ide.energy > energy) {
57  energy = ide.energy;
58  id = ide.trackID;
59  }
60  }
61  return id;
62 
63  } // function GCNFeatureUtils::GetTrackIDFromHit
64 
65  pair<unsigned int, float> GCNFeatureUtils::GetClosestApproach(const SpacePoint sp,
66  const MCParticle p) const {
67 
68  // Spacepoint 3D position
69  geoalgo::Point_t spPos(sp.XYZ()[0], sp.XYZ()[1], sp.XYZ()[2]);
70 
71  // Loop over trajectory segments and find point of closest approach
72  unsigned int id = std::numeric_limits<unsigned int>::max();
74  simb::MCTrajectory traj = p.Trajectory();
75  for (size_t it = 1; it < traj.size(); ++it) {
76  geoalgo::Point_t p1(TVector3(traj.Position(it-1).Vect()));
77  geoalgo::Point_t p2(TVector3(traj.Position(it).Vect()));
78  geoalgo::LineSegment_t ls(p1, p2);
79  float distTmp = fGeoAlgo.SqDist(spPos, ls);
80  if (distTmp < dist) {
81  dist = distTmp;
82  id = it-1;
83  }
84  }
85  return std::make_pair(id, dist);
86 
87  } // function GCNFeatureUtils::GetClosestApproach
88 
89  // Get the number of neighbours within rangeCut cm of this space point
91  const float rangeCut, const std::string &spLabel) const{
92 
93  return GetAllNeighbours(evt,rangeCut,spLabel).at(sp.ID());
94 
95  }
96 
97  // Get the number of neighbours within rangeCut cm of this space point
99  const float rangeCut, const std::vector<art::Ptr<recob::SpacePoint>> &sps) const{
100 
101  return GetAllNeighbours(evt,rangeCut,sps).at(sp.ID());
102 
103  }
104 
105  // Get a map of the number of neighbours for each space point ID. Using this function is much less wasteful
106  // than repeated calls to the above function
107  std::map<int,unsigned int> GCNFeatureUtils::GetAllNeighbours(art::Event const &evt, const float rangeCut, const std::string &spLabel) const{
108 
109  // Get the space points from the event and make the map
110  vector<Ptr<SpacePoint>> spacePoints;
111  auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
112  if (spacePointHandle) {
113  art::fill_ptr_vector(spacePoints, spacePointHandle);
114  }
115  return GetAllNeighbours(evt, rangeCut, spacePoints);
116  }
117 
118  std::map<int,unsigned int> GCNFeatureUtils::GetAllNeighbours(art::Event const &evt, const float rangeCut, const std::vector<art::Ptr<recob::SpacePoint>> &sps) const{
119 
120  map<int,unsigned int> neighbourMap;
121 
122  for(const Ptr<SpacePoint> sp0 : sps){
123  // We want an entry even if it ends up being zero
124  neighbourMap[sp0->ID()] = 0;
125 
126  for(const Ptr<SpacePoint> sp1 : sps){
127 
128  if(sp0->ID() == sp1->ID()) continue;
129 
130  // For some reason we have to use arrays
131  const double *p0 = sp0->XYZ();
132  const double *p1 = sp1->XYZ();
133 
134  // Get the distance between the points
135  const float dx = p1[0] - p0[0];
136  const float dy = p1[1] - p0[1];
137  const float dz = p1[2] - p0[2];
138  const float dist = sqrt(dx*dx + dy*dy + dz*dz);
139 
140  if(dist < rangeCut){
141  ++neighbourMap[sp0->ID()];
142  }
143  }
144  }
145  return neighbourMap;
146  } // function GetAllNeighbours
147 
148  // Sometimes we might want to know the number of neighbours within various radii
149  std::vector<std::map<int,unsigned int>> GCNFeatureUtils::GetNeighboursForRadii(art::Event const &evt, const std::vector<float>& rangeCuts, const std::string &spLabel) const{
150  // Get the space points from the event and make the map
151  vector<Ptr<SpacePoint>> allSpacePoints;
152  auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
153  if (spacePointHandle) {
154  art::fill_ptr_vector(allSpacePoints, spacePointHandle);
155  }
156  return GetNeighboursForRadii(evt,rangeCuts,allSpacePoints);
157  }
158 
159  std::vector<std::map<int,unsigned int>> GCNFeatureUtils::GetNeighboursForRadii(art::Event const &evt,
160  const std::vector<float>& rangeCuts, const std::map<unsigned int,art::Ptr<recob::SpacePoint>> &sps) const{
161  std::vector<art::Ptr<recob::SpacePoint>> vec;
162  for(auto m : sps){
163  vec.push_back(m.second);
164  }
165  return GetNeighboursForRadii(evt, rangeCuts, vec);
166  }
167 
168  std::vector<std::map<int,unsigned int>> GCNFeatureUtils::GetNeighboursForRadii(art::Event const &evt,
169  const std::vector<float>& rangeCuts, const std::vector<art::Ptr<recob::SpacePoint>> &sps) const{
170  std::vector<std::map<int,unsigned int>> result;
171  // Initialise a map for each range cut value
172  for(unsigned int m = 0; m < rangeCuts.size(); ++m){
173  result.push_back(map<int,unsigned int>());
174  }
175  for(const Ptr<SpacePoint> sp0 : sps){
176  // We want an entry even if it ends up being zero
177  for(auto &m : result){
178  m[sp0->ID()] = 0;
179  }
180  for(const Ptr<SpacePoint> sp1 : sps){
181  if(sp0->ID() == sp1->ID()) continue;
182  // For some reason we have to use arrays
183  const double *p0 = sp0->XYZ();
184  const double *p1 = sp1->XYZ();
185  // Get the distance between the points
186  const float dx = p1[0] - p0[0];
187  const float dy = p1[1] - p0[1];
188  const float dz = p1[2] - p0[2];
189  const float dist = sqrt(dx*dx + dy*dy + dz*dz);
190  // Fill the maps if we satify the criteria
191  for(unsigned int r = 0; r < rangeCuts.size(); ++r){
192  if(dist < rangeCuts[r]){
193  result[r][sp0->ID()] = result[r][sp0->ID()] + 1;
194  }
195  }
196  }
197  }
198  return result;
199  }
200 
202  const std::string &spLabel) const{
203  std::vector<art::Ptr<recob::SpacePoint>> allSpacePoints;
204  auto spacePointHandle = evt.getHandle<std::vector<recob::SpacePoint>>(spLabel);
205  if (spacePointHandle) {
206  art::fill_ptr_vector(allSpacePoints, spacePointHandle);
207  }
208  return GetNearestNeighbours(evt,allSpacePoints);
209  }
210 
212  const std::vector<art::Ptr<recob::SpacePoint>> &sps) const{
213  std::map<int,int> closestID;
214  std::map<int,float> closestDistance;
215  // Now loop over all of the space points and simultaneously fill our maps
216  // Get the space points from the event and make the map
217  for(const Ptr<SpacePoint> sp0 : sps){
218  // We want an entry even if it ends up being zero
219  closestID[sp0->ID()] = 0;
220  closestDistance[sp0->ID()] = 99999;
221  for(const Ptr<SpacePoint> sp1 : sps){
222  if(sp0->ID() == sp1->ID()) continue;
223  // For some reason we have to use arrays
224  const double *p0 = sp0->XYZ();
225  const double *p1 = sp1->XYZ();
226  // Get the distance between the points
227  const float dx = p1[0] - p0[0];
228  const float dy = p1[1] - p0[1];
229  const float dz = p1[2] - p0[2];
230  const float dist = sqrt(dx*dx + dy*dy + dz*dz);
231  if(dist < closestDistance[sp0->ID()]){
232  closestDistance[sp0->ID()] = dist;
233  closestID[sp0->ID()] = sp1->ID();
234  }
235  }
236  }
237  return closestID;
238  }
239 
240  // Get the two nearest neighbours to use for calcuation of angles between them and the node in question
241  std::map<int,std::pair<int,int>> GCNFeatureUtils::GetTwoNearestNeighbours(art::Event const &evt,
242  const std::string &spLabel) const{
243  std::vector<art::Ptr<recob::SpacePoint>> allSpacePoints;
244  auto spacePointHandle = evt.getHandle<std::vector<recob::SpacePoint>>(spLabel);
245  if (spacePointHandle) {
246  art::fill_ptr_vector(allSpacePoints, spacePointHandle);
247  }
248  return GetTwoNearestNeighbours(evt,allSpacePoints);
249  }
250 
251  std::map<int,std::pair<int,int>> GCNFeatureUtils::GetTwoNearestNeighbours(art::Event const &evt,
252  const std::map<unsigned int, art::Ptr<recob::SpacePoint>> &sps) const{
253 
254  vector<Ptr<SpacePoint>> vec;
255  for(auto m : sps){
256  vec.push_back(m.second);
257  }
258  return GetTwoNearestNeighbours(evt,vec);
259  }
260 
261  std::map<int,std::pair<int,int>> GCNFeatureUtils::GetTwoNearestNeighbours(art::Event const &evt,
262  const std::vector<art::Ptr<recob::SpacePoint>> &sps) const{
263  std::map<int,int> closestID;
264  std::map<int,int> secondID;
265  // Now loop over all of the space points and simultaneously fill our maps
266  // Get the space points from the event and make the map
267  for(const Ptr<SpacePoint> sp0 : sps){
268  // We want an entry even if it ends up being zero
269  int thisSP = sp0->ID();
270  int closest = -1;
271  int second = -1;
272  float closestDist = 99999;
273  float secondDist = 99999;
274  for(const Ptr<SpacePoint> sp1 : sps){
275  if(thisSP == sp1->ID()) continue;
276  // For some reason we have to use arrays
277  const double *p0 = sp0->XYZ();
278  const double *p1 = sp1->XYZ();
279  // Get the distance between the points
280  const float dx = p1[0] - p0[0];
281  const float dy = p1[1] - p0[1];
282  const float dz = p1[2] - p0[2];
283  const float dist = sqrt(dx*dx + dy*dy + dz*dz);
284  if(dist < closestDist){
285  secondDist = closestDist;
286  closestDist = dist;
287  second = closest;
288  closest = sp1->ID();
289  }
290  else if(dist < secondDist){
291  secondDist = dist;
292  second = sp1->ID();
293  }
294  }
295  closestID.insert(std::make_pair(thisSP,closest));
296  secondID.insert(std::make_pair(thisSP,second));
297  }
298  map<int,pair<int,int>> finalMap;
299  for(unsigned int m = 0; m < closestID.size(); ++m){
300  finalMap[m] = std::make_pair(closestID[m],secondID[m]);
301  }
302  return finalMap;
303  }
304 
305  // Get the angle and the dot product between the vector from the base node to its neighbours
307  const SpacePoint &n1, const SpacePoint &n2, float &dotProduct, float &angle) const{
308  TVector3 basePos(baseNode.XYZ());
309  TVector3 neighbour1Pos(n1.XYZ());
310  TVector3 neighbour2Pos(n2.XYZ());
311  TVector3 baseToNeighbour1 = neighbour1Pos - basePos;
312  TVector3 baseToNeighbour2 = neighbour2Pos - basePos;
313  dotProduct = baseToNeighbour1.Dot(baseToNeighbour2);
314  angle = baseToNeighbour1.Angle(baseToNeighbour2);
315  return;
316  }
317 
318  // Use the association between space points and hits to return a charge
319  std::map<unsigned int, float> GCNFeatureUtils::GetSpacePointChargeMap(
320  std::vector<art::Ptr<recob::SpacePoint>> const& spacePoints,
321  std::vector<std::vector<art::Ptr<recob::Hit>>> const& sp2Hit) const {
322 
323  map<unsigned int, float> ret;
324 
325  for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
326  float charge = 0.0;
327  for (Ptr<Hit> hit : sp2Hit[spIdx]) {
328  charge += hit->Integral();
329  }
330  ret[spacePoints[spIdx]->ID()] = charge;
331  }
332 
333  return ret;
334 
335  } // function GetSpacePointChargeMap
336 
337  std::map<unsigned int, float> GCNFeatureUtils::GetSpacePointChargeMap(
338  art::Event const &evt, const std::string &spLabel) const {
339 
340  vector<Ptr<SpacePoint>> spacePoints;
341  auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
342  if (!spacePointHandle) {
344  << "Could not find spacepoints with module label "
345  << spLabel << "!";
346  }
347  art::fill_ptr_vector(spacePoints, spacePointHandle);
348  art::FindManyP<Hit> fmp(spacePointHandle, evt, spLabel);
349  vector<vector<Ptr<Hit>>> sp2Hit(spacePoints.size());
350  for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) {
351  sp2Hit[spIdx] = fmp.at(spIdx);
352  } // for spacepoint
353 
354  return GetSpacePointChargeMap(spacePoints, sp2Hit);
355 
356  } // function GetSpacePointChargeMap
357 
358  std::map<unsigned int, float> GCNFeatureUtils::GetSpacePointMeanHitRMSMap(
359  art::Event const &evt, const std::string &spLabel) const {
360 
361  map<unsigned int, float> ret;
362 
363  // Get the hits associated to the space points
364  auto allSP = evt.getValidHandle<vector<SpacePoint>>(spLabel);
365  const art::FindManyP<Hit> sp2Hit(allSP, evt, spLabel);
366 
367  for (size_t itSP = 0; itSP < allSP->size(); ++itSP) {
368  const SpacePoint& sp = allSP->at(itSP);
369  float chargeRMS = 0.0;
370  const vector<Ptr<Hit>> spHits = sp2Hit.at(itSP);
371  unsigned int nHits = spHits.size();
372  for (auto const hit : spHits) {
373  chargeRMS += hit->RMS();
374  }
375  ret[sp.ID()] = chargeRMS / static_cast<float>(nHits);
376  }
377 
378  return ret;
379 
380  } // function GetSpacePointMeanHitRMSMap
381 
382  std::map<unsigned int, int> GCNFeatureUtils::GetTrueG4ID(
383  detinfo::DetectorClocksData const& clockData,
384  std::vector<art::Ptr<recob::SpacePoint>> const& spacePoints,
385  std::vector<std::vector<art::Ptr<recob::Hit>>> const& sp2Hit) const {
386 
388  map<unsigned int, int> ret;
389 
390  for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
391  // Use the backtracker to find the G4 IDs associated with these hits
392  std::map<unsigned int, float> trueParticles;
393  for (art::Ptr<recob::Hit> hit : sp2Hit[spIdx]) {
394  auto ides = bt->HitToTrackIDEs(clockData, hit);
395  for (auto ide : ides) {
396  int id = ide.trackID;
397  if (trueParticles.count(id)) trueParticles[id] += ide.energy;
398  else trueParticles[id] = ide.energy;
399  }
400  }
401  ret[spacePoints[spIdx]->ID()] = std::max_element(trueParticles.begin(), trueParticles.end(), [](const pair<unsigned int, float> &lhs,
402  const pair<unsigned int, float> &rhs) { return lhs.second < rhs.second; })->first;
403  }
404  return ret;
405 
406  } // function GetTrueG4ID
407 
408  std::map<unsigned int, int> GCNFeatureUtils::GetTrueG4ID(
409  detinfo::DetectorClocksData const& clockData,
410  art::Event const &evt, const std::string &spLabel) const {
411 
412  vector<Ptr<SpacePoint>> spacePoints;
413  auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
414  if (!spacePointHandle) {
415 
417  << "Could not find spacepoints with module label "
418  << spLabel << "!";
419  }
420  art::fill_ptr_vector(spacePoints, spacePointHandle);
421  art::FindManyP<Hit> fmp(spacePointHandle, evt, spLabel);
422  vector<vector<Ptr<Hit>>> sp2Hit(spacePoints.size());
423  for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) {
424  sp2Hit[spIdx] = fmp.at(spIdx);
425  } // for spacepoint
426  return GetTrueG4ID(clockData, spacePoints, sp2Hit);
427 
428  } // function GetTrueG4ID
429 
430  std::map<unsigned int, int> GCNFeatureUtils::GetTrueG4IDFromHits(
431  detinfo::DetectorClocksData const& clockData,
432  std::vector<art::Ptr<recob::SpacePoint>> const& spacePoints,
433  std::vector<std::vector<art::Ptr<recob::Hit>>> const& sp2Hit) const {
434 
436  map<unsigned int, int> ret;
437 
438  for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
439  // Use the backtracker to find the G4 IDs associated with these hits
440  std::map<unsigned int, unsigned int> trueParticleHits;
441  for (art::Ptr<recob::Hit> hit : sp2Hit[spIdx]) {
442  auto ides = bt->HitToTrackIDEs(clockData, hit);
443  for (auto ide : ides) {
444  int id = ide.trackID;
445  if (trueParticleHits.count(id)) trueParticleHits[id] += 1;
446  else trueParticleHits[id] = 1;
447  }
448  }
449  ret[spacePoints[spIdx]->ID()] = std::max_element(trueParticleHits.begin(), trueParticleHits.end(), [](const pair<unsigned int, unsigned> &lhs,
450  const pair<unsigned int, unsigned int> &rhs) { return lhs.second < rhs.second; })->first;
451  }
452  return ret;
453 
454  } // function GetTrueG4IDFromHits
455 
456  std::map<unsigned int, int> GCNFeatureUtils::GetTrueG4IDFromHits(
457  detinfo::DetectorClocksData const& clockData,
458  art::Event const &evt, const std::string &spLabel) const {
459 
460  vector<Ptr<SpacePoint>> spacePoints;
461  auto spacePointHandle = evt.getHandle<vector<SpacePoint>>(spLabel);
462  if (!spacePointHandle) {
463 
465  << "Could not find spacepoints with module label "
466  << spLabel << "!";
467  }
468  art::fill_ptr_vector(spacePoints, spacePointHandle);
469  art::FindManyP<Hit> fmp(spacePointHandle, evt, spLabel);
470  vector<vector<Ptr<Hit>>> sp2Hit(spacePoints.size());
471  for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) {
472  sp2Hit[spIdx] = fmp.at(spIdx);
473  } // for spacepoint
474  return GetTrueG4IDFromHits(clockData, spacePoints, sp2Hit);
475 
476  } // function GetTrueG4IDFromHits
477 
478  std::map<unsigned int, int> GCNFeatureUtils::GetTruePDG(
479  detinfo::DetectorClocksData const& clockData,
480  art::Event const& evt, const std::string &spLabel, bool useAbsoluteTrackID, bool useHits) const {
481 
482  std::map<unsigned int, int> idMap;
483  if(useHits) idMap = GetTrueG4IDFromHits(clockData, evt,spLabel);
484  else idMap = GetTrueG4ID(clockData, evt,spLabel);
485 
486  map<unsigned int,int> pdgMap;
487 
489 
490  // Now we need to get the true pdg code for each GEANT track ID in the map
491  for(const pair<unsigned int, int> m : idMap){
492  int pdg = 0;
493  if(m.second == 0) std::cout << "Getting particle with ID " << m.second << " for space point " << m.first << std::endl;
494  else{
495  int trackID = m.second;
496  if(useAbsoluteTrackID || trackID >= 0){
497  trackID = abs(trackID);
498  pdg = pi->TrackIdToParticle_P(trackID)->PdgCode();
499  }
500  else pdg = 11; // Dummy value to flag EM activity
501  }
502  pdgMap.insert(std::make_pair(m.first,pdg));
503  }
504 
505  return pdgMap;
506  } // function GetTrueG4ID
507 
508  std::map<unsigned int, std::vector<float>> GCNFeatureUtils::Get2DFeatures(
509  std::vector<art::Ptr<recob::SpacePoint>> const& spacePoints,
510  std::vector<std::vector<art::Ptr<recob::Hit>>> const& sp2Hit) const {
511 
512  // For each spacepoint, we want a size-nine float vector containing wire,
513  // time and charge for each of the associated hits. These features assume
514  // each spacepoint has a maximum of one hit associated from each plane,
515  // and will throw an exception if it finds a spacepoint derived from
516  // more than one hit on the same plane. Be warned!
517 
518  map<unsigned int, vector<float>> ret;
519  for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
520  vector<float> feat(9, 0.);
521  // Loop over hits
522  for (Ptr<Hit> hit : sp2Hit[spIdx]) {
523  int offset = 3 * hit->WireID().Plane;
524  // Throw an error if this plane's info has already been filled
525  if (feat[offset+2] != 0) {
526  std::ostringstream err;
527  err << "2D features for plane " << hit->WireID().Plane << " have already been "
528  << "filled.";
529  throw std::runtime_error(err.str());
530  }
531  feat[offset] = hit->WireID().Wire;
532  feat[offset+1] = hit->PeakTime();
533  feat[offset+2] = hit->Integral();
534  } // for hit
535  ret[spacePoints[spIdx]->ID()] = feat;
536  } // for spacepoint
537  return ret;
538 
539  } // function GCNFeatureUtils::Get2DFeatures
540 
541  // Convert a pixel map into three 2D GCNGraph objects
542  vector<GCNGraph> GCNFeatureUtils::ExtractGraphsFromPixelMap(const PixelMap &pm, const float chargeThreshold) const{
543 
544  // Each pixel map has three vectors of length (nWires*nTDCs)
545  // Each value is the hit charge, and we will make GCNGraph for each view
546  vector<vector<float>> allViews;
547  allViews.push_back(pm.fPEX);
548  allViews.push_back(pm.fPEY);
549  allViews.push_back(pm.fPEZ);
550 
551  const unsigned int nWires = pm.fNWire;
552  const unsigned int nTDCs = pm.fNTdc;
553 
554  vector<GCNGraph> outputGraphs;
555 
556  for(unsigned int v = 0; v < allViews.size(); ++v){
557 
558  GCNGraph newGraph;
559 
560  for(unsigned int w = 0; w < nWires; ++w){
561 
562  for(unsigned int t = 0; t < nTDCs; ++t){
563 
564  const unsigned int index = w*nTDCs + t;
565  const float charge = allViews[v][index];
566 
567  // If the charge is very small then ignore this pixel
568  if(charge < chargeThreshold) continue;
569 
570  // Put the positons into a vector
571  vector<float> pos = {static_cast<float>(w),static_cast<float>(t)};
572  // and the charge
573  vector<float> features = {charge};
574 
575  GCNGraphNode newNode(pos,features);
576  newGraph.AddNode(newNode);
577  } // loop over TDCs
578  } // loop over wires
579  outputGraphs.push_back(newGraph);
580  } // loop over views
581 
582  return outputGraphs;
583  }
584 
585  std::map<unsigned int,unsigned int> GCNFeatureUtils::Get2DGraphNeighbourMap(const GCNGraph &g, const unsigned int npixel) const{
586 
587  map<unsigned int,unsigned int> neighbourMap;
588 
589  // Loop over the nodes and get the wire and tdc
590  for(unsigned int n1 = 0; n1 < g.GetNumberOfNodes(); ++n1){
591 
592  neighbourMap[n1] = 0;
593  const GCNGraphNode &node1 = g.GetNode(n1);
594  const unsigned int w1 = static_cast<unsigned int>(node1.GetPosition()[0]);
595  const unsigned int t1 = static_cast<unsigned int>(node1.GetPosition()[1]);
596 
597  // Loop over the nodes again to compare w,t coordinates
598  for(unsigned int n2 = 0; n2 < g.GetNumberOfNodes(); ++n2){
599 
600  if(n1 == n2) continue;
601 
602  const GCNGraphNode &node2 = g.GetNode(n2);
603  const unsigned int w2 = static_cast<unsigned int>(node2.GetPosition()[0]);
604  const unsigned int t2 = static_cast<unsigned int>(node2.GetPosition()[1]);
605 
606  // Check if the box around the node for neighbours
607  // In this example npixel = 2.
608  //
609  // |---|---|---|---|---|---|---|
610  // | | | | | | | |
611  // |---|---|---|---|---|---|---|
612  // | | x | x | x | x | x | |
613  // |---|---|---|---|---|---|---|
614  // | | x | x | x | x | x | |
615  // |---|---|---|---|---|---|---|
616  // | | x | x |n1 | x | x | | t
617  // |---|---|---|---|---|---|---|
618  // | | x | x | x | x | x | |
619  // |---|---|---|---|---|---|---|
620  // | | x | x | x | x | x | |
621  // |---|---|---|---|---|---|---|
622  // | | | | | | | |
623  // |---|---|---|---|---|---|---|
624  // w
625 
626  if(w2 <= w1 + npixel && w2 >= w1 - npixel){
627  if(t2 <= t1 + npixel && t2 >= t1 - npixel){
628  neighbourMap[n1]++;
629  }
630  }
631 
632  }
633  }
634 
635  return neighbourMap;
636 
637  }
638 
640  detinfo::DetectorClocksData const& clockData,
641  std::vector<art::Ptr<recob::SpacePoint>> const& spacePoints,
642  std::vector<std::vector<art::Ptr<recob::Hit>>> const& sp2Hit, float distCut,
643  std::vector<std::vector<float>>* dirTruth) const{
644 
645  // Fetch cheat services
648 
649  vector<float> ret(spacePoints.size(), -1);
650  if (dirTruth) {
651  dirTruth->clear();
652  dirTruth->resize(spacePoints.size());
653  }
654  vector<pair<size_t, float>> dist;
655  vector<int> trueIDs(spacePoints.size(), -1);
656  vector<int> closestTrajPoint(spacePoints.size(), -1);
657  // vector<float> dist(spacePoints.size(), -1);
658 
659  // Debug counters!
660  size_t nMismatched(0), nNoHit(0), nGood(0);//, nOutOfRange(0), nHitUsed(0), nGood(0);
661  // Loop over each spacepoint
662  for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
663  // Check index and ID agree
664  if (spIdx != (size_t)spacePoints[spIdx]->ID()) {
665  throw art::Exception(art::errors::LogicError) << "Spacepoint index "
666  << spIdx << " mismatched with spacepoint ID "
667  << spacePoints[spIdx]->ID();
668  }
669  bool done = false;
670  Ptr<SpacePoint> sp = spacePoints[spIdx];
671 
672  // Loop over each hit associated with this spacepoint. At this point we
673  // ask three questions to figure out if this should be considered a "true"
674  // spacepoint:
675  // 1) Do all the hits associated with this spacepoint come from a true
676  // particle (ie. are not noise)?
677  // 2) If yes to 1), do all the hits associated with this spacepoint
678  // derive the majority of their energy from the same true particle?
679  // 3) Does the reconstructed position of the spacepoint fall within
680  // some FHICL-configurable distance to the true particle's
681  // trajectory?
682  // 4) Is the spacepoint's hit already used by a different spacepoint?
683  // If yes to all of these, the spacepoint is true. Otherwise, it's false.
684 
685  int trueParticleID = std::numeric_limits<int>::max();
686  bool firstHit = true;
687  for (art::Ptr<recob::Hit> hit : sp2Hit[spIdx]) {
688  int trueID = GetTrackIDFromHit(clockData, *hit);
689  if (trueID == std::numeric_limits<int>::max()) {
690  ++nNoHit;
691  done = true;
692  break;
693  }
694  // Check whether this hit's true particle matches other hits from this spacepoint
695  if (firstHit) {
696  trueParticleID = trueID;
697  firstHit = false;
698  }
699  else if (trueParticleID != trueID) {
700  // If true IDs don't match, quit
701  ++nMismatched;
702  done = true;
703  break;
704  }
705  } // for hit
706 
707  if (done) continue; // if this isn't a valid spacepoint then stop here
708  if (trueParticleID == std::numeric_limits<int>::max()) {
709  ++nNoHit;
710  continue;
711  }
712 
713  // Get the true particle, and find the closest MC trajectory point to spacepoint
714  MCParticle p = pi->TrackIdToParticle(abs(trueParticleID));
715  pair<unsigned int, float> closest = GetClosestApproach(*sp, p);
716  dist.push_back(std::make_pair(spIdx, closest.second));
717  trueIDs[spIdx] = abs(trueParticleID);
718  closestTrajPoint[spIdx] = closest.first;
719  ret[spIdx] = closest.second;
720  ++nGood;
721  //if (closest.second < distCut) {
722  // ++nGood;
723  //ret[spIdx] = -1;
724  //}
725  //else {
726  // ++nOutOfRange;
727  //}
728 
729  } // for spIdx
730 
731  // Sort spacepoints in ascending order of distance to true trajectory
732  //std::sort(std::begin(dist), std::end(dist),
733  // [](const auto& a, const auto&b) { return a.second < b.second; });
734  //set<size_t> usedHits;
735  // Loop over all truth-matched spacepoints
736  //for (auto d : dist) {
737  // // Check whether any of the hits associated with this spacepoint have already been used
738  // for (Ptr<Hit> hit : sp2Hit[d.first]) {
739  // if (usedHits.count(hit.key())) {
740  // ret[d.first] = false; // Remove this spacepoint from the list of good ones
741  // ++nHitUsed;
742  // break;
743  // }
744  // }
745  // if (!ret[d.first]) continue; // If this spacepoint is bad, move on to the next one
746  // // None of the hits have been used, so we'll let this one through but
747  // // add its hits to the list so no other spacepoints can use them
748  // for (Ptr<Hit> hit : sp2Hit[d.first]) {
749  // usedHits.insert(hit.key());
750  // }
751  // ++nGood;
752 
753  //} // for spacepoint
754 
755  // Loop over spacepoints to set directionality
756  for (size_t spIdx = 0; spIdx < spacePoints.size(); ++spIdx) {
757  // If this was a true node & we're storing direction, store it!
758  if (ret[spIdx] && dirTruth) {
759  MCParticle p = pi->TrackIdToParticle(trueIDs[spIdx]);
760  TVector3 dir = p.Trajectory().Momentum(closestTrajPoint[spIdx]).Vect().Unit();
761  for (size_t it = 0; it < 3; ++it) {
762  dirTruth->at(spIdx).push_back(dir[it]);
763  }
764  }
765 
766  } // for spIdx
767 
768  std::cout << "There were " << nGood << " valid spacepoints, " << nNoHit//<< nHitUsed
769 // << " spacepoints with a hit already used, " << nNoHit
770  << " spacepoints with no associated hits or noise hits, and " << nMismatched
771  << " spacepoints produced from mismatched hits." << std::endl;// and " << nOutOfRange
772 // << " spacepoints reconstructed too far away from the true particle trajectory."
773 // << std::endl;
774 
775  return ret;
776 
777  } // function GCNFeatureUtils::GetNodeGroundTruth
778 
779  std::map<unsigned int, unsigned int> GCNFeatureUtils::GetParticleFlowMap(
780  const std::set<unsigned int>& particles) const {
781 
783  map<unsigned int, unsigned int> ret;
784  for (int p : particles) {
785  if (p == 0) continue; // No parent for the event primary
786  ret[p] = pi->TrackIdToParticle(abs(p)).Mother();
787  }
788  return ret;
789 
790  } // function GCNFeatureUtils::GetParticleFlowMap
791 
793  const cvn::GCNGraph* g) {
794 
796 
797  set<int> trackIDs;
798  for (unsigned int i = 0; i < g->GetNumberOfNodes(); ++i) {
799  trackIDs.insert(g->GetNode(i).GetGroundTruth()[0]);
800  }
801 
802  set<int> allIDs = trackIDs; // Copy original so we can safely modify it
803 
804  vector<ptruth> ret;
805 
806  // Add invisible particles to hierarchy
807  for (int id : trackIDs) {
808  const MCParticle* p = pi->TrackIdToParticle_P(abs(id));
809  while (p->Mother() != 0) {
810  allIDs.insert(abs(p->Mother()));
811  p = pi->TrackIdToParticle_P(abs(p->Mother()));
812  }
813  }
814 
815  for (int id : allIDs) {
816  const MCParticle* p = pi->TrackIdToParticle_P(abs(id));
817 
818  ret.push_back(std::make_tuple(abs(id), p->PdgCode(), p->Mother(),
819  p->P(), p->Vx(), p->Vy(), p->Vz(), p->EndX(), p->EndY(), p->EndZ(),
820  p->Process(), p->EndProcess()));
821  }
822 
823  return ret;
824 
825  } // function GCNFeatureUtils::GetParticleTree
826 
827 } // namespace cvn
Node for GCN.
int PdgCode() const
Definition: MCParticle.h:212
static constexpr double g
Definition: Units.h:144
static QCString result
double EndZ() const
Definition: MCParticle.h:228
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::vector< float > fPEY
Vector of Y PE measurements for pixels.
Definition: PixelMap.h:82
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
const simb::MCParticle * TrackIdToParticle_P(int id) const
static std::vector< ptruth > GetParticleTree(const cvn::GCNGraph *g)
int Mother() const
Definition: MCParticle.h:213
double SqDist(const Line_t &line, const Point_t &pt) const
Definition: GeoAlgo.h:87
unsigned int fNWire
Number of wires, length of pixel map.
Definition: PixelMap.h:78
Representation of a simple 3D line segment Defines a finite 3D straight line by having the start and ...
GCNGraph, basic input for the GCN.
Definition: GCNGraph.h:18
struct vector vector
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.
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
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.
const std::vector< float > GetGroundTruth() const
Get the node truth.
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.
string dir
std::string Process() const
Definition: MCParticle.h:215
PixelMap for CVN.
Utility class for truth labels.
double EndY() const
Definition: MCParticle.h:227
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...
T abs(T value)
std::map< unsigned int, unsigned int > GetParticleFlowMap(const std::set< unsigned int > &particles) const
Get hierarchy map from set of particles.
std::vector< float > fPEZ
Vector of Y PE measurements for pixels.
Definition: PixelMap.h:83
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>
simb::MCParticle TrackIdToParticle(int const id) const
bt
Definition: tracks.py:83
std::string EndProcess() const
Definition: MCParticle.h:216
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.
std::map< int, unsigned int > GetAllNeighbours(art::Event const &evt, const float rangeCut, const std::string &spLabel) const
double P(const int i=0) const
Definition: MCParticle.h:234
std::vector< float > fPEX
Vector of X PE measurements for pixels.
Definition: PixelMap.h:81
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
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.
p
Definition: test.py:223
static int max(int a, int b)
const TLorentzVector & Position(const size_type) const
The accessor methods described above.
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.
void err(const char *fmt,...)
Definition: message.cpp:226
Detector simulation of raw signals on wires.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
double Vx(const int i=0) const
Definition: MCParticle.h:221
Declaration of signal hit object.
const GCNGraphNode & GetNode(const unsigned int index) const
Access nodes.
Definition: GCNGraph.cxx:59
int GetTrackIDFromHit(detinfo::DetectorClocksData const &clockData, const recob::Hit &) const
Get primary true G4 ID for hit.
size_type size() const
Definition: MCTrajectory.h:166
Contains all timing reference information for the detector.
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
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.
PixelMap, basic input to CVN neural net.
Definition: PixelMap.h:22
GCNGraph for GCN.
double Vz(const int i=0) const
Definition: MCParticle.h:223
std::vector< cvn::GCNGraph > ExtractGraphsFromPixelMap(const cvn::PixelMap &pm, const float chargeThreshold) const
Convert a pixel map into three 2D GCNGraph objects.
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.
ID_t ID() const
Definition: SpacePoint.h:75
const std::vector< float > GetPosition() const
Get the node position, features or ground truth.
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
float pi
Definition: units.py:11
const unsigned int GetNumberOfNodes() const
Get the number of nodes.
Definition: GCNGraph.cxx:54
TCEvent evt
Definition: DataStructs.cxx:7
unsigned int fNTdc
Number of tdcs, width of pixel map.
Definition: PixelMap.h:79
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
double EndX() const
Definition: MCParticle.h:226
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
std::map< unsigned int, float > GetSpacePointMeanHitRMSMap(art::Event const &evt, const std::string &spLabel) const
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...
const TLorentzVector & Momentum(const size_type) const
Definition: fwd.h:31
double Vy(const int i=0) const
Definition: MCParticle.h:222
geoalgo::GeoAlgo fGeoAlgo
QTextStream & endl(QTextStream &s)
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.