TrackerRecHit2D/src/VectorHit.cc

0001 #include "DataFormats/TrackerRecHit2D/interface/VectorHit.h"
0002 #include "Geometry/CommonDetUnit/interface/StackGeomDet.h"
0003 #include "CLHEP/Units/PhysicalConstants.h"
0004
0005 VectorHit::VectorHit(const GeomDet& idet,
0006                      const LocalPoint& posLower,
0007                      const LocalVector& dir,
0008                      const AlgebraicSymMatrix44& covMatrix,
0009                      const float chi2,
0010                      OmniClusterRef const& lower,
0011                      OmniClusterRef const& upper,
0012                      const float curvature,
0013                      const float curvatureError,
0014                      const float phi)
0015     : BaseTrackerRecHit(idet, trackerHitRTTI::vector),
0016       thePosition(posLower),
0017       theDirection(dir),
0018       theCovMatrix(covMatrix),
0019       theChi2(chi2),
0020       theLowerCluster(lower),
0021       theUpperCluster(upper),
0022       theCurvature(curvature),
0023       theCurvatureError(curvatureError),
0024       thePhi(phi) {}
0025
0026 VectorHit::VectorHit(const GeomDet& idet,
0027                      const VectorHit2D& vh2Dzx,
0028                      const VectorHit2D& vh2Dzy,
0029                      OmniClusterRef const& lower,
0030                      OmniClusterRef const& upper,
0031                      const float curvature,
0032                      const float curvatureError,
0033                      const float phi)
0034     : BaseTrackerRecHit(idet, trackerHitRTTI::vector),
0035       thePosition(LocalPoint(vh2Dzx.localPosition().x(), vh2Dzy.localPosition().x(), 0.)),
0036       theDirection(LocalVector(vh2Dzx.localDirection().x(), vh2Dzy.localDirection().x(), 1.)),
0037       theLowerCluster(lower),
0038       theUpperCluster(upper),
0039       theCurvature(curvature),
0040       theCurvatureError(curvatureError),
0041       thePhi(phi) {
0042   //building the cov matrix 4x4 starting from the 2x2
0043   const AlgebraicSymMatrix22& covMatZX = vh2Dzx.covMatrix();
0044   const AlgebraicSymMatrix22& covMatZY = vh2Dzy.covMatrix();
0045
0046   theCovMatrix = AlgebraicSymMatrix44();
0047
0048   theCovMatrix[0][0] = covMatZX[0][0];  // var(dx/dz)
0049   theCovMatrix[1][1] = covMatZY[0][0];  // var(dy/dz)
0050   theCovMatrix[2][2] = covMatZX[1][1];  // var(x)
0051   theCovMatrix[3][3] = covMatZY[1][1];  // var(y)
0052   theCovMatrix[0][2] = covMatZX[0][1];  // cov(dx/dz,x)
0053   theCovMatrix[1][3] = covMatZY[0][1];  // cov(dy/dz,y)
0054
0055   theChi2 = vh2Dzx.chi2() + vh2Dzy.chi2();
0056 }
0057
0058 bool VectorHit::sharesInput(const TrackingRecHit* other, SharedInputType what) const {
0059   if (what == all && (geographicalId() != other->geographicalId()))
0060     return false;
0061
0062   if (!sameDetModule(*other))
0063     return false;
0064
0065   if (trackerHitRTTI::isVector(*other)) {
0066     const VectorHit* otherVh = static_cast<const VectorHit*>(other);
0067     return sharesClusters(*otherVh, what);
0068   }
0069
0070   if (what == all)
0071     return false;
0072
0073   // what about multi???
0074   auto const& otherClus = reinterpret_cast<const BaseTrackerRecHit*>(other)->firstClusterRef();
0075
0076   if (lowerClusterRef().id() == otherClus.id() || upperClusterRef().id() == otherClus.id())
0077     return (otherClus == lowerClusterRef()) || (otherClus == upperClusterRef());
0078   else {
0079     bool lowerOverlap = false;
0080     bool upperOverlap = false;
0081     if (geographicalId() == other->geographicalId()) {
0082       lowerOverlap = otherClus.stripOverlap(lowerClusterRef());
0083       upperOverlap = otherClus.stripOverlap(upperClusterRef());
0084     }
0085     return (lowerOverlap || upperOverlap);
0086   }
0087 }
0088
0089 bool VectorHit::sharesClusters(VectorHit const& other, SharedInputType what) const {
0090   if (this->lowerClusterRef().id() == other.lowerClusterRef().id() ||
0091       this->upperClusterRef().id() == other.upperClusterRef().id()) {
0092     const bool lowerIdentity = this->lowerClusterRef() == other.lowerClusterRef();
0093     const bool upperIdentity = this->upperClusterRef() == other.upperClusterRef();
0094     return (what == TrackingRecHit::all) ? (lowerIdentity && upperIdentity) : (lowerIdentity || upperIdentity);
0095   } else {
0096     const bool sameDetId = (geographicalId() == other.geographicalId());
0097     bool lowerOverlap = (sameDetId) ? other.lowerClusterRef().stripOverlap(this->lowerClusterRef()) : false;
0098     bool upperOverlap = (sameDetId) ? other.upperClusterRef().stripOverlap(this->upperClusterRef()) : false;
0099     return (what == TrackingRecHit::all) ? (lowerOverlap && upperOverlap) : (lowerOverlap || upperOverlap);
0100   }
0101 }
0102
0103 void VectorHit::getKfComponents4D(KfComponentsHolder& holder) const {
0104   AlgebraicVector4& pars = holder.params<theDimension>();
0105   pars[0] = theDirection.x();
0106   pars[1] = theDirection.y();
0107   pars[2] = thePosition.x();
0108   pars[3] = thePosition.y();
0109
0110   holder.errors<theDimension>() = theCovMatrix;
0111
0112   ProjectMatrix<double, 5, theDimension>& pf = holder.projFunc<theDimension>();
0113   for (int i = 0; i < 4; ++i)
0114     pf.index[i] = i + 1;
0115
0116   holder.measuredParams<theDimension>() = AlgebraicVector4(&holder.tsosLocalParameters().At(1), theDimension);
0117   holder.measuredErrors<theDimension>() = holder.tsosLocalErrors().Sub<AlgebraicSymMatrix44>(1, 1);
0118 }
0119
0120 Global3DPoint VectorHit::lowerGlobalPos() const {
0121   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0122   const PixelGeomDetUnit* geomDetLower = static_cast<const PixelGeomDetUnit*>(stackDet->lowerDet());
0123   return phase2clusterGlobalPos(geomDetLower, lowerCluster());
0124 }
0125
0126 Global3DPoint VectorHit::upperGlobalPos() const {
0127   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0128   const PixelGeomDetUnit* geomDetUpper = static_cast<const PixelGeomDetUnit*>(stackDet->upperDet());
0129   return phase2clusterGlobalPos(geomDetUpper, upperCluster());
0130 }
0131
0132 Global3DPoint VectorHit::phase2clusterGlobalPos(const PixelGeomDetUnit* geomDet, ClusterRef cluster) {
0133   const PixelTopology* topo = &geomDet->specificTopology();
0134   float ix = cluster->center();
0135   float iy = cluster->column() + 0.5f;                   // halfway the column
0136   LocalPoint lp(topo->localX(ix), topo->localY(iy), 0);  // x, y, z
0137   Global3DPoint gp = geomDet->surface().toGlobal(lp);
0138   return gp;
0139 }
0140
0141 GlobalError VectorHit::lowerGlobalPosErr() const {
0142   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0143   const PixelGeomDetUnit* geomDetLower = static_cast<const PixelGeomDetUnit*>(stackDet->lowerDet());
0144   return phase2clusterGlobalPosErr(geomDetLower);
0145 }
0146
0147 GlobalError VectorHit::upperGlobalPosErr() const {
0148   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0149   const PixelGeomDetUnit* geomDetUpper = static_cast<const PixelGeomDetUnit*>(stackDet->upperDet());
0150   return phase2clusterGlobalPosErr(geomDetUpper);
0151 }
0152
0153 GlobalError VectorHit::phase2clusterGlobalPosErr(const PixelGeomDetUnit* geomDet) {
0154   const PixelTopology* topo = &geomDet->specificTopology();
0155   float pitchX = topo->pitch().first;
0156   float pitchY = topo->pitch().second;
0157   constexpr float invTwelve = 1. / 12;
0158   LocalError le(pow(pitchX, 2) * invTwelve, 0, pow(pitchY, 2) * invTwelve);  // e2_xx, e2_xy, e2_yy
0159   GlobalError ge(ErrorFrameTransformer().transform(le, geomDet->surface()));
0160   return ge;
0161 }
0162
0163 Global3DVector VectorHit::globalDirectionVH() const {
0164   Local3DVector theLocalDelta =
0165       LocalVector(theDirection.x() * theDirection.z(), theDirection.y() * theDirection.z(), theDirection.z());
0166   Global3DVector g = det()->surface().toGlobal(theLocalDelta);
0167   return g;
0168 }
0169
0170 Global3DVector VectorHit::globalDirection() const { return (det()->surface().toGlobal(localDirection())); }
0171
0172 float VectorHit::theta() const { return globalDirection().theta(); }
0173
0174 float VectorHit::transverseMomentum(float magField) const {
0175   return magField * (CLHEP::c_light * 1e-5F) / theCurvature;
0176   // pT [GeV] ~ 0.3 * B[T] * R [m], curvature is in cms, using precise value from speed of light
0177   // because curvature is a signed quantity this transverse momentum is also signed
0178 }
0179 float VectorHit::momentum(float magField) const { return transverseMomentum(magField) / (1. * sin(theta())); }
0180
0181 LocalError VectorHit::localPositionError() const {
0182   return LocalError(theCovMatrix[2][2], theCovMatrix[2][3], theCovMatrix[3][3]);
0183 }
0184
0185 LocalError VectorHit::localDirectionError() const {
0186   return LocalError(theCovMatrix[0][0], theCovMatrix[0][1], theCovMatrix[1][1]);
0187 }
0188
0189 const AlgebraicSymMatrix44& VectorHit::covMatrix() const { return theCovMatrix; }
0190
0191 std::ostream& operator<<(std::ostream& os, const VectorHit& vh) {
0192   os << " VectorHit create in the DetId#: " << vh.geographicalId() << "\n"
0193      << " Vectorhit local position      : " << vh.localPosition() << "\n"
0194      << " Vectorhit local direction     : " << vh.localDirection() << "\n"
0195      << " Vectorhit global direction    : " << vh.globalDirection() << "\n"
0196      << " Lower cluster global position : " << vh.lowerGlobalPos() << "\n"
0197      << " Upper cluster global position : " << vh.upperGlobalPos();
0198
0199   return os;
0200 }
0201
0202 /// Access to component RecHits (if any)
0203 std::vector<const TrackingRecHit*> VectorHit::recHits() const { return {}; }
0204
0205 /// Non-const access to component RecHits (if any)
0206 std::vector<TrackingRecHit*> VectorHit::recHits() { return {}; }