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::HepMC3Product>(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<HepMC3Product> srcHandle;
0123   event.getByToken(srcToken_, srcHandle);
0124
0125   const HepMC3::GenEventData* genEventData = srcHandle->GetEvent();
0126   std::unique_ptr<HepMC3::GenEvent> genEvent = std::make_unique<HepMC3::GenEvent>();
0127   genEvent->read_data(*genEventData);
0128
0129   if (!rivetAnalysis_ || !rivetAnalysis_->hasProjection("FS")) {
0130     rivetAnalysis_ = new Rivet::RivetAnalysis(pset_);
0131     analysisHandler_ = std::make_unique<Rivet::AnalysisHandler>();
0132
0133     analysisHandler_->setCheckBeams(false);
0134     analysisHandler_->addAnalysis(rivetAnalysis_);
0135   }
0136
0137   analysisHandler_->analyze(const_cast<GenEvent&>(*genEvent));
0138
0139   // Convert into edm objects
0140   // Prompt neutrinos
0141   for (auto const& p : rivetAnalysis_->neutrinos()) {
0142     neutrinos->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 1, true));
0143   }
0144   std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
0145
0146   // Photons
0147   for (auto const& p : rivetAnalysis_->photons()) {
0148     photons->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 1, true));
0149   }
0150   std::sort(photons->begin(), photons->end(), GreaterByPt<reco::Candidate>());
0151
0152   // Prompt leptons
0153   int iConstituent = -1;
0154   int iTag = -1;
0155   for (auto const& lepton : rivetAnalysis_->leptons()) {
0156     reco::GenJet lepJet;
0157     lepJet.setP4(p4(lepton));
0158     lepJet.setVertex(genVertex_);
0159     lepJet.setPdgId(lepton.pid());
0160     lepJet.setCharge(lepton.charge());
0161
0162     for (auto const& p : lepton.constituents()) {
0163       // ghost taus (momentum scaled with 1e-20 in RivetAnalysis.h already)
0164       if (p.abspid() == 15) {
0165         tags->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 2, true));
0166         lepJet.addDaughter(edm::refToPtr(reco::GenParticleRef(tagsRefHandle, ++iTag)));
0167       }
0168       // electrons, muons, photons
0169       else {
0170         consts->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pid(), 1, true));
0171         lepJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, ++iConstituent)));
0172       }
0173     }
0174
0175     leptons->push_back(lepJet);
0176   }
0177   std::sort(leptons->begin(), leptons->end(), GreaterByPt<reco::GenJet>());
0178
0179   // Jets with constituents and tag particles
0180   for (const auto& jet : rivetAnalysis_->jets()) {
0181     addGenJet(jet, jets, consts, constsRefHandle, iConstituent, tags, tagsRefHandle, iTag);
0182   }
0183   for (const auto& jet : rivetAnalysis_->fatjets()) {
0184     addGenJet(jet, fatjets, consts, constsRefHandle, iConstituent, tags, tagsRefHandle, iTag);
0185   }
0186
0187   // MET
0188   reco::Candidate::LorentzVector metP4(rivetAnalysis_->met().x(),
0189                                        rivetAnalysis_->met().y(),
0190                                        0.,
0191                                        sqrt(pow(rivetAnalysis_->met().x(), 2) + pow(rivetAnalysis_->met().y(), 2)));
0192   mets->push_back(reco::MET(metP4, genVertex_));
0193
0194   event.put(std::move(neutrinos), "neutrinos");
0195   event.put(std::move(photons), "photons");
0196   event.put(std::move(leptons), "leptons");
0197   event.put(std::move(jets), "jets");
0198   event.put(std::move(fatjets), "fatjets");
0199   event.put(std::move(consts), "consts");
0200   event.put(std::move(tags), "tags");
0201   event.put(std::move(mets), "mets");
0202 }
0203
0204 DEFINE_FWK_MODULE(ParticleLevelProducer);