LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ShowerTrackPCADirection_tool.cc
Go to the documentation of this file.
1 //############################################################################
2 //### Name: ShowerTrackPCADirection ###
3 //### Author: Dominic Barker ###
4 //### Date: 13.05.19 ###
5 //### Description: Tool for finding the shower direction using PCA ###
6 //### methods using the initial track hits ###
7 //############################################################################
8 
9 //Framework Includes
11 
12 //LArSoft Includes
19 
20 //Root Includes
21 #include "TPrincipal.h"
22 
23 namespace ShowerRecoTools {
24 
26 
27  public:
29 
30  //Generic Direction Finder
31  int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
32  art::Event& Event,
33  reco::shower::ShowerElementHolder& ShowerEleHolder) override;
34 
35  private:
37  const detinfo::DetectorPropertiesData& detProp,
38  std::vector<art::Ptr<recob::SpacePoint>>& spacePoints_pfp,
39  const art::FindManyP<recob::Hit>& fmh,
40  geo::Point_t& ShowerCentre);
41 
42  //fcl
45  int fVerbose;
46  bool fChargeWeighted; //Should we charge weight the PCA.
47  unsigned int fMinPCAPoints; //Number of spacepoints needed to do the analysis.
48 
52  };
53 
55  : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
56  , fPFParticleLabel(pset.get<art::InputTag>("PFParticleLabel"))
57  , fHitModuleLabel(pset.get<art::InputTag>("HitModuleLabel"))
58  , fVerbose(pset.get<int>("Verbose"))
59  , fChargeWeighted(pset.get<bool>("ChargeWeighted"))
60  , fMinPCAPoints(pset.get<unsigned int>("MinPCAPoints"))
61  , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
62  , fInitialTrackSpacePointsInputLabel(pset.get<std::string>("InitialTrackSpacePointsInputLabel"))
63  , fShowerDirectionOutputLabel(pset.get<std::string>("ShowerDirectionOutputLabel"))
64  {}
65 
67  art::Event& Event,
68  reco::shower::ShowerElementHolder& ShowerEleHolder)
69  {
70 
71  if (!ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
72  if (fVerbose)
73  mf::LogError("ShowerTrackPCA") << "Start Position not set. Stopping" << std::endl;
74  ;
75  return 1;
76  }
77  if (!ShowerEleHolder.CheckElement(fInitialTrackSpacePointsInputLabel)) {
78  if (fVerbose)
79  mf::LogError("ShowerTrackPCA") << "TrackSpacePoints not set, returning " << std::endl;
80  return 1;
81  }
82 
83  //Get the spacepoints handle and the hit assoication
84  auto const spHandle = Event.getValidHandle<std::vector<recob::SpacePoint>>(fPFParticleLabel);
85 
86  const art::FindManyP<recob::Hit>& fmh =
87  ShowerEleHolder.GetFindManyP<recob::Hit>(spHandle, Event, fPFParticleLabel);
88 
89  std::vector<art::Ptr<recob::SpacePoint>> trackSpacePoints;
90  ShowerEleHolder.GetElement(fInitialTrackSpacePointsInputLabel, trackSpacePoints);
91 
92  //We cannot progress with no spacepoints.
93  if (trackSpacePoints.size() < fMinPCAPoints) {
94  if (fVerbose)
95  mf::LogError("ShowerTrackPCA") << "Not enough spacepoints for PCA, returning " << std::endl;
96  return 1;
97  }
98 
99  auto const clockData =
101  auto const detProp =
103 
104  //Find the PCA vector
105  geo::Point_t trackCentre;
106  auto Eigenvector = ShowerPCAVector(clockData, detProp, trackSpacePoints, fmh, trackCentre);
107 
108  //Get the General direction as the vector between the start position and the centre
109  geo::Point_t StartPositionVec = {-999, -999, -999};
110  ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, StartPositionVec);
111  auto const GeneralDir = (trackCentre - StartPositionVec).Unit();
112 
113  //Dot product
114  double DotProduct = Eigenvector.Dot(GeneralDir);
115 
116  //If the dotproduct is negative the Direction needs Flipping
117  if (DotProduct < 0) { Eigenvector *= -1.; }
118 
119  geo::Vector_t EigenvectorErr = {-999, -999, -999};
120  ShowerEleHolder.SetElement(Eigenvector, EigenvectorErr, fShowerDirectionOutputLabel);
121 
122  return 0;
123  }
124 
125  //Function to calculate the shower direction using a charge weight 3D PCA calculation.
127  const detinfo::DetectorClocksData& clockData,
128  const detinfo::DetectorPropertiesData& detProp,
130  const art::FindManyP<recob::Hit>& fmh,
131  geo::Point_t& ShowerCentre)
132  {
133  //Initialise the the PCA.
134  TPrincipal pca(3, "");
135 
136  float TotalCharge = 0;
137 
138  //Get the Shower Centre
139  ShowerCentre =
140  IShowerTool::GetLArPandoraShowerAlg().ShowerCentre(clockData, detProp, sps, fmh, TotalCharge);
141 
142  //Normalise the spacepoints, charge weight and add to the PCA.
143  for (auto& sp : sps) {
144 
145  float wht = 1;
146 
147  //Normalise the spacepoint position.
148  auto const sp_position = sp->position() - ShowerCentre;
149 
150  if (fChargeWeighted) {
151 
152  //Get the charge.
153  float Charge = IShowerTool::GetLArPandoraShowerAlg().SpacePointCharge(sp, fmh);
154 
155  //Get the time of the spacepoint
156  float Time = IShowerTool::GetLArPandoraShowerAlg().SpacePointTime(sp, fmh);
157 
158  //Correct for the lifetime at the moment.
159  Charge *= std::exp((sampling_rate(clockData) * Time) / (detProp.ElectronLifetime() * 1e3));
160 
161  //Charge Weight
162  wht *= std::sqrt(Charge / TotalCharge);
163  }
164 
165  double sp_coord[3];
166  sp_coord[0] = sp_position.X() * wht;
167  sp_coord[1] = sp_position.Y() * wht;
168  sp_coord[2] = sp_position.Z() * wht;
169 
170  //Add to the PCA
171  pca.AddRow(sp_coord);
172  }
173 
174  //Evaluate the PCA
175  pca.MakePrincipals();
176 
177  //Get the Eigenvectors.
178  const TMatrixD* Eigenvectors = pca.GetEigenVectors();
179 
180  return {(*Eigenvectors)[0][0], (*Eigenvectors)[1][0], (*Eigenvectors)[2][0]};
181  }
182 }
183 
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
ShowerTrackPCADirection(const fhicl::ParameterSet &pset)
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
Declaration of signal hit object.
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
STL namespace.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double SpacePointTime(art::Ptr< recob::SpacePoint > const &sp, art::FindManyP< recob::Hit > const &fmh) const
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
parameter set interface
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
geo::Vector_t ShowerPCAVector(const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, std::vector< art::Ptr< recob::SpacePoint >> &spacePoints_pfp, const art::FindManyP< recob::Hit > &fmh, geo::Point_t &ShowerCentre)
Definition: MVAAlg.h:12
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.
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
const art::FindManyP< T1 > & GetFindManyP(const art::ValidHandle< std::vector< T2 >> &handle, const art::Event &evt, const art::InputTag &moduleTag)