File indexing completed on 2023-03-17 13:02:21
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
0078 void init() override {
0079
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
0084 FinalState fs(particle_cut);
0085
0086 declare(fs, "FS");
0087
0088
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
0103
0104
0105 DressedLeptons dressed_leptons(
0106 prompt_photons, prompt_leptons, _lepConeSize, lepton_cut, true);
0107 if (not _usePromptFinalStates)
0108 dressed_leptons = DressedLeptons(photons, charged_leptons, _lepConeSize, lepton_cut, true);
0109 declare(dressed_leptons, "DressedLeptons");
0110
0111 declare(photons, "Photons");
0112
0113
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
0123 declare(FastJets(fsForJets, FastJets::ANTIKT, _fatJetConeSize), "FatJets");
0124
0125
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
0137 declare(MissingMomentum(fs), "MET");
0138 };
0139
0140
0141 void analyze(const Event& event) override {
0142 _jets.clear();
0143 _fatjets.clear();
0144 _leptons.clear();
0145 _photons.clear();
0146 _neutrinos.clear();
0147
0148
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
0155 Particles promptleptons = apply<PromptFinalState>(event, "PromptLeptons").particles();
0156 for (auto& lepton : _leptons) {
0157 const auto& cl = lepton.constituents().front();
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
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
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
0224 void finalize() override{};
0225
0226 std::string status() const override { return "VALIDATED"; }
0227 };
0228
0229 }
0230
0231 #endif