Merge pull request #18402 from intrepid42/ParticleLevelProducer_80X

ParticleLevelProducer based on Rivet

GeneratorInterface/RivetInterface/interface/ParticleLevelProducer.h

@@ -0,0 +1,44 @@
+#ifndef GeneratorInterface_RivetInterface_ParticleLevelProducer_H
+#define GeneratorInterface_RivetInterface_ParticleLevelProducer_H
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/one/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "DataFormats/JetReco/interface/GenJetCollection.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
+
+#include "Rivet/AnalysisHandler.hh"
+#include "GeneratorInterface/RivetInterface/interface/RivetAnalysis.h"
+
+class ParticleLevelProducer : public edm::one::EDProducer<edm::one::SharedResources>
+{
+public:
+  ParticleLevelProducer(const edm::ParameterSet& pset);
+  virtual ~ParticleLevelProducer() {}
+  void produce(edm::Event& event, const edm::EventSetup& eventSetup) override;
+
+private:
+  void addGenJet(Rivet::Jet jet, std::unique_ptr<reco::GenJetCollection> &jets,
+                 std::unique_ptr<reco::GenParticleCollection> &consts, edm::RefProd<reco::GenParticleCollection>& constsRefHandle, int &iConstituent,
+                 std::unique_ptr<reco::GenParticleCollection> &tags, edm::RefProd<reco::GenParticleCollection>& tagsRefHandle, int &iTag);
+
+  template<typename T> reco::Candidate::LorentzVector p4(const T& p) const
+  {
+    return reco::Candidate::LorentzVector(p.px(), p.py(), p.pz(), p.energy());
+  }
+
+  const edm::EDGetTokenT<edm::HepMCProduct> srcToken_;
+
+  reco::Particle::Point genVertex_;
+
+  Rivet::RivetAnalysis* rivetAnalysis_;
+	Rivet::AnalysisHandler analysisHandler_;
+
+};
+
+#endif

GeneratorInterface/RivetInterface/interface/RivetAnalysis.h

