CSCGeometryAsChambers

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Geometry/CSCGeometry/interface/CSCGeometry.h"
#include "Geometry/CSCGeometry/interface/CSCLayer.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include <string>
#include <iomanip>  // for setw() etc.
#include <vector>

class CSCGeometryAsChambers : public edm::one::EDAnalyzer<> {
public:
  explicit CSCGeometryAsChambers(const edm::ParameterSet&);
  ~CSCGeometryAsChambers() override = default;

  void beginJob() override {}
  void analyze(edm::Event const&, edm::EventSetup const&) override;
  void endJob() override {}

  const std::string& myName() { return myName_; }

private:
  const int dashedLineWidth_;
  const std::string dashedLine_;
  const std::string myName_;
  const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> tokGeom_;
};

CSCGeometryAsChambers::CSCGeometryAsChambers(const edm::ParameterSet& iConfig)
    : dashedLineWidth_(132),
      dashedLine_(std::string(dashedLineWidth_, '-')),
      myName_("CSCGeometryAsChambers"),
      tokGeom_(esConsumes()) {}

void CSCGeometryAsChambers::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  std::cout << myName() << ": Analyzer..." << std::endl;
  std::cout << "start " << dashedLine_ << std::endl;

  const edm::ESHandle<CSCGeometry>& pDD = iSetup.getHandle(tokGeom_);
  std::cout << " Geometry node for CSCGeom is  " << &(*pDD) << std::endl;
  std::cout << " I have " << pDD->dets().size() << " detectors" << std::endl;
  std::cout << " I have " << pDD->detTypes().size() << " types"
            << "\n"
            << std::endl;
  std::cout << " I have " << pDD->detUnits().size() << " detUnits" << std::endl;
  std::cout << " I have " << pDD->dets().size() << " dets" << std::endl;
  std::cout << " I have " << pDD->layers().size() << " layers" << std::endl;
  std::cout << " I have " << pDD->chambers().size() << " chambers" << std::endl;

  std::cout << myName() << ": Begin iteration over geometry..." << std::endl;

  int icount = 0;

  const CSCGeometry::ChamberContainer& vc = pDD->chambers();
  std::cout << "No. of chambers stored = " << vc.size() << std::endl;

  std::cout << "\n  #     id(dec)      id(oct)    labels      length       width      thickness   "
               "  g(x=0)   g(y=0)   g(z=0)  g(z=-1)  g(z=+1) "
               "  phi(0)"
            << std::endl;
  std::cout << dashedLine_ << std::endl;

  for (auto chamber : vc) {
    if (chamber) {
      ++icount;

      // What's its DetId?

      CSCDetId detId = chamber->id();
      int id = detId();  // or detId.rawId()

      //	std::cout << "GeomDetUnit is of type " << detId.det() << " and raw id = " << id << std::endl;

      // There's going to be a lot of messing with field width (and precision) so
      // save input values...
      int iw = std::cout.width();      // save current width
      int ip = std::cout.precision();  // save current precision

      //	std::cout << "\n" << "Parameters of chamber# " <<
      std::cout << std::setw(4) << icount << std::setw(12) << id << std::oct << std::setw(12) << id << std::dec
                << std::setw(iw) << "   E" << detId.endcap() << " S" << detId.station() << " R" << detId.ring() << " C"
                << std::setw(2) << detId.chamber() << std::setw(iw);

      // What's its surface?
      // The surface knows how to transform local <-> global

      const Surface& bSurface = chamber->surface();

      //	std::cout << " length=" << bSurface.bounds().length() <<
      //       	             ", width=" << bSurface.bounds().width() <<
      //                     ", thickness=" << bSurface.bounds().thickness() << std::endl;

      std::cout << std::setw(12) << std::setprecision(8) << bSurface.bounds().length() << std::setw(12)
                << std::setprecision(8) << bSurface.bounds().width() << std::setw(12) << std::setprecision(6)
                << bSurface.bounds().thickness();

      // Check global coordinates of centre of CSCChamber, and how
      // local z direction relates to global z direction

      LocalPoint lCentre(0., 0., 0.);
      GlobalPoint gCentre = bSurface.toGlobal(lCentre);

      LocalPoint lCentre1(0., 0., -1.);
      GlobalPoint gCentre1 = bSurface.toGlobal(lCentre1);

      LocalPoint lCentre2(0., 0., 1.);
      GlobalPoint gCentre2 = bSurface.toGlobal(lCentre2);

      //	std::cout << "local(0,0,-1) = global " << gCentre1 << std::endl;
      //	std::cout << "local(0,0)    = global " << gCentre  << std::endl;
      //	std::cout << "local(0,0,+1) = global " << gCentre2 << std::endl;

      double gx = gCentre.x();
      double gy = gCentre.y();
      double gz = gCentre.z();
      double gz1 = gCentre1.z();
      double gz2 = gCentre2.z();
      if (fabs(gx) < 1.e-06)
        gx = 0.;
      if (fabs(gy) < 1.e-06)
        gy = 0.;
      if (fabs(gz) < 1.e-06)
        gz = 0.;
      if (fabs(gz1) < 1.e-06)
        gz1 = 0.;
      if (fabs(gz2) < 1.e-06)
        gz2 = 0.;

      int now = 9;
      int nop = 5;
      std::cout << std::setw(now) << std::setprecision(nop) << gx << std::setw(now) << std::setprecision(nop) << gy
                << std::setw(now) << std::setprecision(nop) << gz << std::setw(now) << std::setprecision(nop) << gz1
                << std::setw(now) << std::setprecision(nop) << gz2;

      // Global Phi of centre of CSCChamber

      //@@ CARE The following attempted conversion to degrees can be easily
      // subverted by GeometryVector/Phi.h enforcing its range convention!
      // Either a) use a separate local double before scaling...
      //        double cphi = gCentre.phi();
      //        double cphiDeg = cphi * radToDeg;
      // or b) use Phi's degree conversion...
      double cphiDeg = gCentre.phi().degrees();

      // I want to display in range 0 to 360

      // Handle some occasional ugly precision problems around zero
      if (fabs(cphiDeg) < 1.e-06) {
        cphiDeg = 0.;
      } else if (cphiDeg < 0.) {
        cphiDeg += 360.;
      } else if (cphiDeg >= 360.) {
        std::cout << "WARNING: resetting phi= " << cphiDeg << " to zero." << std::endl;
        cphiDeg = 0.;
      }

      //        int iphiDeg = static_cast<int>( cphiDeg );
      //	std::cout << "phi(0,0,0) = " << iphiDeg << " degrees" << std::endl;

      now = 9;
      nop = 4;
      std::cout << std::setw(now) << std::setprecision(nop) << cphiDeg << std::endl;

      // Reset the values we changed
      std::cout << std::setprecision(ip) << std::setw(iw);

    } else {
      std::cout << "*** DANGER WILL ROBINSON! *** Stored a null CSCChamber*?  " << std::endl;
    }
  }

  std::cout << dashedLine_ << " end" << std::endl;
}

//define this as a plug-in
DEFINE_FWK_MODULE(CSCGeometryAsChambers);