LArSoft  v09_90_00
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 
17 
19 
20 #include <memory>
21 
22 namespace lar_pandora {
23 
25  public:
26  explicit LArPandoraShowerCreation(fhicl::ParameterSet const& pset);
27 
32 
33  void produce(art::Event& evt) override;
34 
35  private:
43  recob::Shower BuildShower(const int id,
44  const lar_content::LArShowerPCA& larShowerPCA,
45  const pandora::CartesianVector& vertexPosition) const;
46 
52  recob::PCAxis BuildPCAxis(const lar_content::LArShowerPCA& larShowerPCA) const;
53 
54  std::string m_pfParticleLabel;
56  };
57 
59 
60 } // namespace lar_pandora
61 
62 //------------------------------------------------------------------------------------------------------------------------------------------
63 // implementation follows
64 
68 
70 
72 
74 
76 
82 
84 
86 
88 
89 #include <iostream>
90 
91 namespace lar_pandora {
92 
94  : EDProducer{pset}
95  , m_pfParticleLabel(pset.get<std::string>("PFParticleLabel"))
96  , m_useAllParticles(pset.get<bool>("UseAllParticles", false))
97  {
98  produces<std::vector<recob::Shower>>();
99  produces<std::vector<recob::PCAxis>>();
100  produces<art::Assns<recob::PFParticle, recob::Shower>>();
101  produces<art::Assns<recob::PFParticle, recob::PCAxis>>();
102  produces<art::Assns<recob::Shower, recob::Hit>>();
103  produces<art::Assns<recob::Shower, recob::PCAxis>>();
104  }
105 
106  //------------------------------------------------------------------------------------------------------------------------------------------
107 
109  {
110  std::unique_ptr<std::vector<recob::Shower>> outputShowers(new std::vector<recob::Shower>);
111  std::unique_ptr<std::vector<recob::PCAxis>> outputPCAxes(new std::vector<recob::PCAxis>);
112  std::unique_ptr<art::Assns<recob::PFParticle, recob::Shower>> outputParticlesToShowers(
114  std::unique_ptr<art::Assns<recob::PFParticle, recob::PCAxis>> outputParticlesToPCAxes(
116  std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> outputShowersToHits(
118  std::unique_ptr<art::Assns<recob::Shower, recob::PCAxis>> outputShowersToPCAxes(
120 
121  const art::PtrMaker<recob::Shower> makeShowerPtr(evt);
122  const art::PtrMaker<recob::PCAxis> makePCAxisPtr(evt);
123 
124  int showerCounter(0);
125 
126  // Organise inputs
127  PFParticleVector pfParticleVector, extraPfParticleVector;
128  PFParticlesToSpacePoints pfParticlesToSpacePoints;
129  PFParticlesToClusters pfParticlesToClusters;
131  evt, m_pfParticleLabel, pfParticleVector, pfParticlesToSpacePoints);
133  evt, m_pfParticleLabel, extraPfParticleVector, pfParticlesToClusters);
134 
135  VertexVector vertexVector;
136  PFParticlesToVertices pfParticlesToVertices;
137  LArPandoraHelper::CollectVertices(evt, m_pfParticleLabel, vertexVector, pfParticlesToVertices);
138 
139  for (const art::Ptr<recob::PFParticle> pPFParticle : pfParticleVector) {
140  // Select shower-like pfparticles
141  if (!m_useAllParticles && !LArPandoraHelper::IsShower(pPFParticle)) continue;
142 
143  // Obtain associated spacepoints
144  PFParticlesToSpacePoints::const_iterator particleToSpacePointIter(
145  pfParticlesToSpacePoints.find(pPFParticle));
146 
147  if (pfParticlesToSpacePoints.end() == particleToSpacePointIter) {
148  mf::LogDebug("LArPandoraShowerCreation") << "No spacepoints associated to particle ";
149  continue;
150  }
151 
152  // Obtain associated clusters
153  PFParticlesToClusters::const_iterator particleToClustersIter(
154  pfParticlesToClusters.find(pPFParticle));
155 
156  if (pfParticlesToClusters.end() == particleToClustersIter) {
157  mf::LogDebug("LArPandoraShowerCreation") << "No clusters associated to particle ";
158  continue;
159  }
160 
161  // Obtain associated vertex
162  PFParticlesToVertices::const_iterator particleToVertexIter(
163  pfParticlesToVertices.find(pPFParticle));
164 
165  if ((pfParticlesToVertices.end() == particleToVertexIter) ||
166  (1 != particleToVertexIter->second.size())) {
167  mf::LogDebug("LArPandoraShowerCreation") << "Unexpected number of vertices for particle ";
168  continue;
169  }
170 
171  // Copy information into expected pandora form
172  pandora::CartesianPointVector cartesianPointVector;
173  for (const art::Ptr<recob::SpacePoint> spacePoint : particleToSpacePointIter->second)
174  cartesianPointVector.emplace_back(pandora::CartesianVector(
175  spacePoint->XYZ()[0], spacePoint->XYZ()[1], spacePoint->XYZ()[2]));
176 
177  double vertexXYZ[3] = {0., 0., 0.};
178  particleToVertexIter->second.front()->XYZ(vertexXYZ);
179  const pandora::CartesianVector vertexPosition(vertexXYZ[0], vertexXYZ[1], vertexXYZ[2]);
180 
181  // Call pandora "fast" shower fitter
182  try {
183  // Access centroid of shower via this method
184  const lar_content::LArShowerPCA initialLArShowerPCA(
185  lar_content::LArPfoHelper::GetPrincipalComponents(cartesianPointVector, vertexPosition));
186 
187  // Ensure successful creation of all structures before placing results in output containers, remaking LArShowerPCA with updated vertex
188  const pandora::CartesianVector& centroid(initialLArShowerPCA.GetCentroid());
189  const pandora::CartesianVector& primaryAxis(initialLArShowerPCA.GetPrimaryAxis());
190  const pandora::CartesianVector& secondaryAxis(initialLArShowerPCA.GetSecondaryAxis());
191  const pandora::CartesianVector& tertiaryAxis(initialLArShowerPCA.GetTertiaryAxis());
192  const pandora::CartesianVector& eigenvalues(initialLArShowerPCA.GetEigenValues());
193 
194  // Project the PFParticle vertex onto the PCA axis
195  const pandora::CartesianVector projectedVertexPosition(
196  centroid -
197  primaryAxis.GetUnitVector() * (centroid - vertexPosition).GetDotProduct(primaryAxis));
198 
199  // By convention, principal axis should always point away from vertex
200  const float testProjection(primaryAxis.GetDotProduct(projectedVertexPosition - centroid));
201  const float directionScaleFactor(
202  (testProjection > std::numeric_limits<float>::epsilon()) ? -1.f : 1.f);
203 
204  const lar_content::LArShowerPCA larShowerPCA(centroid,
205  primaryAxis * directionScaleFactor,
206  secondaryAxis * directionScaleFactor,
207  tertiaryAxis * directionScaleFactor,
208  eigenvalues);
210  showerCounter++, larShowerPCA, projectedVertexPosition));
211  const recob::PCAxis pcAxis(LArPandoraShowerCreation::BuildPCAxis(larShowerPCA));
212  outputShowers->emplace_back(shower);
213  outputPCAxes->emplace_back(pcAxis);
214  }
215  catch (const pandora::StatusCodeException&) {
216  mf::LogDebug("LArPandoraShowerCreation") << "Unable to extract shower pca";
217  continue;
218  }
219 
220  // Output objects
221  art::Ptr<recob::Shower> pShower(makeShowerPtr(outputShowers->size() - 1));
222  art::Ptr<recob::PCAxis> pPCAxis(makePCAxisPtr(outputPCAxes->size() - 1));
223 
224  HitVector hitsInParticle;
226  evt, m_pfParticleLabel, particleToClustersIter->second, hitsInParticle);
227 
228  // Output associations, after output objects are in place
229  util::CreateAssn(evt, pShower, pPFParticle, *(outputParticlesToShowers.get()));
230  util::CreateAssn(evt, pPCAxis, pPFParticle, *(outputParticlesToPCAxes.get()));
231  util::CreateAssn(evt, *(outputShowers.get()), hitsInParticle, *(outputShowersToHits.get()));
232  util::CreateAssn(evt, pPCAxis, pShower, *(outputShowersToPCAxes.get()));
233  }
234 
235  mf::LogDebug("LArPandora") << " Number of new showers: " << outputShowers->size()
236  << std::endl;
237  mf::LogDebug("LArPandora") << " Number of new pcaxes: " << outputPCAxes->size() << std::endl;
238 
239  evt.put(std::move(outputShowers));
240  evt.put(std::move(outputPCAxes));
241  evt.put(std::move(outputParticlesToShowers));
242  evt.put(std::move(outputParticlesToPCAxes));
243  evt.put(std::move(outputShowersToHits));
244  evt.put(std::move(outputShowersToPCAxes));
245  }
246 
247  //------------------------------------------------------------------------------------------------------------------------------------------
248 
250  const int id,
251  const lar_content::LArShowerPCA& larShowerPCA,
252  const pandora::CartesianVector& vertexPosition) const
253  {
254  const pandora::CartesianVector& showerLength(larShowerPCA.GetAxisLengths());
255  const pandora::CartesianVector& showerDirection(larShowerPCA.GetPrimaryAxis());
256 
257  const float length(showerLength.GetX());
258  const float openingAngle(
259  larShowerPCA.GetPrimaryLength() > 0.f ?
260  std::atan(larShowerPCA.GetSecondaryLength() / larShowerPCA.GetPrimaryLength()) :
261  0.f);
262  const TVector3 direction(
263  showerDirection.GetX(), showerDirection.GetY(), showerDirection.GetZ());
264  const TVector3 vertex(vertexPosition.GetX(), vertexPosition.GetY(), vertexPosition.GetZ());
265 
266  // TODO
267  const TVector3 directionErr;
268  const TVector3 vertexErr;
269  const std::vector<double> totalEnergyErr;
270  const std::vector<double> dEdx;
271  const std::vector<double> dEdxErr;
272  const std::vector<double> totalEnergy;
273  const int bestplane(0);
274 
275  return recob::Shower(direction,
276  directionErr,
277  vertex,
278  vertexErr,
279  totalEnergy,
280  totalEnergyErr,
281  dEdx,
282  dEdxErr,
283  bestplane,
284  id,
285  length,
286  openingAngle);
287  }
288 
289  //------------------------------------------------------------------------------------------------------------------------------------------
290 
292  const lar_content::LArShowerPCA& larShowerPCA) const
293  {
294  const pandora::CartesianVector& showerCentroid(larShowerPCA.GetCentroid());
295  const pandora::CartesianVector& showerDirection(larShowerPCA.GetPrimaryAxis());
296  const pandora::CartesianVector& showerSecondaryVector(larShowerPCA.GetSecondaryAxis());
297  const pandora::CartesianVector& showerTertiaryVector(larShowerPCA.GetTertiaryAxis());
298  const pandora::CartesianVector& showerEigenValues(larShowerPCA.GetEigenValues());
299 
300  const bool svdOK(true);
301  const double eigenValues[3] = {
302  showerEigenValues.GetX(),
303  showerEigenValues.GetY(),
304  showerEigenValues.GetZ()};
305  const double avePosition[3] = {showerCentroid.GetX(),
306  showerCentroid.GetY(),
307  showerCentroid.GetZ()};
308 
309  std::vector<std::vector<double>> eigenVecs = {
311  {showerDirection.GetX(), showerDirection.GetY(), showerDirection.GetZ()},
312  {showerSecondaryVector.GetX(), showerSecondaryVector.GetY(), showerSecondaryVector.GetZ()},
313  {showerTertiaryVector.GetX(), showerTertiaryVector.GetY(), showerTertiaryVector.GetZ()}};
314 
315  // TODO
316  const int numHitsUsed(100);
317  const double aveHitDoca(0.);
318  const size_t iD(util::kBogusI);
319 
320  return recob::PCAxis(svdOK, numHitsUsed, eigenValues, eigenVecs, avePosition, aveHitDoca, iD);
321  }
322 
323 } // 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:2675
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