LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ShowerPCADirection_tool.cc
Go to the documentation of this file.
1 //############################################################################
2 //### Name: ShowerPCADirection ###
3 //### Author: Dominic Barker and Ed Tyley ###
4 //### Date: 13.05.19 ###
5 //### Description: Tool for finding the shower direction using PCA ###
6 //### methods. Derived from LArPandoraModularShowers Method. ###
7 //############################################################################
8 
9 //Framework Includes
11 
12 //LArSoft Includes
22 
23 //C++ Includes
24 #include <Eigen/Dense>
25 
26 namespace ShowerRecoTools {
27 
29 
30  public:
32 
33  //Calculate the direction of the shower.
34  int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
35  art::Event& Event,
36  reco::shower::ShowerElementHolder& ShowerEleHolder) override;
37 
38  private:
39  void InitialiseProducers() override;
40 
41  //Function to add the assoctions
43  art::Event& Event,
44  reco::shower::ShowerElementHolder& ShowerEleHolder) override;
45 
46  // Define standard art tool interface
48  const detinfo::DetectorClocksData& clockData,
49  const detinfo::DetectorPropertiesData& detProp,
50  const std::vector<art::Ptr<recob::SpacePoint>>& spacePoints_pfp,
51  const art::FindManyP<recob::Hit>& fmh,
52  geo::Point_t& ShowerCentre);
53 
55 
56  //fcl
58  int fVerbose;
59  unsigned int
60  fNSegments; //Used in the RMS gradient. How many segments should we split the shower into.
61  bool fUseStartPosition; //If we use the start position the drection of the
62  //PCA vector is decided as (Shower Centre - Shower Start Position).
63  bool fChargeWeighted; //Should the PCA axis be charge weighted.
64 
68  std::string fShowerPCAOutputLabel;
69  };
70 
72  : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
73  , fPFParticleLabel(pset.get<art::InputTag>("PFParticleLabel"))
74  , fVerbose(pset.get<int>("Verbose"))
75  , fNSegments(pset.get<unsigned int>("NSegments"))
76  , fUseStartPosition(pset.get<bool>("UseStartPosition"))
77  , fChargeWeighted(pset.get<bool>("ChargeWeighted"))
78  , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
79  , fShowerDirectionOutputLabel(pset.get<std::string>("ShowerDirectionOutputLabel"))
80  , fShowerCentreOutputLabel(pset.get<std::string>("ShowerCentreOutputLabel"))
81  , fShowerPCAOutputLabel(pset.get<std::string>("ShowerPCAOutputLabel"))
82  {}
83 
85  {
86  InitialiseProduct<std::vector<recob::PCAxis>>(fShowerPCAOutputLabel);
87  InitialiseProduct<art::Assns<recob::Shower, recob::PCAxis>>("ShowerPCAxisAssn");
88  InitialiseProduct<art::Assns<recob::PFParticle, recob::PCAxis>>("PFParticlePCAxisAssn");
89  }
90 
92  art::Event& Event,
93  reco::shower::ShowerElementHolder& ShowerEleHolder)
94  {
95 
96  // Get the assocated pfParicle vertex PFParticles
97  auto const pfpHandle = Event.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleLabel);
98 
100  ShowerEleHolder.GetFindManyP<recob::SpacePoint>(pfpHandle, Event, fPFParticleLabel);
101 
102  //Get the spacepoints handle and the hit assoication
103  auto const spHandle = Event.getValidHandle<std::vector<recob::SpacePoint>>(fPFParticleLabel);
104 
105  const art::FindManyP<recob::Hit>& fmh =
106  ShowerEleHolder.GetFindManyP<recob::Hit>(spHandle, Event, fPFParticleLabel);
107 
108  //SpacePoints
109  std::vector<art::Ptr<recob::SpacePoint>> spacePoints_pfp = fmspp.at(pfparticle.key());
110 
111  //We cannot progress with no spacepoints.
112  if (spacePoints_pfp.size() < 3) {
113  mf::LogWarning("ShowerPCADirection")
114  << spacePoints_pfp.size() << " spacepoints in shower, not calculating direction"
115  << std::endl;
116  return 1;
117  }
118 
119  auto const clockData =
121  auto const detProp =
123 
124  //Find the PCA vector
125  geo::Point_t ShowerCentre;
126  recob::PCAxis PCA = CalculateShowerPCA(clockData, detProp, spacePoints_pfp, fmh, ShowerCentre);
127  auto PCADirection = GetPCAxisVector(PCA);
128 
129  //Save the shower the center for downstream tools
130  geo::Point_t ShowerCentreErr = {-999, -999, -999};
131  ShowerEleHolder.SetElement(ShowerCentre, ShowerCentreErr, fShowerCentreOutputLabel);
132  ShowerEleHolder.SetElement(PCA, fShowerPCAOutputLabel);
133 
134  //Check if we are pointing the correct direction or not, First try the start position
135  if (fUseStartPosition) {
136  if (!ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
137  if (fVerbose)
138  mf::LogError("ShowerPCADirection")
139  << "fUseStartPosition is set but ShowerStartPosition is not set. Bailing" << std::endl;
140  return 1;
141  }
142  //Get the General direction as the vector between the start position and the centre
143  geo::Point_t StartPositionVec = {-999, -999, -999};
144  ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, StartPositionVec);
145 
146  // Calculate the general direction of the shower
147  auto const GeneralDir = (ShowerCentre - StartPositionVec).Unit();
148 
149  //Calculate the dot product between eigenvector and general direction
150  double DotProduct = PCADirection.Dot(GeneralDir);
151 
152  //If the dotproduct is negative the Direction needs Flipping
153  if (DotProduct < 0) { PCADirection *= -1; }
154 
155  //To do
156  geo::Vector_t PCADirectionErr = {-999, -999, -999};
157  ShowerEleHolder.SetElement(PCADirection, PCADirectionErr, fShowerDirectionOutputLabel);
158  return 0;
159  }
160 
161  //Otherwise Check against the RMS of the shower. Method adapated from EMShower Thanks Mike.
163  spacePoints_pfp, ShowerCentre, PCADirection, fNSegments);
164 
165  // If the alg fails to calculate the gradient it will return 0. In this case do nothing
166  // If the gradient is negative, flip the direction of the shower
167  if (RMSGradient < -std::numeric_limits<double>::epsilon()) { PCADirection *= -1; }
168 
169  //To do
170  geo::Vector_t PCADirectionErr = {-999, -999, -999};
171 
172  ShowerEleHolder.SetElement(PCADirection, PCADirectionErr, fShowerDirectionOutputLabel);
173  return 0;
174  }
175 
176  //Function to calculate the shower direction using a charge weight 3D PCA calculation.
178  const detinfo::DetectorClocksData& clockData,
179  const detinfo::DetectorPropertiesData& detProp,
181  const art::FindManyP<recob::Hit>& fmh,
182  geo::Point_t& ShowerCentre)
183  {
184 
185  float TotalCharge = 0;
186  float sumWeights = 0;
187  float xx = 0;
188  float yy = 0;
189  float zz = 0;
190  float xy = 0;
191  float xz = 0;
192  float yz = 0;
193 
194  //Get the Shower Centre
195  if (fChargeWeighted) {
197  clockData, detProp, sps, fmh, TotalCharge);
198  }
199  else {
201  }
202 
203  //Normalise the spacepoints, charge weight and add to the PCA.
204  for (auto& sp : sps) {
205 
206  float wht = 1;
207 
208  //Normalise the spacepoint position.
209  auto const sp_position = sp->position() - ShowerCentre;
210 
211  if (fChargeWeighted) {
212 
213  //Get the charge.
214  float Charge = IShowerTool::GetLArPandoraShowerAlg().SpacePointCharge(sp, fmh);
215 
216  //Get the time of the spacepoint
217  float Time = IShowerTool::GetLArPandoraShowerAlg().SpacePointTime(sp, fmh);
218 
219  //Correct for the lifetime at the moment.
220  Charge *= std::exp((sampling_rate(clockData) * Time) / (detProp.ElectronLifetime() * 1e3));
221 
222  //Charge Weight
223  wht *= std::sqrt(Charge / TotalCharge);
224  }
225 
226  xx += sp_position.X() * sp_position.X() * wht;
227  yy += sp_position.Y() * sp_position.Y() * wht;
228  zz += sp_position.Z() * sp_position.Z() * wht;
229  xy += sp_position.X() * sp_position.Y() * wht;
230  xz += sp_position.X() * sp_position.Z() * wht;
231  yz += sp_position.Y() * sp_position.Z() * wht;
232  sumWeights += wht;
233  }
234 
235  // Using Eigen package
236  Eigen::Matrix3f matrix;
237 
238  // Construct covariance matrix
239  matrix << xx, xy, xz, xy, yy, yz, xz, yz, zz;
240 
241  // Normalise from the sum of weights
242  matrix /= sumWeights;
243 
244  // Run the PCA
245  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigenMatrix(matrix);
246 
247  Eigen::Vector3f eigenValuesVector = eigenMatrix.eigenvalues();
248  Eigen::Matrix3f eigenVectorsMatrix = eigenMatrix.eigenvectors();
249 
250  // Put in the required form for a recob::PCAxis
251  const bool svdOk = true; //TODO: Should probably think about this a bit more
252  const int nHits = sps.size();
253  // For some reason eigen sorts the eigenvalues from smallest to largest, reverse it
254  const double eigenValues[3] = {
255  eigenValuesVector(2), eigenValuesVector(1), eigenValuesVector(0)};
256  std::vector<std::vector<double>> eigenVectors = {
257  {eigenVectorsMatrix(0, 2), eigenVectorsMatrix(1, 2), eigenVectorsMatrix(2, 2)},
258  {eigenVectorsMatrix(0, 1), eigenVectorsMatrix(1, 1), eigenVectorsMatrix(2, 1)},
259  {eigenVectorsMatrix(0, 0), eigenVectorsMatrix(1, 0), eigenVectorsMatrix(2, 0)}};
260  const double avePos[3] = {ShowerCentre.X(), ShowerCentre.Y(), ShowerCentre.Z()};
261 
262  return recob::PCAxis(svdOk, nHits, eigenValues, eigenVectors, avePos);
263  }
264 
266  {
267  //Get the Eigenvectors.
268  std::vector<double> EigenVector = PCAxis.getEigenVectors()[0];
269  return {EigenVector[0], EigenVector[1], EigenVector[2]};
270  }
271 
273  art::Event& Event,
274  reco::shower::ShowerElementHolder& ShowerEleHolder)
275  {
276 
277  //First check the element has been set
278  if (!ShowerEleHolder.CheckElement(fShowerPCAOutputLabel)) {
279  if (fVerbose) mf::LogError("ShowerPCADirection: Add Assns") << "PCA not set." << std::endl;
280  return 1;
281  }
282 
284 
285  const art::Ptr<recob::PCAxis> pcaPtr =
286  GetProducedElementPtr<recob::PCAxis>(fShowerPCAOutputLabel, ShowerEleHolder, ptrSize - 1);
287  const art::Ptr<recob::Shower> showerPtr =
288  GetProducedElementPtr<recob::Shower>("shower", ShowerEleHolder);
289 
290  AddSingle<art::Assns<recob::Shower, recob::PCAxis>>(showerPtr, pcaPtr, "ShowerPCAxisAssn");
291  AddSingle<art::Assns<recob::PFParticle, recob::PCAxis>>(pfpPtr, pcaPtr, "PFParticlePCAxisAssn");
292 
293  return 0;
294  }
295 }
const EigenVectors & getEigenVectors() const
Definition: PCAxis.h:74
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
Double_t xx
Definition: macro.C:12
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
geo::Vector_t GetPCAxisVector(recob::PCAxis &PCAxis)
int AddAssociations(const art::Ptr< recob::PFParticle > &pfpPtr, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
Declaration of signal hit object.
int GetVectorPtrSize(std::string Name)
Definition: IShowerTool.h:158
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
STL namespace.
Double_t zz
Definition: plot.C:276
recob::PCAxis CalculateShowerPCA(const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::SpacePoint >> &spacePoints_pfp, const art::FindManyP< recob::Hit > &fmh, geo::Point_t &ShowerCentre)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double SpacePointTime(art::Ptr< recob::SpacePoint > const &sp, art::FindManyP< recob::Hit > const &fmh) const
ShowerPCADirection(const fhicl::ParameterSet &pset)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
parameter set interface
key_type key() const noexcept
Definition: Ptr.h:166
bool CheckElement(const std::string &Name) const
int GetElement(const std::string &Name, T &Element) const
const shower::LArPandoraShowerAlg & GetLArPandoraShowerAlg() const
Definition: IShowerTool.h:82
double SpacePointCharge(art::Ptr< recob::SpacePoint > const &sp, art::FindManyP< recob::Hit > const &fmh) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
Contains all timing reference information for the detector.
geo::Point_t ShowerCentre(std::vector< art::Ptr< recob::SpacePoint >> const &showersps) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MVAAlg.h:12
double RMSShowerGradient(std::vector< art::Ptr< recob::SpacePoint >> &sps, const geo::Point_t &ShowerCentre, const geo::Vector_t &Direction, const unsigned int nSegments) const
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
const art::FindManyP< T1 > & GetFindManyP(const art::ValidHandle< std::vector< T2 >> &handle, const art::Event &evt, const art::InputTag &moduleTag)