@@ -0,0 +1,169 @@
+#ifndef GeneratorInterface_RivetInterface_RivetAnalysis_H
+#define GeneratorInterface_RivetInterface_RivetAnalysis_H
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "Rivet/Analysis.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Particle.hh"
+#include "Rivet/Particle.fhh"
+#include "Rivet/Event.hh"
+#include "Rivet/Projections/FastJets.hh"
+#include "Rivet/Projections/JetAlg.hh"
+#include "Rivet/Projections/ChargedLeptons.hh"
+#include "Rivet/Projections/PromptFinalState.hh"
+#include "Rivet/Projections/DressedLeptons.hh"
+#include "Rivet/Projections/VetoedFinalState.hh"
+#include "Rivet/Projections/IdentifiedFinalState.hh"
+#include "Rivet/Projections/MissingMomentum.hh"
+#include "Rivet/Tools/RivetHepMC.hh"
+
+namespace Rivet {
+
+  class RivetAnalysis : public Analysis {
+
+    public:
+      std::vector<DressedLepton> leptons() const {return _leptons;}
+      ParticleVector photons() const {return _photons;}
+      ParticleVector neutrinos() const {return _neutrinos;}
+      Jets jets() const {return _jets;}
+      Jets fatjets() const {return _fatjets;}
+      Vector3 met() const {return _met;}
+
+    private:
+      bool _usePromptFinalStates;
+      bool _excludePromptLeptonsFromJetClustering;
+      bool _excludeNeutrinosFromJetClustering;
+
+      double _particleMinPt, _particleMaxEta;
+      double _lepConeSize, _lepMinPt, _lepMaxEta;
+      double _jetConeSize, _jetMinPt, _jetMaxEta;
+      double _fatJetConeSize, _fatJetMinPt, _fatJetMaxEta;
+
+      std::vector<DressedLepton> _leptons;
+      ParticleVector _photons, _neutrinos;
+      Jets _jets, _fatjets;
+      Vector3 _met;
+
+    public:
+      RivetAnalysis(const edm::ParameterSet& pset) : Analysis("RivetAnalysis"),
+      _usePromptFinalStates(pset.getParameter<bool>("usePromptFinalStates")),
+      _excludePromptLeptonsFromJetClustering(pset.getParameter<bool>("excludePromptLeptonsFromJetClustering")),
+      _excludeNeutrinosFromJetClustering(pset.getParameter<bool>("excludeNeutrinosFromJetClustering")),
+
+      _particleMinPt  (pset.getParameter<double>("particleMinPt")),
+      _particleMaxEta (pset.getParameter<double>("particleMaxEta")),
+
+      _lepConeSize (pset.getParameter<double>("lepConeSize")),
+      _lepMinPt    (pset.getParameter<double>("lepMinPt")),
+      _lepMaxEta   (pset.getParameter<double>("lepMaxEta")),
+
+      _jetConeSize (pset.getParameter<double>("jetConeSize")),
+      _jetMinPt    (pset.getParameter<double>("jetMinPt")),
+      _jetMaxEta   (pset.getParameter<double>("jetMaxEta")),
+
+      _fatJetConeSize (pset.getParameter<double>("fatJetConeSize")),
+      _fatJetMinPt    (pset.getParameter<double>("fatJetMinPt")),
+      _fatJetMaxEta   (pset.getParameter<double>("fatJetMaxEta"))
+      {
+      }
+
+      // Initialize Rivet projections
+      void init() {
+        // Cuts
+        Cut particle_cut = (Cuts::abseta < _particleMaxEta) and (Cuts::pT > _particleMinPt*GeV);
+        Cut lepton_cut   = (Cuts::abseta < _lepMaxEta)      and (Cuts::pT > _lepMinPt*GeV);
+
+        // Generic final state
+        FinalState fs(particle_cut);
+
+        // Dressed leptons
+        ChargedLeptons charged_leptons(fs);
+        IdentifiedFinalState photons(fs);
+        photons.acceptIdPair(PID::PHOTON);
+
+        PromptFinalState prompt_leptons(charged_leptons);
+        prompt_leptons.acceptMuonDecays(true);
+        prompt_leptons.acceptTauDecays(true);
+
+        PromptFinalState prompt_photons(photons);
+        prompt_photons.acceptMuonDecays(true);
+        prompt_photons.acceptTauDecays(true);
+
+        // useDecayPhotons=true allows for photons with tau ancestor,
+        // photons from hadrons are vetoed by the PromptFinalState;
+        // will be default DressedLeptons behaviour for Rivet >= 2.5.4
+        DressedLeptons dressed_leptons(prompt_photons, prompt_leptons, _lepConeSize, 
+                       lepton_cut, /*cluster*/ true, /*useDecayPhotons*/ true);
+        if (not _usePromptFinalStates)
+          dressed_leptons = DressedLeptons(photons, charged_leptons, _lepConeSize, 
+                            lepton_cut, /*cluster*/ true, /*useDecayPhotons*/ true);
+        addProjection(dressed_leptons, "DressedLeptons");
+
+        // Photons
+        if (_usePromptFinalStates) {
+          // We remove the photons used up for lepton dressing in this case
+          VetoedFinalState vetoed_prompt_photons(prompt_photons);
+          vetoed_prompt_photons.addVetoOnThisFinalState(dressed_leptons);
+          addProjection(vetoed_prompt_photons, "Photons");
+        }
+        else
+          addProjection(photons, "Photons");
+
+        // Jets
+        VetoedFinalState fsForJets(fs);
+        if (_usePromptFinalStates and _excludePromptLeptonsFromJetClustering)
+          fsForJets.addVetoOnThisFinalState(dressed_leptons);
+        JetAlg::InvisiblesStrategy invisiblesStrategy = JetAlg::DECAY_INVISIBLES;
+        if (_excludeNeutrinosFromJetClustering)
+          invisiblesStrategy = JetAlg::NO_INVISIBLES;
+        addProjection(FastJets(fsForJets, FastJets::ANTIKT, _jetConeSize,
+                               JetAlg::ALL_MUONS, invisiblesStrategy), "Jets");
+
+        // FatJets
+        addProjection(FastJets(fsForJets, FastJets::ANTIKT, _fatJetConeSize), "FatJets");
+
+        // Neutrinos
+        IdentifiedFinalState neutrinos(fs);
+        neutrinos.acceptNeutrinos();
+        if (_usePromptFinalStates) {
+          PromptFinalState prompt_neutrinos(neutrinos);
+          prompt_neutrinos.acceptMuonDecays(true);
+          prompt_neutrinos.acceptTauDecays(true);
+          addProjection(prompt_neutrinos, "Neutrinos");
+        }
+        else
+          addProjection(neutrinos, "Neutrinos");
+
+        // MET
+        addProjection(MissingMomentum(fs), "MET");
+      };
+
+      // Apply Rivet projections
+      void analyze(const Event& event) {
+        _jets.clear();
+        _fatjets.clear();
+        _leptons.clear();
+        _photons.clear();
+        _neutrinos.clear();
+
+        // Get analysis objects from projections
+        Cut jet_cut    = (Cuts::abseta < _jetMaxEta)    and (Cuts::pT > _jetMinPt*GeV);
+        Cut fatjet_cut = (Cuts::abseta < _fatJetMaxEta) and (Cuts::pT > _fatJetMinPt*GeV);
+
+        _leptons   = applyProjection<DressedLeptons>(event, "DressedLeptons").dressedLeptons();
+        _jets      = applyProjection<FastJets>(event, "Jets").jetsByPt(jet_cut);
+        _fatjets   = applyProjection<FastJets>(event, "FatJets").jetsByPt(fatjet_cut);
+        _photons   = applyProjection<FinalState>(event, "Photons").particlesByPt();
+        _neutrinos = applyProjection<FinalState>(event, "Neutrinos").particlesByPt();
+        _met       = -applyProjection<MissingMomentum>(event, "MET").vectorEt();
+      };
+
+      // Do nothing here
+      void finalize() {};
+
+  };
+
+}
+
+#endif

