GEMGeometry/test/ME0GeometryAnalyzer.cc

0001 /** Derived from DTGeometryAnalyzer by Nicola Amapane
0002  *
0003  *  \author M. Maggi - INFN Bari
0004  */
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0013 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0014 #include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
0015 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0016
0017 #include "DataFormats/Math/interface/deltaPhi.h"
0018
0019 #include <memory>
0020 #include <fstream>
0021 #include <string>
0022 #include <cmath>
0023 #include <vector>
0024 #include <iomanip>
0025 #include <set>
0026
0027 class ME0GeometryAnalyzer : public edm::one::EDAnalyzer<> {
0028 public:
0029   ME0GeometryAnalyzer(const edm::ParameterSet& pset);
0030
0031   ~ME0GeometryAnalyzer() override;
0032
0033   void beginJob() override {}
0034   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0035   void endJob() override {}
0036
0037 private:
0038   const std::string& myName() { return myName_; }
0039
0040   const int dashedLineWidth_;
0041   const std::string dashedLine_;
0042   const std::string myName_;
0043   const edm::ESGetToken<ME0Geometry, MuonGeometryRecord> tokGeom_;
0044   std::ofstream ofos;
0045 };
0046
0047 using namespace std;
0048
0049 ME0GeometryAnalyzer::ME0GeometryAnalyzer(const edm::ParameterSet& /*iConfig*/)
0050     : dashedLineWidth_(104),
0051       dashedLine_(string(dashedLineWidth_, '-')),
0052       myName_("ME0GeometryAnalyzer"),
0053       tokGeom_{esConsumes<ME0Geometry, MuonGeometryRecord>(edm::ESInputTag{})} {
0054   ofos.open("ME0testOutput.out");
0055   ofos << "======================== Opening output file" << endl;
0056 }
0057
0058 ME0GeometryAnalyzer::~ME0GeometryAnalyzer() {
0059   ofos.close();
0060   ofos << "======================== Closing output file" << endl;
0061 }
0062
0063 void ME0GeometryAnalyzer::analyze(const edm::Event& /*iEvent*/, const edm::EventSetup& iSetup) {
0064   const auto& pDD = iSetup.getData(tokGeom_);
0065
0066   ofos << myName() << ": Analyzer..." << endl;
0067   ofos << "start " << dashedLine_ << endl;
0068
0069   ofos << " Geometry node for ME0Geom is  " << &pDD << endl;
0070   ofos << " detTypes       \t" << pDD.detTypes().size() << endl;
0071   ofos << " GeomDetUnit       \t" << pDD.detUnits().size() << endl;
0072   ofos << " GeomDet           \t" << pDD.dets().size() << endl;
0073   ofos << " GeomDetUnit DetIds\t" << pDD.detUnitIds().size() << endl;
0074   ofos << " eta partitions \t" << pDD.etaPartitions().size() << endl;
0075   ofos << " layers         \t" << pDD.layers().size() << endl;
0076   ofos << " chambers       \t" << pDD.chambers().size() << endl;
0077   // ofos << " regions        \t"              <<pDD.regions().size() << endl;
0078
0079   // checking uniqueness of roll detIds
0080   bool flagNonUniqueRollID = false;
0081   bool flagNonUniqueRollRawID = false;
0082   int nstrips = 0;
0083   int npads = 0;
0084   for (auto roll1 : pDD.etaPartitions()) {
0085     nstrips += roll1->nstrips();
0086     npads += roll1->npads();
0087     for (auto roll2 : pDD.etaPartitions()) {
0088       if (roll1 != roll2) {
0089         if (roll1->id() == roll2->id())
0090           flagNonUniqueRollID = true;
0091         if (roll1->id().rawId() == roll2->id().rawId())
0092           flagNonUniqueRollRawID = true;
0093       }
0094     }
0095   }
0096   if (flagNonUniqueRollID or flagNonUniqueRollRawID)
0097     ofos << " -- WARNING: non unique roll Ids!!!" << endl;
0098
0099   // checking uniqueness of layer detIds
0100   bool flagNonUniqueLaID = false;
0101   bool flagNonUniqueLaRawID = false;
0102   for (auto la1 : pDD.layers()) {
0103     for (auto la2 : pDD.layers()) {
0104       if (la1 != la2) {
0105         if (la1->id() == la2->id())
0106           flagNonUniqueLaID = true;
0107         if (la1->id().rawId() == la2->id().rawId())
0108           flagNonUniqueLaRawID = true;
0109       }
0110     }
0111   }
0112   if (flagNonUniqueLaID or flagNonUniqueLaRawID)
0113     ofos << " -- WARNING: non unique layer Ids!!!" << endl;
0114
0115   // checking uniqueness of chamber detIds
0116   bool flagNonUniqueChID = false;
0117   bool flagNonUniqueChRawID = false;
0118   for (auto ch1 : pDD.chambers()) {
0119     for (auto ch2 : pDD.chambers()) {
0120       if (ch1 != ch2) {
0121         if (ch1->id() == ch2->id())
0122           flagNonUniqueChID = true;
0123         if (ch1->id().rawId() == ch2->id().rawId())
0124           flagNonUniqueChRawID = true;
0125       }
0126     }
0127   }
0128   if (flagNonUniqueChID or flagNonUniqueChRawID)
0129     ofos << " -- WARNING: non unique chamber Ids!!!" << endl;
0130
0131   // print out number of strips and pads
0132   ofos << " total number of strips\t" << nstrips << endl;
0133   ofos << " total number of pads  \t" << npads << endl;
0134
0135   ofos << myName() << ": Begin iteration over geometry..." << endl;
0136   ofos << "iter " << dashedLine_ << endl;
0137
0138   ofos << myName() << "Begin ME0Geometry TEST" << endl;
0139
0140   /*
0141    * possible checklist for an eta partition:
0142    *   base_bottom, base_top, height, strips, pads
0143    *   cx, cy, cz, ceta, cphi
0144    *   tx, ty, tz, teta, tphi
0145    *   bx, by, bz, beta, bphi
0146    *   pitch center, pitch bottom, pitch top
0147    *   deta, dphi
0148    *   gap thicess
0149    *   sum of all dx + gap = chamber height
0150    */
0151
0152   int i = 1;
0153   for (auto ch : pDD.chambers()) {
0154     ME0DetId chId(ch->id());
0155     int nLayers(ch->nLayers());
0156     ofos << "\tME0Chamber " << i << ", ME0DetId = " << chId.rawId() << ", " << chId << " has " << nLayers << " layers."
0157          << endl;
0158     int j = 1;
0159     for (auto la : ch->layers()) {
0160       ME0DetId laId(la->id());
0161       int nRolls(la->nEtaPartitions());
0162       ofos << "\t\tME0Layer " << j << ", ME0DetId = " << laId.rawId() << ", " << laId << " has " << nRolls
0163            << " eta partitions." << endl;
0164
0165       int k = 1;
0166       auto& rolls(la->etaPartitions());
0167       for (auto roll : rolls) {
0168         // for (auto roll : pDD.etaPartitions()){
0169         ME0DetId rId(roll->id());
0170         ofos << "\t\t\tME0EtaPartition " << k << " , ME0DetId = " << rId.rawId() << ", " << rId << endl;
0171
0172         const BoundPlane& bSurface(roll->surface());
0173         const StripTopology* topology(&(roll->specificTopology()));
0174
0175         // base_bottom, base_top, height, strips, pads (all half length)
0176         auto& parameters(roll->specs()->parameters());
0177         float bottomEdge(parameters[0]);
0178         float topEdge(parameters[1]);
0179         float height(parameters[2]);
0180         float nStrips(parameters[3]);
0181         float nPads(parameters[4]);
0182
0183         LocalPoint lCentre(0., 0., 0.);
0184         GlobalPoint gCentre(bSurface.toGlobal(lCentre));
0185
0186         LocalPoint lTop(0., height, 0.);
0187         GlobalPoint gTop(bSurface.toGlobal(lTop));
0188
0189         LocalPoint lBottom(0., -height, 0.);
0190         GlobalPoint gBottom(bSurface.toGlobal(lBottom));
0191
0192         //   gx, gy, gz, geta, gphi (center)
0193         double cx(gCentre.x());
0194         double cy(gCentre.y());
0195         double cz(gCentre.z());
0196         double ceta(gCentre.eta());
0197         int cphi(static_cast<int>(gCentre.phi().degrees()));
0198         if (cphi < 0)
0199           cphi += 360;
0200
0201         double tx(gTop.x());
0202         double ty(gTop.y());
0203         double tz(gTop.z());
0204         double teta(gTop.eta());
0205         int tphi(static_cast<int>(gTop.phi().degrees()));
0206         if (tphi < 0)
0207           tphi += 360;
0208
0209         double bx(gBottom.x());
0210         double by(gBottom.y());
0211         double bz(gBottom.z());
0212         double beta(gBottom.eta());
0213         int bphi(static_cast<int>(gBottom.phi().degrees()));
0214         if (bphi < 0)
0215           bphi += 360;
0216
0217         // pitch bottom, pitch top, pitch centre
0218         float pitch(roll->pitch());
0219         float topPitch(roll->localPitch(lTop));
0220         float bottomPitch(roll->localPitch(lBottom));
0221
0222         // Type - should be GHA0[1-nRolls]
0223         string type(roll->type().name());
0224
0225         // print info about edges
0226         LocalPoint lEdge1(topology->localPosition(0.));
0227         LocalPoint lEdgeN(topology->localPosition((float)nStrips));
0228
0229         double cstrip1(roll->toGlobal(lEdge1).phi().degrees());
0230         double cstripN(roll->toGlobal(lEdgeN).phi().degrees());
0231         double dphi(cstripN - cstrip1);
0232         if (dphi < 0.)
0233           dphi += 360.;
0234
0235         double deta(abs(beta - teta));
0236         const bool printDetails(true);
0237         if (printDetails)
0238           ofos << "\t\t\t\tType: " << type << endl
0239                << "\t\t\t\tDimensions[cm]: b = " << bottomEdge * 2 << ", B = " << topEdge * 2 << ", H  = " << height * 2
0240                << endl
0241                << "\t\t\t\tnStrips = " << nStrips << ", nPads =  " << nPads << endl
0242                << "\t\t\t\ttop(x,y,z)[cm] = (" << tx << ", " << ty << ", " << tz << "), top (eta,phi) = (" << teta
0243                << ", " << tphi << ")" << endl
0244                << "\t\t\t\tcenter(x,y,z) = (" << cx << ", " << cy << ", " << cz << "), center(eta,phi) = (" << ceta
0245                << ", " << cphi << ")" << endl
0246                << "\t\t\t\tbottom(x,y,z) = (" << bx << ", " << by << ", " << bz << "), bottom(eta,phi) = (" << beta
0247                << ", " << bphi << ")" << endl
0248                << "\t\t\t\tpitch (top,center,bottom) = " << topPitch << " " << pitch << " " << bottomPitch
0249                << ", dEta = " << deta << ", dPhi = " << dphi << endl
0250                << "\t\t\t\tlocal pos at strip 1 " << lEdge1 << " strip N " << lEdgeN << endl;
0251         ++k;
0252       }
0253       ++j;
0254     }
0255     ++i;
0256   }
0257   ofos << dashedLine_ << " end" << endl;
0258 }
0259
0260 //define this as a plug-in
0261 #include "FWCore/Framework/interface/MakerMacros.h"
0262 DEFINE_FWK_MODULE(ME0GeometryAnalyzer);