HiJetAlgos/plugins/ParticleTowerProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoHI/HiJetAlgos
0004 // Class:      ParticleTowerProducer
0005 //
0006 /**\class ParticleTowerProducer RecoHI/HiJetAlgos/plugins/ParticleTowerProducer.cc
0007  Description: [one line class summary]
0008  Implementation:
0009      [Notes on implementation]
0010 */
0011 //
0012 // Original Author:  Yetkin Yilmaz,32 4-A08,+41227673039,
0013 //         Created:  Thu Jan 20 19:53:58 CET 2011
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   // ----------member data ---------------------------
0061
0062   edm::EDGetTokenT<typename edm::View<T> > src_;
0063   const bool useHF_;
0064
0065   // tower edges from fast sim, used starting at index 30 for the HF
0066   static constexpr int ietaMax = 42;
0067 };
0068 //
0069 // constructors and destructor
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 // member functions
0080 //
0081
0082 // ------------ method called to produce the data  ------------
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);  // totally dummy id
0116
0117     // currently sets et = pt, mass to zero
0118     // pt, eta, phi, mass
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 // Taken from FastSimulation/CalorimeterProperties/src/HCALProperties.cc
0128 // Note this returns an abs(ieta)
0129 template <class T>
0130 int ParticleTowerProducer<T>::eta2ieta(double eta) const {
0131   // binary search in the array of towers eta edges
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 // define this as a plug-in
0211 DEFINE_FWK_MODULE(PFTowers);
0212 DEFINE_FWK_MODULE(PackedPFTowers);