LArSoft  v08_29_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 
25 #include <algorithm>
26 #include <cmath>
27 #include <iostream>
28 #include <iterator>
29 
41 
42 #include "TVector3.h"
43 
44 namespace cosmic
45 {
46  class CosmicPCAxisTagger;
47  class SpacePoint;
48 }
49 
50 
52 {
53 public:
54  explicit CosmicPCAxisTagger(fhicl::ParameterSet const & p);
55 
56  void produce(art::Event & e) override;
57 
58 private:
59  typedef std::vector<reco::ClusterHit2D> Hit2DVector;
60 
61  void RecobToClusterHits(const art::PtrVector<recob::Hit>& recobHitVec,
62  Hit2DVector& clusterHitVec,
63  reco::Hit2DListPtr& hit2DListPtr) const;
64 
66  std::string fPCAxisModuleLabel;
67 
70 
74 };
75 
76 
78  EDProducer{p},
79  fDetector(lar::providerFrom<detinfo::DetectorPropertiesService>()),
80  fPcaAlg(p.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
81 // :
82 // Initialize member data here.
83 {
84 
86  fDetector = lar::providerFrom<detinfo::DetectorPropertiesService>();
88 
89  fDetHalfHeight = geo->DetHalfHeight();
90  fDetWidth = 2.*geo->DetHalfWidth();
91  fDetLength = geo->DetLength();
92 
93  float fSamplingRate = fDetector->SamplingRate();
94 
95  fPFParticleModuleLabel = p.get<std::string >("PFParticleModuleLabel");
96  fPCAxisModuleLabel = p.get< std::string >("PCAxisModuleLabel");
97 
98  fTPCXBoundary = p.get< float >("TPCXBoundary", 5);
99  fTPCYBoundary = p.get< float >("TPCYBoundary", 5);
100  fTPCZBoundary = p.get< float >("TPCZBoundary", 5);
101 
102  const double driftVelocity = fDetector->DriftVelocity( fDetector->Efield(), fDetector->Temperature() ); // cm/us
103 
104  //std::cerr << "Drift velocity is " << driftVelocity << " cm/us. Sampling rate is: "<< fSamplingRate << " detector width: " << 2*geo->DetHalfWidth() << std::endl;
105  fDetectorWidthTicks = 2*geo->DetHalfWidth()/(driftVelocity*fSamplingRate/1000); // ~3200 for uB
106 
107  // Call appropriate Produces<>() functions here.
108  produces< std::vector<anab::CosmicTag> >();
109  produces< art::Assns<recob::PFParticle, anab::CosmicTag> >();
110  produces< art::Assns<recob::PCAxis, anab::CosmicTag> >();
111 }
112 
114 {
115  // Instatiate the output
116  std::unique_ptr< std::vector< anab::CosmicTag > > cosmicTagPFParticleVector( new std::vector<anab::CosmicTag> );
117  std::unique_ptr< art::Assns<recob::PCAxis, anab::CosmicTag > > assnOutCosmicTagPCAxis( new art::Assns<recob::PCAxis, anab::CosmicTag>);
118  std::unique_ptr< art::Assns<recob::PFParticle, anab::CosmicTag > > assnOutCosmicTagPFParticle( new art::Assns<recob::PFParticle, anab::CosmicTag>);
119 
120  // Recover handle for PFParticles
122  evt.getByLabel( fPFParticleModuleLabel, pfParticleHandle);
123 
124  if (!pfParticleHandle.isValid())
125  {
126  evt.put( std::move(cosmicTagPFParticleVector) );
127  evt.put( std::move(assnOutCosmicTagPFParticle) );
128  evt.put( std::move(assnOutCosmicTagPCAxis) );
129  return;
130  }
131 
132  // We need a handle to the collection of clusters in the data store so we can
133  // handle associations to hits.
135  evt.getByLabel(fPFParticleModuleLabel, clusterHandle);
136 
137  // Recover the handle for the PCAxes
139  evt.getByLabel( fPCAxisModuleLabel, pcaxisHandle);
140 
141  if (!pcaxisHandle.isValid() || !clusterHandle.isValid())
142  {
143  evt.put( std::move(cosmicTagPFParticleVector) );
144  evt.put( std::move(assnOutCosmicTagPFParticle) );
145  evt.put( std::move(assnOutCosmicTagPCAxis) );
146  return;
147  }
148 
149  // Recover the list of associated PCA axes
150  art::FindManyP<recob::PCAxis> pfPartToPCAxisAssns(pfParticleHandle, evt, fPCAxisModuleLabel);
151 
152  // Add the relations to recover associations cluster hits
153  art::FindManyP<recob::SpacePoint> spacePointAssnVec(pfParticleHandle, evt, fPFParticleModuleLabel);
154 
155  // Recover the collection of associations between PFParticles and clusters, this will
156  // be the mechanism by which we actually deal with clusters
157  art::FindManyP<recob::Cluster> clusterAssns(pfParticleHandle, evt, fPFParticleModuleLabel);
158 
159  // Likewise, recover the collection of associations to hits
160  art::FindManyP<recob::Hit> clusterHitAssns(clusterHandle, evt, fPFParticleModuleLabel);
161 
162  // The outer loop is going to be over PFParticles
163  for(size_t pfPartIdx = 0; pfPartIdx != pfParticleHandle->size(); pfPartIdx++)
164  {
165  art::Ptr<recob::PFParticle> pfParticle(pfParticleHandle, pfPartIdx);
166 
167  // Require the PFParticle be a primary because attached hits will be daughter particles
168  //if (!pfParticle->IsPrimary()) continue;
169  // Oops! Don't do this unless traversing the heirarchy!
170 
171  // Recover the PCAxis vector
172  std::vector<art::Ptr<recob::PCAxis> > pcAxisVec = pfPartToPCAxisAssns.at(pfPartIdx);
173 
174  // Is there an axis associated to this PFParticle?
175  if (pcAxisVec.empty()) continue;
176 
177  // *****************************************************************************************
178  // For what follows below we want the "best" PCAxis object only. However, it can be that
179  // there are two PCAxes for a PFParticle (depending on source) where the "first" axis will
180  // be the "better" one that we want (this statement by fiat, it is defined that way in the
181  // axis producer module).
182  if (pcAxisVec.size() > 1 && pcAxisVec.front()->getID() > pcAxisVec.back()->getID()) std::reverse(pcAxisVec.begin(), pcAxisVec.end());
183  // We need to confirm this!!
184  // *****************************************************************************************
185 
186  // Recover the axis
187  const art::Ptr<recob::PCAxis>& pcAxis = pcAxisVec.front();
188 
189  // Start the tagging process...
190  int isCosmic = 0;
192 
193  // There are two sections to the tagging, in the first we are going to check for hits that are
194  // "out of time" and for this we only need the hit vector. If no hits are out of time then
195  // we need to do a more thorough check of the positions of the hits.
196  // If we find hits that are out of time we'll set the "end points" of our trajectory to
197  // a scale factor past the principle eigen value
198  double eigenVal0 = sqrt(pcAxis->getEigenValues()[0]);
199  double maxArcLen = 3. * eigenVal0;
200 
201  // Recover PCA end points
202  TVector3 vertexPosition(pcAxis->getAvePosition()[0],pcAxis->getAvePosition()[1],pcAxis->getAvePosition()[2]);
203  TVector3 vertexDirection(pcAxis->getEigenVectors()[0][0],pcAxis->getEigenVectors()[0][1],pcAxis->getEigenVectors()[0][2]);
204 
205  TVector3 pcAxisStart = vertexPosition - maxArcLen * vertexDirection;
206  TVector3 pcAxisEnd = vertexPosition + maxArcLen * vertexDirection;
207 
208  // "Track" end points in easily readable form
209  float trackEndPt1_X = pcAxisStart [0];
210  float trackEndPt1_Y = pcAxisStart [1];
211  float trackEndPt1_Z = pcAxisStart [2];
212  float trackEndPt2_X = pcAxisEnd[0];
213  float trackEndPt2_Y = pcAxisEnd[1];
214  float trackEndPt2_Z = pcAxisEnd[2];
215 
216  // Now we get the 2D clusters associated to this PFParticle
217  std::vector<art::Ptr<recob::Cluster> > clusterVec = clusterAssns.at(pfParticle.key());
218 
219  bool dumpMe(false);
220 
221  // Once we have the clusters then we can loop over them to find the associated hits
222  for(const auto& cluster : clusterVec)
223  {
224  // Recover the 2D hits associated to a given cluster
225  std::vector<art::Ptr<recob::Hit> > hitVec = clusterHitAssns.at(cluster->ID());
226 
227  // Once we have the hits the first thing we should do is to check if any are "out of time"
228  // If there are out of time hits then we are going to reject the cluster so no need to do
229  // any further processing.
231  // Check that all hits on particle are "in time"
233  for(const auto& hit : hitVec)
234  {
235  if (dumpMe)
236  {
237  std::cout << "***>> Hit key: " << hit.key() << ", peak - RMS: " << hit->PeakTimeMinusRMS() << ", peak + RMS: " << hit->PeakTimePlusRMS() << ", det width: " << fDetectorWidthTicks << std::endl;
238  }
239  //if( hit->StartTick() < fDetectorWidthTicks || hit->EndTick() > 2.*fDetectorWidthTicks)
240  if( hit->PeakTimeMinusRMS() < fDetectorWidthTicks || hit->PeakTimePlusRMS() > 2.*fDetectorWidthTicks)
241  {
242  isCosmic = 1;
244  break; // If one hit is out of time it must be a cosmic ray
245  }
246  }
247  }
248 
249  // Recover the space points associated to this PFParticle.
250  std::vector<art::Ptr<recob::SpacePoint> > spacePointVec = spacePointAssnVec.at(pfParticle.key());
251 
253  // Now check the TPC boundaries:
255  if(isCosmic==0 && !spacePointVec.empty())
256  {
257  // Do a check on the transverse components of the PCA axes to make sure we are looking at long straight
258  // tracks and not the kind of events we might want to keep
259  double transRMS = sqrt(std::pow(pcAxis->getEigenValues()[1],2) + std::pow(pcAxis->getEigenValues()[1],2));
260 
261  //if (eigenVal0 > 10. && transRMS < 10. && transRMS / eigenVal0 < 0.1)
262  if (eigenVal0 > 0. && transRMS > 0.)
263  {
264  // The idea is to find the maximum extents of this PFParticle using the PCA axis which we
265  // can then use to determine proximity to a TPC boundary.
266  // We implement this by recovering the 3D Space Points and then make a pass through them to
267  // find the space points at the extremes of the distance along the principle axis.
268  // We'll loop through the space points looking for those which have the largest arc lengths along
269  // the principle axis. Set up to do that
270  double arcLengthToFirstHit(9999.);
271  double arcLengthToLastHit(-9999.);
272 
273  for(const auto spacePoint : spacePointVec)
274  {
275  TVector3 spacePointPos(spacePoint->XYZ()[0],spacePoint->XYZ()[1],spacePoint->XYZ()[2]);
276  TVector3 deltaPos = spacePointPos - vertexPosition;
277  double arcLenToHit = deltaPos.Dot(vertexDirection);
278 
279  if (arcLenToHit < arcLengthToFirstHit)
280  {
281  arcLengthToFirstHit = arcLenToHit;
282  pcAxisStart = spacePointPos;
283  }
284 
285  if (arcLenToHit > arcLengthToLastHit)
286  {
287  arcLengthToLastHit = arcLenToHit;
288  pcAxisEnd = spacePointPos;
289  }
290  }
291 
292  // "Track" end points in easily readable form
293  trackEndPt1_X = pcAxisStart [0];
294  trackEndPt1_Y = pcAxisStart [1];
295  trackEndPt1_Z = pcAxisStart [2];
296  trackEndPt2_X = pcAxisEnd[0];
297  trackEndPt2_Y = pcAxisEnd[1];
298  trackEndPt2_Z = pcAxisEnd[2];
299 
300  // In below we check entry and exit points. Note that a special case of a particle entering
301  // and exiting the same surface is considered to be running parallel to the surface and NOT
302  // entering and exiting.
303  // Also, in what follows we make no assumptions on which end point is the "start" or
304  // "end" of the track being considered.
305  bool nBdX[] = {false,false};
306  bool nBdY[] = {false,false};
307  bool nBdZ[] = {false,false};
308 
309  // Check x extents - note that uboone coordinaes system has x=0 at edge
310  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
311  // Also note that, in theory, any cosmic ray entering or exiting the X surfaces will have presumably
312  // been removed already by the checking of "out of time" hits... but this will at least label
313  // neutrino interaction tracks which exit through the X surfaces of the TPC
314  if (fDetWidth - trackEndPt1_X < fTPCXBoundary || trackEndPt1_X < fTPCXBoundary) nBdX[0] = true;
315  if (fDetWidth - trackEndPt2_X < fTPCXBoundary || trackEndPt2_X < fTPCXBoundary) nBdX[1] = true;
316 
317  // Check y extents (note coordinate system change)
318  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
319  if (fDetHalfHeight - trackEndPt1_Y < fTPCYBoundary || fDetHalfHeight + trackEndPt1_Y < fTPCYBoundary) nBdY[0] = true; // one end of track exits out top
320  if (fDetHalfHeight - trackEndPt2_Y < fTPCYBoundary || fDetHalfHeight + trackEndPt2_Y < fTPCYBoundary) nBdY[1] = true; // one end of track exist out bottom
321 
322  // Check z extents
323  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
324  if (fDetLength - trackEndPt1_Z < fTPCZBoundary || trackEndPt1_Z < fTPCZBoundary ) nBdZ[0] = true;
325  if (fDetLength - trackEndPt2_Z < fTPCZBoundary || trackEndPt2_Z < fTPCZBoundary ) nBdZ[1] = true;
326 
327  // Endpoints exiting?
328  bool exitEnd1 = nBdX[0] || nBdY[0]; // end point 1 enters/exits top/bottom or x sides
329  bool exitEnd2 = nBdX[1] || nBdY[1]; // end point 2 enters/exits top/bottom or x sides
330  bool exitEndZ1 = exitEnd1 && nBdZ[1]; // end point 1 enters/exits top/bottom and exits/enters z
331  bool exitEndZ2 = exitEnd1 && nBdZ[0]; // end point 2 enters/exits top/bottom and exits/enters z
332 
333  // This should check for the case of a track which is both entering and exiting
334  // but we consider entering and exiting the z boundaries to be a special case (should it be?)
335  if((exitEnd1 && exitEnd2) || exitEndZ1 || exitEndZ2)
336  {
337  isCosmic = 2;
338  if (nBdX[0] && nBdX[1]) tag_id = anab::CosmicTagID_t::kGeometry_XX;
339  else if (nBdY[0] && nBdY[1]) tag_id = anab::CosmicTagID_t::kGeometry_YY;
340  else if ((nBdX[0] || nBdX[1]) && (nBdY[0] || nBdY[1])) tag_id = anab::CosmicTagID_t::kGeometry_XY;
341  else if ((nBdX[0] || nBdX[1]) && (nBdZ[0] || nBdZ[1])) tag_id = anab::CosmicTagID_t::kGeometry_XZ;
342  else tag_id = anab::CosmicTagID_t::kGeometry_YZ;
343  }
344  // This is the special case of track which appears to enter/exit z boundaries
345  else if ( nBdZ[0] && nBdZ[1])
346  {
347  isCosmic = 3;
349  }
350  // This looks for track which enters/exits a boundary but has other endpoint in TPC
351  else if ( (nBdX[0] || nBdY[0] || nBdZ[0]) != (nBdX[1] || nBdY[1] || nBdZ[1]))
352  {
353  isCosmic = 4 ;
354  if (nBdX[0] || nBdX[1]) tag_id = anab::CosmicTagID_t::kGeometry_X;
355  else if (nBdY[0] || nBdY[1]) tag_id = anab::CosmicTagID_t::kGeometry_Y;
356  else if (nBdZ[0] || nBdZ[1]) tag_id = anab::CosmicTagID_t::kGeometry_Z;
357  }
358  }
359  }
360 
361  std::vector<float> endPt1;
362  std::vector<float> endPt2;
363  endPt1.push_back( trackEndPt1_X );
364  endPt1.push_back( trackEndPt1_Y );
365  endPt1.push_back( trackEndPt1_Z );
366  endPt2.push_back( trackEndPt2_X );
367  endPt2.push_back( trackEndPt2_Y );
368  endPt2.push_back( trackEndPt2_Z );
369 
370  float cosmicScore = isCosmic > 0 ? 1. : 0.;
371 
372  // Handle special cases
373  if (isCosmic == 3) cosmicScore = 0.4; // Enter/Exit at opposite Z boundaries
374  else if (isCosmic == 4) cosmicScore = 0.5; // Enter or Exit but not both
375 
376  // Create the tag object for this PFParticle and make the corresponding association
377  cosmicTagPFParticleVector->emplace_back( endPt1, endPt2, cosmicScore, tag_id);
378 
379  util::CreateAssn(*this, evt, *cosmicTagPFParticleVector, pfParticle, *assnOutCosmicTagPFParticle );
380 
381  // Loop through the tracks resulting from this PFParticle and mark them
382  for(const auto& axis : pcAxisVec)
383  {
384  util::CreateAssn(*this, evt, *cosmicTagPFParticleVector, axis, *assnOutCosmicTagPCAxis );
385  }
386  }
387 
388  evt.put( std::move(cosmicTagPFParticleVector) );
389  evt.put( std::move(assnOutCosmicTagPFParticle) );
390  evt.put( std::move(assnOutCosmicTagPCAxis) );
391 
392  return;
393 
394 } // end of produce
396 
397 //------------------------------------------------------------------------------------------------------------------------------------------
399  Hit2DVector& clusterHitVec,
400  reco::Hit2DListPtr& hit2DListPtr) const
401 {
407  // We'll need the offsets for each plane
408  std::map<geo::View_t, double> viewOffsetMap;
409 
410  viewOffsetMap[geo::kU] = fDetector->GetXTicksOffset(geo::kU, 0, 0)-fDetector->TriggerOffset();
411  viewOffsetMap[geo::kV] = fDetector->GetXTicksOffset(geo::kV, 0, 0)-fDetector->TriggerOffset();
412  viewOffsetMap[geo::kW] = fDetector->GetXTicksOffset(geo::kW, 0, 0)-fDetector->TriggerOffset();
413 
414  // Reserve memory for the hit vector - this is the container of the actual hits
415  clusterHitVec.reserve(recobHitVec.size());
416 
417  // Reserve memory for the list of pointers to these hits - what we'll operate on
418 
419  // Cycle through the recob hits to build ClusterHit2D objects and insert
420  // them into the map
421  for(const auto& recobHit : recobHitVec)
422  {
423  const geo::WireID& hitWireID(recobHit->WireID());
424 
425  double hitPeakTime(recobHit->PeakTime() - viewOffsetMap[recobHit->View()]);
426  double xPosition(fDetector->ConvertTicksToX(recobHit->PeakTime(), hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat));
427 
428  clusterHitVec.emplace_back(reco::ClusterHit2D(0, 0., 0., xPosition, hitPeakTime, hitWireID, recobHit.get()));
429  hit2DListPtr.emplace_back(&clusterHitVec.back());
430  }
431 
432  return;
433 }
434 
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
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
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:7
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.
art framework interface to geometry description
virtual double GetXTicksOffset(int p, int t, int c) const =0