test/stubs/CSCGeometryOfWires.cc

0001 // LAST UPDATED 16.09.2009 ptc
0002
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0008 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0009 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0010 #include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
0011 #include "DataFormats/GeometryVector/interface/Pi.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0013
0014 #include <string>
0015 #include <iomanip>  // for setw() etc.
0016 #include <vector>
0017
0018 class CSCGeometryOfWires : public edm::one::EDAnalyzer<> {
0019 public:
0020   explicit CSCGeometryOfWires(const edm::ParameterSet&);
0021   ~CSCGeometryOfWires() override = default;
0022
0023   void beginJob() override {}
0024   void analyze(edm::Event const&, edm::EventSetup const&) override;
0025   void endJob() override {}
0026
0027   const std::string& myName() { return myName_; }
0028
0029 private:
0030   const int dashedLineWidth_;
0031   const std::string dashedLine_;
0032   const std::string myName_;
0033   const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> tokGeom_;
0034 };
0035
0036 CSCGeometryOfWires::CSCGeometryOfWires(const edm::ParameterSet& iConfig)
0037     : dashedLineWidth_(101),
0038       dashedLine_(std::string(dashedLineWidth_, '-')),
0039       myName_("CSCGeometryOfWires"),
0040       tokGeom_(esConsumes()) {}
0041
0042 void CSCGeometryOfWires::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0043   const double dPi = Geom::pi();
0044   const double radToDeg = 180. / dPi;
0045
0046   std::cout << myName() << ": Analyzer..." << std::endl;
0047   std::cout << "start " << dashedLine_ << std::endl;
0048
0049   const edm::ESHandle<CSCGeometry>& pDD = iSetup.getHandle(tokGeom_);
0050   std::cout << " Geometry node for CSCGeometry is  " << &(*pDD) << std::endl;
0051   std::cout << " I have " << pDD->detTypes().size() << " detTypes" << std::endl;
0052   std::cout << " I have " << pDD->detUnits().size() << " detUnits" << std::endl;
0053   std::cout << " I have " << pDD->dets().size() << " dets" << std::endl;
0054   std::cout << " I have " << pDD->layers().size() << " layers" << std::endl;
0055   std::cout << " I have " << pDD->chambers().size() << " chambers" << std::endl;
0056
0057   std::cout << myName() << ": Begin iteration over geometry..." << std::endl;
0058   std::cout << "iter " << dashedLine_ << std::endl;
0059
0060   int icount = 0;
0061
0062   // Check the DetUnits
0063   for (const auto& it : pDD->detUnits()) {
0064     // Do we really have a CSC layer?
0065
0066     auto layer = dynamic_cast<CSCLayer const*>(it);
0067
0068     if (layer) {
0069       ++icount;
0070
0071       // What's its DetId?
0072
0073       DetId detId = layer->geographicalId();
0074       int id = detId();  // or detId.rawId()
0075
0076       //    std::cout << "GeomDetUnit is of type " << detId.det() << " and raw id = " << id << std::endl;
0077
0078       std::cout << "\n"
0079                 << "Parameters of layer# " << icount << " id= " << id << " = " << std::oct << id << std::dec
0080                 << " (octal) "
0081                 << "   E" << CSCDetId::endcap(id) << " S" << CSCDetId::station(id) << " R" << CSCDetId::ring(id) << " C"
0082                 << CSCDetId::chamber(id) << " L" << CSCDetId::layer(id) << " are:" << std::endl;
0083
0084       const CSCLayerGeometry* geom = layer->geometry();
0085       std::cout << *geom;
0086
0087       const CSCStripTopology* cst = geom->topology();
0088       std::cout << "\n" << *cst;
0089
0090       // What's its surface?
0091       // The surface knows how to transform local <-> global
0092
0093       const Surface& bSurface = layer->surface();
0094
0095       // Check global coordinates of centre of CSCLayer, and how
0096       // local z direction relates to global z direction
0097
0098       LocalPoint lCentre(0., 0., 0.);
0099       GlobalPoint gCentre = bSurface.toGlobal(lCentre);
0100
0101       LocalPoint lCentre1(0., 0., -1.);
0102       GlobalPoint gCentre1 = bSurface.toGlobal(lCentre1);
0103
0104       LocalPoint lCentre2(0., 0., 1.);
0105       GlobalPoint gCentre2 = bSurface.toGlobal(lCentre2);
0106
0107       double gx = gCentre.x();
0108       double gy = gCentre.y();
0109       double gz = gCentre.z();
0110       double gz1 = gCentre1.z();
0111       double gz2 = gCentre2.z();
0112       if (fabs(gx) < 1.e-06)
0113         gx = 0.;
0114       if (fabs(gy) < 1.e-06)
0115         gy = 0.;
0116       if (fabs(gz) < 1.e-06)
0117         gz = 0.;
0118       if (fabs(gz1) < 1.e-06)
0119         gz1 = 0.;
0120       if (fabs(gz2) < 1.e-06)
0121         gz2 = 0.;
0122
0123       // subverted by GeometryVector/Phi.h enforcing its range convention!
0124       // Either a) use a separate local double before scaling...
0125       //        double cphi = gCentre.phi();
0126       //        double cphiDeg = cphi * radToDeg;
0127       // or b) use Phi's degree conversion...
0128       double cphiDeg = gCentre.phi().degrees();
0129
0130       // I want to display in range 0 to 360
0131
0132       // Handle some occasional ugly precision problems around zero
0133       if (fabs(cphiDeg) < 1.e-06) {
0134         cphiDeg = 0.;
0135       } else if (cphiDeg < 0.) {
0136         cphiDeg += 360.;
0137       } else if (cphiDeg >= 360.) {
0138         std::cout << "WARNING: resetting phi= " << cphiDeg << " to zero." << std::endl;
0139         cphiDeg = 0.;
0140       }
0141
0142       //        int iphiDeg = static_cast<int>( cphiDeg );
0143       //    std::cout << "phi(0,0,0) = " << iphiDeg << " degrees" << std::endl;
0144
0145       int nStrips = geom->numberOfStrips();
0146       //        std::cout << std::setw( 4 ) << nStrips;
0147
0148       double cstrip1 = layer->centerOfStrip(1).phi();
0149       double cstripN = layer->centerOfStrip(nStrips).phi();
0150
0151       double phiwid = geom->stripPhiPitch();
0152       double stripwid = geom->stripPitch();
0153       double stripoff = geom->stripOffset();
0154       double phidif = fabs(cstrip1 - cstripN);
0155
0156       // May have one strip at 180-epsilon and other at -180+epsilon
0157       // If so the raw difference is 360-(phi extent of chamber)
0158       // Want to reset that to (phi extent of chamber):
0159       if (phidif > dPi)
0160         phidif = fabs(phidif - 2. * dPi);
0161       double phiwid_check = phidif / double(nStrips - 1);
0162
0163       // Clean up some stupid floating decimal aesthetics
0164       cstrip1 = cstrip1 * radToDeg;
0165       if (fabs(cstrip1) < 1.e-06)
0166         cstrip1 = 0.;
0167       else if (cstrip1 < 0.)
0168         cstrip1 += 360.;
0169
0170       cstripN = cstripN * radToDeg;
0171       if (fabs(cstripN) < 1.e-06)
0172         cstripN = 0.;
0173       else if (cstripN < 0.)
0174         cstripN += 360.;
0175
0176       if (fabs(stripoff) < 1.e-06)
0177         stripoff = 0.;
0178
0179       // Layer geometry:  layer corner phi's... OF ALUMINUM FRAME
0180
0181       std::array<const float, 4> const& parameters = geom->parameters();
0182       // these parameters are half-lengths, due to GEANT
0183       float hBottomEdge = parameters[0];
0184       float hTopEdge = parameters[1];
0185       float hThickness = parameters[2];
0186       float hApothem = parameters[3];
0187
0188       // first the face nearest the interaction
0189       // get the other face by using positive hThickness
0190       LocalPoint upperRightLocal(hTopEdge, hApothem, -hThickness);
0191       LocalPoint upperLeftLocal(-hTopEdge, hApothem, -hThickness);
0192       LocalPoint lowerRightLocal(hBottomEdge, -hApothem, -hThickness);
0193       LocalPoint lowerLeftLocal(-hBottomEdge, -hApothem, -hThickness);
0194
0195       LocalPoint upperEdgeOnY(0., hApothem);
0196       LocalPoint lowerEdgeOnY(0., -hApothem);
0197       LocalPoint leftEdgeOnX(-(hTopEdge + hBottomEdge) / 2., 0.);
0198       LocalPoint rightEdgeOnX((hTopEdge + hBottomEdge) / 2., 0.);
0199
0200       GlobalPoint upperRightGlobal = bSurface.toGlobal(upperRightLocal);
0201       GlobalPoint upperLeftGlobal = bSurface.toGlobal(upperLeftLocal);
0202       GlobalPoint lowerRightGlobal = bSurface.toGlobal(lowerRightLocal);
0203       GlobalPoint lowerLeftGlobal = bSurface.toGlobal(lowerLeftLocal);
0204
0205       float uRGp = upperRightGlobal.phi().degrees();
0206       float uLGp = upperLeftGlobal.phi().degrees();
0207       float lRGp = lowerRightGlobal.phi().degrees();
0208       float lLGp = lowerLeftGlobal.phi().degrees();
0209
0210       std::cout << "\nCHAMBER FRAME corners in local coordinates: \n UR " << upperRightLocal << "\n UL "
0211                 << upperLeftLocal << "\n LR " << lowerRightLocal << "\n LL " << lowerLeftLocal << std::endl;
0212
0213       std::cout << "CHAMBER FRAME corners in global coords: \n UR " << upperRightGlobal << "\n UL " << upperLeftGlobal
0214                 << "\n LR " << lowerRightGlobal << "\n LL " << lowerLeftGlobal << "\n   phi: UR " << uRGp << " UL "
0215                 << uLGp << " LR " << lRGp << " LL " << lLGp << std::endl;
0216
0217       float ang = nStrips * phiwid;
0218       float ctoi = cst->centreToIntersection();
0219
0220       std::cout << "\nSTRIP PLANE:" << std::endl;
0221       std::cout << "============" << std::endl;
0222       std::cout << "centreToIntersection, R = " << ctoi << std::endl;
0223       std::cout << "local y of centre of strip plane, yOff = " << cst->yCentreOfStripPlane() << std::endl;
0224       std::cout << "originToIntersection, R-yOff = " << cst->originToIntersection() << std::endl;
0225       std::cout << "extent of strip plane in local y = " << cst->yExtentOfStripPlane() << std::endl;
0226       std::cout << "no. of strips = " << nStrips << std::endl;
0227
0228       // Local coordinates of ends of strips
0229       MeasurementPoint mp1m(0.5, -0.5);                                // strip 1, -y
0230       MeasurementPoint mp1p(0.5, 0.5);                                 // strip 1, +y
0231       MeasurementPoint mp2m(static_cast<float>(nStrips) - 0.5, -0.5);  // strip N, -y
0232       MeasurementPoint mp2p(static_cast<float>(nStrips) - 0.5, 0.5);   // strip N, +y
0233       LocalPoint lp1m = geom->topology()->localPosition(mp1m);
0234       LocalPoint lp1p = geom->topology()->localPosition(mp1p);
0235       LocalPoint lp2m = geom->topology()->localPosition(mp2m);
0236       LocalPoint lp2p = geom->topology()->localPosition(mp2p);
0237
0238       std::cout << "Strip 1 local coords " << lp1m << " " << lp1p << std::endl;
0239       std::cout << "Strip " << nStrips << " local coords " << lp2m << " " << lp2p << std::endl;
0240
0241       std::cout << "angular width of one strip = " << phiwid << " rads " << std::endl;
0242       std::cout << "angle subtended by layer, A = nstrips x stripPhiPitch = " << ang << " rads = " << ang * radToDeg
0243                 << " deg" << std::endl;
0244       std::cout << "phi (clockwise from local y axis) of one edge = " << cst->phiOfOneEdge() << " rads " << std::endl;
0245
0246       std::cout << "phi width check: (centre strip N - centre strip 1)/(nstrips-1) = " << phiwid_check << std::endl;
0247
0248       std::cout << "Check how well approximations work: " << std::endl;
0249       std::cout << "[T+B = " << hTopEdge + hBottomEdge << "] compared with [2R*tan(A/2) = " << 2. * ctoi * tan(ang / 2.)
0250                 << "]" << std::endl;
0251       std::cout << "[T-B = " << hTopEdge - hBottomEdge
0252                 << "] compared with [2a*tan(A/2) = " << 2. * hApothem * tan(ang / 2.) << "]" << std::endl;
0253       std::cout << "[R = " << ctoi
0254                 << "] compared with [0.5*(T+B)/tan(A/2) = " << 0.5 * (hTopEdge + hBottomEdge) / tan(ang / 2.) << "]"
0255                 << std::endl;
0256
0257       std::cout << "Approximations to where strips intersect: " << std::endl;
0258       std::cout << "RST: match y=0, oi = (hT+hB)/2 / tan(A/2) = "
0259                 << 0.5 * (hTopEdge + hBottomEdge) / tan(0.5 * nStrips * phiwid) << std::endl;
0260       std::cout << "RST: match top, oi = hT / tan(A/2) = " << hTopEdge / tan(0.5 * nStrips * phiwid) << std::endl;
0261       std::cout << "TST: oi = hA*(hT+hB)/(hT-hB) = " << hApothem * (hTopEdge + hBottomEdge) / (hTopEdge - hBottomEdge)
0262                 << std::endl;
0263
0264       std::cout << "\nstrip offset = " << stripoff << std::endl;
0265
0266       std::cout << "\nlocal(0,0,-1) = global " << gCentre1 << std::endl;
0267       std::cout << "local(0,0)    = global " << gCentre << std::endl;
0268       std::cout << "local(0,0,+1) = global " << gCentre2 << std::endl;
0269
0270       // CSCLG::stripAngle(int strip)
0271       std::cout << "CSCLG Angle of strip 1 = " << geom->stripAngle(1) * radToDeg << " deg " << std::endl;
0272       std::cout << "CSCLG Angle of strip " << nStrips / 2 << " = " << geom->stripAngle(nStrips / 2) * radToDeg
0273                 << " deg " << std::endl;
0274       std::cout << "CSCLG Angle of strip " << nStrips << " = " << geom->stripAngle(nStrips) * radToDeg << " deg "
0275                 << std::endl;
0276
0277       // CSCST::stripAngle(float strip) Yes this one's float the CSCLG is int
0278       std::cout << "CSCST Angle of centre of strip 1 = " << cst->stripAngle(0.5) * radToDeg << " deg " << std::endl;
0279       std::cout << "CSCST Angle of centre of strip " << nStrips / 2 << " = "
0280                 << cst->stripAngle(nStrips / 2 - 0.5) * radToDeg << " deg " << std::endl;
0281       std::cout << "CSCST Angle of centre of strip " << nStrips << " = " << cst->stripAngle(nStrips - 0.5) * radToDeg
0282                 << " deg " << std::endl;
0283
0284       std::cout << "Local x of strip 1 on x axis = " << geom->xOfStrip(1, 0.) << std::endl;
0285       std::cout << "Local x of strip " << nStrips / 2 << " on x axis = " << geom->xOfStrip(nStrips / 2, 0.)
0286                 << std::endl;
0287       std::cout << "Local x of strip " << nStrips << " on x axis = " << geom->xOfStrip(nStrips, 0.) << std::endl;
0288
0289       std::cout << "Local x of strip 1 at upper edge = " << geom->xOfStrip(1, hApothem) << std::endl;
0290       std::cout << "Local x of strip " << nStrips / 2 << " at upper edge = " << geom->xOfStrip(nStrips / 2, hApothem)
0291                 << std::endl;
0292       std::cout << "Local x of strip " << nStrips << " at upper edge = " << geom->xOfStrip(nStrips, hApothem)
0293                 << std::endl;
0294
0295       std::cout << "Local x of strip 1 at lower edge = " << geom->xOfStrip(1, -hApothem) << std::endl;
0296       std::cout << "Local x of strip " << nStrips / 2 << " at lower edge = " << geom->xOfStrip(nStrips / 2, -hApothem)
0297                 << std::endl;
0298       std::cout << "Local x of strip " << nStrips << " at lower edge = " << geom->xOfStrip(nStrips, -hApothem)
0299                 << std::endl;
0300
0301       std::cout << "Strip width           = " << stripwid << std::endl;
0302       std::cout << "Strip pitch at middle = " << geom->stripPitch() << std::endl;
0303       std::cout << "Strip pitch at (0,0)  = " << geom->stripPitch(lCentre) << std::endl;
0304
0305       std::cout << "Strip pitch at UR     = " << geom->stripPitch(upperRightLocal) << std::endl;
0306       std::cout << "Strip pitch at UL     = " << geom->stripPitch(upperLeftLocal) << std::endl;
0307       std::cout << "Strip pitch at LL     = " << geom->stripPitch(lowerLeftLocal) << std::endl;
0308       std::cout << "Strip pitch at LR     = " << geom->stripPitch(lowerRightLocal) << std::endl;
0309
0310       std::cout << "Strip pitch at upper edge on y axis = " << geom->stripPitch(upperEdgeOnY) << std::endl;
0311       std::cout << "Strip pitch at lower edge on y axis = " << geom->stripPitch(lowerEdgeOnY) << std::endl;
0312       std::cout << "Strip pitch at left edge on x axis  = " << geom->stripPitch(leftEdgeOnX) << std::endl;
0313       std::cout << "Strip pitch at right edge on x axis = " << geom->stripPitch(rightEdgeOnX) << std::endl;
0314
0315       // Check input to nearestStrip()
0316       std::cout << "Strip units for (0,0) =                " << cst->strip(lCentre) << std::endl;
0317       std::cout << "Strip units for upper edge on y axis = " << cst->strip(upperEdgeOnY) << std::endl;
0318       std::cout << "Strip units for lower edge on y axis = " << cst->strip(lowerEdgeOnY) << std::endl;
0319       std::cout << "Strip units for left edge on x axis  = " << cst->strip(leftEdgeOnX) << std::endl;
0320       std::cout << "Strip units for right edge on x axis = " << cst->strip(rightEdgeOnX) << std::endl;
0321
0322       std::cout << "Nearest strip to (0,0) =                " << geom->nearestStrip(lCentre) << std::endl;
0323       std::cout << "Nearest strip to upper edge on y axis = " << geom->nearestStrip(upperEdgeOnY) << std::endl;
0324       std::cout << "Nearest strip to lower edge on y axis = " << geom->nearestStrip(lowerEdgeOnY) << std::endl;
0325       std::cout << "Nearest strip to left edge on x axis  = " << geom->nearestStrip(leftEdgeOnX) << std::endl;
0326       std::cout << "Nearest strip to right edge on x axis = " << geom->nearestStrip(rightEdgeOnX) << std::endl;
0327
0328       int iNUR = geom->nearestStrip(upperRightLocal);
0329       int iNUL = geom->nearestStrip(upperLeftLocal);
0330       int iNLR = geom->nearestStrip(lowerRightLocal);
0331       int iNLL = geom->nearestStrip(lowerLeftLocal);
0332
0333       int jNUR = geom->nearestWire(upperRightLocal);
0334       int jNUL = geom->nearestWire(upperLeftLocal);
0335       int jNLR = geom->nearestWire(lowerRightLocal);
0336       int jNLL = geom->nearestWire(lowerLeftLocal);
0337
0338       std::cout << "Calculated no. of strips across top    = " << (iNUR - iNUL + 1) << std::endl;
0339       std::cout << "Calculated no. of strips across bottom = " << (iNLR - iNLL + 1) << std::endl;
0340
0341       std::cout << "Nearest strip, wire to UR = " << iNUR << ", " << jNUR << std::endl;
0342       std::cout << "Nearest strip, wire to UL = " << iNUL << ", " << jNUL << std::endl;
0343       std::cout << "Nearest strip, wire to LR = " << iNLR << ", " << jNLR << std::endl;
0344       std::cout << "Nearest strip, wire to LL = " << iNLL << ", " << jNLL << std::endl;
0345
0346       std::cout << "yOfWire(" << jNUR << " , +hTopEdge ) = " << geom->yOfWire(static_cast<float>(jNUR), hTopEdge)
0347                 << std::endl;
0348       std::cout << "yOfWire(" << jNUL << " , -hTopEdge ) = " << geom->yOfWire(static_cast<float>(jNUL), -hTopEdge)
0349                 << std::endl;
0350       std::cout << "yOfWire(" << jNLR << " , hBottomEdge ) = " << geom->yOfWire(static_cast<float>(jNLR), hBottomEdge)
0351                 << std::endl;
0352       std::cout << "yOfWire(" << jNLL << " , -hBottomEdge ) = " << geom->yOfWire(static_cast<float>(jNLL), -hBottomEdge)
0353                 << std::endl;
0354
0355       std::cout << "Examine global phi of strip at top and bottom of chamber frame:" << std::endl;
0356       float phi_1_c = (layer->centerOfStrip(1)).phi();
0357       float phi_n_c = (layer->centerOfStrip(nStrips)).phi();
0358       float phi_c_c = (layer->centerOfStrip(nStrips / 2)).phi();
0359       float x_1_t = geom->xOfStrip(1, hApothem);             // x of strip 1 at top edge
0360       float x_1_b = geom->xOfStrip(1, -hApothem);            // x of strip 1 at bottom edge
0361       float x_n_t = geom->xOfStrip(nStrips, hApothem);       // x of strip n at top edge
0362       float x_n_b = geom->xOfStrip(nStrips, -hApothem);      // x of strip n at bottom edge
0363       float x_c_t = geom->xOfStrip(nStrips / 2, hApothem);   // x of strip n/2 at top edge
0364       float x_c_b = geom->xOfStrip(nStrips / 2, -hApothem);  // x of strip n/2 at bottom edge
0365       GlobalPoint g_1_t = layer->toGlobal(LocalPoint(x_1_t, hApothem, 0.));
0366       GlobalPoint g_1_b = layer->toGlobal(LocalPoint(x_1_b, -hApothem, 0.));
0367       GlobalPoint g_n_t = layer->toGlobal(LocalPoint(x_n_t, hApothem, 0.));
0368       GlobalPoint g_n_b = layer->toGlobal(LocalPoint(x_n_b, -hApothem, 0.));
0369       GlobalPoint g_c_t = layer->toGlobal(LocalPoint(x_c_t, 0., 0.));
0370       GlobalPoint g_c_b = layer->toGlobal(LocalPoint(x_c_b, 0., 0.));
0371       float phi_1_t = g_1_t.phi();
0372       float phi_1_b = g_1_b.phi();
0373       float phi_n_t = g_n_t.phi();
0374       float phi_n_b = g_n_b.phi();
0375       float phi_c_t = g_c_t.phi();
0376       float phi_c_b = g_c_b.phi();
0377       std::cout << " strip  1 top: " << phi_1_t << " centre: " << phi_1_c << " bottom: " << phi_1_b
0378                 << " top-bottom: " << phi_1_t - phi_1_b << std::endl;
0379       std::cout << " strip " << nStrips / 2 << " top: " << phi_c_t << " centre: " << phi_c_c << " bottom: " << phi_c_b
0380                 << " top-bottom: " << phi_c_t - phi_c_b << std::endl;
0381       std::cout << " strip " << nStrips << " top: " << phi_n_t << " centre: " << phi_n_c << " bottom: " << phi_n_b
0382                 << " top-bottom: " << phi_n_t - phi_n_b << std::endl;
0383
0384       int nwg = geom->numberOfWireGroups();
0385       int nw = geom->numberOfWires();
0386
0387       std::cout << "\nWIRE PLANE:" << std::endl;
0388       std::cout << "===========" << std::endl;
0389       std::cout << "wireSpacing = " << geom->wirePitch() << " cm " << std::endl;
0390       std::cout << "wireAngle = " << geom->wireAngle() << " rads " << std::endl;
0391       std::cout << "no. of wires = " << nw << std::endl;
0392       std::cout << "no. of wire groups = " << nwg << std::endl;
0393       std::cout << "no. of wires in wg 1 = " << geom->numberOfWiresPerGroup(1) << std::endl;
0394       std::cout << "no. of wires in wg " << nwg << " = " << geom->numberOfWiresPerGroup(nwg) << std::endl;
0395       std::cout << "wire group containing wire 1 = " << geom->wireGroup(1) << std::endl;
0396       std::cout << "wire group containing wire " << nw << " = " << geom->wireGroup(nw) << std::endl;
0397       std::cout << "length of wg 1 = " << geom->lengthOfWireGroup(1) << std::endl;
0398       std::cout << "length of wg " << nwg << " = " << geom->lengthOfWireGroup(nwg) << std::endl;
0399       std::cout << "middle wire of wg 1 = " << geom->middleWireOfGroup(1) << std::endl;
0400       std::cout << "middle wire of wg " << nwg << " = " << geom->middleWireOfGroup(nwg) << std::endl;
0401       std::cout << "y of wg 1 at x=0 is " << geom->yOfWireGroup(1) << std::endl;
0402       std::cout << "y of wg " << nwg << " at x=0 is " << geom->yOfWireGroup(nwg) << std::endl;
0403       std::cout << "centre of wg 1 is " << geom->localCenterOfWireGroup(1) << std::endl;
0404       std::cout << "centre of wg " << nwg << " is " << geom->localCenterOfWireGroup(nwg) << std::endl;
0405       std::cout << "y of wire 1 at x=0 = " << geom->wireTopology()->yOfWire(1., 0.) << std::endl;
0406       std::cout << "y of wire " << nw << " at x=0 = " << geom->wireTopology()->yOfWire(static_cast<float>(nw), 0.)
0407                 << std::endl;
0408
0409       std::cout << "narrow width of wire plane = " << geom->wireTopology()->narrowWidthOfPlane() << std::endl;
0410       std::cout << "wide width   of wire plane = " << geom->wireTopology()->wideWidthOfPlane() << std::endl;
0411       std::cout << "y extent     of wire plane = " << geom->wireTopology()->lengthOfPlane() << std::endl;
0412       std::cout << "perp extent  of wire plane = " << geom->wireTopology()->extentOfWirePlane() << std::endl;
0413
0414       std::pair<LocalPoint, LocalPoint> lw1 = geom->wireTopology()->wireEnds(geom->middleWireOfGroup(1));
0415       std::cout << "local coords of ends of central wire of wire group 1 = "
0416                 << "(" << lw1.first.x() << ", " << lw1.first.y() << "), "
0417                 << "(" << lw1.second.x() << ", " << lw1.second.y() << ")" << std::endl;
0418
0419       std::pair<LocalPoint, LocalPoint> lw2 = geom->wireTopology()->wireEnds(geom->middleWireOfGroup(nwg));
0420       std::cout << "local coords of ends of central wire of wire group " << nwg << " = "
0421                 << "(" << lw2.first.x() << ", " << lw2.first.y() << "), "
0422                 << "(" << lw2.second.x() << ", " << lw2.second.y() << ")" << std::endl;
0423
0424       // Check idToDetUnit
0425       const GeomDetUnit* gdu = pDD->idToDetUnit(detId);
0426       assert(gdu == layer);
0427       // Check idToDet
0428       const GeomDet* gd = pDD->idToDet(detId);
0429       assert(gd == layer);
0430     } else {
0431       std::cout << "Could not dynamic_cast to a CSCLayer* " << std::endl;
0432     }
0433   }
0434   std::cout << dashedLine_ << " end" << std::endl;
0435 }
0436
0437 //define this as a plug-in
0438 DEFINE_FWK_MODULE(CSCGeometryOfWires);