Public Member Functions | Private Types | Private Attributes | List of all members
cosmic::CosmicPCAxisTagger Class Reference
Inheritance diagram for cosmic::CosmicPCAxisTagger:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 CosmicPCAxisTagger (fhicl::ParameterSet const &p)
 
void produce (art::Event &e) override
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Types

typedef std::vector< reco::ClusterHit2DHit2DVector
 

Private Attributes

std::string fPFParticleModuleLabel
 
std::string fPCAxisModuleLabel
 
lar_cluster3d::PrincipalComponentsAlg fPcaAlg
 Principal Components algorithm. More...
 
int fDetectorWidthTicks
 
float fTPCXBoundary
 
float fTPCYBoundary
 
float fTPCZBoundary
 
float fDetHalfHeight
 
float fDetWidth
 
float fDetLength
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 50 of file CosmicPCAxisTagger_module.cc.

Member Typedef Documentation

Definition at line 57 of file CosmicPCAxisTagger_module.cc.

Constructor & Destructor Documentation

cosmic::CosmicPCAxisTagger::CosmicPCAxisTagger ( fhicl::ParameterSet const &  p)
explicit

Definition at line 69 of file CosmicPCAxisTagger_module.cc.

70  : EDProducer{p}, fPcaAlg(p.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
71 {
73 
75  fDetWidth = 2. * geo->DetHalfWidth();
76  fDetLength = geo->DetLength();
77 
78  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
79  float fSamplingRate = sampling_rate(clock_data);
80 
81  fPFParticleModuleLabel = p.get<std::string>("PFParticleModuleLabel");
82  fPCAxisModuleLabel = p.get<std::string>("PCAxisModuleLabel");
83 
84  fTPCXBoundary = p.get<float>("TPCXBoundary", 5);
85  fTPCYBoundary = p.get<float>("TPCYBoundary", 5);
86  fTPCZBoundary = p.get<float>("TPCZBoundary", 5);
87 
88  auto const detector =
90  const double driftVelocity =
91  detector.DriftVelocity(detector.Efield(), detector.Temperature()); // cm/us
92 
94  2 * geo->DetHalfWidth() / (driftVelocity * fSamplingRate / 1000); // ~3200 for uB
95 
96  produces<std::vector<anab::CosmicTag>>();
97  produces<art::Assns<recob::PFParticle, anab::CosmicTag>>();
98  produces<art::Assns<recob::PCAxis, anab::CosmicTag>>();
99 }
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
std::string string
Definition: nybbler.cc:12
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
p
Definition: test.py:223
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
lar_cluster3d::PrincipalComponentsAlg fPcaAlg
Principal Components algorithm.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
LArSoft geometry interface.
Definition: ChannelGeo.h:16

Member Function Documentation

void cosmic::CosmicPCAxisTagger::produce ( art::Event e)
overridevirtual

Implements art::EDProducer.

Definition at line 102 of file CosmicPCAxisTagger_module.cc.

103 {
104  // Instatiate the output
105  std::unique_ptr<std::vector<anab::CosmicTag>> cosmicTagPFParticleVector(
106  new std::vector<anab::CosmicTag>);
107  std::unique_ptr<art::Assns<recob::PCAxis, anab::CosmicTag>> assnOutCosmicTagPCAxis(
109  std::unique_ptr<art::Assns<recob::PFParticle, anab::CosmicTag>> assnOutCosmicTagPFParticle(
111 
112  // Recover handle for PFParticles
114  evt.getByLabel(fPFParticleModuleLabel, pfParticleHandle);
115 
116  if (!pfParticleHandle.isValid()) {
117  evt.put(std::move(cosmicTagPFParticleVector));
118  evt.put(std::move(assnOutCosmicTagPFParticle));
119  evt.put(std::move(assnOutCosmicTagPCAxis));
120  return;
121  }
122 
123  // We need a handle to the collection of clusters in the data store so we can
124  // handle associations to hits.
126  evt.getByLabel(fPFParticleModuleLabel, clusterHandle);
127 
128  // Recover the handle for the PCAxes
130  evt.getByLabel(fPCAxisModuleLabel, pcaxisHandle);
131 
132  if (!pcaxisHandle.isValid() || !clusterHandle.isValid()) {
133  evt.put(std::move(cosmicTagPFParticleVector));
134  evt.put(std::move(assnOutCosmicTagPFParticle));
135  evt.put(std::move(assnOutCosmicTagPCAxis));
136  return;
137  }
138 
139  // Recover the list of associated PCA axes
140  art::FindManyP<recob::PCAxis> pfPartToPCAxisAssns(pfParticleHandle, evt, fPCAxisModuleLabel);
141 
142  // Add the relations to recover associations cluster hits
143  art::FindManyP<recob::SpacePoint> spacePointAssnVec(
144  pfParticleHandle, evt, fPFParticleModuleLabel);
145 
146  // Recover the collection of associations between PFParticles and clusters, this will
147  // be the mechanism by which we actually deal with clusters
148  art::FindManyP<recob::Cluster> clusterAssns(pfParticleHandle, evt, fPFParticleModuleLabel);
149 
150  // Likewise, recover the collection of associations to hits
151  art::FindManyP<recob::Hit> clusterHitAssns(clusterHandle, evt, fPFParticleModuleLabel);
152 
153  // The outer loop is going to be over PFParticles
154  for (size_t pfPartIdx = 0; pfPartIdx != pfParticleHandle->size(); pfPartIdx++) {
155  art::Ptr<recob::PFParticle> pfParticle(pfParticleHandle, pfPartIdx);
156 
157  // Recover the PCAxis vector
158  std::vector<art::Ptr<recob::PCAxis>> pcAxisVec = pfPartToPCAxisAssns.at(pfPartIdx);
159 
160  // Is there an axis associated to this PFParticle?
161  if (pcAxisVec.empty()) continue;
162 
163  // *****************************************************************************************
164  // For what follows below we want the "best" PCAxis object only. However, it can be that
165  // there are two PCAxes for a PFParticle (depending on source) where the "first" axis will
166  // be the "better" one that we want (this statement by fiat, it is defined that way in the
167  // axis producer module).
168  if (pcAxisVec.size() > 1 && pcAxisVec.front()->getID() > pcAxisVec.back()->getID())
169  std::reverse(pcAxisVec.begin(), pcAxisVec.end());
170  // We need to confirm this!!
171  // *****************************************************************************************
172 
173  // Recover the axis
174  const art::Ptr<recob::PCAxis>& pcAxis = pcAxisVec.front();
175 
176  // Start the tagging process...
177  int isCosmic = 0;
179 
180  // There are two sections to the tagging, in the first we are going to check for hits that are
181  // "out of time" and for this we only need the hit vector. If no hits are out of time then
182  // we need to do a more thorough check of the positions of the hits.
183  // If we find hits that are out of time we'll set the "end points" of our trajectory to
184  // a scale factor past the principle eigen value
185  double eigenVal0 = sqrt(pcAxis->getEigenValues()[0]);
186  double maxArcLen = 3. * eigenVal0;
187 
188  // Recover PCA end points
189  TVector3 vertexPosition(
190  pcAxis->getAvePosition()[0], pcAxis->getAvePosition()[1], pcAxis->getAvePosition()[2]);
191  TVector3 vertexDirection(pcAxis->getEigenVectors()[0][0],
192  pcAxis->getEigenVectors()[0][1],
193  pcAxis->getEigenVectors()[0][2]);
194 
195  TVector3 pcAxisStart = vertexPosition - maxArcLen * vertexDirection;
196  TVector3 pcAxisEnd = vertexPosition + maxArcLen * vertexDirection;
197 
198  // "Track" end points in easily readable form
199  float trackEndPt1_X = pcAxisStart[0];
200  float trackEndPt1_Y = pcAxisStart[1];
201  float trackEndPt1_Z = pcAxisStart[2];
202  float trackEndPt2_X = pcAxisEnd[0];
203  float trackEndPt2_Y = pcAxisEnd[1];
204  float trackEndPt2_Z = pcAxisEnd[2];
205 
206  // Now we get the 2D clusters associated to this PFParticle
207  std::vector<art::Ptr<recob::Cluster>> clusterVec = clusterAssns.at(pfParticle.key());
208 
209  bool dumpMe(false);
210 
211  // Once we have the clusters then we can loop over them to find the associated hits
212  for (const auto& cluster : clusterVec) {
213  // Recover the 2D hits associated to a given cluster
214  std::vector<art::Ptr<recob::Hit>> hitVec = clusterHitAssns.at(cluster->ID());
215 
216  // Once we have the hits the first thing we should do is to check if any are "out of time"
217  // If there are out of time hits then we are going to reject the cluster so no need to do
218  // any further processing.
219  /////////////////////////////////////
220  // Check that all hits on particle are "in time"
221  /////////////////////////////////////
222  for (const auto& hit : hitVec) {
223  if (dumpMe) {
224  std::cout << "***>> Hit key: " << hit.key() << ", peak - RMS: " << hit->PeakTimeMinusRMS()
225  << ", peak + RMS: " << hit->PeakTimePlusRMS()
226  << ", det width: " << fDetectorWidthTicks << std::endl;
227  }
228  if (hit->PeakTimeMinusRMS() < fDetectorWidthTicks ||
229  hit->PeakTimePlusRMS() > 2. * fDetectorWidthTicks) {
230  isCosmic = 1;
232  break; // If one hit is out of time it must be a cosmic ray
233  }
234  }
235  }
236 
237  // Recover the space points associated to this PFParticle.
238  std::vector<art::Ptr<recob::SpacePoint>> spacePointVec = spacePointAssnVec.at(pfParticle.key());
239 
240  /////////////////////////////////
241  // Now check the TPC boundaries:
242  /////////////////////////////////
243  if (isCosmic == 0 && !spacePointVec.empty()) {
244  // Do a check on the transverse components of the PCA axes to make sure we are looking at long straight
245  // tracks and not the kind of events we might want to keep
246  double transRMS =
247  sqrt(std::pow(pcAxis->getEigenValues()[1], 2) + std::pow(pcAxis->getEigenValues()[1], 2));
248 
249  if (eigenVal0 > 0. && transRMS > 0.) {
250  // The idea is to find the maximum extents of this PFParticle using the PCA axis which we
251  // can then use to determine proximity to a TPC boundary.
252  // We implement this by recovering the 3D Space Points and then make a pass through them to
253  // find the space points at the extremes of the distance along the principle axis.
254  // We'll loop through the space points looking for those which have the largest arc lengths along
255  // the principle axis. Set up to do that
256  double arcLengthToFirstHit(9999.);
257  double arcLengthToLastHit(-9999.);
258 
259  for (const auto spacePoint : spacePointVec) {
260  TVector3 spacePointPos(spacePoint->XYZ()[0], spacePoint->XYZ()[1], spacePoint->XYZ()[2]);
261  TVector3 deltaPos = spacePointPos - vertexPosition;
262  double arcLenToHit = deltaPos.Dot(vertexDirection);
263 
264  if (arcLenToHit < arcLengthToFirstHit) {
265  arcLengthToFirstHit = arcLenToHit;
266  pcAxisStart = spacePointPos;
267  }
268 
269  if (arcLenToHit > arcLengthToLastHit) {
270  arcLengthToLastHit = arcLenToHit;
271  pcAxisEnd = spacePointPos;
272  }
273  }
274 
275  // "Track" end points in easily readable form
276  trackEndPt1_X = pcAxisStart[0];
277  trackEndPt1_Y = pcAxisStart[1];
278  trackEndPt1_Z = pcAxisStart[2];
279  trackEndPt2_X = pcAxisEnd[0];
280  trackEndPt2_Y = pcAxisEnd[1];
281  trackEndPt2_Z = pcAxisEnd[2];
282 
283  // In below we check entry and exit points. Note that a special case of a particle entering
284  // and exiting the same surface is considered to be running parallel to the surface and NOT
285  // entering and exiting.
286  // Also, in what follows we make no assumptions on which end point is the "start" or
287  // "end" of the track being considered.
288  bool nBdX[] = {false, false};
289  bool nBdY[] = {false, false};
290  bool nBdZ[] = {false, false};
291 
292  // Check x extents - note that uboone coordinaes system has x=0 at edge
293  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
294  // Also note that, in theory, any cosmic ray entering or exiting the X surfaces will have presumably
295  // been removed already by the checking of "out of time" hits... but this will at least label
296  // neutrino interaction tracks which exit through the X surfaces of the TPC
297  if (fDetWidth - trackEndPt1_X < fTPCXBoundary || trackEndPt1_X < fTPCXBoundary)
298  nBdX[0] = true;
299  if (fDetWidth - trackEndPt2_X < fTPCXBoundary || trackEndPt2_X < fTPCXBoundary)
300  nBdX[1] = true;
301 
302  // Check y extents (note coordinate system change)
303  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
304  if (fDetHalfHeight - trackEndPt1_Y < fTPCYBoundary ||
305  fDetHalfHeight + trackEndPt1_Y < fTPCYBoundary)
306  nBdY[0] = true; // one end of track exits out top
307  if (fDetHalfHeight - trackEndPt2_Y < fTPCYBoundary ||
308  fDetHalfHeight + trackEndPt2_Y < fTPCYBoundary)
309  nBdY[1] = true; // one end of track exist out bottom
310 
311  // Check z extents
312  // Note this counts the case where the track enters and exits the same surface as a "1", not a "2"
313  if (fDetLength - trackEndPt1_Z < fTPCZBoundary || trackEndPt1_Z < fTPCZBoundary)
314  nBdZ[0] = true;
315  if (fDetLength - trackEndPt2_Z < fTPCZBoundary || trackEndPt2_Z < fTPCZBoundary)
316  nBdZ[1] = true;
317 
318  // Endpoints exiting?
319  bool exitEnd1 = nBdX[0] || nBdY[0]; // end point 1 enters/exits top/bottom or x sides
320  bool exitEnd2 = nBdX[1] || nBdY[1]; // end point 2 enters/exits top/bottom or x sides
321  bool exitEndZ1 =
322  exitEnd1 && nBdZ[1]; // end point 1 enters/exits top/bottom and exits/enters z
323  bool exitEndZ2 =
324  exitEnd1 && nBdZ[0]; // end point 2 enters/exits top/bottom and exits/enters z
325 
326  // This should check for the case of a track which is both entering and exiting
327  // but we consider entering and exiting the z boundaries to be a special case (should it be?)
328  if ((exitEnd1 && exitEnd2) || exitEndZ1 || exitEndZ2) {
329  isCosmic = 2;
330  if (nBdX[0] && nBdX[1])
332  else if (nBdY[0] && nBdY[1])
334  else if ((nBdX[0] || nBdX[1]) && (nBdY[0] || nBdY[1]))
336  else if ((nBdX[0] || nBdX[1]) && (nBdZ[0] || nBdZ[1]))
338  else
340  }
341  // This is the special case of track which appears to enter/exit z boundaries
342  else if (nBdZ[0] && nBdZ[1]) {
343  isCosmic = 3;
345  }
346  // This looks for track which enters/exits a boundary but has other endpoint in TPC
347  else if ((nBdX[0] || nBdY[0] || nBdZ[0]) != (nBdX[1] || nBdY[1] || nBdZ[1])) {
348  isCosmic = 4;
349  if (nBdX[0] || nBdX[1])
351  else if (nBdY[0] || nBdY[1])
353  else if (nBdZ[0] || nBdZ[1])
355  }
356  }
357  }
358 
359  std::vector<float> endPt1;
360  std::vector<float> endPt2;
361  endPt1.push_back(trackEndPt1_X);
362  endPt1.push_back(trackEndPt1_Y);
363  endPt1.push_back(trackEndPt1_Z);
364  endPt2.push_back(trackEndPt2_X);
365  endPt2.push_back(trackEndPt2_Y);
366  endPt2.push_back(trackEndPt2_Z);
367 
368  float cosmicScore = isCosmic > 0 ? 1. : 0.;
369 
370  // Handle special cases
371  if (isCosmic == 3)
372  cosmicScore = 0.4; // Enter/Exit at opposite Z boundaries
373  else if (isCosmic == 4)
374  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 
380  *this, evt, *cosmicTagPFParticleVector, pfParticle, *assnOutCosmicTagPFParticle);
381 
382  // Loop through the tracks resulting from this PFParticle and mark them
383  for (const auto& axis : pcAxisVec) {
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 } // end of produce
const EigenVectors & getEigenVectors() const
Definition: PCAxis.h:66
const double * getEigenValues() const
Definition: PCAxis.h:65
enum anab::cosmic_tag_id CosmicTagID_t
constexpr T pow(T x)
Definition: pow.h:72
Cluster finding and building.
bool isValid() const noexcept
Definition: Handle.h:191
def move(depos, offset)
Definition: depos.py:107
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.
Detector simulation of raw signals on wires.
const double * getAvePosition() const
Definition: PCAxis.h:67
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
Definition: qstring.cpp:11649
TCEvent evt
Definition: DataStructs.cxx:7
QTextStream & endl(QTextStream &s)

Member Data Documentation

int cosmic::CosmicPCAxisTagger::fDetectorWidthTicks
private

Definition at line 64 of file CosmicPCAxisTagger_module.cc.

float cosmic::CosmicPCAxisTagger::fDetHalfHeight
private

Definition at line 66 of file CosmicPCAxisTagger_module.cc.

float cosmic::CosmicPCAxisTagger::fDetLength
private

Definition at line 66 of file CosmicPCAxisTagger_module.cc.

float cosmic::CosmicPCAxisTagger::fDetWidth
private

Definition at line 66 of file CosmicPCAxisTagger_module.cc.

lar_cluster3d::PrincipalComponentsAlg cosmic::CosmicPCAxisTagger::fPcaAlg
private

Principal Components algorithm.

Definition at line 62 of file CosmicPCAxisTagger_module.cc.

std::string cosmic::CosmicPCAxisTagger::fPCAxisModuleLabel
private

Definition at line 60 of file CosmicPCAxisTagger_module.cc.

std::string cosmic::CosmicPCAxisTagger::fPFParticleModuleLabel
private

Definition at line 59 of file CosmicPCAxisTagger_module.cc.

float cosmic::CosmicPCAxisTagger::fTPCXBoundary
private

Definition at line 65 of file CosmicPCAxisTagger_module.cc.

float cosmic::CosmicPCAxisTagger::fTPCYBoundary
private

Definition at line 65 of file CosmicPCAxisTagger_module.cc.

float cosmic::CosmicPCAxisTagger::fTPCZBoundary
private

Definition at line 65 of file CosmicPCAxisTagger_module.cc.


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