File indexing completed on 2024-04-06 11:59:07
0001
0002 #include <cmath>
0003 #include <memory>
0004 #include <sstream>
0005 #include <string>
0006 #include <vector>
0007
0008 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0009 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0010 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0011 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0012
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020
0021 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0022 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0023 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0024
0025
0026
0027 class HcalTestThreshold : public edm::one::EDAnalyzer<> {
0028 public:
0029 explicit HcalTestThreshold(edm::ParameterSet const&);
0030 ~HcalTestThreshold() override {}
0031
0032 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0033
0034 private:
0035 void analyze(edm::Event const&, edm::EventSetup const&) override;
0036 double eThreshold(const DetId& id, const CaloGeometry* geo) const;
0037
0038 const int etaMin_, etaMax_, phiValue_;
0039 const std::vector<int> ixNumbers_, iyNumbers_;
0040 const double hitEthrEB_, hitEthrEE0_, hitEthrEE1_;
0041 const double hitEthrEE2_, hitEthrEE3_;
0042 const double hitEthrEELo_, hitEthrEEHi_;
0043
0044 const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> ecalPFRecHitThresholdsToken_;
0045 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0046 };
0047
0048 HcalTestThreshold::HcalTestThreshold(const edm::ParameterSet& iConfig)
0049 : etaMin_(iConfig.getParameter<int>("etaMin")),
0050 etaMax_(iConfig.getParameter<int>("etaMax")),
0051 phiValue_(iConfig.getParameter<int>("phiValue")),
0052 ixNumbers_(iConfig.getParameter<std::vector<int>>("ixEENumbers")),
0053 iyNumbers_(iConfig.getParameter<std::vector<int>>("iyEENumbers")),
0054 hitEthrEB_(iConfig.getParameter<double>("EBHitEnergyThreshold")),
0055 hitEthrEE0_(iConfig.getParameter<double>("EEHitEnergyThreshold0")),
0056 hitEthrEE1_(iConfig.getParameter<double>("EEHitEnergyThreshold1")),
0057 hitEthrEE2_(iConfig.getParameter<double>("EEHitEnergyThreshold2")),
0058 hitEthrEE3_(iConfig.getParameter<double>("EEHitEnergyThreshold3")),
0059 hitEthrEELo_(iConfig.getParameter<double>("EEHitEnergyThresholdLow")),
0060 hitEthrEEHi_(iConfig.getParameter<double>("EEHitEnergyThresholdHigh")),
0061 ecalPFRecHitThresholdsToken_(esConsumes()),
0062 caloGeometryToken_(esConsumes()) {
0063 edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n"
0064 << "\tThreshold for EB " << hitEthrEB_ << "\t for EE " << hitEthrEE0_ << ":"
0065 << hitEthrEE1_ << ":" << hitEthrEE2_ << ":" << hitEthrEE3_ << ":" << hitEthrEELo_
0066 << ":" << hitEthrEEHi_;
0067 edm::LogVerbatim("HcalIsoTrack") << "\tRange in eta for EB " << etaMin_ << ":" << etaMax_ << "\tPhi value "
0068 << phiValue_;
0069 std::ostringstream st1, st2;
0070 st1 << "\t" << ixNumbers_.size() << " EE ix Numbers";
0071 for (const auto& ix : ixNumbers_)
0072 st1 << ": " << ix;
0073 edm::LogVerbatim("HcalIsoTrack") << st1.str();
0074 st2 << "\t" << iyNumbers_.size() << " EE iy Numbers";
0075 for (const auto& iy : iyNumbers_)
0076 st2 << ": " << iy;
0077 edm::LogVerbatim("HcalIsoTrack") << st2.str();
0078 }
0079
0080
0081 void HcalTestThreshold::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0082 edm::LogVerbatim("HcalIsoTrack") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event();
0083
0084 auto const& thresholds = iSetup.getData(ecalPFRecHitThresholdsToken_);
0085 auto const& geo = &iSetup.getData(caloGeometryToken_);
0086
0087
0088 edm::LogVerbatim("HcalIsoTrack") << "\n\nThresholds for EB Towers";
0089 for (int eta = etaMin_; eta < etaMax_; ++eta) {
0090 if (eta != 0) {
0091 EBDetId id(eta, phiValue_, EBDetId::ETAPHIMODE);
0092 edm::LogVerbatim("HcalIsoTrack") << id << " Eta " << (geo->getPosition(id)).eta() << " Thresholds:: old "
0093 << eThreshold(id, geo) << " new " << thresholds[id];
0094 }
0095 }
0096
0097
0098 #ifdef EDM_ML_DEBUG
0099 auto const& validId = geo->getValidDetIds(DetId::Ecal, 2);
0100 edm::LogVerbatim("HcalIsoTrack") << "\n\nList of " << validId.size() << " valid DetId's of EE";
0101 for (const auto& id : validId)
0102 edm::LogVerbatim("HcalIsoTrack") << EEDetId(id) << " SC " << EEDetId(id).isc() << " CR " << EEDetId(id).ic()
0103 << " Eta " << (geo->getPosition(id)).eta();
0104 #endif
0105 edm::LogVerbatim("HcalIsoTrack") << "\n\nThresholds for EE Towers";
0106 for (int zside = 0; zside < 2; ++zside) {
0107 int iz = 2 * zside - 1;
0108 for (const auto& ix : ixNumbers_) {
0109 for (const auto& iy : iyNumbers_) {
0110 EEDetId id(ix, iy, iz, EEDetId::XYMODE);
0111 edm::LogVerbatim("HcalIsoTrack") << id << " Eta " << (geo->getPosition(id)).eta() << " Thresholds:: old "
0112 << eThreshold(id, geo) << " new " << thresholds[id];
0113 }
0114 }
0115 }
0116 }
0117
0118 void HcalTestThreshold::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0119 edm::ParameterSetDescription desc;
0120 std::vector<int> ixNumb = {1, 10, 20, 30, 39};
0121 std::vector<int> iyNumb = {41, 43, 45, 47};
0122 desc.add<int>("etaMin", -85);
0123 desc.add<int>("etaMax", 85);
0124 desc.add<int>("phiValue", 1);
0125 desc.add<std::vector<int>>("ixEENumbers", ixNumb);
0126 desc.add<std::vector<int>>("iyEENumbers", iyNumb);
0127
0128 desc.add<double>("EBHitEnergyThreshold", 0.08);
0129 desc.add<double>("EEHitEnergyThreshold0", 0.30);
0130 desc.add<double>("EEHitEnergyThreshold1", 0.00);
0131 desc.add<double>("EEHitEnergyThreshold2", 0.00);
0132 desc.add<double>("EEHitEnergyThreshold3", 0.00);
0133 desc.add<double>("EEHitEnergyThresholdLow", 0.30);
0134 desc.add<double>("EEHitEnergyThresholdHigh", 0.30);
0135 descriptions.add("hcalTestThreshold", desc);
0136 }
0137
0138 double HcalTestThreshold::eThreshold(const DetId& id, const CaloGeometry* geo) const {
0139 const GlobalPoint& pos = geo->getPosition(id);
0140 double eta = std::abs(pos.eta());
0141 double eThr(hitEthrEB_);
0142 if (id.subdetId() != EcalBarrel) {
0143 eThr = (((eta * hitEthrEE3_ + hitEthrEE2_) * eta + hitEthrEE1_) * eta + hitEthrEE0_);
0144 if (eThr < hitEthrEELo_)
0145 eThr = hitEthrEELo_;
0146 else if (eThr > hitEthrEEHi_)
0147 eThr = hitEthrEEHi_;
0148 }
0149 return eThr;
0150 }
0151
0152
0153 DEFINE_FWK_MODULE(HcalTestThreshold);