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   return (otherClus == lowerClusterRef()) || (otherClus == upperClusterRef());
0076 }
0077
0078 bool VectorHit::sharesClusters(VectorHit const& other, SharedInputType what) const {
0079   bool lower = this->lowerClusterRef() == other.lowerClusterRef();
0080   bool upper = this->upperClusterRef() == other.upperClusterRef();
0081
0082   return (what == TrackingRecHit::all) ? (lower && upper) : (upper || lower);
0083 }
0084
0085 void VectorHit::getKfComponents4D(KfComponentsHolder& holder) const {
0086   AlgebraicVector4& pars = holder.params<theDimension>();
0087   pars[0] = theDirection.x();
0088   pars[1] = theDirection.y();
0089   pars[2] = thePosition.x();
0090   pars[3] = thePosition.y();
0091
0092   holder.errors<theDimension>() = theCovMatrix;
0093
0094   ProjectMatrix<double, 5, theDimension>& pf = holder.projFunc<theDimension>();
0095   for (int i = 0; i < 4; ++i)
0096     pf.index[i] = i + 1;
0097
0098   holder.measuredParams<theDimension>() = AlgebraicVector4(&holder.tsosLocalParameters().At(1), theDimension);
0099   holder.measuredErrors<theDimension>() = holder.tsosLocalErrors().Sub<AlgebraicSymMatrix44>(1, 1);
0100 }
0101
0102 Global3DPoint VectorHit::lowerGlobalPos() const {
0103   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0104   const PixelGeomDetUnit* geomDetLower = static_cast<const PixelGeomDetUnit*>(stackDet->lowerDet());
0105   return phase2clusterGlobalPos(geomDetLower, lowerCluster());
0106 }
0107
0108 Global3DPoint VectorHit::upperGlobalPos() const {
0109   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0110   const PixelGeomDetUnit* geomDetUpper = static_cast<const PixelGeomDetUnit*>(stackDet->upperDet());
0111   return phase2clusterGlobalPos(geomDetUpper, upperCluster());
0112 }
0113
0114 Global3DPoint VectorHit::phase2clusterGlobalPos(const PixelGeomDetUnit* geomDet, ClusterRef cluster) {
0115   const PixelTopology* topo = &geomDet->specificTopology();
0116   float ix = cluster->center();
0117   float iy = cluster->column() + 0.5f;                   // halfway the column
0118   LocalPoint lp(topo->localX(ix), topo->localY(iy), 0);  // x, y, z
0119   Global3DPoint gp = geomDet->surface().toGlobal(lp);
0120   return gp;
0121 }
0122
0123 GlobalError VectorHit::lowerGlobalPosErr() const {
0124   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0125   const PixelGeomDetUnit* geomDetLower = static_cast<const PixelGeomDetUnit*>(stackDet->lowerDet());
0126   return phase2clusterGlobalPosErr(geomDetLower);
0127 }
0128
0129 GlobalError VectorHit::upperGlobalPosErr() const {
0130   const StackGeomDet* stackDet = static_cast<const StackGeomDet*>(det());
0131   const PixelGeomDetUnit* geomDetUpper = static_cast<const PixelGeomDetUnit*>(stackDet->upperDet());
0132   return phase2clusterGlobalPosErr(geomDetUpper);
0133 }
0134
0135 GlobalError VectorHit::phase2clusterGlobalPosErr(const PixelGeomDetUnit* geomDet) {
0136   const PixelTopology* topo = &geomDet->specificTopology();
0137   float pitchX = topo->pitch().first;
0138   float pitchY = topo->pitch().second;
0139   constexpr float invTwelve = 1. / 12;
0140   LocalError le(pow(pitchX, 2) * invTwelve, 0, pow(pitchY, 2) * invTwelve);  // e2_xx, e2_xy, e2_yy
0141   GlobalError ge(ErrorFrameTransformer().transform(le, geomDet->surface()));
0142   return ge;
0143 }
0144
0145 Global3DVector VectorHit::globalDirectionVH() const {
0146   Local3DVector theLocalDelta =
0147       LocalVector(theDirection.x() * theDirection.z(), theDirection.y() * theDirection.z(), theDirection.z());
0148   Global3DVector g = det()->surface().toGlobal(theLocalDelta);
0149   return g;
0150 }
0151
0152 Global3DVector VectorHit::globalDirection() const { return (det()->surface().toGlobal(localDirection())); }
0153
0154 float VectorHit::theta() const { return globalDirection().theta(); }
0155
0156 float VectorHit::transverseMomentum(float magField) const {
0157   return magField * (CLHEP::c_light * 1e-5F) / theCurvature;
0158   // pT [GeV] ~ 0.3 * B[T] * R [m], curvature is in cms, using precise value from speed of light
0159   // because curvature is a signed quantity this transverse momentum is also signed
0160 }
0161 float VectorHit::momentum(float magField) const { return transverseMomentum(magField) / (1. * sin(theta())); }
0162
0163 LocalError VectorHit::localPositionError() const {
0164   return LocalError(theCovMatrix[2][2], theCovMatrix[2][3], theCovMatrix[3][3]);
0165 }
0166
0167 LocalError VectorHit::localDirectionError() const {
0168   return LocalError(theCovMatrix[0][0], theCovMatrix[0][1], theCovMatrix[1][1]);
0169 }
0170
0171 const AlgebraicSymMatrix44& VectorHit::covMatrix() const { return theCovMatrix; }
0172
0173 std::ostream& operator<<(std::ostream& os, const VectorHit& vh) {
0174   os << " VectorHit create in the DetId#: " << vh.geographicalId() << "\n"
0175      << " Vectorhit local position      : " << vh.localPosition() << "\n"
0176      << " Vectorhit local direction     : " << vh.localDirection() << "\n"
0177      << " Vectorhit global direction    : " << vh.globalDirection() << "\n"
0178      << " Lower cluster global position : " << vh.lowerGlobalPos() << "\n"
0179      << " Upper cluster global position : " << vh.upperGlobalPos();
0180
0181   return os;
0182 }
0183
0184 /// Access to component RecHits (if any)
0185 std::vector<const TrackingRecHit*> VectorHit::recHits() const { return {}; }
0186
0187 /// Non-const access to component RecHits (if any)
0188 std::vector<TrackingRecHit*> VectorHit::recHits() { return {}; }