24 #include <Eigen/Dense> 87 InitialiseProduct<art::Assns<recob::Shower, recob::PCAxis>>(
"ShowerPCAxisAssn");
88 InitialiseProduct<art::Assns<recob::PFParticle, recob::PCAxis>>(
"PFParticlePCAxisAssn");
103 auto const spHandle = Event.getValidHandle<std::vector<recob::SpacePoint>>(
fPFParticleLabel);
109 std::vector<art::Ptr<recob::SpacePoint>> spacePoints_pfp = fmspp.at(pfparticle.
key());
112 if (spacePoints_pfp.size() < 3) {
114 << spacePoints_pfp.size() <<
" spacepoints in shower, not calculating direction" 119 auto const clockData =
139 <<
"fUseStartPosition is set but ShowerStartPosition is not set. Bailing" << std::endl;
147 auto const GeneralDir = (ShowerCentre - StartPositionVec).Unit();
150 double DotProduct = PCADirection.Dot(GeneralDir);
153 if (DotProduct < 0) { PCADirection *= -1; }
163 spacePoints_pfp, ShowerCentre, PCADirection,
fNSegments);
167 if (RMSGradient < -std::numeric_limits<double>::epsilon()) { PCADirection *= -1; }
185 float TotalCharge = 0;
186 float sumWeights = 0;
197 clockData, detProp, sps, fmh, TotalCharge);
204 for (
auto& sp : sps) {
209 auto const sp_position = sp->position() - ShowerCentre;
223 wht *= std::sqrt(Charge / TotalCharge);
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;
236 Eigen::Matrix3f matrix;
239 matrix <<
xx, xy, xz, xy, yy, yz, xz, yz,
zz;
242 matrix /= sumWeights;
245 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigenMatrix(matrix);
247 Eigen::Vector3f eigenValuesVector = eigenMatrix.eigenvalues();
248 Eigen::Matrix3f eigenVectorsMatrix = eigenMatrix.eigenvectors();
251 const bool svdOk =
true;
252 const int nHits = sps.size();
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()};
262 return recob::PCAxis(svdOk, nHits, eigenValues, eigenVectors, avePos);
269 return {EigenVector[0], EigenVector[1], EigenVector[2]};
288 GetProducedElementPtr<recob::Shower>(
"shower", ShowerEleHolder);
290 AddSingle<art::Assns<recob::Shower, recob::PCAxis>>(showerPtr, pcaPtr,
"ShowerPCAxisAssn");
291 AddSingle<art::Assns<recob::PFParticle, recob::PCAxis>>(pfpPtr, pcaPtr,
"PFParticlePCAxisAssn");
const EigenVectors & getEigenVectors() const
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Declaration of signal hit object.
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
double ElectronLifetime() const
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.
key_type key() const noexcept
bool CheckElement(const std::string &Name) const
int GetElement(const std::string &Name, T &Element) const
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.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
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
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
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)