GEM uniform angular strip topology

Fireworks/Geometry/src/FWRecoGeometryESProducer.cc

@@ -33,6 +33,7 @@
 #include "Geometry/CommonTopologies/interface/StripTopology.h"
 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
+#include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
 
 #include "TNamed.h"
 #include "FWCore/ParameterSet/interface/FileInPath.h"

Geometry/CommonTopologies/interface/GEMStripTopology.h

@@ -0,0 +1,73 @@
+#ifndef Geometry_CommonTopologies_GEMStripTopology_H
+#define Geometry_CommonTopologies_GEMStripTopology_H
+
+/** \class GEMStripTopology
+ *  based on CSCRadialStripTopology and TrapezoidalStripTopology
+ *  \author Hyunyong Kim - TAMU
+ */
+
+#include "Geometry/CommonTopologies/interface/StripTopology.h"
+
+class GEMStripTopology final : public StripTopology {
+public:
+  GEMStripTopology(int ns, float aw, float dh, float r0);
+  GEMStripTopology(int ns, float aw, float dh, float r0, float yAx);
+  ~GEMStripTopology() override {}
+
+  using StripTopology::localPosition;
+  LocalPoint localPosition(float strip) const override;
+
+  LocalPoint localPosition(const MeasurementPoint&) const override;
+
+  using StripTopology::localError;
+  LocalError localError(float strip, float stripErr2) const override;
+  LocalError localError(const MeasurementPoint&, const MeasurementError&) const override;
+
+  float strip(const LocalPoint&) const override;
+
+  int nearestStrip(const LocalPoint&) const;
+
+  MeasurementPoint measurementPosition(const LocalPoint&) const override;
+
+  MeasurementError measurementError(const LocalPoint&, const LocalError&) const override;
+
+  int channel(const LocalPoint&) const override;
+
+  float phiPitch(void) const { return angularWidth(); }
+
+  float pitch() const override;
+
+  float localPitch(const LocalPoint&) const override;
+
+  float stripAngle(float strip) const override;
+
+  int nstrips() const override { return numberOfStrips_; }
+
+  float stripLength() const override { return detHeight_; }
+
+  float localStripLength(const LocalPoint&) const override;
+
+  float angularWidth() const { return angularWidth_; }
+  float detHeight() const { return detHeight_; }
+  float yExtentOfStripPlane() const { return detHeight_; }
+  float centreToIntersection() const { return centreToIntersection_; }
+  float radius() const { return centreToIntersection_; }
+  float originToIntersection() const { return (centreToIntersection_ - yCentre_); }
+
+  float yDistanceToIntersection(float y) const;
+  float xOfStrip(int strip, float y) const;
+  float yAxisOrientation() const { return yAxisOrientation_; }
+  float phiOfOneEdge() const { return phiOfOneEdge_; }
+  float yCentreOfStripPlane() const { return yCentre_; }
+
+private:
+  int numberOfStrips_;          // total no. of strips in plane of strips
+  float angularWidth_;          // angle subtended by each strip = phi pitch
+  float detHeight_;             // length of long symmetry axis = twice the apothem of the enclosing trapezoid
+  float centreToIntersection_;  // distance centre of detector face to intersection of edge strips (projected)
+  float phiOfOneEdge_;          // local 'phi' of one edge of plane of strips (I choose it negative!)
+  float yAxisOrientation_;      // 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.
+};
+
+#endif

Geometry/CommonTopologies/src/GEMStripTopology.cc

