CommonTopologies/src/TkRadialStripTopology.cc

0001 #include <Geometry/CommonTopologies/interface/TkRadialStripTopology.h>
0002 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0003
0004 #include <iostream>
0005 #include <cmath>
0006 #include <algorithm>
0007 #include <cassert>
0008
0009 #include <vdt/vdtMath.h>
0010
0011 #ifdef MATH_STS
0012 #include <iostream>
0013 #endif
0014 namespace {
0015
0016 #ifdef MATH_STS
0017   struct Stat {
0018     Stat(const char* in) : name(in){};
0019     ~Stat() {
0020       edm::LogVerbatim("CommonTopologies")
0021           << name << ": atan0 calls tot/large/over1: " << natan << "/" << nlarge << "/" << over1;
0022     }
0023
0024     void add(float t) {
0025       auto at = std::abs(t);
0026       ++natan;
0027       if (at > 0.40f)
0028         ++nlarge;
0029       if (at > 1.0)
0030         ++over1;
0031     }
0032     const char* name;
0033     long long natan = 0;
0034     long long nlarge = 0;
0035     long long over1 = 0;
0036   };
0037
0038   Stat statM("mpos");
0039   Stat statS("span");
0040 #endif
0041
0042   // valid for |x| < 0.15  (better then 10^-9
0043   template <typename T>
0044   inline T tan15(T x) {
0045     return x * (T(1) + (x * x) * (T(0.33331906795501708984375) + (x * x) * T(0.135160386562347412109375)));
0046   }
0047
0048   // valid for z < pi/8
0049   //  x * (1 + x*x * (-0.33322894573211669921875 + x*x * (0.1967026889324188232421875 + x*x * (-0.11053790152072906494140625))))  // .1e-7 by Sollya
0050   inline float atan0(float t) {
0051     auto z = t;
0052     // if( t > 0.4142135623730950f ) // * tan pi/8
0053     // z = (t-1.0f)/(t+1.0f);
0054     float z2 = z * z;
0055     float ret =
0056         (((8.05374449538e-2f * z2 - 1.38776856032E-1f) * z2 + 1.99777106478E-1f) * z2 - 3.33329491539E-1f) * z2 * z + z;
0057     // if( t > 0.4142135623730950f ) ret +=0.7853981633974483096f;
0058     return ret;
0059   }
0060
0061   inline float atanClip(float t) {
0062     constexpr float tanPi8 = 0.4142135623730950;
0063     constexpr float pio8 = 3.141592653589793238 / 8;
0064     float at = std::abs(t);
0065     return std::copysign((at < tanPi8) ? atan0(at) : pio8, t);
0066   }
0067
0068 }  // namespace
0069
0070 TkRadialStripTopology::TkRadialStripTopology(int ns, float aw, float dh, float r, int yAx, float yMid)
0071     : theNumberOfStrips(ns),
0072       theAngularWidth(aw),
0073       theAWidthInverse(1.f / aw),
0074       theTanAW(std::tan(aw)),
0075       theDetHeight(dh),
0076       theCentreToIntersection(r),
0077       theYAxisOrientation(yAx),
0078       yCentre(yMid),
0079       theRadialSigma(std::pow(dh, 2.f) * (1.f / 12.f)) {
0080   // Angular offset of extreme edge of detector, so that angle is
0081   // zero for a strip lying along local y axis = long symmetry axis of plane of strips
0082   thePhiOfOneEdge = -(0.5 * theNumberOfStrips) * theAngularWidth;  // always negative!
0083   theTanOfOneEdge = std::tan(std::abs(thePhiOfOneEdge));
0084   assert(std::abs(thePhiOfOneEdge) < 0.15);  //
0085
0086   LogTrace("TkRadialStripTopology") << "TkRadialStripTopology: constructed with"
0087                                     << " strips = " << ns << " width = " << aw << " rad "
0088                                     << " det_height = " << dh << " ctoi = " << r << " phi_edge = " << thePhiOfOneEdge
0089                                     << " rad "
0090                                     << " y_ax_ori = " << theYAxisOrientation << " y_det_centre = " << yCentre << "\n";
0091 }
0092
0093 int TkRadialStripTopology::channel(const LocalPoint& lp) const {
0094   return std::min(int(strip(lp)), theNumberOfStrips - 1);
0095 }
0096
0097 int TkRadialStripTopology::nearestStrip(const LocalPoint& lp) const {
0098   return std::min(nstrips(), static_cast<int>(std::max(float(0), strip(lp))) + 1);
0099 }
0100
0101 float TkRadialStripTopology::yDistanceToIntersection(float y) const {
0102   return yAxisOrientation() * y + originToIntersection();
0103 }
0104
0105 float TkRadialStripTopology::localStripLength(const LocalPoint& lp) const {
0106   return detHeight() * std::sqrt(1.f + std::pow(lp.x() / yDistanceToIntersection(lp.y()), 2.f));
0107 }
0108
0109 float TkRadialStripTopology::xOfStrip(int strip, float y) const {
0110   return yAxisOrientation() * yDistanceToIntersection(y) * std::tan(stripAngle(static_cast<float>(strip) - 0.5f));
0111 }
0112
0113 float TkRadialStripTopology::strip(const LocalPoint& lp) const {
0114   // phi is measured from y axis --> sign of angle is sign of x * yAxisOrientation --> use atan2(x,y), not atan2(y,x)
0115   const float phi = atanClip(lp.x() / yDistanceToIntersection(lp.y()));
0116   const float aStrip = (phi - phiOfOneEdge()) * theAWidthInverse;
0117   return std::max(float(0), std::min((float)nstrips(), aStrip));
0118 }
0119
0120 float TkRadialStripTopology::coveredStrips(const LocalPoint& lp1, const LocalPoint& lp2) const {
0121   // http://en.wikipedia.org/wiki/List_of_trigonometric_identities#Angle_sum_and_difference_identities
0122   // atan(a)-atan(b) = atan( (a-b)/(1+a*b) )
0123   // avoid divisions
0124   // float t1 = lp1.x()/yDistanceToIntersection( lp1.y() );
0125   // float t2 = lp2.x()/yDistanceToIntersection( lp2.y() );
0126   // float t = (t1-t2)/(1.+t1*t2);
0127   auto y1 = yDistanceToIntersection(lp1.y());
0128   auto y2 = yDistanceToIntersection(lp2.y());
0129   auto x1 = lp1.x();
0130   auto x2 = lp2.x();
0131
0132   auto t = (y2 * x1 - y1 * x2) / (y1 * y2 + x1 * x2);
0133
0134 #ifdef MATH_STS
0135   statS.add(t);
0136 #endif
0137   //  edm::LogVerbatim("CommonTopologies") << "atans " << atanClip(t) << " " << std::atan2(lp1.x(), yDistanceToIntersection(lp1.y())) - std::atan2(lp2.x(),yDistanceToIntersection(lp2.y()));
0138   // clip???
0139   return atanClip(t) * theAWidthInverse;
0140   //   return (measurementPosition(lp1)-measurementPosition(lp2)).x();
0141 }
0142
0143 LocalPoint TkRadialStripTopology::localPosition(float strip) const {
0144   return LocalPoint(yAxisOrientation() * originToIntersection() * tan15(stripAngle(strip)), 0);
0145 }
0146
0147 LocalPoint TkRadialStripTopology::localPosition(const MeasurementPoint& mp) const {
0148   const float  // y = (L/cos(phi))*mp.y()*cos(phi)
0149       y(mp.y() * detHeight() + yCentreOfStripPlane()),
0150       x(yAxisOrientation() * yDistanceToIntersection(y) * tan15(stripAngle(mp.x())));
0151   return LocalPoint(x, y);
0152 }
0153
0154 MeasurementPoint TkRadialStripTopology::measurementPosition(const LocalPoint& lp) const {
0155   // phi is [pi/2 - conventional local phi], use atan2(x,y) rather than atan2(y,x)
0156   // clip   ( at pi/8 or detedge+tollerance?)
0157   float t = lp.x() / yDistanceToIntersection(lp.y());
0158 #ifdef MATH_STS
0159   statM.add(t);
0160 #endif
0161   const float phi = atanClip(t);
0162   return MeasurementPoint((phi - phiOfOneEdge()) * theAWidthInverse, (lp.y() - yCentreOfStripPlane()) / detHeight());
0163 }
0164
0165 LocalError TkRadialStripTopology::localError(float strip, float stripErr2) const {
0166   double phi = stripAngle(strip);
0167
0168   const double t1(tan15(phi)),  // std::tan(phif)), // (vdt::fast_tanf(phif)),
0169       t2(t1 * t1),
0170       // s1(std::sin(phi)), c1(std::cos(phi)),
0171       // cs(s1*c1), s2(s1*s1), c2(1-s2), // rotation matrix
0172
0173       tt(stripErr2 * std::pow(centreToIntersection() * angularWidth(), 2.f)),  // tangential sigma^2   *c2
0174       rr(theRadialSigma),  // radial sigma^2( uniform prob density along strip)  *c2
0175
0176       xx(tt + t2 * rr), yy(t2 * tt + rr), xy(t1 * (rr - tt));
0177
0178   return LocalError(xx, xy, yy);
0179 }
0180
0181 LocalError TkRadialStripTopology::localError(const MeasurementPoint& mp, const MeasurementError& me) const {
0182   const double phi(stripAngle(mp.x())), s1(std::sin(phi)), c1(std::cos(phi)), cs(s1 * c1), s2(s1 * s1),
0183       c2(1 - s2),  // rotation matrix
0184
0185       T(angularWidth() * (centreToIntersection() + yAxisOrientation() * mp.y() * detHeight()) /
0186         c1),                // tangential measurement unit (local pitch)
0187       R(detHeight() / c1),  // radial measurement unit (strip length)
0188       tt(me.uu() * T * T),  // tangential sigma^2
0189       rr(me.vv() * R * R),  // radial sigma^2
0190       tr(me.uv() * T * R),
0191
0192       xx(c2 * tt + 2 * cs * tr + s2 * rr), yy(s2 * tt - 2 * cs * tr + c2 * rr), xy(cs * (rr - tt) + tr * (c2 - s2));
0193
0194   return LocalError(xx, xy, yy);
0195 }
0196
0197 MeasurementError TkRadialStripTopology::measurementError(const LocalPoint& p, const LocalError& e) const {
0198   const double yHitToInter(yDistanceToIntersection(p.y())),
0199       t(yAxisOrientation() * p.x() / yHitToInter),  // tan(strip angle)
0200       cs(t / (1 + t * t)), s2(t * cs), c2(1 - s2),  // rotation matrix
0201
0202       T2(1. / (std::pow(angularWidth(), 2.f) *
0203                (std::pow(p.x(), 2.f) + std::pow(yHitToInter, 2)))),  // 1./tangential measurement unit (local pitch) ^2
0204       R2(c2 / std::pow(detHeight(), 2.f)),                           // 1./ radial measurement unit (strip length) ^2
0205
0206       uu((c2 * e.xx() - 2 * cs * e.xy() + s2 * e.yy()) * T2), vv((s2 * e.xx() + 2 * cs * e.xy() + c2 * e.yy()) * R2),
0207       uv((cs * (e.xx() - e.yy()) + e.xy() * (c2 - s2)) * std::sqrt(T2 * R2));
0208
0209   return MeasurementError(uu, uv, vv);
0210 }
0211
0212 // The local pitch is the local x width of the strip at the local (x,y)
0213 float TkRadialStripTopology::localPitch(const LocalPoint& lp) const {
0214   // this should be ~ y*(tan(phi+aw)-tan(phi)) = -x + y*(tan(aw)+tan(phi))/(1.f-tan(aw)*tan(phi)) tan(phi)=x/y
0215   float y = yDistanceToIntersection(lp.y());
0216   float x = std::abs(lp.x());
0217   return y * (y * theTanAW + x) / (y - theTanAW * x) - x;
0218 }
0219
0220 /* old version
0221 float TkRadialStripTopology::localPitch(const LocalPoint& lp) const {
0222   // this should be ~ y*(tan(phi+aw)-tan(phi)) = -tan(phi) + (tan(aw)+tan(phi))/(1.f-tan(aw)*tan(phi))
0223   const int istrip = std::min(nstrips(), static_cast<int>(strip(lp)) + 1); // which strip number
0224   float fangle = stripAngle(static_cast<float>(istrip) - 0.5); // angle of strip centre
0225   float p =
0226     yDistanceToIntersection( lp.y() ) * std::sin(angularWidth()) /
0227     std::pow( std::cos(fangle-0.5f*angularWidth()), 2.f);
0228
0229   float theTanAW = std::tan(theAngularWidth);
0230   float y =  yDistanceToIntersection( lp.y() );
0231   float x = std::abs(lp.x());
0232   float myP = y*(y*theTanAW+x)/(y-theTanAW*x)-x; // (y*theTanAW+x)/(1.f-theTanAW*x/y)-x;
0233   edm::LogVerbatim("CommonTopologies") << "localPitch " << p << " " << myP;
0234
0235   return p;
0236
0237 }
0238 */