test/stubs/CSCGeometryAnalyzer.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0006 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0007 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0008 #include "DataFormats/GeometryVector/interface/Pi.h"
0009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0010
0011 #include <string>
0012 #include <iomanip>  // for setw() etc.
0013 #include <vector>
0014
0015 class CSCGeometryAnalyzer : public edm::one::EDAnalyzer<> {
0016 public:
0017   explicit CSCGeometryAnalyzer(const edm::ParameterSet&);
0018   ~CSCGeometryAnalyzer() override;
0019
0020   void beginJob() override {}
0021   void analyze(edm::Event const&, edm::EventSetup const&) override;
0022   void endJob() override {}
0023
0024   const std::string& myName() { return myName_; }
0025
0026 private:
0027   const int dashedLineWidth_;
0028   const std::string dashedLine_;
0029   const std::string myName_;
0030   const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> ddToken_;
0031 };
0032
0033 CSCGeometryAnalyzer::CSCGeometryAnalyzer(const edm::ParameterSet& iConfig)
0034     : dashedLineWidth_(194),
0035       dashedLine_(std::string(dashedLineWidth_, '-')),
0036       myName_("CSCGeometryAnalyzer"),
0037       ddToken_(esConsumes<CSCGeometry, MuonGeometryRecord>()) {}
0038
0039 CSCGeometryAnalyzer::~CSCGeometryAnalyzer() {}
0040
0041 void CSCGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0042   const double dPi = Geom::pi();
0043   const double radToDeg = 180. / dPi;
0044
0045   std::cout << myName() << ": Analyzer..." << std::endl;
0046   std::cout << "start " << dashedLine_ << std::endl;
0047   std::cout << "pi = " << dPi << ", radToDeg = " << radToDeg << std::endl;
0048
0049   const edm::ESTransientHandle<CSCGeometry> pDD = iSetup.getTransientHandle(ddToken_);
0050   std::cout << " Geometry node for CSCGeom 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 << " Ganged strips? " << pDD->gangedStrips() << std::endl;
0058   std::cout << " Wires only?    " << pDD->wiresOnly() << std::endl;
0059   std::cout << " Real wire geometry? " << pDD->realWireGeometry() << std::endl;
0060   std::cout << " Use offsets to coi? " << pDD->centreTIOffsets() << std::endl;
0061
0062   std::cout << myName() << ": Begin iteration over geometry..." << std::endl;
0063   std::cout << "iter " << dashedLine_ << std::endl;
0064
0065   std::cout << "\n  #     id(dec)      id(oct)                   "
0066                "  g(x=0)   g(y=0)   g(z=0)  g(z=-1)  g(z=+1)  Ns "
0067                "  phi(0)  phi(s1)  phi(sN)    dphi    dphi'      ds     off"
0068                "       uR       uL       lR       lL"
0069             << std::endl;
0070
0071   int icount = 0;
0072
0073   // Check the DetUnits
0074   for (const auto& it : pDD->detUnits()) {
0075     // Do we really have a CSC layer?
0076
0077     auto layer = dynamic_cast<CSCLayer const*>(it);
0078
0079     if (layer) {
0080       ++icount;
0081
0082       // What's its DetId?
0083
0084       DetId detId = layer->geographicalId();
0085       int id = detId();  // or detId.rawId()
0086
0087       //    std::cout << "GeomDetUnit is of type " << detId.det() << " and raw id = " << id << std::endl;
0088
0089       // There's going to be a lot of messing with field width (and precision) so
0090       // save input values...
0091       int iw = std::cout.width();      // save current width
0092       int ip = std::cout.precision();  // save current precision
0093
0094       //    std::cout << "\n" << "Parameters of layer# " <<
0095       std::cout << std::setw(4) << icount << std::setw(12) << id << std::oct << std::setw(12) << id << std::dec
0096                 << std::setw(iw) << "   E" << CSCDetId::endcap(id) << " S" << CSCDetId::station(id) << " R"
0097                 << CSCDetId::ring(id) << " C" << std::setw(2) << CSCDetId::chamber(id) << std::setw(iw) << " L"
0098                 << CSCDetId::layer(id);
0099       // CSCDetId::layer(id)  << " are:" << std::endl;
0100
0101       // Check global coordinates of centre of CSCLayer, and how
0102       // local z direction relates to global z direction
0103
0104       LocalPoint lCentre(0., 0., 0.);
0105       GlobalPoint gCentre = layer->toGlobal(lCentre);
0106
0107       LocalPoint lCentre1(0., 0., -1.);
0108       GlobalPoint gCentre1 = layer->toGlobal(lCentre1);
0109
0110       LocalPoint lCentre2(0., 0., 1.);
0111       GlobalPoint gCentre2 = layer->toGlobal(lCentre2);
0112
0113       //        std::cout << "\nlocal(0,0,-1) = global " << gCentre1 << std::endl;
0114       //        std::cout <<   "local(0,0)    = global " << gCentre  << std::endl;
0115       //        std::cout <<   "local(0,0,+1) = global " << gCentre2 << std::endl;
0116
0117       double gx = gCentre.x();
0118       double gy = gCentre.y();
0119       double gz = gCentre.z();
0120       double gz1 = gCentre1.z();
0121       double gz2 = gCentre2.z();
0122       if (fabs(gx) < 1.e-06)
0123         gx = 0.;
0124       if (fabs(gy) < 1.e-06)
0125         gy = 0.;
0126       if (fabs(gz) < 1.e-06)
0127         gz = 0.;
0128       if (fabs(gz1) < 1.e-06)
0129         gz1 = 0.;
0130       if (fabs(gz2) < 1.e-06)
0131         gz2 = 0.;
0132
0133       int now = 9;
0134       int nop = 5;
0135       std::cout << std::setw(now) << std::setprecision(nop) << gx << std::setw(now) << std::setprecision(nop) << gy
0136                 << std::setw(now) << std::setprecision(nop) << gz << std::setw(now) << std::setprecision(nop) << gz1
0137                 << std::setw(now) << std::setprecision(nop) << gz2;
0138
0139       // Global Phi of centre of CSCLayer
0140
0141       //@@ CARE The following attempted conversion to degrees can be easily
0142       // subverted by GeometryVector/Phi.h enforcing its range convention!
0143       // Either a) use a separate local double before scaling...
0144       //        double cphi = gCentre.phi();
0145       //        double cphiDeg = cphi * radToDeg;
0146       // or b) use Phi's degree conversion...
0147       double cphiDeg = gCentre.phi().degrees();
0148
0149       // I want to display in range 0 to 360
0150       if (fabs(cphiDeg) < 1.e-06) {
0151         cphiDeg = 0.;
0152       } else if (cphiDeg < 0.) {
0153         cphiDeg += 360.;
0154       } else if (cphiDeg >= 360.) {
0155         std::cout << "WARNING: resetting phi= " << cphiDeg << " to zero." << std::endl;
0156         cphiDeg = 0.;
0157       }
0158
0159       // Handle some occasional ugly precision problems around zero
0160       if (fabs(cphiDeg) < 1.e-06)
0161         cphiDeg = 0.;
0162       //        int iphiDeg = static_cast<int>( cphiDeg );
0163       //    std::cout << "phi(0,0,0) = " << iphiDeg << " degrees" << std::endl;
0164
0165       int nStrips = layer->geometry()->numberOfStrips();
0166       std::cout << std::setw(4) << nStrips;
0167
0168       double cstrip1 = layer->centerOfStrip(1).phi();
0169       double cstripN = layer->centerOfStrip(nStrips).phi();
0170
0171       double phiwid = layer->geometry()->stripPhiPitch();
0172       double stripwid = layer->geometry()->stripPitch();
0173       double stripoff = layer->geometry()->stripOffset();
0174       double phidif = fabs(cstrip1 - cstripN);
0175
0176       // May have one strip at 180-epsilon and other at -180+epsilon
0177       // If so the raw difference is 360-(phi extent of chamber)
0178       // Want to reset that to (phi extent of chamber):
0179       if (phidif > dPi)
0180         phidif = fabs(phidif - 2. * dPi);
0181       double phiwid_check = phidif / double(nStrips - 1);
0182
0183       // Clean up some stupid floating decimal aesthetics
0184       cstrip1 = cstrip1 * radToDeg;
0185       if (fabs(cstrip1) < 1.e-06)
0186         cstrip1 = 0.;
0187       else if (cstrip1 < 0.)
0188         cstrip1 += 360.;
0189
0190       cstripN = cstripN * radToDeg;
0191       if (fabs(cstripN) < 1.e-06)
0192         cstripN = 0.;
0193       else if (cstripN < 0.)
0194         cstripN += 360.;
0195
0196       if (fabs(stripoff) < 1.e-06)
0197         stripoff = 0.;
0198
0199       now = 9;
0200       nop = 4;
0201       std::cout << std::setw(now) << std::setprecision(nop) << cphiDeg << std::setw(now) << std::setprecision(nop)
0202                 << cstrip1 << std::setw(now) << std::setprecision(nop) << cstripN << std::setw(now)
0203                 << std::setprecision(nop) << phiwid << std::setw(now) << std::setprecision(nop) << phiwid_check
0204                 << std::setw(now) << std::setprecision(nop) << stripwid << std::setw(now) << std::setprecision(nop)
0205                 << stripoff;
0206       //          << std::setw(8) << (layer->getOwner())->sector() ; //@@ No sector yet!
0207
0208       // Layer geometry:  layer corner phi's...
0209
0210       std::array<const float, 4> const& parameters = layer->geometry()->parameters();
0211       // these parameters are half-lengths, due to GEANT
0212       float hBottomEdge = parameters[0];
0213       float hTopEdge = parameters[1];
0214       float hThickness = parameters[2];
0215       float hApothem = parameters[3];
0216
0217       //@@ CHECK
0218       //    std::cout << "hB=" << hBottomEdge << " hT=" << hTopEdge << " hTh=" << hThickness <<
0219       //      " hAp=" << hApothem << std::endl;
0220
0221       // first the face nearest the interaction
0222       // get the other face by using positive hThickness
0223       LocalPoint upperRightLocal(hTopEdge, hApothem, -hThickness);
0224       LocalPoint upperLeftLocal(-hTopEdge, hApothem, -hThickness);
0225       LocalPoint lowerRightLocal(hBottomEdge, -hApothem, -hThickness);
0226       LocalPoint lowerLeftLocal(-hBottomEdge, -hApothem, -hThickness);
0227
0228       GlobalPoint upperRightGlobal = layer->toGlobal(upperRightLocal);
0229       GlobalPoint upperLeftGlobal = layer->toGlobal(upperLeftLocal);
0230       GlobalPoint lowerRightGlobal = layer->toGlobal(lowerRightLocal);
0231       GlobalPoint lowerLeftGlobal = layer->toGlobal(lowerLeftLocal);
0232
0233       float uRGp = upperRightGlobal.phi().degrees();
0234       float uLGp = upperLeftGlobal.phi().degrees();
0235       float lRGp = lowerRightGlobal.phi().degrees();
0236       float lLGp = lowerLeftGlobal.phi().degrees();
0237
0238       now = 9;
0239       std::cout << std::setw(now) << uRGp << std::setw(now) << uLGp << std::setw(now) << lRGp << std::setw(now) << lLGp
0240                 << std::endl;
0241
0242       // Reset the values we changed
0243       std::cout << std::setprecision(ip) << std::setw(iw);
0244
0245       // Check idToDetUnit
0246       const GeomDetUnit* gdu = pDD->idToDetUnit(detId);
0247       assert(gdu == layer);
0248       // Check idToDet
0249       const GeomDet* gd = pDD->idToDet(detId);
0250       assert(gd == layer);
0251     } else {
0252       std::cout << "Could not dynamic_cast to a CSCLayer* " << std::endl;
0253     }
0254   }
0255   std::cout << dashedLine_ << " end" << std::endl;
0256 }
0257
0258 //define this as a plug-in
0259 DEFINE_FWK_MODULE(CSCGeometryAnalyzer);