DTRecHit/src/DTRecSegment4D.cc

0001 /** \file
0002  *
0003  * \author Stefano Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>
0004  * \author Riccardo Bellan - INFN TO <riccardo.bellan@cern.ch>
0005  */
0006
0007 /* This Class Header */
0008 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
0009
0010 /* Collaborating Class Header */
0011 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 /* C++ Headers */
0014
0015 DTRecSegment4D::DTRecSegment4D(const DTChamberRecSegment2D& phiSeg,
0016                                const DTSLRecSegment2D& zedSeg,
0017                                const LocalPoint& posZInCh,
0018                                const LocalVector& dirZInCh)
0019     : RecSegment(phiSeg.chamberId()), theProjection(full), thePhiSeg(phiSeg), theZedSeg(zedSeg), theDimension(4) {
0020   // Check consistency of 2 sub-segments
0021   if (DTChamberId(phiSeg.geographicalId().rawId()) != DTChamberId(zedSeg.geographicalId().rawId()))
0022     throw cms::Exception("DTRecSegment4D")
0023         << "the z Segment and the phi segment have different chamber id" << std::endl;
0024
0025   // The position of 2D segments are defined in the SL frame: I must first
0026   // extrapolate that position at the Chamber reference plane
0027   LocalPoint posZAt0 = posZInCh + dirZInCh * (-posZInCh.z()) / cos(dirZInCh.theta());
0028
0029   thePosition = LocalPoint(phiSeg.localPosition().x(), posZAt0.y(), 0.);
0030   LocalVector dirPhiInCh = phiSeg.localDirection();
0031
0032   // given the actual definition of chamber refFrame, (with z poiniting to IP),
0033   // the zed component of direction is negative.
0034   theDirection = LocalVector(dirPhiInCh.x() / fabs(dirPhiInCh.z()), dirZInCh.y() / fabs(dirZInCh.z()), -1.);
0035   theDirection = theDirection.unit();
0036
0037   // set cov matrix
0038   theCovMatrix = AlgebraicSymMatrix(4);
0039   theCovMatrix[0][0] = phiSeg.covMatrix()[0][0];  //sigma (dx/dz)
0040   theCovMatrix[0][2] = phiSeg.covMatrix()[0][1];  //cov(dx/dz,x)
0041   theCovMatrix[2][2] = phiSeg.covMatrix()[1][1];  //sigma (x)
0042   setCovMatrixForZed(posZInCh);
0043 }
0044
0045 DTRecSegment4D::DTRecSegment4D(const DTChamberRecSegment2D& phiSeg)
0046     : RecSegment(phiSeg.chamberId()),
0047       theProjection(phi),
0048       thePhiSeg(phiSeg),
0049       theZedSeg(DTSLRecSegment2D()),
0050       theDimension(2) {
0051   thePosition = thePhiSeg.localPosition();
0052
0053   theDirection = thePhiSeg.localDirection();
0054
0055   // set cov matrix
0056   theCovMatrix = AlgebraicSymMatrix(4);
0057   theCovMatrix[0][0] = phiSeg.covMatrix()[0][0];  //sigma (dx/dz)
0058   theCovMatrix[0][2] = phiSeg.covMatrix()[0][1];  //cov(dx/dz,x)
0059   theCovMatrix[2][2] = phiSeg.covMatrix()[1][1];  //sigma (x)
0060 }
0061
0062 DTRecSegment4D::DTRecSegment4D(const DTSLRecSegment2D& zedSeg, const LocalPoint& posZInCh, const LocalVector& dirZInCh)
0063     : RecSegment(zedSeg.superLayerId().chamberId()),
0064       theProjection(Z),
0065       thePhiSeg(DTChamberRecSegment2D()),
0066       theZedSeg(zedSeg),
0067       theDimension(2) {
0068   LocalPoint posZAt0 = posZInCh + dirZInCh * (-posZInCh.z() / cos(dirZInCh.theta()));
0069
0070   thePosition = posZAt0;
0071   theDirection = dirZInCh;
0072
0073   // set cov matrix
0074   theCovMatrix = AlgebraicSymMatrix(4);
0075   setCovMatrixForZed(posZInCh);
0076 }
0077
0078 DTRecSegment4D::~DTRecSegment4D() {}
0079
0080 AlgebraicVector DTRecSegment4D::parameters() const {
0081   if (dimension() == 4) {
0082     // (dx/dz,dy/dz,x,y)
0083     AlgebraicVector result(4);
0084     result[2] = thePosition.x();
0085     result[3] = thePosition.y();
0086     result[0] = theDirection.x() / theDirection.z();
0087     result[1] = theDirection.y() / theDirection.z();
0088     return result;
0089   }
0090
0091   AlgebraicVector result(2);
0092   if (theProjection == phi) {
0093     // (dx/dz,x)
0094     result[1] = thePosition.x();
0095     result[0] = theDirection.x() / theDirection.z();
0096   } else if (theProjection == Z) {
0097     // (dy/dz,y) (note we are in the chamber r.f.)
0098     result[1] = thePosition.y();
0099     result[0] = theDirection.y() / theDirection.z();
0100   }
0101   return result;
0102 }
0103
0104 AlgebraicSymMatrix DTRecSegment4D::parametersError() const {
0105   if (dimension() == 4) {
0106     return theCovMatrix;
0107   }
0108
0109   AlgebraicSymMatrix result(2);
0110   if (theProjection == phi) {
0111     result[0][0] = theCovMatrix[0][0];  //S(dx/dz)
0112     result[0][1] = theCovMatrix[0][2];  //Cov(dx/dz,x)
0113     result[1][1] = theCovMatrix[2][2];  //S(x)
0114   } else if (theProjection == Z) {
0115     result[0][0] = theCovMatrix[1][1];  //S(dy/dz)
0116     result[0][1] = theCovMatrix[1][3];  //Cov(dy/dz,y)
0117     result[1][1] = theCovMatrix[3][3];  //S(y)
0118   }
0119   return result;
0120 }
0121
0122 //These methods are only used to initialize the const static values
0123 // used by projectionMatrix().
0124 static AlgebraicMatrix initThe4DProjectionMatrix() {
0125   AlgebraicMatrix the4DProjectionMatrix(4, 5, 0);
0126   the4DProjectionMatrix[0][1] = 1;
0127   the4DProjectionMatrix[1][2] = 1;
0128   the4DProjectionMatrix[2][3] = 1;
0129   the4DProjectionMatrix[3][4] = 1;
0130   return the4DProjectionMatrix;
0131 }
0132 static const AlgebraicMatrix the4DProjectionMatrix{initThe4DProjectionMatrix()};
0133
0134 static AlgebraicMatrix initThe2DPhiProjMatrix() {
0135   AlgebraicMatrix the2DPhiProjMatrix(2, 5, 0);
0136   the2DPhiProjMatrix[0][1] = 1;
0137   the2DPhiProjMatrix[1][3] = 1;
0138   return the2DPhiProjMatrix;
0139 }
0140 static const AlgebraicMatrix the2DPhiProjMatrix{initThe2DPhiProjMatrix()};
0141
0142 static AlgebraicMatrix initThe2DZProjMatrix() {
0143   AlgebraicMatrix the2DZProjMatrix(2, 5, 0);
0144   the2DZProjMatrix[0][2] = 1;
0145   the2DZProjMatrix[1][4] = 1;
0146   return the2DZProjMatrix;
0147 }
0148 static const AlgebraicMatrix the2DZProjMatrix{initThe2DZProjMatrix()};
0149
0150 AlgebraicMatrix DTRecSegment4D::projectionMatrix() const {
0151   if (dimension() == 4) {
0152     return the4DProjectionMatrix;
0153   } else if (theProjection == phi) {
0154     return the2DPhiProjMatrix;
0155   } else if (theProjection == Z) {
0156     return the2DZProjMatrix;
0157   } else {
0158     return AlgebraicMatrix();
0159   }
0160 }
0161
0162 LocalError DTRecSegment4D::localPositionError() const {
0163   return LocalError(theCovMatrix[2][2], theCovMatrix[2][3], theCovMatrix[3][3]);
0164 }
0165
0166 LocalError DTRecSegment4D::localDirectionError() const {
0167   return LocalError(theCovMatrix[0][0], theCovMatrix[0][1], theCovMatrix[1][1]);
0168 }
0169
0170 double DTRecSegment4D::chi2() const {
0171   double result = 0;
0172   if (hasPhi())
0173     result += thePhiSeg.chi2();
0174   if (hasZed())
0175     result += theZedSeg.chi2();
0176   return result;
0177 }
0178
0179 int DTRecSegment4D::degreesOfFreedom() const {
0180   int result = 0;
0181   if (hasPhi())
0182     result += thePhiSeg.degreesOfFreedom();
0183   if (hasZed())
0184     result += theZedSeg.degreesOfFreedom();
0185   return result;
0186 }
0187
0188 void DTRecSegment4D::setCovMatrixForZed(const LocalPoint& posZInCh) {
0189   // Warning!!! the covariance matrix for Theta SL segment is defined in the SL
0190   // reference frame, here that in the Chamber ref frame must be used.
0191   // For direction, no problem, but the position is extrapolated, so we must
0192   // propagate the error properly.
0193
0194   // many thanks to Paolo Ronchese for the help in deriving the formulas!
0195
0196   // y=m*z+q in SL frame
0197   // y=m'*z+q' in CH frame
0198
0199   // var(m') = var(m)
0200   theCovMatrix[1][1] = theZedSeg.parametersError()[0][0];  //sigma (dy/dz)
0201
0202   // cov(m',q') = DeltaZ*Var(m) + Cov(m,q)
0203   theCovMatrix[1][3] =
0204       posZInCh.z() * theZedSeg.parametersError()[0][0] + theZedSeg.parametersError()[0][1];  //cov(dy/dz,y)
0205
0206   // Var(q') = DeltaZ^2*Var(m) + Var(q) + 2*DeltaZ*Cov(m,q)
0207   // cout << "Var(q') = DeltaZ^2*Var(m) + Var(q) + 2*DeltaZ*Cov(m,q)" << endl;
0208   // cout << "Var(q')= " << posZInCh.z()*posZInCh.z() << "*" <<
0209   //   theZedSeg.parametersError()[0][0] << " + " <<
0210   //   theZedSeg.parametersError()[1][1] << " + " <<
0211   //   2*posZInCh.z() << "*" << theZedSeg.parametersError()[0][1] ;
0212   theCovMatrix[3][3] = 2. * (posZInCh.z() * posZInCh.z()) * theZedSeg.parametersError()[0][0] +
0213                        theZedSeg.parametersError()[1][1] + 2. * posZInCh.z() * theZedSeg.parametersError()[0][1];
0214   // cout << " = " << theCovMatrix[3][3] << endl;
0215 }
0216
0217 std::ostream& operator<<(std::ostream& os, const DTRecSegment4D& seg) {
0218   os << "Pos " << seg.localPosition() << " Dir: " << seg.localDirection() << " dim: " << seg.dimension()
0219      << " chi2/ndof: " << seg.chi2() << "/" << seg.degreesOfFreedom() << " :";
0220   if (seg.hasPhi())
0221     os << seg.phiSegment()->recHits().size();
0222   else
0223     os << 0;
0224   os << ":";
0225   if (seg.hasZed())
0226     os << seg.zSegment()->recHits().size();
0227   else
0228     os << 0;
0229   return os;
0230 }
0231
0232 /// Access to component RecHits (if any)
0233 std::vector<const TrackingRecHit*> DTRecSegment4D::recHits() const {
0234   std::vector<const TrackingRecHit*> pointersOfRecHits;
0235
0236   if (hasPhi())
0237     pointersOfRecHits.push_back(phiSegment());
0238   if (hasZed())
0239     pointersOfRecHits.push_back(zSegment());
0240
0241   return pointersOfRecHits;
0242 }
0243
0244 /// Non-const access to component RecHits (if any)
0245 std::vector<TrackingRecHit*> DTRecSegment4D::recHits() {
0246   std::vector<TrackingRecHit*> pointersOfRecHits;
0247
0248   if (hasPhi())
0249     pointersOfRecHits.push_back(phiSegment());
0250   if (hasZed())
0251     pointersOfRecHits.push_back(zSegment());
0252
0253   return pointersOfRecHits;
0254 }
0255
0256 DTChamberId DTRecSegment4D::chamberId() const { return DTChamberId(geographicalId()); }