Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-08 08:16:24

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 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0006 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0007 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0008 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0009 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0010 #include "Geometry/HcalTowerAlgo/interface/HcalFlexiHardcodeGeometryLoader.h"
0011 #include "Geometry/HcalTowerAlgo/interface/HcalHardcodeGeometryLoader.h"
0012 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0013 #include "DataFormats/HcalDetId/interface/HcalTrigTowerDetId.h"
0014 #include <iostream>
0015 #include <string>
0016 
0017 class HcalGeometryTester : public edm::one::EDAnalyzer<> {
0018 public:
0019   explicit HcalGeometryTester(const edm::ParameterSet&);
0020   ~HcalGeometryTester(void) override {}
0021 
0022   void beginJob() override {}
0023   void analyze(edm::Event const&, edm::EventSetup const&) override;
0024   void endJob() override {}
0025 
0026 private:
0027   void testValidDetIds(CaloSubdetectorGeometry* geom,
0028                        const HcalTopology& topo,
0029                        DetId::Detector det,
0030                        int subdet,
0031                        const std::string& label);
0032   void testClosestCells(CaloSubdetectorGeometry* geom, const HcalTopology& top);
0033   void testClosestCell(const HcalDetId& detId, CaloSubdetectorGeometry* geom);
0034   void testTriggerGeometry(const HcalTopology& topology);
0035   void testFlexiValidDetIds(CaloSubdetectorGeometry* geom,
0036                             const HcalTopology& topology,
0037                             DetId::Detector det,
0038                             int subdet,
0039                             const std::string& label,
0040                             std::vector<int>& dins);
0041   void testFlexiGeomHF(CaloSubdetectorGeometry* geom);
0042 
0043   bool useOld_;
0044   edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
0045   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0046 };
0047 
0048 HcalGeometryTester::HcalGeometryTester(const edm::ParameterSet& iConfig) {
0049   useOld_ = iConfig.getParameter<bool>("UseOldLoader");
0050   tok_ddrec_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>();
0051   tok_htopo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0052 }
0053 
0054 void HcalGeometryTester::analyze(const edm::Event& /*iEvent*/, const edm::EventSetup& iSetup) {
0055   const HcalDDDRecConstants hcons = iSetup.getData(tok_ddrec_);
0056   const HcalTopology topology = iSetup.getData(tok_htopo_);
0057   CaloSubdetectorGeometry* geom(nullptr);
0058   if (useOld_) {
0059     HcalHardcodeGeometryLoader m_loader;
0060     geom = m_loader.load(topology);
0061   } else {
0062     HcalFlexiHardcodeGeometryLoader m_loader;
0063     geom = m_loader.load(topology, hcons);
0064   }
0065 
0066   testValidDetIds(geom, topology, DetId::Hcal, HcalBarrel, " BARREL ");
0067   testValidDetIds(geom, topology, DetId::Hcal, HcalEndcap, " ENDCAP ");
0068   testValidDetIds(geom, topology, DetId::Hcal, HcalOuter, " OUTER ");
0069   testValidDetIds(geom, topology, DetId::Hcal, HcalForward, " FORWARD ");
0070 
0071   testTriggerGeometry(topology);
0072 
0073   testClosestCells(geom, topology);
0074 
0075   edm::LogVerbatim("HCalGeom") << "HcalGeometryTester::Test SLHC Hcal Geometry";
0076   std::vector<int> dins;
0077 
0078   testFlexiValidDetIds(geom, topology, DetId::Hcal, HcalBarrel, " BARREL ", dins);
0079   testFlexiValidDetIds(geom, topology, DetId::Hcal, HcalEndcap, " ENDCAP ", dins);
0080   testFlexiValidDetIds(geom, topology, DetId::Hcal, HcalOuter, " OUTER ", dins);
0081   testFlexiValidDetIds(geom, topology, DetId::Hcal, HcalForward, " FORWARD ", dins);
0082 
0083   testFlexiGeomHF(geom);
0084 }
0085 
0086 void HcalGeometryTester::testValidDetIds(CaloSubdetectorGeometry* caloGeom,
0087                                          const HcalTopology& topology,
0088                                          DetId::Detector det,
0089                                          int subdet,
0090                                          const std::string& label) {
0091   std::stringstream s;
0092   s << label << " : \n";
0093   const std::vector<DetId>& idshb = caloGeom->getValidDetIds(det, subdet);
0094 
0095   int counter = 0;
0096   for (std::vector<DetId>::const_iterator i = idshb.begin(); i != idshb.end(); i++, ++counter) {
0097     HcalDetId hid = (*i);
0098     s << counter << ": din " << topology.detId2denseId(*i) << ":" << hid;
0099     auto cell = caloGeom->getGeometry(*i);
0100     s << *cell << std::endl;
0101   }
0102 
0103   s << "=== Total " << counter << " cells in " << label << std::endl;
0104   edm::LogVerbatim("HCalGeom") << s.str();
0105 }
0106 
0107 void HcalGeometryTester::testClosestCells(CaloSubdetectorGeometry* g, const HcalTopology& topology) {
0108   // make sure each cell is its own closest cell
0109   HcalDetId barrelDet1(HcalBarrel, 1, 1, 1);
0110   HcalDetId barrelDet2(HcalBarrel, 16, 50, 1);
0111   HcalDetId barrelDet3(HcalBarrel, -6, 20, 3);
0112   HcalDetId barrelDet4(HcalBarrel, -12, 33, 4);
0113   HcalDetId endcapDet1(HcalEndcap, -17, 72, 1);
0114   HcalDetId endcapDet2(HcalEndcap, 29, 35, 1);
0115   HcalDetId endcapDet3(HcalEndcap, 20, 51, 3);
0116   HcalDetId endcapDet4(HcalEndcap, -24, 63, 5);
0117   HcalDetId forwardDet1(HcalForward, 30, 71, 1);
0118   HcalDetId forwardDet3(HcalForward, -40, 71, 1);
0119 
0120   if (topology.valid(barrelDet1))
0121     testClosestCell(barrelDet1, g);
0122   if (topology.valid(barrelDet2))
0123     testClosestCell(barrelDet2, g);
0124   if (topology.valid(barrelDet3))
0125     testClosestCell(barrelDet3, g);
0126   if (topology.valid(barrelDet4))
0127     testClosestCell(barrelDet4, g);
0128   if (topology.valid(endcapDet1))
0129     testClosestCell(endcapDet1, g);
0130   if (topology.valid(endcapDet2))
0131     testClosestCell(endcapDet2, g);
0132   if (topology.valid(endcapDet3))
0133     testClosestCell(endcapDet3, g);
0134   if (topology.valid(endcapDet4))
0135     testClosestCell(endcapDet4, g);
0136   if (topology.valid(forwardDet1))
0137     testClosestCell(forwardDet1, g);
0138   if (topology.valid(forwardDet3))
0139     testClosestCell(forwardDet3, g);
0140 
0141   const std::vector<DetId>& idsb = g->getValidDetIds(DetId::Hcal, HcalBarrel);
0142   for (auto id : idsb) {
0143     testClosestCell(HcalDetId(id), g);
0144   }
0145 
0146   const std::vector<DetId>& idse = g->getValidDetIds(DetId::Hcal, HcalEndcap);
0147   for (auto id : idse) {
0148     testClosestCell(HcalDetId(id), g);
0149   }
0150 }
0151 
0152 void HcalGeometryTester::testClosestCell(const HcalDetId& detId, CaloSubdetectorGeometry* geom) {
0153   auto cell = geom->getGeometry(detId);
0154   HcalDetId closest = geom->getClosestCell(cell->getPosition());
0155   edm::LogVerbatim("HCalGeom") << "i/p " << detId << " position " << cell->getPosition() << " closest " << closest;
0156 
0157   if (closest != detId) {
0158     edm::LogVerbatim("HCalGeom") << "HcalGeometryTester::Mismatch.  Original HCAL cell is " << detId
0159                                  << " while nearest is " << closest << " ***ERROR***";
0160   }
0161 }
0162 
0163 void HcalGeometryTester::testTriggerGeometry(const HcalTopology& topology) {
0164   HcalTrigTowerGeometry trigTowers(&topology);
0165   edm::LogVerbatim("HCalGeom") << "HCAL trigger tower eta bounds:";
0166   for (int ieta = 1; ieta <= 32; ++ieta) {
0167     double eta1, eta2;
0168     trigTowers.towerEtaBounds(ieta, 0, eta1, eta2);
0169     edm::LogVerbatim("HCalGeom") << "[" << ieta << "] " << eta1 << " " << eta2;
0170   }
0171 
0172   // now test some cell mappings
0173   HcalDetId barrelDet(HcalBarrel, 1, 1, 1);
0174   HcalDetId endcapDet(HcalEndcap, 29, 1, 1);
0175   HcalDetId forwardDet1(HcalForward, 29, 71, 1);
0176   HcalDetId forwardDet2(HcalForward, 29, 71, 2);
0177   HcalDetId forwardDet3(HcalForward, 40, 71, 1);
0178 
0179   using TowerDets = std::vector<HcalTrigTowerDetId>;
0180   if (topology.valid(barrelDet)) {
0181     TowerDets barrelTowers = trigTowers.towerIds(barrelDet);
0182     edm::LogVerbatim("HCalGeom") << "Trigger Tower Size: Barrel " << barrelTowers.size();
0183     assert(barrelTowers.size() == 1);
0184     edm::LogVerbatim("HCalGeom") << "Tower[0] " << barrelTowers[0];
0185   }
0186   if (topology.valid(endcapDet)) {
0187     TowerDets endcapTowers = trigTowers.towerIds(endcapDet);
0188     edm::LogVerbatim("HCalGeom") << "Trigger Tower Size: Endcap " << endcapTowers.size();
0189     assert(!endcapTowers.empty());
0190     for (unsigned int k = 0; k < endcapTowers.size(); ++k)
0191       edm::LogVerbatim("HCalGeom") << "Tower[" << k << "] " << endcapTowers[k];
0192   }
0193   if (topology.valid(forwardDet1)) {
0194     TowerDets forwardTowers1 = trigTowers.towerIds(forwardDet1);
0195     edm::LogVerbatim("HCalGeom") << "Trigger Tower Size: Forward1 " << forwardTowers1.size();
0196     assert(!forwardTowers1.empty());
0197     for (unsigned int k = 0; k < forwardTowers1.size(); ++k)
0198       edm::LogVerbatim("HCalGeom") << "Tower[" << k << "] " << forwardTowers1[k];
0199   }
0200   if (topology.valid(forwardDet2)) {
0201     TowerDets forwardTowers2 = trigTowers.towerIds(forwardDet2);
0202     edm::LogVerbatim("HCalGeom") << "Trigger Tower Size: Forward2 " << forwardTowers2.size();
0203     assert(!forwardTowers2.empty());
0204     for (unsigned int k = 0; k < forwardTowers2.size(); ++k)
0205       edm::LogVerbatim("HCalGeom") << "Tower[" << k << "] " << forwardTowers2[k];
0206   }
0207   if (topology.valid(forwardDet3)) {
0208     TowerDets forwardTowers3 = trigTowers.towerIds(forwardDet3);
0209     edm::LogVerbatim("HCalGeom") << "Trigger Tower Size: Forward3 " << forwardTowers3.size();
0210     assert(!forwardTowers3.empty());
0211     for (unsigned int k = 0; k < forwardTowers3.size(); ++k)
0212       edm::LogVerbatim("HCalGeom") << "Tower[" << k << "] " << forwardTowers3[k];
0213   }
0214 }
0215 
0216 void HcalGeometryTester::testFlexiValidDetIds(CaloSubdetectorGeometry* caloGeom,
0217                                               const HcalTopology& topology,
0218                                               DetId::Detector det,
0219                                               int subdet,
0220                                               const std::string& label,
0221                                               std::vector<int>& dins) {
0222   std::stringstream s;
0223   s << label << " : " << std::endl;
0224   const std::vector<DetId>& idshb = caloGeom->getValidDetIds(det, subdet);
0225 
0226   int counter = 0;
0227   for (std::vector<DetId>::const_iterator i = idshb.begin(); i != idshb.end(); i++, ++counter) {
0228     HcalDetId hid = (*i);
0229     s << counter << ": din " << topology.detId2denseId(*i) << ":" << hid;
0230     dins.emplace_back(topology.detId2denseId(*i));
0231 
0232     auto cell = caloGeom->getGeometry(*i);
0233     s << *cell << std::endl;
0234   }
0235 
0236   std::sort(dins.begin(), dins.end());
0237   s << "=== Total " << counter << " cells in " << label << std::endl;
0238 
0239   counter = 0;
0240   for (std::vector<int>::const_iterator i = dins.begin(); i != dins.end(); ++i, ++counter) {
0241     HcalDetId hid = (topology.denseId2detId(*i));
0242     HcalDetId ihid = (topology.denseId2detId(dins[counter]));
0243     s << counter << ": din " << (*i) << " :" << hid << " == " << ihid << std::endl;
0244   }
0245   edm::LogVerbatim("HCalGeom") << s.str();
0246 }
0247 
0248 void HcalGeometryTester::testFlexiGeomHF(CaloSubdetectorGeometry* caloGeom) {
0249   std::stringstream s;
0250   s << "Test HF Geometry : " << std::endl;
0251   for (int ieta = 29; ieta <= 41; ++ieta) {
0252     HcalDetId cell3(HcalForward, ieta, 3, 1);
0253     auto cellGeometry3 = caloGeom->getGeometry(cell3);
0254     if (cellGeometry3) {
0255       s << "cell geometry iphi=3 -> ieta=" << ieta << " eta " << cellGeometry3->getPosition().eta() << "+-"
0256         << std::abs(cellGeometry3->getCorners()[0].eta() - cellGeometry3->getCorners()[2].eta()) / 2 << " phi "
0257         << cellGeometry3->getPosition().phi() / 3.1415 * 180 << "+-"
0258         << std::abs(cellGeometry3->getCorners()[0].phi() - cellGeometry3->getCorners()[2].phi()) / 3.1415 * 180 / 2.
0259         << std::endl;
0260     }
0261   }
0262   edm::LogVerbatim("HCalGeom") << s.str();
0263 }
0264 
0265 DEFINE_FWK_MODULE(HcalGeometryTester);