@@ -0,0 +1,134 @@
+/**GEMStripTopology
+ * based on CSCRadialStripTopology and TrapezoidalStripTopology
+ * \author Hyunyong Kim - TAMU
+ */
+#include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+
+GEMStripTopology::GEMStripTopology(int ns, float aw, float dh, float r0)
+    : numberOfStrips_(ns), angularWidth_(aw), detHeight_(dh), centreToIntersection_(r0) {
+  assert(angularWidth_ != 0);
+  assert(detHeight_ != 0);
+  yAxisOrientation_ = 1;
+  phiOfOneEdge_ = -(0.5 * numberOfStrips_) * angularWidth_ * yAxisOrientation_;
+  yCentre_ = 0;
+  LogTrace("GEMStripTopology") << "Constructing GEMStripTopology with"
+                               << " nstrips = " << ns << " angular width = " << aw << " det. height = " << dh
+                               << " r0 = " << r0 << "\n";
+}
+
+GEMStripTopology::GEMStripTopology(int ns, float aw, float dh, float r0, float yAx)
+    : numberOfStrips_(ns), angularWidth_(aw), detHeight_(dh), centreToIntersection_(r0), yAxisOrientation_(yAx) {
+  assert(angularWidth_ != 0);
+  assert(detHeight_ != 0);
+  phiOfOneEdge_ = -(0.5 * numberOfStrips_) * angularWidth_ * yAxisOrientation_;
+  yCentre_ = 0;
+  LogTrace("GEMStripTopology") << "Constructing GEMStripTopology with"
+                               << " nstrips = " << ns << " angular width = " << aw << " det. height = " << dh
+                               << " r0 = " << r0 << " yAxOrientation = " << yAx << "\n";
+}
+
+LocalPoint GEMStripTopology::localPosition(float strip) const {
+  return LocalPoint(yAxisOrientation() * originToIntersection() * tan(stripAngle(strip)), 0);
+}
+
+LocalPoint GEMStripTopology::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);
+}
+
+LocalError GEMStripTopology::localError(float strip, float stripErr2) const {
+  const double phi(stripAngle(strip)), t1(std::tan(phi)), t2(t1 * t1),
+      tt(stripErr2 * std::pow(centreToIntersection() * angularWidth(), 2)),  // tangential sigma^2   *c2
+      rr(std::pow(detHeight(), 2) * (1.f / 12.f)),  // 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 GEMStripTopology::localError(const MeasurementPoint& mp, const MeasurementError& me) const {
+  const double phi(stripAngle(mp.x())), s1(std::sin(phi)), c1(std::cos(phi));
+  assert(c1 != 0);
+  const double 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);
+}
+
+float GEMStripTopology::strip(const LocalPoint& lp) const {
+  const float  // phi is measured from y axis --> sign of angle is sign of x * yAxisOrientation --> use atan2(x,y), not atan2(y,x)
+      phi(std::atan2(lp.x(), yDistanceToIntersection(lp.y()))),
+      aStrip((phi - yAxisOrientation() * phiOfOneEdge()) / angularWidth());
+  return std::max(float(0), std::min((float)nstrips(), aStrip));
+}
+
+int GEMStripTopology::nearestStrip(const LocalPoint& lp) const {
+  return std::min(nstrips(), static_cast<int>(std::max(float(0), strip(lp))) + 1);
+}
+
+MeasurementPoint GEMStripTopology::measurementPosition(const LocalPoint& lp) const {
+  const float  // phi is [pi/2 - conventional local phi], use atan2(x,y) rather than atan2(y,x)
+      phi(yAxisOrientation() * std::atan2(lp.x(), yDistanceToIntersection(lp.y())));
+  return MeasurementPoint(yAxisOrientation() * (phi - phiOfOneEdge()) / angularWidth(),
+                          (lp.y() - yCentreOfStripPlane()) / detHeight());
+}
+
+MeasurementError GEMStripTopology::measurementError(const LocalPoint& p, const LocalError& e) const {
+  const double yHitToInter(yDistanceToIntersection(p.y()));
+  assert(yHitToInter != 0);
+  const double 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) *
+               (std::pow(p.x(), 2) + std::pow(yHitToInter, 2)))),  // 1./tangential measurement unit (local pitch) ^2
+      R2(c2 / std::pow(detHeight(), 2)),                           // 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);
+}
+
+int GEMStripTopology::channel(const LocalPoint& lp) const { return std::min(int(strip(lp)), numberOfStrips_ - 1); }
+
+float GEMStripTopology::pitch() const { return localPitch(LocalPoint(0, 0)); }
+
+float GEMStripTopology::localPitch(const LocalPoint& lp) const {
+  const int istrip = std::min(nstrips(), static_cast<int>(strip(lp)) + 1);  // which strip number
+  const float fangle = stripAngle(static_cast<float>(istrip) - 0.5);        // angle of strip centre
+  assert(std::cos(fangle - 0.5f * angularWidth()) != 0);
+  return yDistanceToIntersection(lp.y()) * std::sin(angularWidth()) /
+         std::pow(std::cos(fangle - 0.5f * angularWidth()), 2);
+}
+
+float GEMStripTopology::stripAngle(float strip) const {
+  return phiOfOneEdge() + yAxisOrientation() * strip * angularWidth();
+}
+
+float GEMStripTopology::localStripLength(const LocalPoint& lp) const {
+  assert(yDistanceToIntersection(lp.y()) != 0);
+  return detHeight() * std::sqrt(1.f + std::pow(lp.x() / yDistanceToIntersection(lp.y()), 2));
+}
+
+float GEMStripTopology::yDistanceToIntersection(float y) const {
+  return yAxisOrientation() * y + originToIntersection();
+}
+
+float GEMStripTopology::xOfStrip(int strip, float y) const {
+  return yAxisOrientation() * yDistanceToIntersection(y) * std::tan(stripAngle(static_cast<float>(strip) - 0.5));
+}

Geometry/GEMGeometry/src/GEMEtaPartition.cc

@@ -1,6 +1,6 @@
 #include "Geometry/GEMGeometry/interface/GEMEtaPartition.h"
 #include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h"
