Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-12 05:44:31

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