CSCGeometry/src/CSCLayerGeometry.cc

0001 // This is CSCLayerGeometry.cc
0002
0003 #include <Geometry/CSCGeometry/interface/CSCGeometry.h>
0004 #include <Geometry/CSCGeometry/interface/CSCLayerGeometry.h>
0005
0006 #include <Geometry/CSCGeometry/src/CSCUngangedStripTopology.h>
0007 #include <Geometry/CSCGeometry/src/CSCGangedStripTopology.h>
0008 #include <Geometry/CSCGeometry/interface/CSCWireGroupPackage.h>
0009
0010 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0011
0012 #include "DataFormats/Math/interface/GeantUnits.h"
0013
0014 #include <algorithm>
0015 #include <iostream>
0016
0017 using namespace geant_units;
0018 using namespace geant_units::operators;
0019
0020 // Complicated initialization list since the chamber Bounds passed in must have its thickness reset to that of the layer
0021 // Note for half-widths, t + b = 2w and the TPB only has accessors for t and w so b = 2w - t
0022
0023 CSCLayerGeometry::CSCLayerGeometry(const CSCGeometry* geom,
0024                                    int iChamberType,
0025                                    const TrapezoidalPlaneBounds& bounds,
0026                                    int nstrips,
0027                                    float stripOffset,
0028                                    float stripPhiPitch,
0029                                    float whereStripsMeet,
0030                                    float extentOfStripPlane,
0031                                    float yCentreOfStripPlane,
0032                                    const CSCWireGroupPackage& wg,
0033                                    float wireAngleInDegrees,
0034                                    double yOfFirstWire,
0035                                    float hThickness)
0036     : TrapezoidalPlaneBounds(
0037           bounds.widthAtHalfLength() - bounds.width() / 2., bounds.width() / 2., bounds.length() / 2., hThickness),
0038       theWireTopology(nullptr),
0039       theStripTopology(nullptr),
0040       hBottomEdge(bounds.widthAtHalfLength() - bounds.width() / 2.),
0041       hTopEdge(bounds.width() / 2.),
0042       apothem(bounds.length() / 2.),
0043       myName("CSCLayerGeometry"),
0044       chamberType(iChamberType) {
0045   LogTrace("CSCLayerGeometry|CSC") << myName << ": being constructed, this=" << this;
0046
0047   // Ganged strips in ME1A?
0048   bool gangedME1A = (iChamberType == 1 && geom->gangedStrips());
0049
0050   CSCStripTopology* aStripTopology = new CSCUngangedStripTopology(
0051       nstrips, stripPhiPitch, extentOfStripPlane, whereStripsMeet, stripOffset, yCentreOfStripPlane);
0052
0053   if (gangedME1A) {
0054     theStripTopology = new CSCGangedStripTopology(*aStripTopology, 16);
0055     delete aStripTopology;
0056   } else {
0057     theStripTopology = aStripTopology;
0058   }
0059
0060   if (!geom->realWireGeometry()) {
0061     // Approximate ORCA_8_8_0 and earlier calculated geometry...
0062     float wangler = convertDegToRad(wireAngleInDegrees);
0063     float wireCos = cos(wangler);
0064     float wireSin = sin(wangler);
0065     float y2 = apothem * wireCos + hBottomEdge * fabs(wireSin);
0066     float wireSpacing = convertMmToCm(wg.wireSpacing);
0067     float wireOffset = -y2 + wireSpacing / 2.;
0068     yOfFirstWire = wireOffset / wireCos;
0069   }
0070
0071   theWireTopology = new CSCWireTopology(wg, yOfFirstWire, wireAngleInDegrees);
0072   LogTrace("CSCLayerGeometry|CSC") << myName << ": constructed: " << *this;
0073 }
0074
0075 CSCLayerGeometry::CSCLayerGeometry(const CSCLayerGeometry& melg)
0076     : TrapezoidalPlaneBounds(melg.hBottomEdge, melg.hTopEdge, melg.apothem, 0.5 * melg.thickness()),
0077       theWireTopology(nullptr),
0078       theStripTopology(nullptr),
0079       hBottomEdge(melg.hBottomEdge),
0080       hTopEdge(melg.hTopEdge),
0081       apothem(melg.apothem),
0082       chamberType(melg.chamberType) {
0083   // CSCStripTopology is abstract, so need clone()
0084   if (melg.theStripTopology)
0085     theStripTopology = melg.theStripTopology->clone();
0086   // CSCWireTopology is concrete, so direct copy
0087   if (melg.theWireTopology)
0088     theWireTopology = new CSCWireTopology(*(melg.theWireTopology));
0089 }
0090
0091 CSCLayerGeometry& CSCLayerGeometry::operator=(const CSCLayerGeometry& melg) {
0092   if (&melg != this) {
0093     delete theStripTopology;
0094     if (melg.theStripTopology)
0095       theStripTopology = melg.theStripTopology->clone();
0096     else
0097       theStripTopology = nullptr;
0098
0099     delete theWireTopology;
0100     if (melg.theWireTopology)
0101       theWireTopology = new CSCWireTopology(*(melg.theWireTopology));
0102     else
0103       theWireTopology = nullptr;
0104
0105     hBottomEdge = melg.hBottomEdge;
0106     hTopEdge = melg.hTopEdge;
0107     apothem = melg.apothem;
0108   }
0109   return *this;
0110 }
0111
0112 CSCLayerGeometry::~CSCLayerGeometry() {
0113   LogTrace("CSCLayerGeometry|CSC") << myName << ": being destroyed, this=" << this
0114                                    << "\nDeleting theStripTopology=" << theStripTopology
0115                                    << " and theWireTopology=" << theWireTopology;
0116   delete theStripTopology;
0117   delete theWireTopology;
0118 }
0119
0120 LocalPoint CSCLayerGeometry::stripWireIntersection(int strip, float wire) const {
0121   // This allows _float_ wire no. so that we can calculate the
0122   // intersection of a strip with the mid point of a wire group
0123   // containing an even no. of wires (which is not an actual wire),
0124   // as well as for a group containing an odd no. of wires.
0125
0126   // Equation of wire and strip as straight lines in local xy
0127   // y = mx + c where m = tan(angle w.r.t. x axis)
0128   // At the intersection x = -(cs-cw)/(ms-mw)
0129   // At y=0, 0 = ms * xOfStrip(strip) + cs => cs = -ms*xOfStrip
0130   // At x=0, yOfWire(wire) = 0 + cw => cw = yOfWire
0131
0132   float ms = tan(stripAngle(strip));
0133   float mw = tan(wireAngle());
0134   float xs = xOfStrip(strip);
0135   float xi = (ms * xs + yOfWire(wire)) / (ms - mw);
0136   float yi = ms * (xi - xs);
0137
0138   return LocalPoint(xi, yi);
0139 }
0140
0141 LocalPoint CSCLayerGeometry::stripWireGroupIntersection(int strip, int wireGroup) const {
0142   // middleWire is only an actual wire for a group with an odd no. of wires
0143   float middleWire = middleWireOfGroup(wireGroup);
0144   return stripWireIntersection(strip, middleWire);
0145 }
0146
0147 float CSCLayerGeometry::stripAngle(int strip) const {
0148   // Cleverly subtly change meaning of stripAngle once more.
0149   // In TrapezoidalStripTopology it is angle measured
0150   // counter-clockwise from y axis.
0151   // In APTST and RST it is angle measured
0152   // clockwise from y axis.
0153   // Output of this function is measured counter-clockwise
0154   // from x-axis, so it is a conventional 2-dim azimuthal angle
0155   // in the (x,y) local coordinates
0156
0157   // We want angle at centre of strip (strip N covers
0158   // *float* range N-1 to N-epsilon)
0159
0160   return 0.5_pi - theStripTopology->stripAngle(strip - 0.5);
0161 }
0162
0163 LocalPoint CSCLayerGeometry::localCenterOfWireGroup(int wireGroup) const {
0164   // It can use CSCWireTopology::yOfWireGroup for y,
0165   // But x requires mixing with 'extent' of wire plane
0166
0167   // If the wires are NOT tilted, default to simple calculation...
0168   if (fabs(wireAngle()) < 1.E-6) {
0169     float y = yOfWireGroup(wireGroup);
0170     return LocalPoint(0., y);
0171   } else {
0172     // w is "wire" at the center of the wire group
0173     float w = middleWireOfGroup(wireGroup);
0174     std::vector<float> store = theWireTopology->wireValues(w);
0175     return LocalPoint(store[0], store[1]);
0176   }
0177 }
0178
0179 float CSCLayerGeometry::lengthOfWireGroup(int wireGroup) const {
0180   // Return length of 'wire' in the middle of the wire group
0181   float w = middleWireOfGroup(wireGroup);
0182   std::vector<float> store = theWireTopology->wireValues(w);
0183   return store[2];
0184 }
0185
0186 std::pair<LocalPoint, float> CSCLayerGeometry::possibleRecHitPosition(float s, int w1, int w2) const {
0187   LocalPoint sw1 = intersectionOfStripAndWire(s, w1);
0188   LocalPoint sw2 = intersectionOfStripAndWire(s, w2);
0189
0190   // Average the two points
0191   LocalPoint midpt((sw1.x() + sw2.x()) / 2., (sw1.y() + sw2.y()) / 2);
0192
0193   // Length of strip crossing this group of wires
0194   float length = sqrt((sw1.x() - sw2.x()) * (sw1.x() - sw2.x()) + (sw1.y() - sw2.y()) * (sw1.y() - sw2.y()));
0195
0196   return std::pair<LocalPoint, float>(midpt, length);
0197 }
0198
0199 LocalPoint CSCLayerGeometry::intersectionOfStripAndWire(float s, int w) const {
0200   std::pair<float, float> pw = theWireTopology->equationOfWire(static_cast<float>(w));
0201   std::pair<float, float> ps = theStripTopology->equationOfStrip(s);
0202   LocalPoint sw = intersectionOfTwoLines(ps, pw);
0203
0204   // If point falls outside wire plane, at extremes in local y,
0205   // replace its y by that of appropriate edge of wire plane
0206   if (!(theWireTopology->insideYOfWirePlane(sw.y()))) {
0207     float y = theWireTopology->restrictToYOfWirePlane(sw.y());
0208     // and adjust x to match new y
0209     float x = sw.x() + (y - sw.y()) * tan(theStripTopology->stripAngle(s));
0210     sw = LocalPoint(x, y);
0211   }
0212
0213   return sw;
0214 }
0215
0216 LocalPoint CSCLayerGeometry::intersectionOfTwoLines(std::pair<float, float> p1, std::pair<float, float> p2) const {
0217   // Calculate the point of intersection of two straight lines (in 2-dim)
0218   // input arguments are pair(m1,c1) and pair(m2,c2) where m=slope, c=intercept (y=mx+c)
0219   // BEWARE! Do not call with m1 = m2 ! No trapping !
0220
0221   float m1 = p1.first;
0222   float c1 = p1.second;
0223   float m2 = p2.first;
0224   float c2 = p2.second;
0225   float x = (c2 - c1) / (m1 - m2);
0226   float y = (m1 * c2 - m2 * c1) / (m1 - m2);
0227   return LocalPoint(x, y);
0228 }
0229
0230 LocalError CSCLayerGeometry::localError(int strip, float sigmaStrip, float sigmaWire) const {
0231   // Input sigmas are expected to be in _distance units_
0232   // - uncertainty in strip measurement (typically from Gatti fit, value is in local x units)
0233   // - uncertainty in wire measurement (along direction perpendicular to wires)
0234
0235   float wangle = this->wireAngle();
0236   float strangle = this->stripAngle(strip);
0237
0238   float sinAngdif = sin(strangle - wangle);
0239   float sinAngdif2 = sinAngdif * sinAngdif;
0240
0241   float du = sigmaStrip / sin(strangle);  // sigmaStrip is just x-component of strip error
0242   float dv = sigmaWire;
0243
0244   // The notation is
0245   // wsins = wire resol  * sin(strip angle)
0246   // wcoss = wire resol  * cos(strip angle)
0247   // ssinw = strip resol * sin(wire angle)
0248   // scosw = strip resol * cos(wire angle)
0249
0250   float wsins = dv * sin(strangle);
0251   float wcoss = dv * cos(strangle);
0252   float ssinw = du * sin(wangle);
0253   float scosw = du * cos(wangle);
0254
0255   float dx2 = (scosw * scosw + wcoss * wcoss) / sinAngdif2;
0256   float dy2 = (ssinw * ssinw + wsins * wsins) / sinAngdif2;
0257   float dxy = (scosw * ssinw + wcoss * wsins) / sinAngdif2;
0258
0259   return LocalError(dx2, dxy, dy2);
0260 }
0261
0262 bool CSCLayerGeometry::inside(const Local3DPoint& lp) const {
0263   bool result = false;
0264   const float epsilon = 1.e-06;
0265   if (fabs(lp.z()) < thickness() / 2.) {  // thickness of TPB is that of gas layer
0266     std::pair<float, float> ylims = yLimitsOfStripPlane();
0267     if ((lp.y() > ylims.first) && (lp.y() < ylims.second)) {
0268       // 'strip' returns float value between 0. and float(Nstrips) and value outside
0269       // is set to 0. or float(Nstrips)... add a conservative precision of 'epsilon'
0270       if ((theStripTopology->strip(lp) > epsilon) && (theStripTopology->strip(lp) < (numberOfStrips() - epsilon)))
0271         result = true;
0272     }
0273   }
0274   return result;
0275 }
0276
0277 bool CSCLayerGeometry::inside(const Local2DPoint& lp) const {
0278   LocalPoint lp2(lp.x(), lp.y(), 0.);
0279   return inside(lp2);
0280 }
0281
0282 bool CSCLayerGeometry::inside(const Local3DPoint& lp, const LocalError& le, float scale) const {
0283   // Effectively consider that the LocalError components extend the area which is acceptable.
0284   // Form a little box centered on the point, with x, y diameters defined by the errors
0285   // and require that ALL four corners of the box fall outside the strip region for failure
0286
0287   // Note that LocalError is 2-dim x,y and doesn't supply a z error
0288   float deltaX = scale * sqrt(le.xx());
0289   float deltaY = scale * sqrt(le.yy());
0290
0291   LocalPoint lp1(lp.x() - deltaX, lp.y() - deltaY, lp.z());
0292   LocalPoint lp2(lp.x() - deltaX, lp.y() + deltaY, lp.z());
0293   LocalPoint lp3(lp.x() + deltaX, lp.y() + deltaY, lp.z());
0294   LocalPoint lp4(lp.x() + deltaX, lp.y() - deltaY, lp.z());
0295
0296   return (inside(lp1) || inside(lp2) || inside(lp3) || inside(lp4));
0297 }
0298
0299 void CSCLayerGeometry::setTopology(CSCStripTopology* newTopology) {
0300   delete theStripTopology;
0301   theStripTopology = newTopology;
0302 }
0303
0304 std::ostream& operator<<(std::ostream& stream, const CSCLayerGeometry& lg) {
0305   stream << "LayerGeometry " << std::endl
0306          << "------------- " << std::endl
0307          << "numberOfStrips               " << lg.numberOfStrips() << std::endl
0308          << "numberOfWires                " << lg.numberOfWires() << std::endl
0309          << "numberOfWireGroups           " << lg.numberOfWireGroups() << std::endl
0310          << "wireAngle  (rad)             " << lg.wireAngle()
0311          << std::endl
0312          //         << "wireAngle  (deg)      " << lg.theWireAngle << std::endl
0313          //         << "sin(wireAngle)        " << lg.theWireSin << std::endl
0314          //         << "cos(wireAngle)        " << lg.theWireCos << std::endl
0315          << "wirePitch                    " << lg.wirePitch() << std::endl
0316          << "stripPitch                   " << lg.stripPitch()
0317          << std::endl
0318          //         << "numberOfWiresPerGroup " << lg.theNumberOfWiresPerGroup << std::endl
0319          //         << "numberOfWiresInLastGroup " << lg.theNumberOfWiresInLastGroup << std::endl
0320          //         << "wireOffset            " << lg.theWireOffset << std::endl
0321          //         << "whereStripsMeet       " << lg.whereStripsMeet << std::endl;
0322          << "hBottomEdge                  " << lg.hBottomEdge << std::endl
0323          << "hTopEdge                     " << lg.hTopEdge << std::endl
0324          << "apothem                      " << lg.apothem << std::endl
0325          << "length (should be 2xapothem) " << lg.length() << std::endl
0326          << "thickness                    " << lg.thickness() << std::endl;
0327   return stream;
0328 }