HcalCalibAlgos/test/HcalTestThreshold.cc

0001 // system include files
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 //#define EDM_ML_DEBUG
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 // ------------ method called when starting to processes a run  ------------
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   // First EB
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   // Next EE
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   // energy thershold for ECAL (from Egamma group)
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 //define this as a plug-in
0153 DEFINE_FWK_MODULE(HcalTestThreshold);