RivetInterface/interface/RivetAnalysis.h

0001 #ifndef GeneratorInterface_RivetInterface_RivetAnalysis_H
0002 #define GeneratorInterface_RivetInterface_RivetAnalysis_H
0003
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005
0006 #include "Rivet/Analysis.hh"
0007 #include "Rivet/Projections/FinalState.hh"
0008 #include "Rivet/Particle.hh"
0009 #include "Rivet/Particle.fhh"
0010 #include "Rivet/Event.hh"
0011 #include "Rivet/Projections/FastJets.hh"
0012 #include "Rivet/Projections/ChargedLeptons.hh"
0013 #include "Rivet/Projections/PromptFinalState.hh"
0014 #include "Rivet/Projections/LeptonFinder.hh"
0015 #include "Rivet/Projections/VetoedFinalState.hh"
0016 #include "Rivet/Projections/IdentifiedFinalState.hh"
0017 #include "Rivet/Projections/MissingMomentum.hh"
0018 #include "Rivet/Tools/RivetHepMC.hh"
0019
0020 namespace Rivet {
0021
0022   class RivetAnalysis : public Analysis {
0023   public:
0024     Particles leptons() const { return _leptons; }
0025     Particles photons() const { return _photons; }
0026     Particles neutrinos() const { return _neutrinos; }
0027     Jets jets() const { return _jets; }
0028     Jets fatjets() const { return _fatjets; }
0029     Vector3 met() const { return _met; }
0030
0031   private:
0032     bool _usePromptFinalStates;
0033     bool _excludePromptLeptonsFromJetClustering;
0034     bool _excludeNeutrinosFromJetClustering;
0035     bool _doJetClustering;
0036
0037     double _particleMinPt, _particleMaxEta;
0038     double _lepConeSize, _lepMinPt, _lepMaxEta;
0039     double _jetConeSize, _jetMinPt, _jetMaxEta;
0040     double _fatJetConeSize, _fatJetMinPt, _fatJetMaxEta;
0041
0042     double _phoMinPt, _phoMaxEta, _phoIsoConeSize, _phoMaxRelIso;
0043
0044     Particles _leptons, _photons, _neutrinos;
0045     Jets _jets, _fatjets;
0046     Vector3 _met;
0047
0048   public:
0049     RivetAnalysis(const edm::ParameterSet& pset)
0050         : Analysis("RivetAnalysis"),
0051           _usePromptFinalStates(pset.getParameter<bool>("usePromptFinalStates")),
0052           _excludePromptLeptonsFromJetClustering(pset.getParameter<bool>("excludePromptLeptonsFromJetClustering")),
0053           _excludeNeutrinosFromJetClustering(pset.getParameter<bool>("excludeNeutrinosFromJetClustering")),
0054           _doJetClustering(pset.getParameter<bool>("doJetClustering")),
0055
0056           _particleMinPt(pset.getParameter<double>("particleMinPt")),
0057           _particleMaxEta(pset.getParameter<double>("particleMaxEta")),
0058
0059           _lepConeSize(pset.getParameter<double>("lepConeSize")),
0060           _lepMinPt(pset.getParameter<double>("lepMinPt")),
0061           _lepMaxEta(pset.getParameter<double>("lepMaxEta")),
0062
0063           _jetConeSize(pset.getParameter<double>("jetConeSize")),
0064           _jetMinPt(pset.getParameter<double>("jetMinPt")),
0065           _jetMaxEta(pset.getParameter<double>("jetMaxEta")),
0066
0067           _fatJetConeSize(pset.getParameter<double>("fatJetConeSize")),
0068           _fatJetMinPt(pset.getParameter<double>("fatJetMinPt")),
0069           _fatJetMaxEta(pset.getParameter<double>("fatJetMaxEta")),
0070
0071           _phoMinPt(pset.getParameter<double>("phoMinPt")),
0072           _phoMaxEta(pset.getParameter<double>("phoMaxEta")),
0073           _phoIsoConeSize(pset.getParameter<double>("phoIsoConeSize")),
0074           _phoMaxRelIso(pset.getParameter<double>("phoMaxRelIso")) {}
0075
0076     // Initialize Rivet projections
0077     void init() override {
0078       // Cuts
0079       Cut particle_cut = (Cuts::abseta < _particleMaxEta) and (Cuts::pT > _particleMinPt * GeV);
0080       Cut lepton_cut = (Cuts::abseta < _lepMaxEta) and (Cuts::pT > _lepMinPt * GeV);
0081
0082       // Generic final state
0083       FinalState fs(particle_cut);
0084
0085       declare(fs, "FS");
0086
0087       // Dressed leptons
0088       ChargedLeptons charged_leptons(fs);
0089       IdentifiedFinalState photons(fs);
0090       photons.acceptIdPair(PID::PHOTON);
0091
0092       PromptFinalState prompt_leptons(charged_leptons);
0093       prompt_leptons.acceptMuonDecays(true);
0094       prompt_leptons.acceptTauDecays(true);
0095       declare(prompt_leptons, "PromptLeptons");
0096
0097       PromptFinalState prompt_photons(photons);
0098       prompt_photons.acceptMuonDecays(true);
0099       prompt_photons.acceptTauDecays(true);
0100
0101       LeptonFinder prompt_dressed_leptons(prompt_leptons, prompt_photons, _lepConeSize, lepton_cut);
0102       LeptonFinder dressed_leptons(charged_leptons, photons, _lepConeSize, lepton_cut);
0103       declare(_usePromptFinalStates ? prompt_dressed_leptons : dressed_leptons, "DressedLeptons");
0104
0105       declare(photons, "Photons");
0106
0107       // Jets
0108       VetoedFinalState fsForJets(fs);
0109       if (_usePromptFinalStates and _excludePromptLeptonsFromJetClustering)
0110         fsForJets.addVetoOnThisFinalState(prompt_dressed_leptons);
0111       JetInvisibles invisiblesStrategy = JetInvisibles::DECAY;
0112       if (_excludeNeutrinosFromJetClustering)
0113         invisiblesStrategy = JetInvisibles::NONE;
0114       declare(FastJets(fsForJets, JetAlg::ANTIKT, _jetConeSize, JetMuons::ALL, invisiblesStrategy), "Jets");
0115
0116       // FatJets
0117       declare(FastJets(fsForJets, JetAlg::ANTIKT, _fatJetConeSize, JetMuons::ALL, invisiblesStrategy), "FatJets");
0118
0119       // Neutrinos
0120       IdentifiedFinalState neutrinos(fs);
0121       neutrinos.acceptNeutrinos();
0122       if (_usePromptFinalStates) {
0123         PromptFinalState prompt_neutrinos(neutrinos);
0124         prompt_neutrinos.acceptMuonDecays(true);
0125         prompt_neutrinos.acceptTauDecays(true);
0126         declare(prompt_neutrinos, "Neutrinos");
0127       } else
0128         declare(neutrinos, "Neutrinos");
0129
0130       // MET
0131       declare(MissingMomentum(fs), "MET");
0132     };
0133
0134     // Apply Rivet projections
0135     void analyze(const Event& event) override {
0136       _jets.clear();
0137       _fatjets.clear();
0138       _leptons.clear();
0139       _photons.clear();
0140       _neutrinos.clear();
0141
0142       // Get analysis objects from projections
0143       Cut jet_cut = (Cuts::abseta < _jetMaxEta) and (Cuts::pT > _jetMinPt * GeV);
0144       Cut fatjet_cut = (Cuts::abseta < _fatJetMaxEta) and (Cuts::pT > _fatJetMinPt * GeV);
0145
0146       _leptons = apply<LeptonFinder>(event, "DressedLeptons").particlesByPt();
0147
0148       // search tau ancestors
0149       Particles promptleptons = apply<PromptFinalState>(event, "PromptLeptons").particles();
0150       for (auto& lepton : _leptons) {
0151         const auto& cl = lepton.constituents().front();  // this has no ancestors anymore :(
0152         for (auto const& pl : promptleptons) {
0153           if (cl.momentum() == pl.momentum()) {
0154             for (auto& p : pl.ancestors()) {
0155               if (p.abspid() == 15) {
0156                 p.setMomentum(p.momentum() * 1e-20);
0157                 lepton.addConstituent(p, false);
0158               }
0159             }
0160             break;
0161           }
0162         }
0163       }
0164
0165       // Photons
0166       Particles fsparticles = apply<FinalState>(event, "FS").particles();
0167
0168       for (auto& photon : apply<FinalState>(event, "Photons").particlesByPt()) {
0169         if (photon.pt() < _phoMinPt)
0170           continue;
0171
0172         if (abs(photon.eta()) > _phoMaxEta)
0173           continue;
0174
0175         double photonptsum = 0;
0176
0177         for (auto& fsparticle : fsparticles) {
0178           if (deltaR(fsparticle, photon) == 0)
0179             continue;
0180
0181           if (deltaR(fsparticle, photon) > _phoIsoConeSize)
0182             continue;
0183
0184           photonptsum += fsparticle.pt();
0185         }
0186
0187         if (photonptsum / photon.pt() > _phoMaxRelIso)
0188           continue;
0189
0190         _photons.push_back(photon);
0191       }
0192
0193       if (_doJetClustering) {
0194         _jets = apply<FastJets>(event, "Jets").jetsByPt(jet_cut);
0195         _fatjets = apply<FastJets>(event, "FatJets").jetsByPt(fatjet_cut);
0196       }
0197       _neutrinos = apply<FinalState>(event, "Neutrinos").particlesByPt();
0198       _met = apply<MissingMomentum>(event, "MET").missingMomentum().p3();
0199
0200       // check for leptonic decays of tag particles
0201       for (auto& jet : _jets) {
0202         for (auto& pt : jet.tags()) {
0203           for (auto& p : pt.children(isLepton)) {
0204             pt.addConstituent(p, false);
0205           }
0206         }
0207       }
0208       for (auto& jet : _fatjets) {
0209         for (auto& pt : jet.tags()) {
0210           for (auto& p : pt.children(isLepton)) {
0211             pt.addConstituent(p, false);
0212           }
0213         }
0214       }
0215     };
0216
0217     // Do nothing here
0218     void finalize() override {}
0219
0220     std::string status() const override { return "VALIDATED"; }
0221   };
0222
0223 }  // namespace Rivet
0224
0225 #endif