CommonTopologies/src/TrapezoidalStripTopology.cc

0001 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include <iostream>
0005 #include <cmath>
0006 #include <algorithm>
0007
0008 //#define EDM_ML_DEBUG
0009
0010 TrapezoidalStripTopology::TrapezoidalStripTopology(int ns, float p, float l, float r0)
0011     : theNumberOfStrips(ns), thePitch(p), theDistToBeam(r0), theDetHeight(l) {
0012   theOffset = -theNumberOfStrips / 2. * thePitch;
0013   theYAxOr = 1;
0014 #ifdef EDM_ML_DEBUG
0015   edm::LogVerbatim("CommonTopologies") << "Constructing TrapezoidalStripTopology with nstrips = " << ns
0016                                        << " pitch = " << p << " length = " << l << " r0 =" << r0;
0017 #endif
0018 }
0019
0020 TrapezoidalStripTopology::TrapezoidalStripTopology(int ns, float p, float l, float r0, int yAx)
0021     : theNumberOfStrips(ns), thePitch(p), theDistToBeam(r0), theDetHeight(l), theYAxOr(yAx) {
0022   theOffset = -theNumberOfStrips / 2. * thePitch;
0023 #ifdef EDM_ML_DEBUG
0024   edm::LogVerbatim("CommonTopologies") << "Constructing TrapezoidalStripTopology with nstrips = " << ns
0025                                        << " pitch = " << p << " length = " << l << " r0 =" << r0
0026                                        << " yAxOrientation =" << yAx;
0027 #endif
0028 }
0029
0030 LocalPoint TrapezoidalStripTopology::localPosition(float strip) const {
0031   return LocalPoint(strip * thePitch + theOffset, 0.0);
0032 }
0033
0034 LocalPoint TrapezoidalStripTopology::localPosition(const MeasurementPoint& mp) const {
0035   float y = mp.y() * theDetHeight;
0036   float x = (mp.x() * thePitch + theOffset) * (theYAxOr * y + theDistToBeam) / theDistToBeam;
0037   return LocalPoint(x, y);
0038 }
0039
0040 LocalError TrapezoidalStripTopology::localError(float strip, float stripErr2) const {
0041   float lt, lc2, ls2, lslc;
0042   float localL2, localP2;
0043   float sl2, sp2;
0044   // angle from strip to local frame (see CMS TN / 95-170)
0045   lt = -(strip * thePitch + theOffset) * theYAxOr / theDistToBeam;
0046   lc2 = 1.f / (1. + lt * lt);
0047   lslc = lt * lc2;
0048   ls2 = 1.f - lc2;
0049   localL2 = theDetHeight * theDetHeight / lc2;
0050   localP2 = thePitch * thePitch * lc2;
0051   sl2 = localL2 / 12.;
0052   sp2 = stripErr2 * localP2;
0053   return LocalError(lc2 * sp2 + ls2 * sl2, lslc * (sp2 - sl2), ls2 * sp2 + lc2 * sl2);
0054 }
0055
0056 LocalError TrapezoidalStripTopology::localError(const MeasurementPoint& mp, const MeasurementError& merr) const {
0057   float lt, lc2, ls2, lslc;
0058   float localL, localP;
0059   float sl2, sp2, spl;
0060   // angle from strip to local frame (see CMS TN / 95-170)
0061   lt = -(mp.x() * thePitch + theOffset) * theYAxOr / theDistToBeam;
0062   lc2 = 1. / (1. + lt * lt);
0063   lslc = lt * lc2;
0064   ls2 = 1.f - lc2;
0065   localL = theDetHeight / std::sqrt(lc2);
0066   localP = localPitch(localPosition(mp));
0067   sp2 = merr.uu() * localP * localP;
0068   sl2 = merr.vv() * localL * localL;
0069   spl = merr.uv() * localP * localL;
0070   return LocalError(lc2 * sp2 + ls2 * sl2 - 2 * lslc * spl,
0071                     lslc * (sp2 - sl2) + (lc2 - ls2) * spl,
0072                     ls2 * sp2 + lc2 * sl2 + 2 * lslc * spl);
0073 }
0074
0075 float TrapezoidalStripTopology::strip(const LocalPoint& lp) const {
0076   float aStrip = ((lp.x() * theDistToBeam / (theYAxOr * lp.y() + theDistToBeam)) - theOffset) / thePitch;
0077   if (aStrip < 0)
0078     aStrip = 0;
0079   else if (aStrip > theNumberOfStrips)
0080     aStrip = theNumberOfStrips;
0081   return aStrip;
0082 }
0083
0084 MeasurementPoint TrapezoidalStripTopology::measurementPosition(const LocalPoint& lp) const {
0085   return MeasurementPoint(((lp.x() * theDistToBeam / (theYAxOr * lp.y() + theDistToBeam)) - theOffset) / thePitch,
0086                           lp.y() / theDetHeight);
0087 }
0088
0089 MeasurementError TrapezoidalStripTopology::measurementError(const LocalPoint& lp, const LocalError& lerr) const {
0090   float lt, lc2, ls2, lslc;
0091   float localL, localP;
0092   float sl2, sp2, spl;
0093   lt = -lp.x() / (theYAxOr * lp.y() + theDistToBeam) * theYAxOr;
0094   lc2 = 1. / (1. + lt * lt);
0095   lslc = lt * lc2;
0096   ls2 = 1. - lc2;
0097   localL = theDetHeight / std::sqrt(lc2);
0098   localP = localPitch(lp);
0099   sp2 = lc2 * lerr.xx() + ls2 * lerr.yy() + 2 * lslc * lerr.xy();
0100   sl2 = ls2 * lerr.xx() + lc2 * lerr.yy() - 2 * lslc * lerr.xy();
0101   spl = lslc * (lerr.yy() - lerr.xx()) + (lc2 - ls2) * lerr.xy();
0102   return MeasurementError(sp2 / (localP * localP), spl / (localP * localL), sl2 / (localL * localL));
0103 }
0104
0105 int TrapezoidalStripTopology::channel(const LocalPoint& lp) const {
0106   return std::min(int(strip(lp)), theNumberOfStrips - 1);
0107 }
0108
0109 float TrapezoidalStripTopology::pitch() const { return thePitch; }
0110
0111 float TrapezoidalStripTopology::localPitch(const LocalPoint& lp) const {
0112   float x = lp.x();
0113   float y = theYAxOr * lp.y() + theDistToBeam;
0114   return thePitch * y / (theDistToBeam * std::sqrt(1.f + x * x / (y * y)));
0115 }
0116
0117 float TrapezoidalStripTopology::stripAngle(float strip) const {
0118   return std::atan(-(strip * thePitch + theOffset) * theYAxOr / theDistToBeam);
0119 }
0120
0121 int TrapezoidalStripTopology::nstrips() const { return theNumberOfStrips; }
0122
0123 float TrapezoidalStripTopology::shiftOffset(float pitch_fraction) {
0124   theOffset += thePitch * pitch_fraction;
0125   return theOffset;
0126 }
0127
0128 float TrapezoidalStripTopology::localStripLength(const LocalPoint& lp) const {
0129   float ltan = -lp.x() / (theYAxOr * lp.y() + theDistToBeam) * theYAxOr;
0130   float localL = theDetHeight * std::sqrt(1.f + ltan * ltan);
0131   //  float lcos2 = 1.f/(1.f+ltan*ltan);
0132   //  float localL = theDetHeight / std::sqrt(lcos2);
0133
0134   return localL;
0135 }