Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:48:50

0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004 
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 
0014 #include "DataFormats/DetId/interface/DetId.h"
0015 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0016 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0017 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0018 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0019 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0020 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0021 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0022 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0023 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0026 #include "CoralBase/Exception.h"
0027 
0028 class HGCalTestRecHitTool : public edm::one::EDAnalyzer<> {
0029 public:
0030   explicit HGCalTestRecHitTool(const edm::ParameterSet&);
0031   ~HGCalTestRecHitTool() override;
0032 
0033   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0034 
0035 private:
0036   std::vector<double> retrieveLayerPositions(unsigned int);
0037   template <typename GEOM>
0038   void check_geom(const GEOM* geom) const;
0039   template <typename DDD>
0040   void check_ddd(const DDD* ddd) const;
0041   const HGCalDDDConstants* get_ddd(const DetId& detid) const;
0042   const HcalDDDRecConstants* get_ddd(const HcalDetId& detid) const;
0043   double getLayerZ(DetId const&) const;
0044   double getLayerZ(int type, int layer) const;
0045   GlobalPoint getPosition(DetId const&) const;
0046 
0047   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0048   const CaloGeometry* geom_;
0049   const int mode_;
0050   int layerEE_, layerFH_, layerBH_;
0051   int layerEE1000_, layerFH1000_, layerBH1000_;
0052 };
0053 
0054 HGCalTestRecHitTool::HGCalTestRecHitTool(const edm::ParameterSet& iC)
0055     : geomToken_{esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})},
0056       geom_(nullptr),
0057       mode_(iC.getParameter<int>("Mode")) {
0058   layerEE_ = layerFH_ = layerBH_ = 0;
0059   layerEE1000_ = layerFH1000_ = layerBH1000_ = 0;
0060 }
0061 
0062 HGCalTestRecHitTool::~HGCalTestRecHitTool() {}
0063 
0064 void HGCalTestRecHitTool::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0065   if (auto pG = iSetup.getHandle(geomToken_)) {
0066     geom_ = pG.product();
0067     auto geomEE = ((mode_ == 0) ? static_cast<const HGCalGeometry*>(
0068                                       geom_->getSubdetectorGeometry(DetId::Forward, ForwardSubdetector::HGCEE))
0069                                 : static_cast<const HGCalGeometry*>(
0070                                       geom_->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty)));
0071     layerEE_ = (geomEE->topology().dddConstants()).layers(true);
0072     layerEE1000_ = (geomEE->topology().dddConstants()).getLayer(10000., true);
0073     auto geomFH = ((mode_ == 0) ? static_cast<const HGCalGeometry*>(
0074                                       geom_->getSubdetectorGeometry(DetId::Forward, ForwardSubdetector::HGCHEF))
0075                                 : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0076                                       DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty)));
0077     layerFH_ = (geomFH->topology().dddConstants()).layers(true);
0078     layerFH1000_ = (geomFH->topology().dddConstants()).getLayer(10000., true);
0079     if (mode_ == 0) {
0080       auto geomBH =
0081           static_cast<const HcalGeometry*>(geom_->getSubdetectorGeometry(DetId::Hcal, HcalSubdetector::HcalEndcap));
0082       layerBH_ = (geomBH->topology().dddConstants())->getMaxDepth(1);
0083     } else {
0084       auto geomBH = static_cast<const HGCalGeometry*>(
0085           geom_->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0086       layerBH_ = (geomBH->topology().dddConstants()).layers(true);
0087       layerBH1000_ = (geomBH->topology().dddConstants()).getLayer(10000., true);
0088     }
0089     edm::LogVerbatim("HGCalGeom") << "Layers " << layerEE_ << ":" << layerFH_ << ":" << layerBH_
0090                                   << "\nLayer # at 1000 cm " << layerEE1000_ << ":" << layerFH1000_ << ":"
0091                                   << layerBH1000_;
0092     for (int layer = 1; layer <= layerEE_; ++layer)
0093       edm::LogVerbatim("HGCalGeom") << "EE Layer " << layer << " Wafers "
0094                                     << (geomEE->topology().dddConstants()).wafers(layer, 0) << ":"
0095                                     << (geomEE->topology().dddConstants()).wafers(layer, 1) << ":"
0096                                     << (geomEE->topology().dddConstants()).wafers(layer, 2) << std::endl;
0097     for (int layer = 1; layer <= layerFH_; ++layer)
0098       edm::LogVerbatim("HGCalGeom") << "FH Layer " << layer << " Wafers "
0099                                     << (geomFH->topology().dddConstants()).wafers(layer, 0) << ":"
0100                                     << (geomFH->topology().dddConstants()).wafers(layer, 1) << ":"
0101                                     << (geomFH->topology().dddConstants()).wafers(layer, 2) << std::endl;
0102     int nlayer = ((mode_ == 0) ? (layerEE_ + layerFH_ + layerBH_) : (layerEE_ + layerFH_));
0103     retrieveLayerPositions(nlayer);
0104   } else {
0105     edm::LogVerbatim("HGCalGeom") << "Cannot get valid CaloGeometry Object" << std::endl;
0106   }
0107 }
0108 
0109 template <typename GEOM>
0110 void HGCalTestRecHitTool::check_geom(const GEOM* geom) const {
0111   if (nullptr == geom) {
0112     throw cms::Exception("HGCalTestRecHitTools") << "Geometry not provided yet to HGCalTestRecHitTools!";
0113   }
0114 }
0115 
0116 template <typename DDD>
0117 void HGCalTestRecHitTool::check_ddd(const DDD* ddd) const {
0118   if (nullptr == ddd) {
0119     throw cms::Exception("HGCalTestRecHitTools") << "DDDConstants not accessibl to HGCalTestRecHitTools!";
0120   }
0121 }
0122 
0123 std::vector<double> HGCalTestRecHitTool::retrieveLayerPositions(unsigned layers) {
0124   DetId id;
0125   std::vector<double> layerPositions(layers);
0126   for (int ilayer = 1; ilayer <= (int)(layers); ++ilayer) {
0127     int lay(ilayer), type(0);
0128     if (ilayer <= layerEE_) {
0129       id = ((mode_ == 0) ? static_cast<DetId>(HGCalDetId(ForwardSubdetector::HGCEE, 1, lay, 1, 2, 1))
0130                          : static_cast<DetId>(HGCSiliconDetId(DetId::HGCalEE, 1, 0, lay, 3, 3, 1, 1)));
0131     } else if (ilayer > layerEE_ && ilayer <= (layerEE_ + layerFH_)) {
0132       lay = ilayer - layerEE_;
0133       id = ((mode_ == 0) ? static_cast<DetId>(HGCalDetId(ForwardSubdetector::HGCHEF, 1, lay, 1, 2, 1))
0134                          : static_cast<DetId>(HGCSiliconDetId(DetId::HGCalHSi, 1, 0, lay, 3, 3, 1, 1)));
0135       type = 1;
0136     } else {
0137       lay = ilayer - (layerEE_ + layerFH_);
0138       id = HcalDetId(HcalSubdetector::HcalEndcap, 50, 100, lay);
0139       type = 2;
0140     }
0141     const GlobalPoint pos = getPosition(id);
0142     if (id.det() == DetId::Hcal) {
0143       edm::LogVerbatim("HGCalGeom") << "GEOM  layer " << ilayer << " ID " << HcalDetId(id);
0144     } else if (id.det() == DetId::Forward) {
0145       edm::LogVerbatim("HGCalGeom") << "GEOM  layer " << ilayer << " ID " << HGCalDetId(id);
0146     } else if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) {
0147       edm::LogVerbatim("HGCalGeom") << "GEOM  layer " << ilayer << " ID " << HGCSiliconDetId(id);
0148     } else {
0149       edm::LogVerbatim("HGCalGeom") << "GEOM  layer " << ilayer << " ID " << HGCScintillatorDetId(id);
0150     }
0151     edm::LogVerbatim("HGCalGeom") << " Z " << pos.z() << ":" << getLayerZ(id) << ":" << getLayerZ(type, lay)
0152                                   << std::endl;
0153     layerPositions[ilayer - 1] = pos.z();
0154   }
0155   return layerPositions;
0156 }
0157 
0158 const HcalDDDRecConstants* HGCalTestRecHitTool::get_ddd(const HcalDetId& detid) const {
0159   auto geom = geom_->getSubdetectorGeometry(detid);
0160   check_geom(geom);
0161   const HcalGeometry* hc = static_cast<const HcalGeometry*>(geom);
0162   const HcalDDDRecConstants* ddd = hc->topology().dddConstants();
0163   check_ddd(ddd);
0164   return ddd;
0165 }
0166 
0167 const HGCalDDDConstants* HGCalTestRecHitTool::get_ddd(const DetId& detid) const {
0168   auto geom = geom_->getSubdetectorGeometry(detid);
0169   check_geom(geom);
0170   const HGCalGeometry* hg = static_cast<const HGCalGeometry*>(geom);
0171   const HGCalDDDConstants* ddd = &(hg->topology().dddConstants());
0172   check_ddd(ddd);
0173   return ddd;
0174 }
0175 
0176 GlobalPoint HGCalTestRecHitTool::getPosition(const DetId& id) const {
0177   auto geom = geom_->getSubdetectorGeometry(id);
0178   check_geom(geom);
0179   GlobalPoint position;
0180   if (id.det() == DetId::Hcal) {
0181     position = geom->getGeometry(id)->getPosition();
0182   } else {
0183     position = (dynamic_cast<const HGCalGeometry*>(geom))->getPosition(id);
0184   }
0185   return position;
0186 }
0187 
0188 double HGCalTestRecHitTool::getLayerZ(const DetId& id) const {
0189   double zpos(0);
0190   if (id.det() == DetId::Hcal) {
0191     auto geom = geom_->getSubdetectorGeometry(id);
0192     check_geom(geom);
0193     zpos = (static_cast<const HcalGeometry*>(geom))->getGeometry(id)->getPosition().z();
0194   } else {
0195     const HGCalDDDConstants* ddd = get_ddd(id);
0196     int layer = ((id.det() == DetId::Forward) ? HGCalDetId(id).layer()
0197                                               : ((id.det() != DetId::HGCalHSc) ? HGCSiliconDetId(id).layer()
0198                                                                                : HGCScintillatorDetId(id).layer()));
0199     zpos = ddd->waferZ(layer, true);
0200   }
0201   return zpos;
0202 }
0203 
0204 double HGCalTestRecHitTool::getLayerZ(int type, int layer) const {
0205   double zpos(0);
0206   if (type == 2) {
0207     auto geom =
0208         static_cast<const HcalGeometry*>(geom_->getSubdetectorGeometry(DetId::Hcal, HcalSubdetector::HcalEndcap));
0209     check_geom(geom);
0210     auto ddd = (geom->topology().dddConstants());
0211     check_ddd(ddd);
0212     std::pair<int, int> etar = ddd->getEtaRange(1);
0213     zpos = ddd->getRZ(2, etar.second, layer);
0214   } else {
0215     auto geom = ((type == 1) ? ((mode_ == 0) ? static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0216                                                    DetId::Forward, ForwardSubdetector::HGCHEF))
0217                                              : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0218                                                    DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty)))
0219                              : ((mode_ == 0) ? static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0220                                                    DetId::Forward, ForwardSubdetector::HGCEE))
0221                                              : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0222                                                    DetId::HGCalEE, ForwardSubdetector::ForwardEmpty))));
0223     check_geom(geom);
0224     auto ddd = &(geom->topology().dddConstants());
0225     check_ddd(ddd);
0226     zpos = ddd->waferZ(layer, true);
0227   }
0228   return zpos;
0229 }
0230 
0231 //define this as a plug-in
0232 DEFINE_FWK_MODULE(HGCalTestRecHitTool);