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