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