test/stubs/CSCGACwithB.cc

0001 // Modified version of CSCGeometryAsChambers.cc which just adds value of magnetic field
0002 // at centre of each CSC (Nikolai Terentiev needed the values for CSC gain studies.)
0003 // Tim Cox 11.03.2011 - drop ME1/A
0004
0005 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0011 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0012 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0013 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015
0016 #include <string>
0017 #include <iomanip>  // for setw() etc.
0018 #include <vector>
0019
0020 class CSCGACwithB : public edm::one::EDAnalyzer<> {
0021 public:
0022   explicit CSCGACwithB(const edm::ParameterSet&);
0023   ~CSCGACwithB() override = default;
0024
0025   void beginJob() override {}
0026   void analyze(edm::Event const&, edm::EventSetup const&) override;
0027   void endJob() override {}
0028
0029   const std::string& myName() { return myName_; }
0030
0031 private:
0032   const int dashedLineWidth_;
0033   const std::string dashedLine_;
0034   const std::string myName_;
0035   const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> tokGeom_;
0036   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tokField_;
0037 };
0038
0039 CSCGACwithB::CSCGACwithB(const edm::ParameterSet& iConfig)
0040     : dashedLineWidth_(175),
0041       dashedLine_(std::string(dashedLineWidth_, '-')),
0042       myName_("CSCGACwithB"),
0043       tokGeom_(esConsumes()),
0044       tokField_(esConsumes()) {}
0045
0046 void CSCGACwithB::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0047   std::cout << myName() << ": Analyzer..." << std::endl;
0048   std::cout << "start " << dashedLine_ << std::endl;
0049
0050   // Access CSC geometry
0051   const edm::ESHandle<CSCGeometry>& pDD = iSetup.getHandle(tokGeom_);
0052   std::cout << " Geometry node for CSCGeom is  " << &(*pDD) << std::endl;
0053   std::cout << " " << pDD->dets().size() << " detectors" << std::endl;
0054   std::cout << " I have " << pDD->detTypes().size() << " types"
0055             << "\n"
0056             << std::endl;
0057   std::cout << " I have " << pDD->detUnits().size() << " detUnits" << std::endl;
0058   std::cout << " I have " << pDD->dets().size() << " dets" << std::endl;
0059   std::cout << " I have " << pDD->layers().size() << " layers" << std::endl;
0060   std::cout << " I have " << pDD->chambers().size() << " chambers" << std::endl;
0061
0062   // Access magnetic field
0063   const edm::ESHandle<MagneticField>& magneticField = iSetup.getHandle(tokField_);
0064
0065   std::cout << myName() << ": Begin iteration over geometry..." << std::endl;
0066
0067   int icount = 0;  // non-ME1/A chambers
0068   int jcount = 0;  // all chambers
0069
0070   const CSCGeometry::ChamberContainer& vc = pDD->chambers();
0071   std::cout << "No. of chambers stored = " << vc.size() << std::endl;
0072
0073   std::cout << "\n  #     id(dec)      id(oct)    labels      length       width      thickness   "
0074                "  g(x=0)   g(y=0)   g(z=0)  g(z=-1)  g(z=+1)   bx(centre)  by(centre)  bz(centre)"
0075                "  phi(0)"
0076             << std::endl;
0077   std::cout << dashedLine_ << std::endl;
0078
0079   for (auto chamber : vc) {
0080     if (chamber) {
0081       ++jcount;
0082
0083       CSCDetId detId = chamber->id();
0084       int id = detId();  // or detId.rawId()
0085
0086       // *** Nikolai doesn't want ME1/A polluting the values ***
0087
0088       if (detId.ring() == 4)
0089         continue;
0090       ++icount;
0091
0092       // There's going to be a lot of messing with field width (and precision) so
0093       // save input values...
0094       int iw = std::cout.width();      // save current width
0095       int ip = std::cout.precision();  // save current precision
0096
0097       std::cout << std::setw(4) << icount << std::setw(12) << id << std::oct << std::setw(12) << id << std::dec
0098                 << std::setw(iw) << "   E" << detId.endcap() << " S" << detId.station() << " R" << detId.ring() << " C"
0099                 << std::setw(2) << detId.chamber() << std::setw(iw);
0100
0101       // What's its surface?
0102       // The surface knows how to transform local <-> global
0103
0104       const Surface& bSurface = chamber->surface();
0105
0106       //    std::cout << " length=" << bSurface.bounds().length() <<
0107       //                     ", width=" << bSurface.bounds().width() <<
0108       //                     ", thickness=" << bSurface.bounds().thickness() << std::endl;
0109
0110       std::cout << std::setw(12) << std::setprecision(8) << bSurface.bounds().length() << std::setw(12)
0111                 << std::setprecision(8) << bSurface.bounds().width() << std::setw(12) << std::setprecision(6)
0112                 << bSurface.bounds().thickness();
0113
0114       // Check global coordinates of centre of CSCChamber, and how
0115       // local z direction relates to global z direction
0116
0117       LocalPoint lCentre(0., 0., 0.);
0118       GlobalPoint gCentre = bSurface.toGlobal(lCentre);
0119
0120       LocalPoint lCentre1(0., 0., -1.);
0121       GlobalPoint gCentre1 = bSurface.toGlobal(lCentre1);
0122
0123       LocalPoint lCentre2(0., 0., 1.);
0124       GlobalPoint gCentre2 = bSurface.toGlobal(lCentre2);
0125
0126       double gx = gCentre.x();
0127       double gy = gCentre.y();
0128       double gz = gCentre.z();
0129       double gz1 = gCentre1.z();
0130       double gz2 = gCentre2.z();
0131       if (fabs(gx) < 1.e-06)
0132         gx = 0.;
0133       if (fabs(gy) < 1.e-06)
0134         gy = 0.;
0135       if (fabs(gz) < 1.e-06)
0136         gz = 0.;
0137       if (fabs(gz1) < 1.e-06)
0138         gz1 = 0.;
0139       if (fabs(gz2) < 1.e-06)
0140         gz2 = 0.;
0141
0142       int now = 9;
0143       int nop = 5;
0144       std::cout << std::setw(now) << std::setprecision(nop) << gx << std::setw(now) << std::setprecision(nop) << gy
0145                 << std::setw(now) << std::setprecision(nop) << gz << std::setw(now) << std::setprecision(nop) << gz1
0146                 << std::setw(now) << std::setprecision(nop) << gz2;
0147
0148       // Magnetic field at centre of CSC (global position = gCentre)
0149       float bx = magneticField->inTesla(gCentre).x();
0150       float by = magneticField->inTesla(gCentre).y();
0151       float bz = magneticField->inTesla(gCentre).z();
0152
0153       now = 12;
0154       std::cout << std::setw(now) << std::setprecision(nop) << bx << std::setw(now) << std::setprecision(nop) << by
0155                 << std::setw(now) << std::setprecision(nop) << bz;
0156
0157       // Global Phi of centre of CSC
0158
0159       //@@ CARE The following attempted conversion to degrees can be easily
0160       // subverted by GeometryVector/Phi.h enforcing its range convention!
0161       // Either a) use a separate local double before scaling...
0162       //        double cphi = gCentre.phi();
0163       //        double cphiDeg = cphi * radToDeg;
0164       // or b) use Phi's degree conversion...
0165       double cphiDeg = gCentre.phi().degrees();
0166
0167       // I want to display in range 0 to 360
0168
0169       // Handle some occasional ugly precision problems around zero
0170       if (fabs(cphiDeg) < 1.e-06) {
0171         cphiDeg = 0.;
0172       } else if (cphiDeg < 0.) {
0173         cphiDeg += 360.;
0174       } else if (cphiDeg >= 360.) {
0175         std::cout << "WARNING: resetting phi= " << cphiDeg << " to zero." << std::endl;
0176         cphiDeg = 0.;
0177       }
0178
0179       //        int iphiDeg = static_cast<int>( cphiDeg );
0180       //    std::cout << "phi(0,0,0) = " << iphiDeg << " degrees" << std::endl;
0181
0182       now = 9;
0183       nop = 4;
0184       std::cout << std::setw(now) << std::setprecision(nop) << cphiDeg << std::endl;
0185
0186       // Reset the values we changed
0187       std::cout << std::setprecision(ip) << std::setw(iw);
0188
0189     } else {
0190       std::cout << "*** DANGER WILL ROBINSON! *** Stored a null CSCChamber*?  " << std::endl;
0191     }
0192   }
0193
0194   std::cout << dashedLine_ << " end" << std::endl;
0195   std::cout << "processed " << jcount << " chambers, and dumped values for " << icount << " chambers." << std::endl;
0196 }
0197
0198 //define this as a plug-in
0199 DEFINE_FWK_MODULE(CSCGACwithB);