-#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
+#include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
 
 GEMEtaPartition::GEMEtaPartition(GEMDetId id, const BoundPlane::BoundPlanePointer& bp, GEMEtaPartitionSpecs* rrs)
     : GeomDet(bp), id_(id), specs_(rrs) {

Geometry/GEMGeometry/src/GEMEtaPartitionSpecs.cc

@@ -1,8 +1,8 @@
 #include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h"
-#include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
-#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
+#include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
 
 using namespace GeomDetEnumerators;
+using namespace angle_units::operators;
 
 GEMEtaPartitionSpecs::GEMEtaPartitionSpecs(SubDetector rss, const std::string& name, const GEMSpecs& pars)
     : GeomDetType(name, rss), _p(pars), _n(name) {
@@ -13,14 +13,17 @@ GEMEtaPartitionSpecs::GEMEtaPartitionSpecs(SubDetector rss, const std::string& n
     float r0 = h * (B + b) / (B - b);
     float striplength = h * 2;
     float strips = _p[3];
-    float pitch = (b + B) / strips;
+    float dphi = _p[5];
+    float phiPitch = dphi / strips;
+
     int nstrip = static_cast<int>(strips);
-    _top = new TrapezoidalStripTopology(nstrip, pitch, striplength, r0);
+    _top = new GEMStripTopology(nstrip, phiPitch, striplength, r0);
 
     float pads = _p[4];
-    float pad_pitch = (b + B) / pads;
+    float padPhiPitch = dphi / pads;
+
     int npad = static_cast<int>(pads);
-    _top_pad = new TrapezoidalStripTopology(npad, pad_pitch, striplength, r0);
+    _top_pad = new GEMStripTopology(npad, padPhiPitch, striplength, r0);
   } else {
     _top = nullptr;
     _top_pad = nullptr;

Geometry/GEMGeometryBuilder/data/GEMSpecs.xml

@@ -9,5 +9,9 @@
       <PartSelector path="//GHA1.."/>
       <Parameter name="nPads" value="192"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v10/GEMSpecs.xml

@@ -25,5 +25,16 @@
       <PartSelector path="//GHA2.."/>
       <Parameter name="nPads" value="384"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE21" eval="true">
+      <PartSelector path="//GHA2.."/>
+      <Parameter name="dPhi" value="20.30*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiME0" eval="true">
+      <PartSelector path="//ME0Box/ME0L/GHA0.."/>
+      <Parameter name="dPhi" value="20.22*deg"/>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v11/GEMSpecs.xml

@@ -25,5 +25,17 @@
       <PartSelector path="//GHA2.."/>
       <Parameter name="nPads" value="384"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE21" eval="true">
+      <PartSelector path="//GHA2.."/>
+      <Parameter name="dPhi" value="20.30*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiME0" eval="true">
+      <PartSelector path="//GHA0.."/>
+      <Parameter name="dPhi" value="20.22*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v12/GEMSpecs.xml

@@ -25,5 +25,17 @@
       <PartSelector path="//GHA2.."/>
       <Parameter name="nPads" value="384"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE21" eval="true">
+      <PartSelector path="//GHA2.."/>
+      <Parameter name="dPhi" value="20.30*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE0" eval="true">
+      <PartSelector path="//GHA0.."/>
+      <Parameter name="dPhi" value="20.22*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v3/GEMSpecs.xml

@@ -9,5 +9,9 @@
       <PartSelector path="//GHA1.."/>
       <Parameter name="nPads" value="96"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v4/GEMSpecs.xml

@@ -9,5 +9,9 @@
       <PartSelector path="//GHA1.."/>
       <Parameter name="nPads" value="192"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v5/GEMSpecs.xml

@@ -17,5 +17,13 @@
       <PartSelector path="//GHA2.."/>
       <Parameter name="nPads" value="384"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE21" eval="true">
+      <PartSelector path="//GHA2.."/>
+      <Parameter name="dPhi" value="20.30*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v7/GEMSpecs.xml

@@ -25,5 +25,17 @@
       <PartSelector path="//GHA2.."/>
       <Parameter name="nPads" value="384"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE21" eval="true">
+      <PartSelector path="//GHA2.."/>
+      <Parameter name="dPhi" value="20.30*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiME0" eval="true">
+      <PartSelector path="//GHA0.."/>
+      <Parameter name="dPhi" value="20.22*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>

Geometry/GEMGeometryBuilder/data/v9/GEMSpecs.xml

@@ -25,5 +25,17 @@
       <PartSelector path="//GHA2.."/>
       <Parameter name="nPads" value="384"/>
     </SpecPar>
+    <SpecPar name="dPhiGE11" eval="true">
+      <PartSelector path="//GHA1.."/>
+      <Parameter name="dPhi" value="10.15*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE21" eval="true">
+      <PartSelector path="//GHA2.."/>
+      <Parameter name="dPhi" value="20.30*deg"/>
+    </SpecPar>
+    <SpecPar name="dPhiGE0" eval="true">
+      <PartSelector path="//GHA0.."/>
+      <Parameter name="dPhi" value="20.22*deg"/>
+    </SpecPar>
   </SpecParSection>
 </DDDefinition>