Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:53

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