HcalDetIdTester

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
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
#include "Geometry/HcalTowerAlgo/interface/HcalFlexiHardcodeGeometryLoader.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include <iostream>

class HcalDetIdTester : public edm::one::EDAnalyzer<> {
public:
  explicit HcalDetIdTester(const edm::ParameterSet&);
  ~HcalDetIdTester(void) override;

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

private:
  bool geomDB_;
  edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
  edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
};

HcalDetIdTester::HcalDetIdTester(const edm::ParameterSet& iConfig) {
  geomDB_ = iConfig.getParameter<bool>("GeometryFromDB");
  tok_ddrec_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>();
  tok_htopo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
  tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
}

HcalDetIdTester::~HcalDetIdTester(void) {}

void HcalDetIdTester::analyze(const edm::Event& /*iEvent*/, const edm::EventSetup& iSetup) {
  const HcalDDDRecConstants hcons = iSetup.getData(tok_ddrec_);
  const HcalTopology topology = iSetup.getData(tok_htopo_);

  CaloSubdetectorGeometry* caloGeom(nullptr);
  if (geomDB_) {
    const CaloGeometry* geo = &iSetup.getData(tok_geom_);
    caloGeom = (CaloSubdetectorGeometry*)(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
  } else {
    HcalFlexiHardcodeGeometryLoader m_loader;
    caloGeom = m_loader.load(topology, hcons);
  }

  int maxDepth = (topology.maxDepthHB() > topology.maxDepthHE()) ? topology.maxDepthHB() : topology.maxDepthHE();

  int nfail0(0);
  for (int det = 1; det <= HcalForward; det++) {
    for (int eta = -HcalDetId::kHcalEtaMask2; eta <= (int)(HcalDetId::kHcalEtaMask2); eta++) {
      for (int depth = 1; depth <= maxDepth; depth++) {
        for (unsigned int phi = 0; phi <= HcalDetId::kHcalPhiMask2; phi++) {
          HcalDetId detId((HcalSubdetector)det, eta, phi, depth);
          if (topology.valid(detId)) {
            auto cell = caloGeom->getGeometry((DetId)(detId));
            if (cell) {
              edm::LogVerbatim("HCalGeom") << detId << " " << cell->getPosition();
            } else {
              edm::LogVerbatim("HCalGeom") << detId << " position not found";
              ++nfail0;
            }
          }
        }
      }
    }
  }

  int nfail1(0);
  const std::vector<DetId>& ids = caloGeom->getValidDetIds();
  for (auto id : ids) {
    if (!topology.valid(id)) {
      edm::LogVerbatim("HCalGeom") << HcalDetId(id) << " declared as invalid == ERROR";
      ++nfail1;
    }
  }

  edm::LogVerbatim("HCalGeom") << "\nNumber of failures:\nTopology certifies but geometry fails " << nfail0
                               << "\nGeometry certifies but Topology fails " << nfail1 << std::endl;
}

DEFINE_FWK_MODULE(HcalDetIdTester);