PileupAlgos/plugins/SoftKillerProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CommonTools/PileupAlgos
0004 // Class:      SoftKillerProducer
0005 //
0006 /**\class SoftKillerProducer SoftKillerProducer.cc CommonTools/PileupAlgos/plugins/SoftKillerProducer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Salvatore Rappoccio
0015 //         Created:  Wed, 22 Oct 2014 15:14:20 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <vector>
0022
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/stream/EDProducer.h"
0026
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "DataFormats/Common/interface/View.h"
0032 #include "DataFormats/Candidate/interface/Candidate.h"
0033 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0034 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0035 #include "DataFormats/Common/interface/ValueMap.h"
0036 #include "fastjet/contrib/SoftKiller.hh"
0037
0038 //
0039 // class declaration
0040 //
0041
0042 class SoftKillerProducer : public edm::stream::EDProducer<> {
0043 public:
0044   typedef math::XYZTLorentzVector LorentzVector;
0045   typedef std::vector<LorentzVector> LorentzVectorCollection;
0046   typedef std::vector<reco::PFCandidate> PFOutputCollection;
0047
0048   explicit SoftKillerProducer(const edm::ParameterSet&);
0049   ~SoftKillerProducer() override;
0050
0051 private:
0052   void produce(edm::Event&, const edm::EventSetup&) override;
0053
0054   edm::EDGetTokenT<reco::CandidateView> tokenPFCandidates_;
0055
0056   double Rho_EtaMax_;
0057   double rParam_;
0058 };
0059
0060 //
0061 // constants, enums and typedefs
0062 //
0063
0064 //
0065 // static data member definitions
0066 //
0067
0068 //
0069 // constructors and destructor
0070 //
0071 SoftKillerProducer::SoftKillerProducer(const edm::ParameterSet& iConfig)
0072     : Rho_EtaMax_(iConfig.getParameter<double>("Rho_EtaMax")), rParam_(iConfig.getParameter<double>("rParam")) {
0073   tokenPFCandidates_ = consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("PFCandidates"));
0074
0075   produces<edm::ValueMap<LorentzVector> >("SoftKillerP4s");
0076   produces<PFOutputCollection>();
0077 }
0078
0079 SoftKillerProducer::~SoftKillerProducer() {}
0080
0081 //
0082 // member functions
0083 //
0084
0085 // ------------ method called to produce the data  ------------
0086 void SoftKillerProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0087   std::unique_ptr<PFOutputCollection> pOutput(new PFOutputCollection);
0088
0089   // get PF Candidates
0090   edm::Handle<reco::CandidateView> pfCandidates;
0091   iEvent.getByToken(tokenPFCandidates_, pfCandidates);
0092
0093   std::vector<fastjet::PseudoJet> fjInputs;
0094   for (auto i = pfCandidates->begin(), ibegin = pfCandidates->begin(), iend = pfCandidates->end(); i != iend; ++i) {
0095     fjInputs.push_back(fastjet::PseudoJet(i->px(), i->py(), i->pz(), i->energy()));
0096     fjInputs.back().set_user_index(i - ibegin);
0097   }
0098
0099   // soft killer:
0100   fastjet::contrib::SoftKiller soft_killer(Rho_EtaMax_, rParam_);
0101
0102   double pt_threshold = 0.;
0103   std::vector<fastjet::PseudoJet> soft_killed_event;
0104   soft_killer.apply(fjInputs, soft_killed_event, pt_threshold);
0105
0106   std::unique_ptr<edm::ValueMap<LorentzVector> > p4SKOut(new edm::ValueMap<LorentzVector>());
0107   LorentzVectorCollection skP4s;
0108
0109   static const reco::PFCandidate dummySinceTranslateIsNotStatic;
0110
0111   // To satisfy the value map, the size of the "killed" collection needs to be the
0112   // same size as the input collection, so if the constituent is killed, just set E = 0
0113   for (auto j = fjInputs.begin(), jend = fjInputs.end(); j != jend; ++j) {
0114     const reco::Candidate& cand = pfCandidates->at(j->user_index());
0115     auto id = dummySinceTranslateIsNotStatic.translatePdgIdToType(cand.pdgId());
0116     const reco::PFCandidate* pPF = dynamic_cast<const reco::PFCandidate*>(&cand);
0117     reco::PFCandidate pCand(pPF ? *pPF : reco::PFCandidate(cand.charge(), cand.p4(), id));
0118     auto val = j->user_index();
0119     auto skmatch = find_if(soft_killed_event.begin(), soft_killed_event.end(), [&val](fastjet::PseudoJet const& i) {
0120       return i.user_index() == val;
0121     });
0122     LorentzVector pVec;
0123     if (skmatch != soft_killed_event.end()) {
0124       pVec.SetPxPyPzE(skmatch->px(), skmatch->py(), skmatch->pz(), skmatch->E());
0125     } else {
0126       pVec.SetPxPyPzE(0., 0., 0., 0.);
0127     }
0128     pCand.setP4(pVec);
0129     skP4s.push_back(pVec);
0130     pOutput->push_back(pCand);
0131   }
0132
0133   //Compute the modified p4s
0134   edm::ValueMap<LorentzVector>::Filler p4SKFiller(*p4SKOut);
0135   p4SKFiller.insert(pfCandidates, skP4s.begin(), skP4s.end());
0136   p4SKFiller.fill();
0137
0138   iEvent.put(std::move(p4SKOut), "SoftKillerP4s");
0139   iEvent.put(std::move(pOutput));
0140 }
0141
0142 //define this as a plug-in
0143 DEFINE_FWK_MODULE(SoftKillerProducer);