HcalTowerAlgo/test/CaloTowerGeometryAnalyzer.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0008 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0009 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0010 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0011 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0012
0013 #include <fstream>
0014 #include <iostream>
0015 #include <cmath>
0016
0017 namespace {
0018   bool AreSame(double a, double b, double epsilon) { return std::fabs(a - b) < epsilon; }
0019
0020   /*
0021   inline double
0022   round( double n, int digits )
0023   {
0024     double mult = pow( 10, digits );
0025     return floor( n * mult ) / mult;
0026   }
0027   */
0028
0029   std::map<int, std::string> hcalmapping = {
0030       {1, "HB"}, {2, "HE"}, {3, "HO"}, {4, "HF"}, {5, "HT"}, {6, "ZDC"}, {0, "Empty"}};
0031 }  // namespace
0032
0033 class CaloTowerGeometryAnalyzer : public edm::one::EDAnalyzer<> {
0034 public:
0035   explicit CaloTowerGeometryAnalyzer(const edm::ParameterSet&);
0036   ~CaloTowerGeometryAnalyzer(void) override;
0037
0038   void beginJob() override {}
0039   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0040   void endJob() override {}
0041
0042 private:
0043   std::string m_fname;
0044   double m_epsilon;
0045   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0046 };
0047
0048 CaloTowerGeometryAnalyzer::CaloTowerGeometryAnalyzer(const edm::ParameterSet& iConfig)
0049     : m_fname("CaloTower.cells"), m_epsilon(0.004) {
0050   m_fname = iConfig.getParameter<std::string>("FileName");
0051   m_epsilon = iConfig.getParameter<double>("Epsilon");
0052   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0053 }
0054
0055 CaloTowerGeometryAnalyzer::~CaloTowerGeometryAnalyzer(void) {}
0056
0057 void CaloTowerGeometryAnalyzer::analyze(const edm::Event& /*iEvent*/, const edm::EventSetup& iSetup) {
0058   std::fstream fAll(std::string(m_fname + ".all").c_str(), std::ios_base::out);
0059   std::fstream f(std::string(m_fname + ".diff").c_str(), std::ios_base::out);
0060
0061   const CaloGeometry* caloGeom = &iSetup.getData(tok_geom_);
0062
0063   const std::vector<DetId>& dha(caloGeom->getSubdetectorGeometry(DetId::Hcal, 1)->getValidDetIds());
0064
0065   const std::vector<DetId>& ids(caloGeom->getValidDetIds());
0066
0067   fAll << std::setw(4) << "iEta" << std::setw(4) << "iPhi" << std::setw(4) << "Z" << std::setw(10) << "Sub:Depth"
0068        << std::setw(6) << "Eta" << std::setw(13) << "Phi" << std::setw(18) << "Corners in Eta\n";
0069
0070   f << std::setw(4) << "iEta" << std::setw(4) << "iPhi" << std::setw(4) << "Z" << std::setw(10) << "Sub:Depth"
0071     << std::setw(6) << "Eta" << std::setw(13) << "Phi" << std::setw(18) << "Corners in Eta\n";
0072
0073   for (auto id : ids) {
0074     const CaloGenericDetId cgid(id);
0075     if (cgid.isCaloTower()) {
0076       const CaloTowerDetId cid(id);
0077       const int ie(cid.ieta());
0078       const int ip(cid.iphi());
0079       const int iz(cid.zside());
0080
0081       fAll << std::setw(4) << ie << std::setw(4) << ip << std::setw(4) << iz << std::setw(6) << "-";
0082
0083       auto cell = caloGeom->getGeometry(id);
0084       assert(cell);
0085       const GlobalPoint& pos = cell->getPosition();
0086       double eta = pos.eta();
0087       double phi = pos.phi();
0088       fAll << std::setw(10) << eta << std::setw(13) << phi;
0089
0090       const CaloCellGeometry::CornersVec& corners(cell->getCorners());
0091       for (unsigned int i(0); i != corners.size(); ++i) {
0092         const GlobalPoint& cpos = corners[i];
0093         double ceta = cpos.eta();
0094
0095         fAll << std::setw(13) << ceta;
0096       }
0097       fAll << "\n";
0098
0099       for (auto ii : dha) {
0100         const HcalDetId hid(ii);
0101         const int iie(hid.ieta());
0102         const int iip(hid.iphi());
0103         const int iiz(hid.zside());
0104         const int iid(hid.depth());
0105
0106         if (ie == iie && ip == iip && iz == iiz) {
0107           fAll << std::setw(4) << iie << std::setw(4) << iip << std::setw(4) << iiz;
0108
0109           std::map<int, std::string>::const_iterator iter = hcalmapping.find(hid.subdet());
0110           if (iter != hcalmapping.end()) {
0111             fAll << std::setw(4) << iter->second.c_str() << ":" << iid;
0112           }
0113
0114           auto hcell = caloGeom->getGeometry(ii);
0115           assert(hcell);
0116           const GlobalPoint& hpos = hcell->getPosition();
0117           double heta = hpos.eta();
0118           double hphi = hpos.phi();
0119
0120           fAll << std::setw(10) << heta << std::setw(13) << hphi;
0121
0122           const CaloCellGeometry::CornersVec& hcorners(hcell->getCorners());
0123           for (unsigned int i(0); i != hcorners.size(); ++i) {
0124             const GlobalPoint& hcpos = hcorners[i];
0125             double hceta = hcpos.eta();
0126
0127             fAll << std::setw(13) << hceta;
0128           }
0129
0130           if (!AreSame(eta, heta, m_epsilon)) {
0131             fAll << "*DIFFER in Eta*";
0132
0133             f << std::setw(4) << ie << std::setw(4) << ip << std::setw(4) << iz << std::setw(6) << "-";
0134
0135             f << std::setw(10) << eta << std::setw(13) << phi;
0136
0137             for (unsigned int i(0); i != corners.size(); ++i) {
0138               const GlobalPoint& cpos = corners[i];
0139               double ceta = cpos.eta();
0140
0141               f << std::setw(9) << ceta;
0142             }
0143             f << "\n";
0144
0145             f << std::setw(4) << iie << std::setw(4) << iip << std::setw(4) << iiz;
0146
0147             if (iter != hcalmapping.end()) {
0148               f << std::setw(4) << iter->second.c_str() << ":" << iid;
0149             }
0150
0151             f << std::setw(10) << heta << std::setw(13) << hphi;
0152
0153             for (unsigned int i(0); i != hcorners.size(); ++i) {
0154               const GlobalPoint& hcpos = hcorners[i];
0155               double hceta = hcpos.eta();
0156
0157               f << std::setw(9) << hceta;
0158             }
0159             f << "\n\n";
0160           }
0161
0162           if (!AreSame(phi, hphi, m_epsilon))
0163             fAll << " *DIFFER in Phi*";
0164           fAll << "\n";
0165         }
0166       }
0167     }
0168   }
0169
0170   fAll.close();
0171   f.close();
0172 }
0173
0174 DEFINE_FWK_MODULE(CaloTowerGeometryAnalyzer);