CSCGeometry/src/CSCWireGeometry.cc

0001 #include "Geometry/CSCGeometry/interface/CSCWireGeometry.h"
0002
0003 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0004
0005 #include <cmath>
0006
0007 LocalPoint CSCWireGeometry::intersection(float m1, float c1, float m2, float c2) const {
0008   // Calculate the point of intersection of two straight lines (in 2-dim)
0009   // BEWARE! Do not call with m1 = m2 ! No trapping !
0010
0011   float x = (c2 - c1) / (m1 - m2);
0012   float y = (m1 * c2 - m2 * c1) / (m1 - m2);
0013   return LocalPoint(x, y);
0014 }
0015
0016 std::vector<float> CSCWireGeometry::wireValues(float wire) const {
0017   // return x and y of mid-point of wire, and length of wire, as 3-dim vector.
0018   // If wire does not intersect active area the returned vector is filled with 0's.
0019
0020   std::pair<LocalPoint, LocalPoint> ends = wireEnds(wire);
0021
0022   std::vector<float> buf(3);  // note all elem init to 0
0023
0024   buf[0] = (ends.first.x() + ends.second.x()) / 2.;  // x is first elem of first & second pairs
0025   buf[1] = (ends.first.y() + ends.second.y()) / 2.;  // y is second elem of first & second pairs
0026   float d2 = (ends.first.x() - ends.second.x()) * (ends.first.x() - ends.second.x()) +
0027              (ends.first.y() - ends.second.y()) * (ends.first.y() - ends.second.y());
0028   buf[2] = sqrt(d2);
0029   return buf;
0030 }
0031
0032 std::pair<LocalPoint, LocalPoint> CSCWireGeometry::wireEnds(float wire) const {
0033   // return local (x, y) of each end of wire.
0034   // If wire does not intersect active area set all values to 0.
0035
0036   const float fprec = 1.E-06;
0037
0038   // slope of wire
0039   float wangle = wireAngle();
0040   float mw = 0;
0041   if (fabs(wangle) > fprec)
0042     mw = tan(wireAngle());
0043
0044   // intercept of wire
0045   float cw = yOfWire(wire);
0046
0047   LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: wire=" << wire << ", wire angle = " << wangle
0048                                   << ", intercept on y axis=" << cw;
0049
0050   // Find extent of wire plane
0051   double ww = wideWidthOfPlane();
0052   double nw = narrowWidthOfPlane();
0053   double len = lengthOfPlane();
0054
0055   // slope & intercept of line defining one non-parallel edge of wire-plane trapezoid
0056   float m1 = 2. * len / (ww - nw);
0057   float c1 = 0.;
0058   if (fabs(wangle) < fprec) {
0059     c1 = yOfFirstWire() - nw * len / (ww - nw);  // non-ME11
0060   } else {
0061     c1 = -len / 2. - nw * len / (ww - nw);  // ME11
0062   }
0063
0064   // slope & intercept of other non-parallel edge of wire-plane trapezoid
0065   float m2 = -m1;
0066   float c2 = c1;
0067
0068   // wire intersects edge 1 at
0069   LocalPoint pw1 = intersection(mw, cw, m1, c1);
0070   // wire intersects edge 2 at
0071   LocalPoint pw2 = intersection(mw, cw, m2, c2);
0072
0073   float x1 = pw1.x();
0074   float y1 = pw1.y();
0075
0076   float x2 = pw2.x();
0077   float y2 = pw2.y();
0078
0079   LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: wire intersects edges of plane at "
0080                                   << "\n  x1=" << x1 << " y1=" << y1 << " x2=" << x2 << " y2=" << y2;
0081
0082   // WIRES ARE NOT TILTED?
0083
0084   if (fabs(wangle) < fprec) {
0085     LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: wires are not tilted ";
0086     return std::pair<LocalPoint, LocalPoint>(LocalPoint(x1, y1), LocalPoint(x2, y2));
0087   }
0088
0089   // WIRES ARE TILTED
0090
0091   // ht and hb will be used to check where wire intersects edges of wire plane
0092   float ht = ww / 2.;
0093   float hb = nw / 2.;
0094   float mt = 0.;         // slope of top edge
0095   float mb = 0.;         //slope of bottom edge
0096   float cb = -len / 2.;  // intercept bottom edge of wire plane
0097   float ct = len / 2.;   // intercept top edge of wire plane
0098
0099   LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: slopes & intercepts "
0100                                   << "\n  mt=" << mt << " ct=" << ct << " mb=" << mb << " cb=" << cb << "\n  m1=" << m1
0101                                   << " c1=" << c1 << " m2=" << m2 << " c2=" << c2 << "\n  mw=" << mw << " cw=" << cw;
0102
0103   // wire intersects top edge at
0104   LocalPoint pwt = intersection(mw, cw, mt, ct);
0105   // wire intersects bottom edge at
0106   LocalPoint pwb = intersection(mw, cw, mb, cb);
0107
0108   // get the local coordinates
0109   float xt = pwt.x();
0110   float yt = pwt.y();
0111
0112   float xb = pwb.x();
0113   float yb = pwb.y();
0114
0115   LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: wire intersects top & bottom of wire plane at "
0116                                   << "\n  xt=" << xt << " yt=" << yt << " xb=" << xb << " yb=" << yb;
0117
0118   float xWireEnd[4], yWireEnd[4];
0119
0120   int i = 0;
0121   if (fabs(x1) >= hb && fabs(x1) <= ht) {
0122     // wire does intersect side edge 1 of wire plane
0123     xWireEnd[i] = x1;
0124     yWireEnd[i] = y1;
0125     i++;
0126   }
0127   if (fabs(xb) <= hb) {
0128     // wire does intersect bottom edge of wire plane
0129     xWireEnd[i] = xb;
0130     yWireEnd[i] = yb;
0131     i++;
0132   }
0133   if (fabs(x2) >= hb && fabs(x2) <= ht) {
0134     // wire does intersect side edge 2 of wire plane
0135     xWireEnd[i] = x2;
0136     yWireEnd[i] = y2;
0137     i++;
0138   }
0139   if (fabs(xt) <= ht) {
0140     // wire does intersect top edge of wire plane
0141     xWireEnd[i] = xt;
0142     yWireEnd[i] = yt;
0143     i++;
0144   }
0145
0146   if (i != 2) {
0147     // the wire does not intersect the wire plane (!)
0148
0149     LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: wire does not intersect wire plane!!";
0150     //     throw cms::Exception("BadCSCGeometry") << "the wire has " << i <<
0151     //       " ends!" << "\n";
0152
0153     return std::pair<LocalPoint, LocalPoint>(LocalPoint(0., 0.), LocalPoint(0., 0.));
0154   }
0155
0156   LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: ME11 wire ends ";
0157   for (int j = 0; j < i; j++) {
0158     LogTrace("CSCWireGeometry|CSC") << "  x = " << xWireEnd[j] << " y = " << yWireEnd[j];
0159   }
0160
0161   return std::pair<LocalPoint, LocalPoint>(LocalPoint(xWireEnd[0], yWireEnd[0]), LocalPoint(xWireEnd[1], yWireEnd[1]));
0162 }
0163
0164 //@@ COULD/SHOULD BE IMPLEMENTED IN Slanted & Nonslanted DERIVED CLASSES
0165
0166 std::pair<float, float> CSCWireGeometry::equationOfWire(float wire) const {
0167   const float fprec = 1.E-06;
0168
0169   // slope of wire
0170   float wangle = wireAngle();
0171   float mw = 0;
0172   if (fabs(wangle) > fprec)
0173     mw = tan(wangle);
0174
0175   // intercept of wire
0176   float cw = yOfWire(wire);
0177
0178   LogTrace("CSCWireGeometry|CSC") << "CSCWireGeometry: wire=" << wire << ", wire angle = " << wangle
0179                                   << ", intercept on y axis=" << cw;
0180
0181   return std::pair<float, float>(mw, cw);
0182 }
0183
0184 //@@ COULD/SHOULD BE IMPLEMENTED IN Slanted & Nonslanted DERIVED CLASSES
0185
0186 std::pair<float, float> CSCWireGeometry::yLimitsOfWirePlane() const {
0187   const float fprec = 0.1;      // wire angle is either 0 or 29 degrees = 0.506 rads
0188   float ylow = yOfFirstWire();  // non-ME11 chambers
0189   float wangle = wireAngle();
0190   if (fabs(wangle) > fprec) {
0191     ylow += tan(std::abs(wangle)) * narrowWidthOfPlane() / 2.;  // correction for ME11
0192   }
0193   float yhigh = ylow + lengthOfPlane();  // add extent of wire plane in y
0194
0195   return std::pair<float, float>(ylow, yhigh);
0196 }