RivetInterface/plugins/ParticleLevelProducer.cc

0001 #include <memory>
0002
0003 #include "GeneratorInterface/RivetInterface/interface/ParticleLevelProducer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "DataFormats/Candidate/interface/Candidate.h"
0007 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0008 #include "DataFormats/Common/interface/RefToPtr.h"
0009 #include "DataFormats/METReco/interface/METFwd.h"
0010 #include "DataFormats/METReco/interface/MET.h"
0011 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0012 #include "CommonTools/Utils/interface/PtComparator.h"
0013
0014 #include "Rivet/Analysis.hh"
0015
0016 using namespace std;
0017 using namespace edm;
0018 using namespace reco;
0019 using namespace Rivet;
0020
0021 ParticleLevelProducer::ParticleLevelProducer(const edm::ParameterSet& pset)
0022     : srcToken_(consumes<edm::HepMCProduct>(pset.getParameter<edm::InputTag>("src"))), pset_(pset) {
0023   usesResource("Rivet");
0024   genVertex_ = reco::Particle::Point(0, 0, 0);
0025
0026   produces<reco::GenParticleCollection>("neutrinos");
0027   produces<reco::GenParticleCollection>("photons");
0028   produces<reco::GenJetCollection>("leptons");
0029   produces<reco::GenJetCollection>("jets");
0030   produces<reco::GenJetCollection>("fatjets");
0031   produces<reco::GenParticleCollection>("consts");
0032   produces<reco::GenParticleCollection>("tags");
0033   produces<reco::METCollection>("mets");
0034 }
0035
0036 void ParticleLevelProducer::addGenJet(Rivet::Jet jet,
0037                                       std::unique_ptr<reco::GenJetCollection>& jets,
0038                                       std::unique_ptr<reco::GenParticleCollection>& consts,
0039                                       edm::RefProd<reco::GenParticleCollection>& constsRefHandle,
0040                                       int& iConstituent,
0041                                       std::unique_ptr<reco::GenParticleCollection>& tags,
0042                                       edm::RefProd<reco::GenParticleCollection>& tagsRefHandle,
0043                                       int& iTag) {
0044   const auto pjet = jet.pseudojet();
0045
0046   reco::GenJet genJet;
0047   genJet.setP4(p4(jet));
0048   genJet.setVertex(genVertex_);
0049   if (jet.bTagged())
0050     genJet.setPdgId(5);
0051   else if (jet.cTagged())
0052     genJet.setPdgId(4);
0053   genJet.setJetArea(pjet.has_area() ? pjet.area() : 0);
0054
0055   for (auto const& p : jet.particles()) {
0056     auto pp4 = p4(p);
0057     bool match = false;
0058     int iMatch = -1;
0059     for (auto const& q : *consts) {
0060       ++iMatch;
0061       if (q.p4() == pp4) {
0062         match = true;
0063         break;
0064       }
0065     }
0066     if (match) {
0067       genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, iMatch)));
0068     } else {
0069       consts->push_back(reco::GenParticle(p.charge(), pp4, genVertex_, p.pid(), 1, true));
0070       genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, ++iConstituent)));
0071     }
0072   }
0073   for (auto const& p : jet.tags()) {
0074     // The tag particles are accessible as jet daughters, so scale down p4 for safety.
0075     // p4 needs to be multiplied by 1e20 for fragmentation analysis.
0076     auto pp4 = p4(p) * 1e-20;
0077     bool match = false;
0078     int iMatch = -1;
0079     for (auto const& q : *tags) {
0080       ++iMatch;
0081       if (q.p4() == pp4) {
0082         match = true;
0083         break;
0084       }
0085     }
0086     if (match) {
0087       genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(tagsRefHandle, iMatch)));
0088     } else {
0089       tags->push_back(reco::GenParticle(p.charge(), p4(p) * 1e-20, genVertex_, p.pid(), 2, true));
0090       genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(tagsRefHandle, ++iTag)));
0091       // Also save lepton+neutrino daughters of tag particles
0092       int iTagMother = iTag;
0093       for (auto const& d : p.constituents()) {
0094         ++iTag;
0095         int d_status = (d.isStable()) ? 1 : 2;
0096         tags->push_back(reco::GenParticle(d.charge(), p4(d) * 1e-20, genVertex_, d.pid(), d_status, true));
0097         tags->at(iTag).addMother(reco::GenParticleRef(tagsRefHandle, iTagMother));
0098         tags->at(iTagMother).addDaughter(reco::GenParticleRef(tagsRefHandle, iTag));
0099         genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(tagsRefHandle, iTag)));
0100       }
0101     }
0102   }
0103
0104   jets->push_back(genJet);
0105 }
0106
0107 void ParticleLevelProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup) {
0108   using namespace Rivet;
0109   typedef reco::Candidate::LorentzVector LorentzVector;
0110
0111   std::unique_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
0112   std::unique_ptr<reco::GenParticleCollection> photons(new reco::GenParticleCollection);
0113   std::unique_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
0114   std::unique_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
0115   std::unique_ptr<reco::GenJetCollection> fatjets(new reco::GenJetCollection);
0116   std::unique_ptr<reco::GenParticleCollection> consts(new reco::GenParticleCollection);
0117   std::unique_ptr<reco::GenParticleCollection> tags(new reco::GenParticleCollection);
0118   std::unique_ptr<reco::METCollection> mets(new reco::METCollection);
0119   auto constsRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("consts");
0120   auto tagsRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("tags");
0121
0122   edm::Handle<HepMCProduct> srcHandle;
0123   event.getByToken(srcToken_, srcHandle);
0124
0125   const HepMC::GenEvent* genEvent = srcHandle->GetEvent();
0126
0127   if (!rivetAnalysis_ || !rivetAnalysis_->hasProjection("FS")) {
0128     rivetAnalysis_ = new Rivet::RivetAnalysis(pset_);
0129     analysisHandler_ = std::make_unique<Rivet::AnalysisHandler>();
0130
0131     analysisHandler_->setIgnoreBeams(true);
0132     analysisHandler_->addAnalysis(rivetAnalysis_);
0133   }
0134
0135   analysisHandler_->analyze(*genEvent);
0136
0137   // Convert into edm objects
0138   // Prompt neutrinos
0139   for (auto const& p : rivetAnalysis_->neutrinos()) {
0140     neutrinos->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 1, true));
0141   }
0142   std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
0143
0144   // Photons
0145   for (auto const& p : rivetAnalysis_->photons()) {
0146     photons->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 1, true));
0147   }
0148   std::sort(photons->begin(), photons->end(), GreaterByPt<reco::Candidate>());
0149
0150   // Prompt leptons
0151   int iConstituent = -1;
0152   int iTag = -1;
0153   for (auto const& lepton : rivetAnalysis_->leptons()) {
0154     reco::GenJet lepJet;
0155     lepJet.setP4(p4(lepton));
0156     lepJet.setVertex(genVertex_);
0157     lepJet.setPdgId(lepton.pid());
0158     lepJet.setCharge(lepton.charge());
0159
0160     for (auto const& p : lepton.constituents()) {
0161       // ghost taus (momentum scaled with 10e-20 in RivetAnalysis.h already)
0162       if (p.abspid() == 15) {
0163         tags->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 2, true));
0164         lepJet.addDaughter(edm::refToPtr(reco::GenParticleRef(tagsRefHandle, ++iTag)));
0165       }
0166       // electrons, muons, photons
0167       else {
0168         consts->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 1, true));
0169         lepJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, ++iConstituent)));
0170       }
0171     }
0172
0173     leptons->push_back(lepJet);
0174   }
0175   std::sort(leptons->begin(), leptons->end(), GreaterByPt<reco::GenJet>());
0176
0177   // Jets with constituents and tag particles
0178   for (const auto& jet : rivetAnalysis_->jets()) {
0179     addGenJet(jet, jets, consts, constsRefHandle, iConstituent, tags, tagsRefHandle, iTag);
0180   }
0181   for (const auto& jet : rivetAnalysis_->fatjets()) {
0182     addGenJet(jet, fatjets, consts, constsRefHandle, iConstituent, tags, tagsRefHandle, iTag);
0183   }
0184
0185   // MET
0186   reco::Candidate::LorentzVector metP4(rivetAnalysis_->met().x(),
0187                                        rivetAnalysis_->met().y(),
0188                                        0.,
0189                                        sqrt(pow(rivetAnalysis_->met().x(), 2) + pow(rivetAnalysis_->met().y(), 2)));
0190   mets->push_back(reco::MET(metP4, genVertex_));
0191
0192   event.put(std::move(neutrinos), "neutrinos");
0193   event.put(std::move(photons), "photons");
0194   event.put(std::move(leptons), "leptons");
0195   event.put(std::move(jets), "jets");
0196   event.put(std::move(fatjets), "fatjets");
0197   event.put(std::move(consts), "consts");
0198   event.put(std::move(tags), "tags");
0199   event.put(std::move(mets), "mets");
0200 }
0201
0202 DEFINE_FWK_MODULE(ParticleLevelProducer);