File indexing completed on 2024-09-08 23:51:48
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
0077 void init() override {
0078
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
0083 FinalState fs(particle_cut);
0084
0085 declare(fs, "FS");
0086
0087
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
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
0117 declare(FastJets(fsForJets, JetAlg::ANTIKT, _fatJetConeSize, JetMuons::ALL, invisiblesStrategy), "FatJets");
0118
0119
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
0131 declare(MissingMomentum(fs), "MET");
0132 };
0133
0134
0135 void analyze(const Event& event) override {
0136 _jets.clear();
0137 _fatjets.clear();
0138 _leptons.clear();
0139 _photons.clear();
0140 _neutrinos.clear();
0141
0142
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
0149 Particles promptleptons = apply<PromptFinalState>(event, "PromptLeptons").particles();
0150 for (auto& lepton : _leptons) {
0151 const auto& cl = lepton.constituents().front();
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
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
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
0218 void finalize() override {}
0219
0220 std::string status() const override { return "VALIDATED"; }
0221 };
0222
0223 }
0224
0225 #endif