30 #include "Eigen/Dense" 31 #include "Eigen/Eigenvalues" 32 #include "Eigen/Geometry" 33 #include "Eigen/Jacobi" 70 float doca3DScl)
const 102 int maxIterations(3);
106 maxRange = sclFctr * 0.5 * (maxRange + aveDoca);
109 int totalRejects(numRejHits);
110 int maxRejects(0.4 * hitPairVector.size());
113 while (maxIterations-- && numRejHits > 0 && totalRejects < maxRejects) {
138 bool skeletonOnly)
const 149 Eigen::Vector3d meanPos(Eigen::Vector3d::Zero());
150 double meanWeightSum(0.);
154 double minimumDeltaPeakSig(0.00001);
160 std::vector<double> hitChiSquareVec;
162 hitChiSquareVec.resize(hitPairVector.size());
164 std::transform(hitPairVector.begin(),
166 hitChiSquareVec.begin(),
167 [](
const auto&
hit) {
return hit->getHitChiSquare(); });
168 std::sort(hitChiSquareVec.begin(), hitChiSquareVec.end());
170 size_t numToKeep = 0.8 * hitChiSquareVec.size();
172 hitChiSquareVec.resize(numToKeep);
174 double aveValue = std::accumulate(hitChiSquareVec.begin(), hitChiSquareVec.end(), double(0.)) /
175 double(hitChiSquareVec.size());
176 double rms = std::sqrt(std::inner_product(hitChiSquareVec.begin(),
177 hitChiSquareVec.end(),
178 hitChiSquareVec.begin(),
181 [aveValue](
const auto&
left,
const auto&
right) {
182 return (
left - aveValue) * (
right - aveValue);
184 double(hitChiSquareVec.size()));
186 minimumDeltaPeakSig = std::max(minimumDeltaPeakSig, aveValue - rms);
190 for (
const auto&
hit : hitPairVector) {
196 double weight = std::max(minimumDeltaPeakSig,
197 double(
hit->getHitChiSquare()));
200 meanPos(0) +=
hit->getPosition()[0] *
weight;
201 meanPos(1) +=
hit->getPosition()[1] *
weight;
202 meanPos(2) +=
hit->getPosition()[2] *
weight;
208 meanPos /= meanWeightSum;
217 double weightSum(0.);
220 for (
const auto&
hit : hitPairVector) {
225 double weight = 1. / std::max(minimumDeltaPeakSig,
226 double(
hit->getHitChiSquare()));
228 double x = (
hit->getPosition()[0] - meanPos(0)) * weight;
229 double y = (
hit->getPosition()[1] - meanPos(1)) * weight;
230 double z = (
hit->getPosition()[2] - meanPos(2)) * weight;
232 weightSum += weight *
weight;
245 sig << xi2, xiyi, xizi, xiyi, yi2, yizi, xizi, yizi, zi2;
247 sig *= 1. / weightSum;
249 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigenMat(sig);
251 if (eigenMat.info() == Eigen::ComputationInfo::Success) {
257 eigenMat.eigenvectors().transpose().cast<
float>();
260 if (std::isnan(recobEigenVals[0])) {
261 std::cout <<
"==> Third eigenvalue returns a nan" << std::endl;
263 recobEigenVals[0] = 0.;
266 recobEigenVecs.row(0) = recobEigenVecs.row(1).cross(recobEigenVecs.row(2));
271 true, numPairsInt, recobEigenVals, recobEigenVecs, meanPos.cast<
float>());
274 mf::LogDebug(
"Cluster3D") <<
"PCA decompose failure, numPairs = " << numPairsInt << std::endl;
284 bool updateAvePos)
const 298 float aveHitDoca(0.);
299 Eigen::Vector3f avePosUpdate(0., 0., 0.);
310 std::vector<Eigen::Vector3f> hitPosVec;
313 for (
const auto& hit3D : hitPairVector) {
315 for (
const auto&
hit : hit3D->getHits()) {
323 auto const wirePosArr = wire_geom.
GetCenter();
325 Eigen::Vector3f wireCenter(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
326 Eigen::Vector3f wireDirVec(
330 float hitPeak(
hit->getHit()->PeakTime());
332 Eigen::Vector3f wirePos(
338 Eigen::Vector3f xAxis(1., 0., 0.);
339 Eigen::Vector3f planeNormal =
340 xAxis.cross(wireDirVec);
342 float docaInPlane(wirePos[0] - avePosition[0]);
343 float arcLenToPlane(0.);
344 float cosAxisToPlaneNormal = axisDirVec.dot(planeNormal);
346 Eigen::Vector3f axisPlaneIntersection = wirePos;
347 Eigen::Vector3f hitPosTVec = wirePos;
349 if (fabs(cosAxisToPlaneNormal) > 0.) {
350 Eigen::Vector3f deltaPos = wirePos - avePosition;
352 arcLenToPlane = deltaPos.dot(planeNormal) / cosAxisToPlaneNormal;
353 axisPlaneIntersection = avePosition + arcLenToPlane * axisDirVec;
354 docaInPlane = wirePos[0] - axisPlaneIntersection[0];
356 Eigen::Vector3f axisToInter = axisPlaneIntersection - wirePos;
357 float arcLenToDoca = axisToInter.dot(wireDirVec);
359 hitPosTVec += arcLenToDoca * wireDirVec;
363 Eigen::Vector3f wVec = avePosition - wirePos;
366 float a(axisDirVec.dot(axisDirVec));
367 float b(axisDirVec.dot(wireDirVec));
368 float c(wireDirVec.dot(wireDirVec));
369 float d(axisDirVec.dot(wVec));
370 float e(wireDirVec.dot(wVec));
372 float den(a * c - b * b);
378 arcLen1 = (b *
e - c *
d) / den;
379 arcLen2 = (a *
e - b *
d) / den;
382 mf::LogDebug(
"Cluster3D") <<
"** Parallel lines, need a solution here" << std::endl;
393 Eigen::Vector3f hitPos = wirePos + arcLen2 * wireDirVec;
394 Eigen::Vector3f axisPos = avePosition + arcLen1 * axisDirVec;
395 float deltaX = hitPos(0) - axisPos(0);
396 float deltaY = hitPos(1) - axisPos(1);
397 float deltaZ = hitPos(2) - axisPos(2);
398 float doca2 = deltaX * deltaX + deltaY * deltaY + deltaZ * deltaZ;
399 float doca = sqrt(doca2);
403 aveHitDoca += fabs(docaInPlane);
418 hit->setDocaToAxis(fabs(docaInPlane));
419 hit->setArcLenToPoca(arcLenToPlane);
423 if (
hit->getStatusBits() & 0x80) {
continue; }
427 avePosUpdate += hitPosTVec;
429 hitPosVec.push_back(hitPosTVec);
436 avePosUpdate /= float(nHits);
439 aveHitDoca /= float(nHits);
441 if (updateAvePos) { avePosition = avePosUpdate; }
444 for (
auto& hitPos : hitPosVec) {
446 float x = hitPos[0] - avePosition[0];
447 float y = hitPos[1] - avePosition[1];
448 float z = hitPos[2] - avePosition[2];
463 sig << xi2, xiyi, xizi, xiyi, yi2, yizi, xizi, yizi, zi2;
465 sig *= 1. / (nHits - 1);
467 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigenMat(sig);
469 if (eigenMat.info() == Eigen::ComputationInfo::Success) {
475 eigenMat.eigenvectors().transpose().cast<
float>();
479 true, nHits, recobEigenVals, recobEigenVecs, avePosition, aveHitDoca);
482 mf::LogDebug(
"Cluster3D") <<
"PCA decompose failure, numPairs = " << nHits << std::endl;
497 Eigen::Vector3f avePosition(
505 for (
const auto* clusterHit3D : hitPairVector) {
507 clusterHit3D->clearStatusBits(0x80);
511 Eigen::Vector3f clusPos(clusterHit3D->getPosition()[0],
512 clusterHit3D->getPosition()[1],
513 clusterHit3D->getPosition()[2]);
516 Eigen::Vector3f clusToHitVec = clusPos - avePosition;
519 float arclenToPoca = clusToHitVec.dot(axisDirVec);
522 Eigen::Vector3f docaPos = avePosition + arclenToPoca * axisDirVec;
525 Eigen::Vector3f docaPosToClusPos = clusPos - docaPos;
526 float docaToAxis = docaPosToClusPos.norm();
528 aveDoca3D += docaToAxis;
531 clusterHit3D->setDocaToAxis(docaToAxis);
532 clusterHit3D->setArclenToPoca(arclenToPoca);
536 aveDoca3D /= float(hitPairVector.size());
551 Eigen::Vector3f avePosition(
559 float aveHitDoca(0.);
562 for (
const auto*
hit : hit2DListPtr) {
571 auto const wirePosArr = wire_geom.
GetCenter();
573 Eigen::Vector3f wireCenter(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
574 Eigen::Vector3f wireDirVec(
578 Eigen::Vector3f wirePos(
hit->getXPosition(), wireCenter[1], wireCenter[2]);
581 Eigen::Vector3f xAxis(1., 0., 0.);
582 Eigen::Vector3f planeNormal =
583 xAxis.cross(wireDirVec);
585 float arcLenToPlane(0.);
586 float docaInPlane(wirePos[0] - avePosition[0]);
587 float cosAxisToPlaneNormal = axisDirVec.dot(planeNormal);
589 Eigen::Vector3f axisPlaneIntersection = wirePos;
592 if (fabs(cosAxisToPlaneNormal) > 0.) {
593 Eigen::Vector3f deltaPos = wirePos - avePosition;
596 std::min(
float(deltaPos.dot(planeNormal) / cosAxisToPlaneNormal), maxArcLen);
597 axisPlaneIntersection = avePosition + arcLenToPlane * axisDirVec;
599 Eigen::Vector3f axisToInter = axisPlaneIntersection - wirePos;
600 float arcLenToDoca = axisToInter.dot(wireDirVec);
603 if (fabs(arcLenToDoca) > wire_geom.
HalfL()) arcLenToDoca = wire_geom.
HalfL();
608 Eigen::Vector3f docaVec = axisPlaneIntersection - (wirePos + arcLenToDoca * wireDirVec);
609 docaInPlane = docaVec.norm();
612 aveHitDoca += fabs(docaInPlane);
615 hit->setDocaToAxis(fabs(docaInPlane));
616 hit->setArcLenToPoca(arcLenToPlane);
620 aveHitDoca /= float(hit2DListPtr.size());
629 float maxDocaAllowed)
const 639 for (
const auto& hit3D : hitPairVector) {
641 for (
const auto&
hit : hit3D->getHits()) {
643 hit->clearStatusBits(0x80);
645 if (
hit->getDocaToAxis() > maxDocaAllowed) {
646 hit->setStatusBit(0x80);
658 float maxDocaAllowed)
const 668 Eigen::Vector3f avePosition(
673 for (
const auto* clusterHit3D : hitPairVector) {
675 clusterHit3D->clearStatusBits(0x80);
679 Eigen::Vector3f clusPos(clusterHit3D->getPosition()[0],
680 clusterHit3D->getPosition()[1],
681 clusterHit3D->getPosition()[2]);
684 Eigen::Vector3f clusToHitVec = clusPos - avePosition;
687 float arclenToPoca = clusToHitVec.dot(axisDirVec);
690 Eigen::Vector3f docaPos = avePosition + arclenToPoca * axisDirVec;
693 Eigen::Vector3f docaPosToClusPos = clusPos - docaPos;
694 float docaToAxis = docaPosToClusPos.norm();
697 clusterHit3D->setDocaToAxis(docaToAxis);
698 clusterHit3D->setArclenToPoca(arclenToPoca);
701 if (clusterHit3D->getDocaToAxis() > maxDocaAllowed) {
702 clusterHit3D->setStatusBit(0x80);
float m_parallel
means lines are parallel
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
void setAveHitDoca(double doca) const
std::list< const reco::ClusterHit2D * > Hit2DListPtr
export some data structure definitions
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
int PCAAnalysis_reject3DOutliers(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca, float aveHitDoca) const
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Utilities related to art service access.
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
PrincipalComponentsAlg(fhicl::ParameterSet const &pset)
Constructor.
Declaration of signal hit object.
int PCAAnalysis_reject2DOutliers(const reco::HitPairListPtr &hitPairVector, float aveHitDoca) const
CryostatID_t Cryostat
Index of cryostat.
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
float getDocaToAxis() const
float getAveHitDoca() const
const EigenValues & getEigenValues() const
void PCAAnalysis_2D(const detinfo::DetectorPropertiesData &detProp, const reco::HitPairListPtr &hitPairVector, reco::PrincipalComponents &pca, bool updateAvePos=false) const
Vector_t Direction() const
Eigen::Vector3f EigenValues
T get(std::string const &key) const
const Eigen::Vector3f & getAvePosition() const
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
std::list< const reco::ClusterHit3D * > HitPairListPtr
PlaneID_t Plane
Index of the plane within its TPC.
Definition of data types for geometry description.
const geo::Geometry * m_geometry
Detector simulation of raw signals on wires.
This header file defines the interface to a principal components analysis designed to be used within ...
Encapsulate the geometry of a wire.
double ConvertTicksToX(double ticks, int p, int t, int c) const
double HalfL() const
Returns half the length of the wire [cm].
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
Eigen::Matrix3f EigenVectors
float getArclenToPoca() const
TPCID_t TPC
Index of the TPC within its cryostat.
void PCAAnalysis_calc2DDocas(const reco::Hit2DListPtr &hit2DVector, const reco::PrincipalComponents &pca) const
art framework interface to geometry description
const EigenVectors & getEigenVectors() const
void PCAAnalysis(const detinfo::DetectorPropertiesData &detProp, const reco::HitPairListPtr &hitPairVector, reco::PrincipalComponents &pca, float doca3DScl=3.) const
Run the Principal Components Analysis.