File indexing completed on 2024-04-06 12:14:25
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();
0074
0075
0076
0077 int iw = std::cout.width();
0078 int ip = std::cout.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
0086
0087
0088 const Surface& bSurface = layer->surface();
0089
0090
0091
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
0125
0126
0127
0128
0129
0130
0131
0132 double cphiDeg = gCentre.phi().degrees();
0133
0134
0135
0136
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
0147
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
0161
0162
0163 if (phidif > dPi)
0164 phidif = fabs(phidif - 2. * dPi);
0165 double phiwid_check = phidif / double(nStrips - 1);
0166
0167
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
0191
0192
0193
0194 std::array<const float, 4> const& parameters = layer->geometry()->parameters();
0195
0196 float hBottomEdge = parameters[0];
0197 float hTopEdge = parameters[1];
0198 float hThickness = parameters[2];
0199 float hApothem = parameters[3];
0200
0201
0202
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
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
0234 DEFINE_FWK_MODULE(CSCGeometryAsLayers);