TkRadialStripTopology.cc

Go to the documentation of this file.

#include <Geometry/CommonTopologies/interface/TkRadialStripTopology.h>
#include <FWCore/MessageLogger/interface/MessageLogger.h>

#include <iostream>
#include <cmath>
#include <algorithm>
#include <cassert>

#include <vdt/vdtMath.h>

#ifdef MATH_STS
#include<iostream>
#endif
namespace {

#ifdef MATH_STS
  struct Stat {
    Stat(const char * in) : name(in){};
    ~Stat() {
      std::cout << name << ": atan0 calls tot/large/over1: " << natan << "/" << nlarge <<"/" << over1 << std::endl; 
    }

    void add(float t) { 
      auto at = std::abs(t);
      ++natan;
      if (at>0.40f) ++nlarge;
      if (at>1.0) ++over1;
    }
    const char * name;
    long long natan=0;
    long long nlarge=0;
    long long over1=0;    
  };

  Stat statM("mpos");
  Stat statS("span");
#endif


   // valid for |x| < 0.15  (better then 10^-9
   inline double tan15(double x) {
      return x * (1 + (x*x) * (0.33331906795501708984375 + (x*x) * 0.135160386562347412109375));
   }

  // valid for z < pi/8
  //  x * (1 + x*x * (-0.33322894573211669921875 + x*x * (0.1967026889324188232421875 + x*x * (-0.11053790152072906494140625))))  // .1e-7 by Sollya
  inline 
  float atan0(float t) {
    auto z=t;
    // if( t > 0.4142135623730950f ) // * tan pi/8 
    // z = (t-1.0f)/(t+1.0f);
    float z2 = z * z;
    float ret=
      ((( 8.05374449538e-2f * z2
      - 1.38776856032E-1f) * z2
    + 1.99777106478E-1f) * z2
       - 3.33329491539E-1f) * z2 * z
      + z;
    // if( t > 0.4142135623730950f ) ret +=0.7853981633974483096f;
    return ret;
  }

  inline 
  float atanClip(float t) {
    constexpr float tanPi8 = 0.4142135623730950;
    constexpr float pio8 = 3.141592653589793238/8;
    float at = std::abs(t);
    return std::copysign( (at< tanPi8 )  ? atan0(at) : pio8, t );
  }

}

TkRadialStripTopology::TkRadialStripTopology(int ns, float aw, float dh, float r, int yAx, float yMid) :
  theNumberOfStrips(ns), theAngularWidth(aw), theAWidthInverse(1.f/aw),theTanAW(std::tan(aw)),
  theDetHeight(dh), theCentreToIntersection(r),
  theYAxisOrientation(yAx), yCentre( yMid),
  theRadialSigma(std::pow(dh, 2.f) * (1.f/12.f)) {   
  // Angular offset of extreme edge of detector, so that angle is
  // zero for a strip lying along local y axis = long symmetry axis of plane of strips
  thePhiOfOneEdge = -(0.5*theNumberOfStrips) * theAngularWidth; // always negative!
  theTanOfOneEdge = std::tan(std::abs(thePhiOfOneEdge));
  assert(std::abs(thePhiOfOneEdge)<0.15); // 

  LogTrace("TkRadialStripTopology") << "TkRadialStripTopology: constructed with"
        << " strips = " << ns
        << " width = " << aw << " rad "
        << " det_height = " << dh
        << " ctoi = " << r 
        << " phi_edge = " << thePhiOfOneEdge << " rad "
        << " y_ax_ori = " << theYAxisOrientation
    << " y_det_centre = " << yCentre 
        << "\n";
}    

int TkRadialStripTopology::channel(const LocalPoint& lp) const { return   std::min( int( strip(lp) ), theNumberOfStrips-1 ) ;}

int TkRadialStripTopology::nearestStrip(const LocalPoint & lp) const {   return std::min( nstrips(), static_cast<int>( std::max(float(0), strip(lp)) ) + 1);}

