LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
LArPandoraShowerCreation_module.cc
Go to the documentation of this file.
1 
10 
11 #include "fhiclcpp/ParameterSet.h"
12 
18 
20 
21 #include <memory>
22 
23 namespace lar_pandora {
24 
26  public:
27  explicit LArPandoraShowerCreation(fhicl::ParameterSet const& pset);
28 
33 
34  void produce(art::Event& evt) override;
35 
36  private:
44  recob::Shower BuildShower(const int id,
45  const lar_content::LArShowerPCA& larShowerPCA,
46  const pandora::CartesianVector& vertexPosition) const;
47 
53  recob::PCAxis BuildPCAxis(const lar_content::LArShowerPCA& larShowerPCA) const;
54 
55  std::string m_pfParticleLabel;
57  };
58 
60 
61 } // namespace lar_pandora
62 
63 //------------------------------------------------------------------------------------------------------------------------------------------
64 // implementation follows
65 
69 
71 
73 
75 
77 
83 
85 
87 
89 
90 #include <iostream>
91 
92 namespace lar_pandora {
93 
95  : EDProducer{pset}
96  , m_pfParticleLabel(pset.get<std::string>("PFParticleLabel"))
97  , m_useAllParticles(pset.get<bool>("UseAllParticles", false))
98  {
99  produces<std::vector<recob::Shower>>();
100  produces<std::vector<recob::PCAxis>>();
101  produces<art::Assns<recob::PFParticle, recob::Shower>>();
102  produces<art::Assns<recob::PFParticle, recob::PCAxis>>();
103  produces<art::Assns<recob::Shower, recob::Hit>>();
104  produces<art::Assns<recob::Shower, recob::PCAxis>>();
105  }
106 
107  //------------------------------------------------------------------------------------------------------------------------------------------
108 
110  {
111  std::unique_ptr<std::vector<recob::Shower>> outputShowers(new std::vector<recob::Shower>);
112  std::unique_ptr<std::vector<recob::PCAxis>> outputPCAxes(new std::vector<recob::PCAxis>);
113  std::unique_ptr<art::Assns<recob::PFParticle, recob::Shower>> outputParticlesToShowers(
115  std::unique_ptr<art::Assns<recob::PFParticle, recob::PCAxis>> outputParticlesToPCAxes(
117  std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> outputShowersToHits(
119  std::unique_ptr<art::Assns<recob::Shower, recob::PCAxis>> outputShowersToPCAxes(
121 
122  const art::PtrMaker<recob::Shower> makeShowerPtr(evt);
123  const art::PtrMaker<recob::PCAxis> makePCAxisPtr(evt);
124 
125  int showerCounter(0);
126 
127  // Organise inputs
128  PFParticleVector pfParticleVector, extraPfParticleVector;
129  PFParticlesToSpacePoints pfParticlesToSpacePoints;
130  PFParticlesToClusters pfParticlesToClusters;
132  evt, m_pfParticleLabel, pfParticleVector, pfParticlesToSpacePoints);
134  evt, m_pfParticleLabel, extraPfParticleVector, pfParticlesToClusters);
135 
136  VertexVector vertexVector;
137  PFParticlesToVertices pfParticlesToVertices;
138  LArPandoraHelper::CollectVertices(evt, m_pfParticleLabel, vertexVector, pfParticlesToVertices);
139 
140  for (const art::Ptr<recob::PFParticle> pPFParticle : pfParticleVector) {
141  // Select shower-like pfparticles
142  if (!m_useAllParticles && !LArPandoraHelper::IsShower(pPFParticle)) continue;
143 
144  // Obtain associated spacepoints
145  PFParticlesToSpacePoints::const_iterator particleToSpacePointIter(
146  pfParticlesToSpacePoints.find(pPFParticle));
147 
148  if (pfParticlesToSpacePoints.end() == particleToSpacePointIter) {
149  mf::LogDebug("LArPandoraShowerCreation") << "No spacepoints associated to particle ";
150  continue;
151  }
152 
153  // Obtain associated clusters
154  PFParticlesToClusters::const_iterator particleToClustersIter(
155  pfParticlesToClusters.find(pPFParticle));
156 
157  if (pfParticlesToClusters.end() == particleToClustersIter) {
158  mf::LogDebug("LArPandoraShowerCreation") << "No clusters associated to particle ";
159  continue;
160  }
161 
162  // Obtain associated vertex
163  PFParticlesToVertices::const_iterator particleToVertexIter(
164  pfParticlesToVertices.find(pPFParticle));
165 
166  if ((pfParticlesToVertices.end() == particleToVertexIter) ||
167  (1 != particleToVertexIter->second.size())) {
168  mf::LogDebug("LArPandoraShowerCreation") << "Unexpected number of vertices for particle ";
169  continue;
170  }
171 
172  // Copy information into expected pandora form
173  pandora::CartesianPointVector cartesianPointVector;
174  for (const art::Ptr<recob::SpacePoint> spacePoint : particleToSpacePointIter->second)
175  cartesianPointVector.emplace_back(pandora::CartesianVector(
176  spacePoint->XYZ()[0], spacePoint->XYZ()[1], spacePoint->XYZ()[2]));
177 
178  double vertexXYZ[3] = {0., 0., 0.};
179  particleToVertexIter->second.front()->XYZ(vertexXYZ);
180  const pandora::CartesianVector vertexPosition(vertexXYZ[0], vertexXYZ[1], vertexXYZ[2]);
181 
182  // Call pandora "fast" shower fitter
183  try {
184  // Access centroid of shower via this method
185  const lar_content::LArShowerPCA initialLArShowerPCA(
186  lar_content::LArPfoHelper::GetPrincipalComponents(cartesianPointVector, vertexPosition));
187 
188  // Ensure successful creation of all structures before placing results in output containers, remaking LArShowerPCA with updated vertex
189  const pandora::CartesianVector& centroid(initialLArShowerPCA.GetCentroid());
190  const pandora::CartesianVector& primaryAxis(initialLArShowerPCA.GetPrimaryAxis());
191  const pandora::CartesianVector& secondaryAxis(initialLArShowerPCA.GetSecondaryAxis());
192  const pandora::CartesianVector& tertiaryAxis(initialLArShowerPCA.GetTertiaryAxis());
193  const pandora::CartesianVector& eigenvalues(initialLArShowerPCA.GetEigenValues());
194 
195  // Project the PFParticle vertex onto the PCA axis
196  const pandora::CartesianVector projectedVertexPosition(
197  centroid -
198  primaryAxis.GetUnitVector() * (centroid - vertexPosition).GetDotProduct(primaryAxis));
199 
200  // By convention, principal axis should always point away from vertex
201  const float testProjection(primaryAxis.GetDotProduct(projectedVertexPosition - centroid));
202  const float directionScaleFactor(
203  (testProjection > std::numeric_limits<float>::epsilon()) ? -1.f : 1.f);
204 
205  const lar_content::LArShowerPCA larShowerPCA(centroid,
206  primaryAxis * directionScaleFactor,
207  secondaryAxis * directionScaleFactor,
208  tertiaryAxis * directionScaleFactor,
209  eigenvalues);
211  showerCounter++, larShowerPCA, projectedVertexPosition));
212  const recob::PCAxis pcAxis(LArPandoraShowerCreation::BuildPCAxis(larShowerPCA));
213  outputShowers->emplace_back(shower);
214  outputPCAxes->emplace_back(pcAxis);
215  }
216  catch (const pandora::StatusCodeException&) {
217  mf::LogDebug("LArPandoraShowerCreation") << "Unable to extract shower pca";
218  continue;
219  }
220 
221  // Output objects
222  art::Ptr<recob::Shower> pShower(makeShowerPtr(outputShowers->size() - 1));
223  art::Ptr<recob::PCAxis> pPCAxis(makePCAxisPtr(outputPCAxes->size() - 1));
224 
225  HitVector hitsInParticle;
227  evt, m_pfParticleLabel, particleToClustersIter->second, hitsInParticle);
228 
229  // Output associations, after output objects are in place
230  util::CreateAssn(evt, pShower, pPFParticle, *(outputParticlesToShowers.get()));
231  util::CreateAssn(evt, pPCAxis, pPFParticle, *(outputParticlesToPCAxes.get()));
232  util::CreateAssn(evt, *(outputShowers.get()), hitsInParticle, *(outputShowersToHits.get()));
233  util::CreateAssn(evt, pPCAxis, pShower, *(outputShowersToPCAxes.get()));
234  }
235 
236  mf::LogDebug("LArPandora") << " Number of new showers: " << outputShowers->size()
237  << std::endl;
238  mf::LogDebug("LArPandora") << " Number of new pcaxes: " << outputPCAxes->size() << std::endl;
239 
240  evt.put(std::move(outputShowers));
241  evt.put(std::move(outputPCAxes));
242  evt.put(std::move(outputParticlesToShowers));
243  evt.put(std::move(outputParticlesToPCAxes));
244  evt.put(std::move(outputShowersToHits));
245  evt.put(std::move(outputShowersToPCAxes));
246  }
247 
248  //------------------------------------------------------------------------------------------------------------------------------------------
249 
251  const int id,
252  const lar_content::LArShowerPCA& larShowerPCA,
253  const pandora::CartesianVector& vertexPosition) const
254  {
255  const pandora::CartesianVector& showerLength(larShowerPCA.GetAxisLengths());
256  const pandora::CartesianVector& showerDirection(larShowerPCA.GetPrimaryAxis());
257 
258  const float length(showerLength.GetX());
259  const float openingAngle(
260  larShowerPCA.GetPrimaryLength() > 0.f ?
261  std::atan(larShowerPCA.GetSecondaryLength() / larShowerPCA.GetPrimaryLength()) :
262  0.f);
263  const TVector3 direction(
264  showerDirection.GetX(), showerDirection.GetY(), showerDirection.GetZ());
265  const TVector3 vertex(vertexPosition.GetX(), vertexPosition.GetY(), vertexPosition.GetZ());
266 
267  // TODO
268  const TVector3 directionErr;
269  const TVector3 vertexErr;
270  const std::vector<double> totalEnergyErr;
271  const std::vector<double> dEdx;
272  const std::vector<double> dEdxErr;
273  const std::vector<double> totalEnergy;
274  const int bestplane(0);
275 
276  return recob::Shower(direction,
277  directionErr,
278  vertex,
279  vertexErr,
280  totalEnergy,
281  totalEnergyErr,
282  dEdx,
283  dEdxErr,
284  bestplane,
285  id,
286  length,
287  openingAngle);
288  }
289 
290  //------------------------------------------------------------------------------------------------------------------------------------------
291 
293  const lar_content::LArShowerPCA& larShowerPCA) const
294  {
295  const pandora::CartesianVector& showerCentroid(larShowerPCA.GetCentroid());
296  const pandora::CartesianVector& showerDirection(larShowerPCA.GetPrimaryAxis());
297  const pandora::CartesianVector& showerSecondaryVector(larShowerPCA.GetSecondaryAxis());
298  const pandora::CartesianVector& showerTertiaryVector(larShowerPCA.GetTertiaryAxis());
299  const pandora::CartesianVector& showerEigenValues(larShowerPCA.GetEigenValues());
300 
301  const bool svdOK(true);
302  const double eigenValues[3] = {
303  showerEigenValues.GetX(),
304  showerEigenValues.GetY(),
305  showerEigenValues.GetZ()};
306  const double avePosition[3] = {showerCentroid.GetX(),
307  showerCentroid.GetY(),
308  showerCentroid.GetZ()};
309 
310  std::vector<std::vector<double>> eigenVecs = {
312  {showerDirection.GetX(), showerDirection.GetY(), showerDirection.GetZ()},
313  {showerSecondaryVector.GetX(), showerSecondaryVector.GetY(), showerSecondaryVector.GetZ()},
314  {showerTertiaryVector.GetX(), showerTertiaryVector.GetY(), showerTertiaryVector.GetZ()}};
315 
316  // TODO
317  const int numHitsUsed(100);
318  const double aveHitDoca(0.);
319  const size_t iD(util::kBogusI);
320 
321  return recob::PCAxis(svdOK, numHitsUsed, eigenValues, eigenVecs, avePosition, aveHitDoca, iD);
322  }
323 
324 } // namespace lar_pandora
static LArShowerPCA GetPrincipalComponents(const pandora::CartesianPointVector &pointVector, const pandora::CartesianVector &vertexPosition)
Perform PCA analysis on a set of 3D points and return results.
Header file for the pfo helper class.
recob::PCAxis BuildPCAxis(const lar_content::LArShowerPCA &larShowerPCA) const
Build a recob::PCAxis object.
LArPandoraShowerCreation & operator=(LArPandoraShowerCreation const &)=delete
const pandora::CartesianVector & GetCentroid() const
Return centroid.
std::map< art::Ptr< recob::PFParticle >, ClusterVector > PFParticlesToClusters
bool m_useAllParticles
Build a recob::Track for every recob::PFParticle.
float GetPrimaryLength() const
Return primary length.
Header file for lar pfo objects.
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
static bool IsShower(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as shower-like.
recob::Shower BuildShower(const int id, const lar_content::LArShowerPCA &larShowerPCA, const pandora::CartesianVector &vertexPosition) const
Build a recob::Shower object.
constexpr int kBogusI
obviously bogus integer value
intermediate_table::const_iterator const_iterator
std::map< art::Ptr< recob::PFParticle >, VertexVector > PFParticlesToVertices
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
TFile f
Definition: plotHisto.C:6
static void GetAssociatedHits(const art::Event &evt, const std::string &label, const std::vector< art::Ptr< T >> &inputVector, HitVector &associatedHits, const pandora::IntVector *const indexVector=nullptr)
Get all hits associated with input clusters.
std::string m_pfParticleLabel
The pf particle label.
LArPandoraShowerCreation(fhicl::ParameterSet const &pset)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
const pandora::CartesianVector & GetEigenValues() const
Return vector of eigenvalues.
std::vector< art::Ptr< recob::PFParticle > > PFParticleVector
static void CollectVertices(const art::Event &evt, const std::string &label, VertexVector &vertexVector, PFParticlesToVertices &particlesToVertices)
Collect the reconstructed PFParticles and associated Vertices from the ART event record.
const Double32_t * XYZ() const
Definition: SpacePoint.h:78
const pandora::CartesianVector & GetAxisLengths() const
Return vector of lengths.
std::map< art::Ptr< recob::PFParticle >, SpacePointVector > PFParticlesToSpacePoints
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2671
static void CollectPFParticles(const art::Event &evt, const std::string &label, PFParticleVector &particleVector)
Collect the reconstructed PFParticles from the ART event record.
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
std::vector< art::Ptr< recob::Hit > > HitVector
Utility object to perform functions of association.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
std::vector< art::Ptr< recob::Vertex > > VertexVector
const pandora::CartesianVector & GetSecondaryAxis() const
Return secondary axis.
const pandora::CartesianVector & GetTertiaryAxis() const
Return tertiary axis.
LArShowerPCA class.
TCEvent evt
Definition: DataStructs.cxx:8
float GetSecondaryLength() const
Return secondary length.
helper function for LArPandoraInterface producer module
art framework interface to geometry description
const pandora::CartesianVector & GetPrimaryAxis() const
Return primary axis.
vertex reconstruction