GeneratorInterface/RivetInterface/plugins/BuildFile.xml

@@ -9,6 +9,8 @@
 <use   name="DQMServices/Core"/>
 <use   name="SimGeneral/HepPDTRecord"/>
 <use   name="DataFormats/HepMCCandidate"/>
+<use   name="DataFormats/JetReco"/>
+<use   name="DataFormats/METReco"/>
 <library   file="*.cc" name="GeneratorInterfaceRivetInterface_plugins">
   <flags   EDM_PLUGIN="1"/>
 </library>

GeneratorInterface/RivetInterface/plugins/ParticleLevelProducer.cc

@@ -0,0 +1,167 @@
+#include "GeneratorInterface/RivetInterface/interface/ParticleLevelProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/Common/interface/RefToPtr.h"
+#include "DataFormats/METReco/interface/METFwd.h"
+#include "DataFormats/METReco/interface/MET.h"
+#include "RecoJets/JetProducers/interface/JetSpecific.h"
+#include "CommonTools/Utils/interface/PtComparator.h"
+
+#include "Rivet/Analysis.hh"
+
+using namespace std;
+using namespace edm;
+using namespace reco;
+using namespace Rivet;
+
+ParticleLevelProducer::ParticleLevelProducer(const edm::ParameterSet& pset):
+  srcToken_(consumes<edm::HepMCProduct>(pset.getParameter<edm::InputTag>("src"))),
+  rivetAnalysis_(new Rivet::RivetAnalysis(pset))
+{
+  usesResource();
+
+  genVertex_ = reco::Particle::Point(0,0,0);
+
+  produces<reco::GenParticleCollection>("neutrinos");
+  produces<reco::GenParticleCollection>("photons");
+  produces<reco::GenJetCollection>("leptons");
+  produces<reco::GenJetCollection>("jets");
+  produces<reco::GenJetCollection>("fatjets");
+  produces<reco::GenParticleCollection>("consts");
+  produces<reco::GenParticleCollection>("tags");
+  produces<reco::METCollection>("mets");
+
+  analysisHandler_.addAnalysis(rivetAnalysis_);
+}
+
+void ParticleLevelProducer::addGenJet(Rivet::Jet jet, std::unique_ptr<reco::GenJetCollection> &jets,
+                                      std::unique_ptr<reco::GenParticleCollection> &consts, edm::RefProd<reco::GenParticleCollection>& constsRefHandle, int &iConstituent,
+                                      std::unique_ptr<reco::GenParticleCollection> &tags, edm::RefProd<reco::GenParticleCollection>& tagsRefHandle, int &iTag)
+{
+  const auto pjet = jet.pseudojet();
+
+  reco::GenJet genJet;
+  genJet.setP4(p4(jet));
+  genJet.setVertex(genVertex_);
+  if ( jet.bTagged() ) genJet.setPdgId(5);
+  genJet.setJetArea(pjet.has_area() ? pjet.area() : 0);
+
+  for ( auto const & p : jet.particles()) {
+    auto pp4 = p4(p);
+    bool match = false; int iMatch = -1;
+    for ( auto const & q : *consts ) {
+      ++iMatch;
+      if (q.p4() == pp4) { match = true; break; }
+    }
+    if (match){
+      genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, iMatch)));
+    }
+    else {
+      consts->push_back(reco::GenParticle(p.charge(), pp4, genVertex_, p.pdgId(), 1, true));
+      genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, ++iConstituent)));
+    }
+  }
+  for ( auto const & p : jet.tags()) {
+    // The tag particles are accessible as jet daughters, so scale down p4 for safety.
+    // p4 needs to be multiplied by 1e20 for fragmentation analysis.
+    auto pp4 = p4(p)*1e-20;
+    bool match = false; int iMatch = -1;
+    for ( auto const & q : *tags ) {
+      ++iMatch;
+      if (q.p4() == pp4) { match = true; break; }
+    }
+    if (match){
+      genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(tagsRefHandle, iMatch)));
+    }
+    else {    
+      tags->push_back(reco::GenParticle(p.charge(), p4(p)*1e-20, genVertex_, p.pdgId(), 2, true));
+      genJet.addDaughter(edm::refToPtr(reco::GenParticleRef(tagsRefHandle, ++iTag)));
+    }
+  }
+
+  jets->push_back(genJet);
+}
+
+void ParticleLevelProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup)
+{
+  using namespace Rivet;
+  typedef reco::Candidate::LorentzVector LorentzVector;
+
+  std::unique_ptr<reco::GenParticleCollection> neutrinos(new reco::GenParticleCollection);
+  std::unique_ptr<reco::GenParticleCollection> photons(new reco::GenParticleCollection);
+  std::unique_ptr<reco::GenJetCollection> leptons(new reco::GenJetCollection);
+  std::unique_ptr<reco::GenJetCollection> jets(new reco::GenJetCollection);
+  std::unique_ptr<reco::GenJetCollection> fatjets(new reco::GenJetCollection);
+  std::unique_ptr<reco::GenParticleCollection> consts(new reco::GenParticleCollection);
+  std::unique_ptr<reco::GenParticleCollection> tags(new reco::GenParticleCollection);
+  std::unique_ptr<reco::METCollection> mets(new reco::METCollection);
+  auto constsRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("consts");
+  auto tagsRefHandle = event.getRefBeforePut<reco::GenParticleCollection>("tags");
+
+  edm::Handle<HepMCProduct> srcHandle;
+  event.getByToken(srcToken_, srcHandle);
+
+  const HepMC::GenEvent* genEvent = srcHandle->GetEvent();
+  analysisHandler_.analyze(*genEvent);
+
+  // Convert into edm objects
+  // Prompt neutrinos
+  for ( auto const & p : rivetAnalysis_->neutrinos() ) {
+    neutrinos->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pdgId(), 1, true));
+  }
+  std::sort(neutrinos->begin(), neutrinos->end(), GreaterByPt<reco::Candidate>());
+
+  // Photons
+  for ( auto const & p : rivetAnalysis_->photons() ) {
+    photons->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pdgId(), 1, true));
+  }
+  std::sort(photons->begin(), photons->end(), GreaterByPt<reco::Candidate>());
+
+  // Prompt leptons
+  int iConstituent = -1;
+  for ( auto const & lepton : rivetAnalysis_->leptons() ) {
+    reco::GenJet lepJet;
+    lepJet.setP4(p4(lepton));
+    lepJet.setVertex(genVertex_);
+    lepJet.setPdgId(lepton.pdgId());
+    lepJet.setCharge(lepton.charge());
+
+    const auto cl = lepton.constituentLepton();
+    consts->push_back(reco::GenParticle(cl.charge(), p4(cl), genVertex_, cl.pdgId(), 1, true));
+    lepJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, ++iConstituent)));
+
+    for ( auto const & p : lepton.constituentPhotons()) {
+      consts->push_back(reco::GenParticle(p.charge(), p4(p), genVertex_, p.pdgId(), 1, true));
+      lepJet.addDaughter(edm::refToPtr(reco::GenParticleRef(constsRefHandle, ++iConstituent)));
+    }
+
+    leptons->push_back(lepJet);
+  }
+  std::sort(leptons->begin(), leptons->end(), GreaterByPt<reco::GenJet>());
+
+  // Jets with constituents and tag particles
+  int iTag = -1;
+  for ( auto jet : rivetAnalysis_->jets() ) {
+    addGenJet(jet, jets, consts, constsRefHandle, iConstituent, tags, tagsRefHandle, iTag);
+  }
+  for ( auto jet : rivetAnalysis_->fatjets() ) {
+    addGenJet(jet, fatjets, consts, constsRefHandle, iConstituent, tags, tagsRefHandle, iTag);
+  }
+
+  // MET
+  reco::Candidate::LorentzVector metP4(rivetAnalysis_->met().x(), rivetAnalysis_->met().y(), 0., sqrt(pow(rivetAnalysis_->met().x(), 2) + pow(rivetAnalysis_->met().y(), 2)));
+  mets->push_back(reco::MET(metP4, genVertex_));
+
+  event.put(std::move(neutrinos), "neutrinos");
+  event.put(std::move(photons), "photons");
+  event.put(std::move(leptons), "leptons");
+  event.put(std::move(jets), "jets");
+  event.put(std::move(fatjets), "fatjets");
+  event.put(std::move(consts), "consts");
+  event.put(std::move(tags), "tags");
+  event.put(std::move(mets), "mets");
+}
+
+DEFINE_FWK_MODULE(ParticleLevelProducer);