float TkRadialStripTopology::yDistanceToIntersection( float y ) const { return   yAxisOrientation()*y + originToIntersection() ;}

float TkRadialStripTopology::localStripLength(const LocalPoint& lp) const {  
  return detHeight() * std::sqrt(1.f + std::pow( lp.x()/yDistanceToIntersection(lp.y()), 2.f) );
}

float TkRadialStripTopology::xOfStrip(int strip, float y) const { 
  return   
    yAxisOrientation() * yDistanceToIntersection( y ) * std::tan( stripAngle(static_cast<float>(strip) - 0.5f ) );
}

float TkRadialStripTopology::strip(const LocalPoint& lp) const {
     // phi is measured from y axis --> sign of angle is sign of x * yAxisOrientation --> use atan2(x,y), not atan2(y,x)
  const float phi =  atanClip(lp.x()/yDistanceToIntersection(lp.y()));
  const float aStrip = ( phi -  phiOfOneEdge() )*theAWidthInverse;
  return  std::max(float(0), std::min( (float)nstrips(), aStrip ));
}

float TkRadialStripTopology::coveredStrips(const LocalPoint& lp1, const LocalPoint& lp2)  const {
  // http://en.wikipedia.org/wiki/List_of_trigonometric_identities#Angle_sum_and_difference_identities
  // atan(a)-atan(b) = atan( (a-b)/(1+a*b) )  
  // avoid divisions
  // float t1 = lp1.x()/yDistanceToIntersection( lp1.y() );
  // float t2 = lp2.x()/yDistanceToIntersection( lp2.y() );
  // float t = (t1-t2)/(1.+t1*t2);
  auto y1 = yDistanceToIntersection( lp1.y() );
  auto y2 = yDistanceToIntersection( lp2.y() );
  auto x1 = lp1.x();
  auto x2 = lp2.x();

  auto t = (y2*x1 -y1*x2)/(y1*y2 + x1*x2);


#ifdef MATH_STS
  statS.add(t);
#endif
//   std::cout << "atans " << atanClip(t)
//                        <<" "<< std::atan2(lp1.x(),yDistanceToIntersection(lp1.y()) ) 
//                               -std::atan2(lp2.x(),yDistanceToIntersection(lp2.y()) ) << std::endl;
  // clip???
  return atanClip(t)*theAWidthInverse;
  //   return (measurementPosition(lp1)-measurementPosition(lp2)).x();
}  

LocalPoint TkRadialStripTopology::localPosition(float strip) const {
  return LocalPoint( yAxisOrientation() * originToIntersection() * std::tan( stripAngle(strip) ), 0 );
}

LocalPoint TkRadialStripTopology::localPosition(const MeasurementPoint& mp) const {
  const float  // y = (L/cos(phi))*mp.y()*cos(phi) 
    y( mp.y()*detHeight()  +  yCentreOfStripPlane() ),
    x( yAxisOrientation() * yDistanceToIntersection( y ) * std::tan ( stripAngle( mp.x() ) ) );
  return LocalPoint( x, y );
}

MeasurementPoint TkRadialStripTopology::measurementPosition(const LocalPoint& lp) const {
  // phi is [pi/2 - conventional local phi], use atan2(x,y) rather than atan2(y,x)
  // clip   ( at pi/8 or detedge+tollerance?)
  float t =  lp.x()/yDistanceToIntersection(lp.y());
#ifdef MATH_STS
  statM.add(t);
#endif
  const float phi = atanClip(t);
  return MeasurementPoint( ( phi-phiOfOneEdge() )*theAWidthInverse,
               ( lp.y() - yCentreOfStripPlane() )        / detHeight() );
}


