CommonTopologies/src/GEMStripTopology.cc

0001 /**GEMStripTopology
0002  * based on CSCRadialStripTopology and TrapezoidalStripTopology
0003  * \author Hyunyong Kim - TAMU
0004  */
0005 #include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include <iostream>
0008 #include <cmath>
0009 #include <algorithm>
0010
0011 GEMStripTopology::GEMStripTopology(int ns, float aw, float dh, float r0)
0012     : numberOfStrips_(ns), angularWidth_(aw), detHeight_(dh), centreToIntersection_(r0) {
0013   assert(angularWidth_ != 0);
0014   assert(detHeight_ != 0);
0015   yAxisOrientation_ = 1;
0016   phiOfOneEdge_ = -(0.5 * numberOfStrips_) * angularWidth_ * yAxisOrientation_;
0017   yCentre_ = 0;
0018   LogTrace("GEMStripTopology") << "Constructing GEMStripTopology with"
0019                                << " nstrips = " << ns << " angular width = " << aw << " det. height = " << dh
0020                                << " r0 = " << r0 << "\n";
0021 }
0022
0023 GEMStripTopology::GEMStripTopology(int ns, float aw, float dh, float r0, float yAx)
0024     : numberOfStrips_(ns), angularWidth_(aw), detHeight_(dh), centreToIntersection_(r0), yAxisOrientation_(yAx) {
0025   assert(angularWidth_ != 0);
0026   assert(detHeight_ != 0);
0027   phiOfOneEdge_ = -(0.5 * numberOfStrips_) * angularWidth_ * yAxisOrientation_;
0028   yCentre_ = 0;
0029   LogTrace("GEMStripTopology") << "Constructing GEMStripTopology with"
0030                                << " nstrips = " << ns << " angular width = " << aw << " det. height = " << dh
0031                                << " r0 = " << r0 << " yAxOrientation = " << yAx << "\n";
0032 }
0033
0034 LocalPoint GEMStripTopology::localPosition(float strip) const {
0035   return LocalPoint(yAxisOrientation() * originToIntersection() * tan(stripAngle(strip)), 0);
0036 }
0037
0038 LocalPoint GEMStripTopology::localPosition(const MeasurementPoint& mp) const {
0039   const float  // y = (L/cos(phi))*mp.y()*cos(phi)
0040       y(mp.y() * detHeight() + yCentreOfStripPlane()),
0041       x(yAxisOrientation() * yDistanceToIntersection(y) * std::tan(stripAngle(mp.x())));
0042   return LocalPoint(x, y);
0043 }
0044
0045 LocalError GEMStripTopology::localError(float strip, float stripErr2) const {
0046   const double phi(stripAngle(strip)), t1(std::tan(phi)), t2(t1 * t1),
0047       tt(stripErr2 * std::pow(centreToIntersection() * angularWidth(), 2)),  // tangential sigma^2   *c2
0048       rr(std::pow(detHeight(), 2) * (1.f / 12.f)),  // radial sigma^2( uniform prob density along strip)  *c2
0049
0050       xx(tt + t2 * rr), yy(t2 * tt + rr), xy(t1 * (rr - tt));
0051
0052   return LocalError(xx, xy, yy);
0053 }
0054
0055 LocalError GEMStripTopology::localError(const MeasurementPoint& mp, const MeasurementError& me) const {
0056   const double phi(stripAngle(mp.x())), s1(std::sin(phi)), c1(std::cos(phi));
0057   assert(c1 != 0);
0058   const double cs(s1 * c1), s2(s1 * s1),
0059       c2(1 - s2),  // rotation matrix
0060
0061       T(angularWidth() * (centreToIntersection() + yAxisOrientation() * mp.y() * detHeight()) /
0062         c1),                // tangential measurement unit (local pitch)
0063       R(detHeight() / c1),  // radial measurement unit (strip length)
0064       tt(me.uu() * T * T),  // tangential sigma^2
0065       rr(me.vv() * R * R),  // radial sigma^2
0066       tr(me.uv() * T * R),
0067
0068       xx(c2 * tt + 2 * cs * tr + s2 * rr), yy(s2 * tt - 2 * cs * tr + c2 * rr), xy(cs * (rr - tt) + tr * (c2 - s2));
0069
0070   return LocalError(xx, xy, yy);
0071 }
0072
0073 float GEMStripTopology::strip(const LocalPoint& lp) const {
0074   const float  // phi is measured from y axis --> sign of angle is sign of x * yAxisOrientation --> use atan2(x,y), not atan2(y,x)
0075       phi(std::atan2(lp.x(), yDistanceToIntersection(lp.y()))),
0076       aStrip((phi - yAxisOrientation() * phiOfOneEdge()) / angularWidth());
0077   return std::max(float(0), std::min((float)nstrips(), aStrip));
0078 }
0079
0080 int GEMStripTopology::nearestStrip(const LocalPoint& lp) const {
0081   return std::min(nstrips(), static_cast<int>(std::max(float(0), strip(lp))) + 1);
0082 }
0083
0084 MeasurementPoint GEMStripTopology::measurementPosition(const LocalPoint& lp) const {
0085   const float  // phi is [pi/2 - conventional local phi], use atan2(x,y) rather than atan2(y,x)
0086       phi(yAxisOrientation() * std::atan2(lp.x(), yDistanceToIntersection(lp.y())));
0087   return MeasurementPoint(yAxisOrientation() * (phi - phiOfOneEdge()) / angularWidth(),
0088                           (lp.y() - yCentreOfStripPlane()) / detHeight());
0089 }
0090
0091 MeasurementError GEMStripTopology::measurementError(const LocalPoint& p, const LocalError& e) const {
0092   const double yHitToInter(yDistanceToIntersection(p.y()));
0093   assert(yHitToInter != 0);
0094   const double t(yAxisOrientation() * p.x() / yHitToInter),  // tan(strip angle)
0095       cs(t / (1 + t * t)), s2(t * cs), c2(1 - s2),           // rotation matrix
0096
0097       T2(1. / (std::pow(angularWidth(), 2) *
0098                (std::pow(p.x(), 2) + std::pow(yHitToInter, 2)))),  // 1./tangential measurement unit (local pitch) ^2
0099       R2(c2 / std::pow(detHeight(), 2)),                           // 1./ radial measurement unit (strip length) ^2
0100
0101       uu((c2 * e.xx() - 2 * cs * e.xy() + s2 * e.yy()) * T2), vv((s2 * e.xx() + 2 * cs * e.xy() + c2 * e.yy()) * R2),
0102       uv((cs * (e.xx() - e.yy()) + e.xy() * (c2 - s2)) * std::sqrt(T2 * R2));
0103
0104   return MeasurementError(uu, uv, vv);
0105 }
0106
0107 int GEMStripTopology::channel(const LocalPoint& lp) const { return std::min(int(strip(lp)), numberOfStrips_ - 1); }
0108
0109 float GEMStripTopology::pitch() const { return localPitch(LocalPoint(0, 0)); }
0110
0111 float GEMStripTopology::localPitch(const LocalPoint& lp) const {
0112   const int istrip = std::min(nstrips(), static_cast<int>(strip(lp)) + 1);  // which strip number
0113   const float fangle = stripAngle(static_cast<float>(istrip) - 0.5);        // angle of strip centre
0114   assert(std::cos(fangle - 0.5f * angularWidth()) != 0);
0115   return yDistanceToIntersection(lp.y()) * std::sin(angularWidth()) /
0116          std::pow(std::cos(fangle - 0.5f * angularWidth()), 2);
0117 }
0118
0119 float GEMStripTopology::stripAngle(float strip) const {
0120   return phiOfOneEdge() + yAxisOrientation() * strip * angularWidth();
0121 }
0122
0123 float GEMStripTopology::localStripLength(const LocalPoint& lp) const {
0124   assert(yDistanceToIntersection(lp.y()) != 0);
0125   return detHeight() * std::sqrt(1.f + std::pow(lp.x() / yDistanceToIntersection(lp.y()), 2));
0126 }
0127
0128 float GEMStripTopology::yDistanceToIntersection(float y) const {
0129   return yAxisOrientation() * y + originToIntersection();
0130 }
0131
0132 float GEMStripTopology::xOfStrip(int strip, float y) const {
0133   return yAxisOrientation() * yDistanceToIntersection(y) * std::tan(stripAngle(static_cast<float>(strip) - 0.5));
0134 }