Merge pull request cms-sw#2296 from VinInn/TrackGeom

Track geom

Author: ktf

DataFormats/GeometrySurface/interface/GloballyPositioned.h

@@ -59,6 +59,14 @@ class GloballyPositioned {
     ToLocal(GloballyPositioned const & frame) :
       thePos(frame.position()), theRot(frame.rotation().transposed()){}
 
+    LocalPoint operator()(const GlobalPoint& gp) const {
+         return toLocal(gp);
+    }
+
+    LocalVector operator()(const GlobalVector& gv) const {
+       	 return	toLocal(gv);
+    }
+
     LocalPoint toLocal( const GlobalPoint& gp) const {
       return LocalPoint( theRot.multiplyInverse( gp.basicVector() -
 			 thePos.basicVector()) 
@@ -69,7 +77,7 @@ class GloballyPositioned {
       return LocalVector(theRot.multiplyInverse(gv.basicVector()));
     } 
 
-  private:
+  // private:
     PositionType  thePos;
     RotationType  theRot;
 

DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h

@@ -17,6 +17,7 @@ class RectangularPlaneBounds GCC11_FINAL : public Bounds {
   /// Construct from  half width (extension in local X),
   /// half length (Y) and half thickness (Z)
   RectangularPlaneBounds( float w, float h, float t);
+  ~RectangularPlaneBounds();
 
   /// Lenght along local Y
   virtual float length()    const { return 2*halfLength;}
@@ -26,11 +27,25 @@ class RectangularPlaneBounds GCC11_FINAL : public Bounds {
   virtual float thickness() const { return 2*halfThickness;}
 
   // basic bounds function
-  virtual bool inside( const Local2DPoint& p) const;
+  virtual bool inside( const Local2DPoint& p) const {
+    return
+     (std::abs(p.x()) < halfWidth) &
+     (std::abs(p.y()) < halfLength);
+  }
 
-  virtual bool inside( const Local3DPoint& p) const;
+  virtual bool inside( const Local3DPoint& p) const {
+    return
+     (std::abs(p.x()) < halfWidth) &
+     (std::abs(p.y()) < halfLength) &
+     (std::abs(p.z()) < halfThickness);
+  }
 
-  virtual bool inside(const Local2DPoint& p, float tollerance) const;
+
+
+  virtual bool inside(const Local2DPoint& p, float tollerance) const {
+    return (std::abs(p.x()) < (halfWidth  + tollerance) ) &
+           (std::abs(p.y()) < (halfLength + tollerance) );
+  }
 
 
   virtual bool inside( const Local3DPoint& p, const LocalError& err,

DataFormats/GeometrySurface/interface/SimpleDiskBounds.h

@@ -20,7 +20,10 @@ class SimpleDiskBounds GCC11_FINAL : public Bounds {
   virtual float width()     const { return 2*theRmax;}
   virtual float thickness() const { return theZmax-theZmin;}
 
-  virtual bool inside( const Local3DPoint& p) const;
+  virtual bool inside( const Local3DPoint& p) const {
+    return  ((p.z() > theZmin) & (p.z() < theZmax)) &&
+    ( (p.perp2() > theRmin*theRmin) & (p.perp2() < theRmax*theRmax) );
+  }
 
   virtual bool inside( const Local3DPoint& p, const LocalError& err, float scale) const;
 
@@ -32,6 +35,9 @@ class SimpleDiskBounds GCC11_FINAL : public Bounds {
   float innerRadius() const {return theRmin;}
   float outerRadius() const {return theRmax;}
 
+  float minZ() const { return theZmin;}
+  float        maxZ() const { return theZmax;}
+
 private:
   float theRmin;
   float theRmax;

DataFormats/GeometrySurface/src/RectangularPlaneBounds.cc

@@ -6,49 +6,23 @@ RectangularPlaneBounds::RectangularPlaneBounds( float w, float h, float t) :
   halfWidth(w), halfLength(h), halfThickness(t) {}
 
 
-bool RectangularPlaneBounds::inside( const Local2DPoint& p) const {
-  return std::abs(p.x()) <= halfWidth && std::abs(p.y()) <= halfLength;
-}
-
-bool RectangularPlaneBounds::inside( const Local3DPoint& p) const {
-  return 
-    std::abs(p.x()) < halfWidth && 
-    std::abs(p.y()) < halfLength &&
-    std::abs(p.z()) < halfThickness;
-}
-
-bool RectangularPlaneBounds::inside(const Local2DPoint& p, float tollerance) const {
-  return std::abs(p.x()) < halfWidth  + tollerance &&
-         std::abs(p.y()) < halfLength + tollerance;
-
-}
+RectangularPlaneBounds::~RectangularPlaneBounds(){}
 
 bool RectangularPlaneBounds::inside(const Local3DPoint& p, const LocalError& err,
 				    float scale) const {
-  if(scale >=0){
-    return 
-      std::abs(p.z()) < halfThickness &&
-      (std::abs(p.x()) < halfWidth  || std::abs(p.x()) < halfWidth  + std::sqrt(err.xx())*scale) &&
-      (std::abs(p.y()) < halfLength || std::abs(p.y()) < halfLength + std::sqrt(err.yy())*scale);
-  }else{
-    return
-      std::abs(p.z()) < halfThickness &&
+  if( (scale >=0) && inside(p) ) return true;
+  return
+      std::abs(p.z()) <  halfThickness &&
       (std::abs(p.x()) < halfWidth  + std::sqrt(err.xx())*scale) &&
       (std::abs(p.y()) < halfLength + std::sqrt(err.yy())*scale);
-  }
 }
 
 bool RectangularPlaneBounds::inside( const Local2DPoint& p, const LocalError& err, 
 				     float scale) const {
-  if(scale >=0){
-    return 
-      (std::abs(p.x()) < halfWidth  || std::abs(p.x()) < halfWidth  + std::sqrt(err.xx())*scale) &&
-      (std::abs(p.y()) < halfLength || std::abs(p.y()) < halfLength + std::sqrt(err.yy())*scale);
-  }else{
-    return 
+  if( (scale >=0) && inside(p) ) return true;
+  return 
       (std::abs(p.x()) < halfWidth  + std::sqrt(err.xx())*scale) &&
       (std::abs(p.y()) < halfLength + std::sqrt(err.yy())*scale);
-  }
 }
 
 

DataFormats/GeometrySurface/src/SimpleDiskBounds.cc

@@ -11,10 +11,6 @@ SimpleDiskBounds::SimpleDiskBounds( float rmin, float rmax, float zmin, float zm
   if ( theZmin > theZmax) std::swap( theZmin, theZmax);
 }
 
-bool SimpleDiskBounds::inside( const Local3DPoint& p) const {
-  return p.z()    > theZmin && p.z()    < theZmax &&
-    p.perp() > theRmin && p.perp() < theRmax;
-}
 
 bool SimpleDiskBounds::inside( const Local2DPoint& p, const LocalError& err) const {
   return Bounds::inside(p,err);

DataFormats/GeometrySurface/src/TrapezoidalPlaneBounds.cc

@@ -12,32 +12,30 @@ TrapezoidalPlaneBounds::TrapezoidalPlaneBounds( float be, float te,
   // angle of the trapezoid for faster inside() implementation.
 
   offset = a * (te+be) / (te-be);  // check sign if te < be !!! 
-  tan_a = te / fabs(offset + a);
+  tan_a = te / std::abs(offset + a);
 }
 
 
 int TrapezoidalPlaneBounds::yAxisOrientation() const {
-  int yAx = 1;
-  if(hbotedge>htopedge) yAx = -1;
-  return yAx;
+  return (hbotedge>htopedge) ? -1 : 1;
 }
 
 bool TrapezoidalPlaneBounds::inside( const Local2DPoint& p) const {
-  return fabs(p.y()) < hapothem && 
-    fabs(p.x())/fabs(p.y()+offset) < tan_a;
+  return (std::abs(p.y()) < hapothem) & 
+    (std::abs(p.x()) < tan_a*std::abs(p.y()+offset));
 }
 
 bool TrapezoidalPlaneBounds::inside( const Local3DPoint& p) const {
-  return fabs(p.y()) < hapothem &&
-    fabs(p.x())/fabs(p.y()+offset) < tan_a &&
-    fabs(p.z()) < hthickness;
+  return ((std::abs(p.y()) < hapothem) & (std::abs(p.z()) < hthickness)) && std::abs(p.x()) < tan_a*std::abs(p.y()+offset);
 }
 
 bool TrapezoidalPlaneBounds::inside( const Local3DPoint& p,
 				     const LocalError& err, float scale) const {
-  TrapezoidalPlaneBounds tmp( hbotedge + sqrt(err.xx())*scale,
-			      htopedge + sqrt(err.xx())*scale,
-			      hapothem + sqrt(err.yy())*scale,
+  if (scale>=0 && inside(p)) return true;
+
+  TrapezoidalPlaneBounds tmp( hbotedge + std::sqrt(err.xx())*scale,
+			      htopedge + std::sqrt(err.xx())*scale,
+			      hapothem + std::sqrt(err.yy())*scale,
 			      hthickness);
   return tmp.inside(p);
 }

DataFormats/GeometryVector/interface/VectorUtil.h

@@ -4,56 +4,22 @@
 
 
 #include "DataFormats/GeometryVector/interface/Pi.h"
-#include "Math/VectorUtil.h" 
+#include "DataFormats/Math/interface/deltaR.h"
 #include <cmath>
 
 
 namespace Geom {
+   using reco::deltaPhi;
+   using reco::deltaR2;
+   using reco::deltaR;
 
-  /**
-     Find aximutal Angle difference between two generic vectors ( v2.Phi() - v1.Phi() ) 
-     The only requirements on the Vector classes is that they implement the Phi() method
-     \param v1  Vector of any type implementing the Phi() operator
-     \param v2  Vector of any type implementing the Phi() operator
-     \return  Phi difference
-     \f[ \Delta \phi = \phi_2 - \phi_1 \f]
-  */
-  inline double deltaBarePhi(double phi1, double phi2) { 
-    double dphi = phi2-phi1; 
-    if ( dphi > M_PI ) {
-      dphi -= 2.0*M_PI;
-    } else if ( dphi <= -M_PI ) {
-      dphi += 2.0*M_PI;
-    }
-    return dphi;
-  }
-  inline double deltaPhi(float phi1, float phi2) { 
-    using ROOT::Math::VectorUtil::Phi_mpi_pi;
-    return deltaBarePhi(Phi_mpi_pi(phi2),Phi_mpi_pi(phi1));
-  }
-  inline double deltaPhi(double phi1, double phi2) { 
-    using ROOT::Math::VectorUtil::Phi_mpi_pi;
-    return deltaBarePhi(Phi_mpi_pi(phi2),Phi_mpi_pi(phi1));
-  }
-  template <class Vector1, class Vector2> 
-  double deltaPhi( const Vector1 & v1, const Vector2 & v2) { 
-    return deltaBarePhi(v1.phi(),v2.phi()); 
-  }
 
 
   /** Definition of ordering of azimuthal angles.
    *  phi1 is less than phi2 if the angle covered by a point going from
    *  phi1 to phi2 in the counterclockwise direction is smaller than pi.
    *  It makes sense only if ALL phis are in a single hemisphere...
    */
-  /*
-  inline bool phiLess( float phi1, float phi2) {
-    float diff = fmod(phi2 - phi1, 2.0*M_PI);
-    // float diff = phi2-phi1; 
-    if ( diff < 0) diff += 2*M_PI;
-    return diff < M_PI;
-  }
-  */
   inline bool phiLess(float phi1, float phi2) {
     return deltaPhi(phi1,phi2)<0;
   }
@@ -62,29 +28,10 @@ namespace Geom {
   }
   template <class Vector1, class Vector2> 
   bool phiLess(const Vector1 & v1, const Vector2 & v2) {
-    return deltaPhi(v1,v2)<0.; 
+    return deltaPhi(v1.phi(),v2.phi())<0.; 
   }
 
 
-  /**
-     Find difference in pseudorapidity (Eta) and Phi betwen two generic vectors
-     The only requirements on the Vector classes is that they implement the Phi() and Eta() method
-     \param v1  Vector 1  
-     \param v2  Vector 2
-     \return   Angle between the two vectors
-     \f[ \Delta R = \sqrt{  ( \Delta \phi )^2 + ( \Delta \eta )^2 } \f]
-  */ 
-  template <class Vector1, class Vector2> 
-  double deltaR2( const Vector1 & v1, const Vector2 & v2) { 
-    double dphi = deltaPhi(v1,v2); 
-    double deta = v2.eta() - v1.eta(); 
-    return dphi*dphi + deta*deta; 
-  }
-  template <class Vector1, class Vector2> 
-  double deltaR( const Vector1 & v1, const Vector2 & v2) { 
-    return std::sqrt( deltaR2(v1,v2));
-  }
-  
 }
 
 #endif

DataFormats/GeometryVector/test/BuildFile.xml

@@ -1,3 +1,4 @@
 <use  name="DataFormats/GeometryVector"/>
-<bin  file="BasicVector_t.cpp" name="DataFormatsBasicVector_t">
-</bin>
+<bin  file="BasicVector_t.cpp" name="DataFormatsBasicVector_t"/>
+<bin  file="phiLess_t.cpp" name="GeomVectPhiLess_t"/>
+

DataFormats/GeometryVector/test/phiLess_t.cpp

@@ -0,0 +1,21 @@
+#include "DataFormats/GeometryVector/interface/VectorUtil.h"
+
+#include <cassert>
+
+
+int main() {
+
+  assert(Geom::phiLess(0.f,2.f));
+  assert(Geom::phiLess(6.f,0.f));
+  assert(Geom::phiLess(3.2f,0.f));
+  assert(Geom::phiLess(3.0f,3.2f));
+
+  assert(Geom::phiLess(-0.3f,0.f));
+  assert(Geom::phiLess(-0.3f,0.1f));
+  assert(Geom::phiLess(-3.0f,0.f));
+  assert(Geom::phiLess(3.0f,-3.0f));
+  assert(Geom::phiLess(0.f,-3.4f));
+
+
+
+}

Geometry/CommonTopologies/BuildFile.xml

@@ -1,6 +1,7 @@
 <use   name="FWCore/MessageLogger"/>
 <use   name="Utilities/General"/>
+<use   name="vdt"/>
 <export>
   <lib   name="1"/>
 </export>
-# <flags   CXXFLAGS="-O3 -fno-math-errno -funsafe-loop-optimizations -ftree-loop-if-convert-stores --param max-completely-peel-times=1"/>
+# <flags   CXXFLAGS="-O3 -ffast-math -funsafe-loop-optimizations -ftree-loop-if-convert-stores"/>

Geometry/CommonTopologies/interface/TkRadialStripTopology.h

@@ -119,7 +119,8 @@ class TkRadialStripTopology GCC11_FINAL : public RadialStripTopology {
    * whereas values 1, 2, ... nstrips correspond to the upper phi edges of
    * the strips.
    */
-  float stripAngle(float strip) const;
+  float stripAngle(float strip) const { return   yAxisOrientation() * (phiOfOneEdge() +  strip * angularWidth()) ;}
+
 
   /** 
    * Total number of strips 
@@ -246,6 +247,8 @@ class TkRadialStripTopology GCC11_FINAL : public RadialStripTopology {
   float theTanOfOneEdge;   // the positive tangent of the above...
   float theYAxisOrientation; // 1 means y axis going from smaller to larger side, -1 means opposite direction
   float yCentre; // Non-zero if offset in local y between midpoint of detector (strip plane) extent and local origin.
+  double theRadialSigma;    // radial sigma^2( uniform prob density along strip)
+
 };
 
 #endif