LocalError TkRadialStripTopology::localError(float strip, float stripErr2) const {
  double phi = stripAngle(strip);

  const double
    t1(tan15(phi)),    // std::tan(phif)), // (vdt::fast_tanf(phif)), 
    t2(t1*t1),
    // s1(std::sin(phi)), c1(std::cos(phi)),
    // cs(s1*c1), s2(s1*s1), c2(1-s2), // rotation matrix

    tt( stripErr2 * std::pow( centreToIntersection()*angularWidth() ,2.f) ), // tangential sigma^2   *c2
    rr( theRadialSigma),                                   // radial sigma^2( uniform prob density along strip)  *c2

    xx( tt + t2*rr  ),
    yy( t2*tt + rr  ),
    xy( t1*( rr - tt ) );
  
  return LocalError( xx, xy, yy );
}

LocalError TkRadialStripTopology::localError(const MeasurementPoint& mp, const MeasurementError& me) const {
  const double
    phi(stripAngle(mp.x())), s1(std::sin(phi)), c1(std::cos(phi)),
    cs(s1*c1), s2(s1*s1), c2(1-s2), // rotation matrix

    T( angularWidth() * ( centreToIntersection() + yAxisOrientation()*mp.y()*detHeight()) / c1 ), // tangential measurement unit (local pitch)
    R( detHeight()/ c1 ),  // radial measurement unit (strip length)
    tt( me.uu() * T*T ),   // tangential sigma^2
    rr( me.vv() * R*R   ), // radial sigma^2
    tr( me.uv() * T*R ),  
    
    xx(  c2*tt  +  2*cs*tr  +  s2*rr      ),
    yy(  s2*tt  -  2*cs*tr  +  c2*rr      ),
    xy( cs*( rr - tt )  +  tr*( c2 - s2 ) );

  return LocalError( xx, xy, yy );
}

MeasurementError TkRadialStripTopology::measurementError(const LocalPoint& p,  const LocalError& e) const {
  const double
    yHitToInter(yDistanceToIntersection(p.y())),
    t(yAxisOrientation() * p.x() / yHitToInter),   // tan(strip angle) 
    cs(t/(1+t*t)), s2(t*cs), c2(1-s2),             // rotation matrix

    T2( 1./(std::pow(angularWidth(),2.f) * ( std::pow(p.x(),2.f) + std::pow(yHitToInter,2)) )), // 1./tangential measurement unit (local pitch) ^2
    R2( c2/std::pow(detHeight(),2.f) ),                                    // 1./ radial measurement unit (strip length) ^2

    uu(       ( c2*e.xx() - 2*cs*e.xy() + s2*e.yy() )   * T2 ),
    vv(       ( s2*e.xx() + 2*cs*e.xy() + c2*e.yy() )   * R2 ),
    uv( ( cs*( e.xx() - e.yy() ) + e.xy()*( c2 - s2 ) )  * std::sqrt (T2*R2) );
  
  return MeasurementError(uu, uv, vv);
}
 

 // The local pitch is the local x width of the strip at the local (x,y)
float TkRadialStripTopology::localPitch(const LocalPoint& lp) const { 
  // 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
  float y =  yDistanceToIntersection( lp.y() );
  float x = std::abs(lp.x());
  return y*(y*theTanAW+x)/(y-theTanAW*x)-x;
}
 
/* old version
float TkRadialStripTopology::localPitch(const LocalPoint& lp) const { 
  // this should be ~ y*(tan(phi+aw)-tan(phi)) = -tan(phi) + (tan(aw)+tan(phi))/(1.f-tan(aw)*tan(phi)) 
  const int istrip = std::min(nstrips(), static_cast<int>(strip(lp)) + 1); // which strip number
  float fangle = stripAngle(static_cast<float>(istrip) - 0.5); // angle of strip centre
  float p =
    yDistanceToIntersection( lp.y() ) * std::sin(angularWidth()) /
    std::pow( std::cos(fangle-0.5f*angularWidth()), 2.f);

  float theTanAW = std::tan(theAngularWidth);
  float y =  yDistanceToIntersection( lp.y() );
  float x = std::abs(lp.x());
  float myP = y*(y*theTanAW+x)/(y-theTanAW*x)-x; // (y*theTanAW+x)/(1.f-theTanAW*x/y)-x;
  std::cout << "localPitch " << p << " " << myP << std::endl;

  return p;

}
*/