File indexing completed on 2024-04-06 12:25:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
0018 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0019 #include "DataFormats/Candidate/interface/Particle.h"
0020 #include "DataFormats/Common/interface/Handle.h"
0021 #include "DataFormats/Common/interface/SortedCollection.h"
0022 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
0028 #include "DataFormats/Math/interface/deltaPhi.h"
0029 #include "DataFormats/DetId/interface/DetId.h"
0030 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0031 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0032 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0033 #include "FWCore/Framework/interface/stream/EDProducer.h"
0034 #include "FWCore/Framework/interface/EventSetup.h"
0035 #include "FWCore/Utilities/interface/EDGetToken.h"
0036 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0037 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0038 #include "RecoHI/HiJetAlgos/plugins/HITowerHelper.h"
0039
0040 #include <cmath>
0041 #include <cstdlib>
0042 #include <memory>
0043 #include <string>
0044 #include <vector>
0045 #include <map>
0046 #include <utility>
0047
0048 template <class T>
0049 class ParticleTowerProducer : public edm::stream::EDProducer<> {
0050 public:
0051 explicit ParticleTowerProducer(const edm::ParameterSet&);
0052 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0053
0054 private:
0055 void produce(edm::Event&, const edm::EventSetup&) override;
0056 int eta2ieta(double eta) const;
0057 int phi2iphi(double phi, int ieta) const;
0058 double ieta2eta(int ieta) const;
0059 double iphi2phi(int iphi, int ieta) const;
0060
0061
0062 edm::EDGetTokenT<typename edm::View<T> > src_;
0063 const bool useHF_;
0064
0065
0066 static constexpr int ietaMax = 42;
0067 };
0068
0069
0070
0071 template <class T>
0072 ParticleTowerProducer<T>::ParticleTowerProducer(const edm::ParameterSet& iConfig)
0073 : src_(consumes<typename edm::View<T> >(iConfig.getParameter<edm::InputTag>("src"))),
0074 useHF_(iConfig.getParameter<bool>("useHF")) {
0075 produces<CaloTowerCollection>();
0076 }
0077
0078
0079
0080
0081
0082
0083 template <class T>
0084 void ParticleTowerProducer<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0085 using namespace edm;
0086
0087 typedef std::pair<int, int> EtaPhi;
0088 typedef std::map<EtaPhi, double> EtaPhiMap;
0089 EtaPhiMap towers_;
0090 towers_.clear();
0091
0092 auto const& inputs = iEvent.get(src_);
0093 for (auto const& particle : inputs) {
0094 double eta = particle.eta();
0095
0096 int ieta = eta2ieta(eta);
0097 int iphi = phi2iphi(particle.phi(), ieta);
0098
0099 if (!useHF_ && abs(ieta) > 29)
0100 continue;
0101
0102 EtaPhi ep(ieta, iphi);
0103 towers_[ep] += particle.et();
0104 }
0105
0106 auto prod = std::make_unique<CaloTowerCollection>();
0107 prod->reserve(towers_.size());
0108 for (auto const& tower : towers_) {
0109 EtaPhi ep = tower.first;
0110 double et = tower.second;
0111
0112 int ieta = ep.first;
0113 int iphi = ep.second;
0114
0115 CaloTowerDetId newTowerId(ieta, iphi);
0116
0117
0118
0119 reco::Particle::PolarLorentzVector p4(et, ieta2eta(ieta), iphi2phi(iphi, ieta), 0.);
0120 GlobalPoint point(p4.x(), p4.y(), p4.z());
0121 prod->emplace_back(newTowerId, p4.e(), 0, 0, 0, 0, p4, point, point);
0122 }
0123
0124 iEvent.put(std::move(prod));
0125 }
0126
0127
0128
0129 template <class T>
0130 int ParticleTowerProducer<T>::eta2ieta(double eta) const {
0131
0132
0133 int ieta = 1;
0134 double xeta = fabs(eta);
0135 while (xeta > hi::etaedge[ieta] && ieta < ietaMax - 1) {
0136 ++ieta;
0137 }
0138
0139 if (eta < 0)
0140 ieta = -ieta;
0141 return ieta;
0142 }
0143
0144 template <class T>
0145 int ParticleTowerProducer<T>::phi2iphi(double phi, int ieta) const {
0146 phi = angle0to2pi::make0To2pi(phi);
0147 int nphi = 72;
0148 int n = 1;
0149 if (abs(ieta) > 20)
0150 n = 2;
0151 if (abs(ieta) >= 40)
0152 n = 4;
0153
0154 int iphi = (int)std::ceil(phi / 2.0 / M_PI * nphi / n);
0155
0156 iphi = n * (iphi - 1) + 1;
0157
0158 return iphi;
0159 }
0160
0161 template <class T>
0162 double ParticleTowerProducer<T>::iphi2phi(int iphi, int ieta) const {
0163 double phi = 0;
0164 int nphi = 72;
0165
0166 int n = 1;
0167 if (abs(ieta) > 20)
0168 n = 2;
0169 if (abs(ieta) >= 40)
0170 n = 4;
0171
0172 int myphi = (iphi - 1) / n + 1;
0173
0174 phi = 2. * M_PI * (myphi - 0.5) / nphi * n;
0175 while (phi > M_PI)
0176 phi -= 2. * M_PI;
0177
0178 return phi;
0179 }
0180
0181 template <class T>
0182 double ParticleTowerProducer<T>::ieta2eta(int ieta) const {
0183 int sign = 1;
0184 if (ieta < 0) {
0185 sign = -1;
0186 ieta = -ieta;
0187 }
0188
0189 double eta = sign * (hi::etaedge[ieta] + hi::etaedge[ieta - 1]) / 2.;
0190 return eta;
0191 }
0192
0193 template <class T>
0194 void ParticleTowerProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0195 edm::ParameterSetDescription desc;
0196 desc.add<bool>("useHF", true);
0197 if (typeid(T) == typeid(reco::PFCandidate)) {
0198 desc.add<edm::InputTag>("src", edm::InputTag("particleFlow"));
0199 descriptions.add("PFTowers", desc);
0200 }
0201 if (typeid(T) == typeid(pat::PackedCandidate)) {
0202 desc.add<edm::InputTag>("src", edm::InputTag("packedPFCandidates"));
0203 descriptions.add("PackedPFTowers", desc);
0204 }
0205 }
0206
0207 using PFTowers = ParticleTowerProducer<reco::PFCandidate>;
0208 using PackedPFTowers = ParticleTowerProducer<pat::PackedCandidate>;
0209
0210
0211 DEFINE_FWK_MODULE(PFTowers);
0212 DEFINE_FWK_MODULE(PackedPFTowers);