LArSoft  v08_20_00
Liquid Argon Software toolkit - http://larsoft.org/
CosmicPCAxisTagger_module.cc
Go to the documentation of this file.
1 // Class: CosmicPCAxisTagger
3 // Module Type: producer
4 // File: CosmicPCAxisTagger_module.cc
5 // This module checks timing and TPC volume boundaries as a
6 // way to tag potential cosmic rays
7 // This particular module uses PFParticles as input and handles
8 // special cases associated with them. Instead of tracks, it uses
9 // PCAxis objects for getting the start/end points of the
10 // candidate CR's
11 // This module started life as CosmicTrackTagger_module, written
12 // by Sarah Lockwitz, and small alterations made to handle the
13 // PFParticle input
14 //
15 // Generated at Wed Sep 17 19:17:00 2014 by Tracy Usher by cloning CosmicTrackTagger
16 // from art v1_02_02.
17 // artmod -e beginJob -e reconfigure -e endJob producer trkf::CosmicPCAxisTagger
19 
24 #include "art_root_io/TFileService.h"
26 
27 #include <iostream>
28 #include <cmath>
29 #include <algorithm>
30 #include <numeric>
31 #include <iterator>
32 
35 
42 
48 
49 #include "TVector3.h"
50 
51 namespace cosmic
52 {
53  class CosmicPCAxisTagger;
54  class SpacePoint;
55  class Track;
56 }
57 
58 
60 {
61 public:
62  explicit CosmicPCAxisTagger(fhicl::ParameterSet const & p);
63 
64  void produce(art::Event & e) override;
65 
66 private:
67  typedef std::vector<reco::ClusterHit2D> Hit2DVector;
68 
69  void RecobToClusterHits(const art::PtrVector<recob::Hit>& recobHitVec,
70  Hit2DVector& clusterHitVec,
71  reco::Hit2DListPtr& hit2DListPtr) const;
72 
74  std::string fPCAxisModuleLabel;
75 
78 
82 };
83 
84 
86  EDProducer{p},
87  fDetector(lar::providerFrom<detinfo::DetectorPropertiesService>()),
88  fPcaAlg(p.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
89 // :
90 // Initialize member data here.
91 {
92 
94  fDetector = lar::providerFrom<detinfo::DetectorPropertiesService>();
96 
97  fDetHalfHeight = geo->DetHalfHeight();
98  fDetWidth = 2.*geo->DetHalfWidth();
99  fDetLength = geo->DetLength();
100 
101  float fSamplingRate = fDetector->SamplingRate();
102 
103  fPFParticleModuleLabel = p.get<std::string >("PFParticleModuleLabel");
104  fPCAxisModuleLabel = p.get< std::string >("PCAxisModuleLabel");
105 
106  fPcaAlg.reconfigure(p.get<fhicl::ParameterSet>("PrincipalComponentsAlg"));
107 
108  fTPCXBoundary = p.get< float >("TPCXBoundary", 5);
109  fTPCYBoundary = p.get< float >("TPCYBoundary", 5);
110  fTPCZBoundary = p.get< float >("TPCZBoundary", 5);
111 
112  const double driftVelocity = fDetector->DriftVelocity( fDetector->Efield(), fDetector->Temperature() ); // cm/us
113 
114  //std::cerr << "Drift velocity is " << driftVelocity << " cm/us. Sampling rate is: "<< fSamplingRate << " detector width: " << 2*geo->DetHalfWidth() << std::endl;
115  fDetectorWidthTicks = 2*geo->DetHalfWidth()/(driftVelocity*fSamplingRate/1000); // ~3200 for uB
116 
117  // Call appropriate Produces<>() functions here.
118  produces< std::vector<anab::CosmicTag> >();
119  produces< art::Assns<recob::PFParticle, anab::CosmicTag> >();
120  produces< art::Assns<recob::PCAxis, anab::CosmicTag> >();
121 }
122 
124 {
125  // Instatiate the output
126  std::unique_ptr< std::vector< anab::CosmicTag > > cosmicTagPFParticleVector( new std::vector<anab::CosmicTag> );
127  std::unique_ptr< art::Assns<recob::PCAxis, anab::CosmicTag > > assnOutCosmicTagPCAxis( new art::Assns<recob::PCAxis, anab::CosmicTag>);
128  std::unique_ptr< art::Assns<recob::PFParticle, anab::CosmicTag > > assnOutCosmicTagPFParticle( new art::Assns<recob::PFParticle, anab::CosmicTag>);
129 
130  // Recover handle for PFParticles
132  evt.getByLabel( fPFParticleModuleLabel, pfParticleHandle);
133 
134  if (!pfParticleHandle.isValid())
135  {
136  evt.put( std::move(cosmicTagPFParticleVector) );
137  evt.put( std::move(assnOutCosmicTagPFParticle) );
138  evt.put( std::move(assnOutCosmicTagPCAxis) );
139  return;
140  }
141 
142  // We need a handle to the collection of clusters in the data store so we can
143  // handle associations to hits.
145  evt.getByLabel(fPFParticleModuleLabel, clusterHandle);
146 
147  // Recover the handle for the PCAxes
149  evt.getByLabel( fPCAxisModuleLabel, pcaxisHandle);
150 
151  if (!pcaxisHandle.isValid() || !clusterHandle.isValid())
152  {
153  evt.put( std::move(cosmicTagPFParticleVector) );
154  evt.put( std::move(assnOutCosmicTagPFParticle) );
155  evt.put( std::move(assnOutCosmicTagPCAxis) );
156  return;
157  }
158 
159  // Recover the list of associated PCA axes
160  art::FindManyP<recob::PCAxis> pfPartToPCAxisAssns(pfParticleHandle, evt, fPCAxisModuleLabel);
161 
162  // Add the relations to recover associations cluster hits
163  art::FindManyP<recob::SpacePoint> spacePointAssnVec(pfParticleHandle, evt, fPFParticleModuleLabel);
164 
165  // Recover the collection of associations between PFParticles and clusters, this will
166  // be the mechanism by which we actually deal with clusters
167  art::FindManyP<recob::Cluster> clusterAssns(pfParticleHandle, evt, fPFParticleModuleLabel);
168 
169  // Likewise, recover the collection of associations to hits
170  art::FindManyP<recob::Hit> clusterHitAssns(clusterHandle, evt, fPFParticleModuleLabel);
171 
172  // The outer loop is going to be over PFParticles
173  for(size_t pfPartIdx = 0; pfPartIdx != pfParticleHandle->size(); pfPartIdx++)
174  {
175  art::Ptr<recob::PFParticle> pfParticle(pfParticleHandle, pfPartIdx);
176 
177  // Require the PFParticle be a primary because attached hits will be daughter particles
178  //if (!pfParticle->IsPrimary()) continue;
179  // Oops! Don't do this unless traversing the heirarchy!
180 
181  // Recover the PCAxis vector
182  std::vector<art::Ptr<recob::PCAxis> > pcAxisVec = pfPartToPCAxisAssns.at(pfPartIdx);
183 
184  // Is there an axis associated to this PFParticle?
185  if (pcAxisVec.empty()) continue;
186 
187  // *****************************************************************************************
188  // For what follows below we want the "best" PCAxis object only. However, it can be that
189  // there are two PCAxes for a PFParticle (depending on source) where the "first" axis will
190  // be the "better" one that we want (this statement by fiat, it is defined that way in the
191  // axis producer module).
192  if (pcAxisVec.size() > 1 && pcAxisVec.front()->getID() > pcAxisVec.back()->getID()) std::reverse(pcAxisVec.begin(), pcAxisVec.end());
193  // We need to confirm this!!
194  // *****************************************************************************************
195 
196  // Recover the axis
197  const art::Ptr<recob::PCAxis>& pcAxis = pcAxisVec.front();
198 
199  // Start the tagging process...
200  int isCosmic = 0;
202 
203  // There are two sections to the tagging, in the first we are going to check for hits that are
204  // "out of time" and for this we only need the hit vector. If no hits are out of time then
205  // we need to do a more thorough check of the positions of the hits.
206  // If we find hits that are out of time we'll set the "end points" of our trajectory to
207  // a scale factor past the principle eigen value
208  double eigenVal0 = sqrt(pcAxis->getEigenValues()[0]);
209  double maxArcLen = 3. * eigenVal0;
210 
211  // Recover PCA end points
212  TVector3 vertexPosition(pcAxis->getAvePosition()[0],pcAxis->getAvePosition()[1],pcAxis->getAvePosition()[2]);
213  TVector3 vertexDirection(pcAxis->getEigenVectors()[0][0],pcAxis->getEigenVectors()[0][1],pcAxis->getEigenVectors()[0][2]);
214 
215  TVector3 pcAxisStart = vertexPosition - maxArcLen * vertexDirection;
216  TVector3 pcAxisEnd = vertexPosition + maxArcLen * vertexDirection;
217 
218  // "Track" end points in easily readable form
219  float trackEndPt1_X = pcAxisStart [0];
220  float trackEndPt1_Y = pcAxisStart [1];
221  float trackEndPt1_Z = pcAxisStart [2];
222  float trackEndPt2_X = pcAxisEnd[0];
223  float trackEndPt2_Y = pcAxisEnd[1];
224  float trackEndPt2_Z = pcAxisEnd[2];
225 
226  // Now we get the 2D clusters associated to this PFParticle
227  std::vector<art::Ptr<recob::Cluster> > clusterVec = clusterAssns.at(pfParticle.key());
228 
229  bool dumpMe(false);
230 
231  // Once we have the clusters then we can loop over them to find the associated hits
232  for(const auto& cluster : clusterVec)
233  {
234  // Recover the 2D hits associated to a given cluster
235  std::vector<art::Ptr<recob::Hit> > hitVec = clusterHitAssns.at(cluster->ID());
236 
237  // Once we have the hits the first thing we should do is to check if any are "out of time"
238  // If there are out of time hits then we are going to reject the cluster so no need to do
239  // any further processing.
241  // Check that all hits on particle are "in time"
243  for(const auto& hit : hitVec)
244  {
245  if (dumpMe)
246  {
247  std::cout << "***>> Hit key: " << hit.key() << ", peak - RMS: " << hit->PeakTimeMinusRMS() << ", peak + RMS: " << hit->PeakTimePlusRMS() << ", det width: " << fDetectorWidthTicks << std::endl;
248  }
249  //if( hit->StartTick() < fDetectorWidthTicks || hit->EndTick() > 2.*fDetectorWidthTicks)
250  if( hit->PeakTimeMinusRMS() < fDetectorWidthTicks || hit->PeakTimePlusRMS() > 2.*fDetectorWidthTicks)
251  {
252  isCosmic = 1;
254  break; // If one hit is out of time it must be a cosmic ray
255  }
256  }
257  }
258 
259  // Recover the space points associated to this PFParticle.
260  std::vector<art::Ptr<recob::SpacePoint> > spacePointVec = spacePointAssnVec.at(pfParticle.key());
261 
263  // Now check the TPC boundaries:
265  if(isCosmic==0 && !spacePointVec.empty())
266  {
267  // Do a check on the transverse components of the PCA axes to make sure we are looking at long straight
268  // tracks and not the kind of events we might want to keep
269  double transRMS = sqrt(std::pow(pcAxis->getEigenValues()[1],2) + std::pow(pcAxis->getEigenValues()[1],2));
270 
271  //if (eigenVal0 > 10. && transRMS < 10. && transRMS / eigenVal0 < 0.1)
272  if (eigenVal0 > 0. && transRMS > 0.)
273  {
274  // The idea is to find the maximum extents of this PFParticle using the PCA axis which we
275  // can then use to determine proximity to a TPC boundary.
276  // We implement this by recovering the 3D Space Points and then make a pass through them to
277  // find the space points at the extremes of the distance along the principle axis.
278  // We'll loop through the space points looking for those which have the largest arc lengths along
279  // the principle axis. Set up to do that
280  double arcLengthToFirstHit(9999.);
281  double arcLengthToLastHit(-9999.);
282 
283  for(const auto spacePoint : spacePointVec)
284  {
285  TVector3 spacePointPos(spacePoint->XYZ()[0],spacePoint->XYZ()[1],spacePoint->XYZ()[2]);
286  TVector3 deltaPos = spacePointPos - vertexPosition;
287  double arcLenToHit = deltaPos.Dot(vertexDirection);
288 
289  if (arcLenToHit < arcLengthToFirstHit)
290  {
291  arcLengthToFirstHit = arcLenToHit;
292  pcAxisStart = spacePointPos;
293  }
294 
295  if (arcLenToHit > arcLengthToLastHit)
296  {
297  arcLengthToLastHit = arcLenToHit;
298  pcAxisEnd = spacePointPos;
299  }
300  }
301 
302  // "Track" end points in easily readable form
303  trackEndPt1_X = pcAxisStart [0];
304  trackEndPt1_Y = pcAxisStart [1];
305  trackEndPt1_Z = pcAxisStart [2];
306  trackEndPt2_X = pcAxisEnd[0];
307  trackEndPt2_Y = pcAxisEnd[1];
308  trackEndPt2_Z = pcAxisEnd[2];
309 
310  // In below we check entry and exit points. Note that a special case of a particle entering
311  // and exiting the same surface is considered to be running parallel to the surface and NOT
312  // entering and exiting.
313  // Also, in what follows we make no assumptions on which end point is the "start" or
314  // "end" of the track being considered.
315  bool nBdX[] = {false,false};
316  bool nBdY[] = {false,false};
317  bool nBdZ[] = {false,false};
318 
319  // Check x extents - note that uboone coordinaes system has x=0 at edge
320  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
321  // Also note that, in theory, any cosmic ray entering or exiting the X surfaces will have presumably
322  // been removed already by the checking of "out of time" hits... but this will at least label
323  // neutrino interaction tracks which exit through the X surfaces of the TPC
324  if (fDetWidth - trackEndPt1_X < fTPCXBoundary || trackEndPt1_X < fTPCXBoundary) nBdX[0] = true;
325  if (fDetWidth - trackEndPt2_X < fTPCXBoundary || trackEndPt2_X < fTPCXBoundary) nBdX[1] = true;
326 
327  // Check y extents (note coordinate system change)
328  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
329  if (fDetHalfHeight - trackEndPt1_Y < fTPCYBoundary || fDetHalfHeight + trackEndPt1_Y < fTPCYBoundary) nBdY[0] = true; // one end of track exits out top
330  if (fDetHalfHeight - trackEndPt2_Y < fTPCYBoundary || fDetHalfHeight + trackEndPt2_Y < fTPCYBoundary) nBdY[1] = true; // one end of track exist out bottom
331 
332  // Check z extents
333  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
334  if (fDetLength - trackEndPt1_Z < fTPCZBoundary || trackEndPt1_Z < fTPCZBoundary ) nBdZ[0] = true;
335  if (fDetLength - trackEndPt2_Z < fTPCZBoundary || trackEndPt2_Z < fTPCZBoundary ) nBdZ[1] = true;
336 
337  // Endpoints exiting?
338  bool exitEnd1 = nBdX[0] || nBdY[0]; // end point 1 enters/exits top/bottom or x sides
339  bool exitEnd2 = nBdX[1] || nBdY[1]; // end point 2 enters/exits top/bottom or x sides
340  bool exitEndZ1 = exitEnd1 && nBdZ[1]; // end point 1 enters/exits top/bottom and exits/enters z
341  bool exitEndZ2 = exitEnd1 && nBdZ[0]; // end point 2 enters/exits top/bottom and exits/enters z
342 
343  // This should check for the case of a track which is both entering and exiting
344  // but we consider entering and exiting the z boundaries to be a special case (should it be?)
345  if((exitEnd1 && exitEnd2) || exitEndZ1 || exitEndZ2)
346  {
347  isCosmic = 2;
348  if (nBdX[0] && nBdX[1]) tag_id = anab::CosmicTagID_t::kGeometry_XX;
349  else if (nBdY[0] && nBdY[1]) tag_id = anab::CosmicTagID_t::kGeometry_YY;
350  else if ((nBdX[0] || nBdX[1]) && (nBdY[0] || nBdY[1])) tag_id = anab::CosmicTagID_t::kGeometry_XY;
351  else if ((nBdX[0] || nBdX[1]) && (nBdZ[0] || nBdZ[1])) tag_id = anab::CosmicTagID_t::kGeometry_XZ;
352  else tag_id = anab::CosmicTagID_t::kGeometry_YZ;
353  }
354  // This is the special case of track which appears to enter/exit z boundaries
355  else if ( nBdZ[0] && nBdZ[1])
356  {
357  isCosmic = 3;
359  }
360  // This looks for track which enters/exits a boundary but has other endpoint in TPC
361  else if ( (nBdX[0] || nBdY[0] || nBdZ[0]) != (nBdX[1] || nBdY[1] || nBdZ[1]))
362  {
363  isCosmic = 4 ;
364  if (nBdX[0] || nBdX[1]) tag_id = anab::CosmicTagID_t::kGeometry_X;
365  else if (nBdY[0] || nBdY[1]) tag_id = anab::CosmicTagID_t::kGeometry_Y;
366  else if (nBdZ[0] || nBdZ[1]) tag_id = anab::CosmicTagID_t::kGeometry_Z;
367  }
368  }
369  }
370 
371  std::vector<float> endPt1;
372  std::vector<float> endPt2;
373  endPt1.push_back( trackEndPt1_X );
374  endPt1.push_back( trackEndPt1_Y );
375  endPt1.push_back( trackEndPt1_Z );
376  endPt2.push_back( trackEndPt2_X );
377  endPt2.push_back( trackEndPt2_Y );
378  endPt2.push_back( trackEndPt2_Z );
379 
380  float cosmicScore = isCosmic > 0 ? 1. : 0.;
381 
382  // Handle special cases
383  if (isCosmic == 3) cosmicScore = 0.4; // Enter/Exit at opposite Z boundaries
384  else if (isCosmic == 4) cosmicScore = 0.5; // Enter or Exit but not both
385 
386  // Create the tag object for this PFParticle and make the corresponding association
387  cosmicTagPFParticleVector->emplace_back( endPt1, endPt2, cosmicScore, tag_id);
388 
389  util::CreateAssn(*this, evt, *cosmicTagPFParticleVector, pfParticle, *assnOutCosmicTagPFParticle );
390 
391  // Loop through the tracks resulting from this PFParticle and mark them
392  for(const auto& axis : pcAxisVec)
393  {
394  util::CreateAssn(*this, evt, *cosmicTagPFParticleVector, axis, *assnOutCosmicTagPCAxis );
395  }
396  }
397 
398  evt.put( std::move(cosmicTagPFParticleVector) );
399  evt.put( std::move(assnOutCosmicTagPFParticle) );
400  evt.put( std::move(assnOutCosmicTagPCAxis) );
401 
402  return;
403 
404 } // end of produce
406 
407 //------------------------------------------------------------------------------------------------------------------------------------------
409  Hit2DVector& clusterHitVec,
410  reco::Hit2DListPtr& hit2DListPtr) const
411 {
417  // We'll need the offsets for each plane
418  std::map<geo::View_t, double> viewOffsetMap;
419 
420  viewOffsetMap[geo::kU] = fDetector->GetXTicksOffset(geo::kU, 0, 0)-fDetector->TriggerOffset();
421  viewOffsetMap[geo::kV] = fDetector->GetXTicksOffset(geo::kV, 0, 0)-fDetector->TriggerOffset();
422  viewOffsetMap[geo::kW] = fDetector->GetXTicksOffset(geo::kW, 0, 0)-fDetector->TriggerOffset();
423 
424  // Reserve memory for the hit vector - this is the container of the actual hits
425  clusterHitVec.reserve(recobHitVec.size());
426 
427  // Reserve memory for the list of pointers to these hits - what we'll operate on
428 
429  // Cycle through the recob hits to build ClusterHit2D objects and insert
430  // them into the map
431  for(const auto& recobHit : recobHitVec)
432  {
433  const geo::WireID& hitWireID(recobHit->WireID());
434 
435  double hitPeakTime(recobHit->PeakTime() - viewOffsetMap[recobHit->View()]);
436  double xPosition(fDetector->ConvertTicksToX(recobHit->PeakTime(), hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat));
437 
438  clusterHitVec.emplace_back(reco::ClusterHit2D(0, 0., 0., xPosition, hitPeakTime, hitWireID, recobHit.get()));
439  hit2DListPtr.emplace_back(&clusterHitVec.back());
440  }
441 
442  return;
443 }
444 
const EigenVectors & getEigenVectors() const
Definition: PCAxis.h:66
std::list< const reco::ClusterHit2D * > Hit2DListPtr
export some data structure definitions
Definition: Cluster3D.h:334
const double * getEigenValues() const
Definition: PCAxis.h:65
virtual int TriggerOffset() const =0
Collect all the geometry header files together.
detinfo::DetectorProperties const * fDetector
Pointer to the detector properties.
void RecobToClusterHits(const art::PtrVector< recob::Hit > &recobHitVec, Hit2DVector &clusterHitVec, reco::Hit2DListPtr &hit2DListPtr) const
Planes which measure V.
Definition: geo_types.h:78
void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
enum anab::cosmic_tag_id CosmicTagID_t
Declaration of signal hit object.
EDProducer()=default
std::vector< reco::ClusterHit2D > Hit2DVector
virtual double SamplingRate() const =0
Returns the period of the TPC readout electronics clock.
Cluster finding and building.
bool isValid() const
Definition: Handle.h:183
Planes which measure U.
Definition: geo_types.h:77
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:437
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:56
CosmicPCAxisTagger(fhicl::ParameterSet const &p)
virtual double Temperature() const =0
key_type key() const
Definition: Ptr.h:238
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
Declaration of cluster object.
void produce(art::Event &e) override
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
This header file defines the interface to a principal components analysis designed to be used within ...
const double * getAvePosition() const
Definition: PCAxis.h:67
virtual double ConvertTicksToX(double ticks, int p, int t, int c) const =0
Utility object to perform functions of association.
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
lar_cluster3d::PrincipalComponentsAlg fPcaAlg
Principal Components algorithm.
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:79
virtual double Efield(unsigned int planegap=0) const =0
Returns the nominal electric field in the specified volume.
TCEvent evt
Definition: DataStructs.cxx:5
Float_t e
Definition: plot.C:34
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:692
Namespace collecting geometry-related classes utilities.
Definition: Track.hh:43
art framework interface to geometry description
virtual double GetXTicksOffset(int p, int t, int c) const =0