SimpleAnalysis.hh

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DirectFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/LeptonFinder.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/Smearing.hh"
#include "Rivet/Projections/TauFinder.hh"
#include "Rivet/Tools/Cutflow.hh"
 
namespace Rivet {
 

  class SimpleAnalysis : public Analysis {
  public:

    enum class IDClass { LOOSE, MEDIUM, TIGHT };

    SimpleAnalysis(const std::string& name)
        : Analysis(name) { }

    virtual ~SimpleAnalysis() { }
 


    void postInit() {
 
      // Electron projections
      FinalState es(_acut && _ecut && Cuts::abspid == PID::ELECTRON);
      LeptonFinder es_direct(_acut && _ecut && Cuts::abspid == PID::ELECTRON, 0.1);
      for (const auto& classkv : _idclasses()) {
        const IDClass& idc = classkv.first;
        const string pre = "Electrons";
        const string aname = _pckey(true, idc, pre);
        if (!hasProjection(aname)) {
          declare(SmearedParticles(es, _effs_electron[idc], _smear_electron), aname);
        }
        const string pname = _pckey(false, idc, pre);
        if (!hasProjection(pname)) {
          declare(SmearedParticles(es_direct, _effs_electron[idc], _smear_electron), pname);
        }
      }
 
      // Muon projections
      FinalState ms(_acut && _mcut && Cuts::abspid == PID::MUON);
      LeptonFinder ms_direct(_acut && _mcut && Cuts::abspid == PID::MUON, 0.1);
      for (const auto& classkv : _idclasses()) {
        const IDClass& idc = classkv.first;
        const string pre = "Muons";
        const string aname = _pckey(true, idc, pre);
        if (!hasProjection(aname)) {
          declare(SmearedParticles(ms, _effs_muon[idc], _smear_muon), aname);
        }
        const string pname = _pckey(false, idc, pre);
        if (!hasProjection(pname)) {
          declare(SmearedParticles(ms_direct, _effs_muon[idc], _smear_muon), pname);
        }
      }
 
      // Photon projections
      FinalState ys(_acut && _phocut && Cuts::abspid == PID::PHOTON);
      DirectFinalState ys_direct(_acut && _phocut && Cuts::abspid == PID::PHOTON);
      for (const auto& classkv : _idclasses()) {
        const IDClass& idc = classkv.first;
        const string pre = "Photons";
        const string aname = _pckey(true, idc, pre);
        if (!hasProjection(aname)) {
          declare(SmearedParticles(ys, _effs_photon[idc], _smear_photon), aname);
        }
        const string pname = _pckey(false, idc, pre);
        if (!hasProjection(pname)) {
          declare(SmearedParticles(ys_direct, _effs_photon[idc], _smear_photon), pname);
        }
      }
 
      // Tau projections
      TauFinder ts(TauDecay::HADRONIC, LeptonOrigin::ANY, _acut && _tcut);
      TauFinder ts_direct(TauDecay::HADRONIC, LeptonOrigin::NODECAY, _acut && _tcut);
      for (const auto& classkv : _idclasses()) {
        const IDClass& idc = classkv.first;
        const string pre = "Taus";
        const string aname = _pckey(true, idc, pre);
        if (!hasProjection(aname)) {
          declare(SmearedParticles(ts, _effs_tau[idc], _smear_tau), aname);
        }
        const string pname = _pckey(false, idc, pre);
        if (!hasProjection(pname)) {
          declare(SmearedParticles(ts_direct, _effs_tau[idc], _smear_tau), pname);
        }
      }
 
      // Track projection
      const FinalState trkfs(_acut && _trkcut && Cuts::abscharge > 0);
      if (!hasProjection("Tracks")) {
        declare(SmearedParticles(trkfs, _eff_trk, _smear_trk), "Tracks");
      }
 
      // Jet projections
      const FinalState allfs(_acut);
      declare(allfs, "AllParticles");
      for (const auto& jditem : _jdefs) {
        const JetScheme& jd = jditem.second;
        FastJets jetfs(allfs, jd.alg, jd.R, jd.muons, jd.invis);
        if (!hasProjection("Jets" + jditem.first)) {
          declare(SmearedJets(jetfs, _smears_jet[jditem.first], _effs_btag[jditem.first]),
                  "Jets" + jditem.first);
        }
      }
 
      // Missing momentum projection
      if (!hasProjection("MET")) {
        declare(SmearedMET(MissingMomentum(allfs), _smearps_met, _smear_met), "MET");
      }
    }
 

    void preAnalyze(const Event& event) {
      _electrons.clear();
      _muons.clear();
      _taus.clear();
      _photons.clear();
      _jets.clear();
      _bjets.clear();
    }

 



    void setElectronReco(const ParticleEffFn& eff,
                         const ParticleSmearFn& smear = PARTICLE_SMEAR_IDENTITY,
                         IDClass idc = IDClass::MEDIUM) {
      _effs_electron[idc] = eff;
      _smear_electron = smear;
    }

    void setMuonReco(const ParticleEffFn& eff,
                     const ParticleSmearFn& smear = PARTICLE_SMEAR_IDENTITY,
                     IDClass idc = IDClass::MEDIUM) {
      _effs_muon[idc] = eff;
      _smear_muon = smear;
    }

    void setPhotonReco(const ParticleEffFn& eff,
                       const ParticleSmearFn& smear = PARTICLE_SMEAR_IDENTITY,
                       IDClass idc = IDClass::MEDIUM) {
      _effs_photon[idc] = eff;
      _smear_photon = smear;
    }

    void setTauReco(const ParticleEffFn& eff,
                    const ParticleSmearFn& smear = PARTICLE_SMEAR_IDENTITY,
                    IDClass idc = IDClass::MEDIUM) {
      _effs_tau[idc] = eff;
      _smear_tau = smear;
    }

    void setTrkReco(const ParticleEffFn& eff, const ParticleSmearFn& smear = PARTICLE_SMEAR_IDENTITY) {
      _eff_trk = eff;
      _smear_trk = smear;
    }

    void setJetReco(const JetSmearFn& smear,
                    const JetEffFn& btageff,
                    const std::string& jetname = "DEFAULT") {
      _smears_jet[jetname] = smear;
      _effs_btag[jetname] = btageff;
    }

    void setMETReco(const METSmearParamsFn& smearps) {
      _smearps_met = smearps;
    }

    void setMETReco(const METSmearParamsFn& smearps, const METSmearFn& smear) {
      _smearps_met = smearps;
      _smear_met = smear;
    }

 

    void setAccCut(const Cut& cut) {
      _acut = cut;
    }
    void setElectronCut(const Cut& cut) {
      _ecut = cut;
    }
    void setMuonCut(const Cut& cut) {
      _mcut = cut;
    }
    void setPhotonCut(const Cut& cut) {
      _phocut = cut;
    }
    void setTauCut(const Cut& cut) {
      _tcut = cut;
    }
    void setTrkCut(const Cut& cut) {
      _trkcut = cut;
    }
    void setJetDef(const JetScheme& jdef, const std::string& jetname = "DEFAULT") {
      _jdefs[jetname] = jdef;
    }
    void setJetCut(const Cut& cut, const std::string& jetname = "DEFAULT") {
      _jcuts[jetname] = cut;
    }

 


    const Particles& electrons(IDClass eclass = IDClass::MEDIUM, bool include_indirect = false) const {
      const std::string kpc = _pckey(include_indirect, eclass);
      if (_electrons.find(kpc) == _electrons.end()) { //< if not cached
        _electrons[kpc] = apply<ParticleFinder>(currentEvent(), "Electrons" + kpc).particlesByPt();
      }
      return _electrons.at(kpc);
    }

    Particles electrons(const Cut& cut,
                        IDClass eclass = IDClass::MEDIUM,
                        bool include_indirect = false) const {
      return select(electrons(eclass, include_indirect), cut);
    }
 

    const Particles& muons(IDClass muclass = IDClass::MEDIUM, bool include_indirect = false) const {
      const std::string kpc = _pckey(include_indirect, muclass);
      if (_muons.find(kpc) == _muons.end()) { //< if not cached
        _muons[kpc] = apply<ParticleFinder>(currentEvent(), "Muons" + kpc).particlesByPt();
      }
      return _muons.at(kpc);
    }

    Particles muons(const Cut& cut, IDClass muclass = IDClass::MEDIUM, bool include_indirect = false) const {
      return select(muons(muclass, include_indirect), cut);
    }
 

    const Particles& photons(IDClass phoclass = IDClass::MEDIUM,
                             bool include_indirect = false,
                             double dRiso = 0.4,
                             double isofrac = 0.1) const {
      const std::string kpc = _pckey(include_indirect, phoclass);
      if (_photons.find(kpc) == _photons.end()) { //< if not cached
        Particles ystmp = apply<ParticleFinder>(currentEvent(), "Photons" + kpc).particlesByPt();
        const Particles& clusters = apply<ParticleFinder>(currentEvent(), "AllParticles").particles();
        _photons[kpc].reserve(ystmp.size());
        for (const Particle& y : ystmp) {
          double sumpt = 0;
          for (const Particle& cl : select(clusters, deltaRLess(y, dRiso))) {
            if (deltaR(cl, y) == 0.01) continue; //< irresolvably close
            sumpt += cl.pT();
          }
          if (sumpt / y.pT() < isofrac) _photons[kpc].push_back(y);
        }
      }
      return _photons.at(kpc);
    }

    Particles photons(const Cut& cut,
                      IDClass phoclass = IDClass::MEDIUM,
                      bool include_indirect = false) const {
      return select(photons(phoclass, include_indirect), cut);
    }
 

    const Particles& taus(IDClass tauclass = IDClass::MEDIUM, bool include_indirect = false) const {
      const std::string kpc = _pckey(include_indirect, tauclass);
      if (_taus.find(kpc) == _taus.end()) { //< if not cached
        _taus[kpc] = apply<ParticleFinder>(currentEvent(), "Taus" + kpc).particlesByPt();
      }
      return _taus.at(kpc);
    }

    Particles taus(const Cut& cut, IDClass tauclass = IDClass::MEDIUM, bool include_indirect = false) const {
      return select(taus(tauclass, include_indirect), cut);
    }
 

    const Particles tracks() const {
      return apply<ParticleFinder>(currentEvent(), "Tracks").particlesByPt();
    }

    Particles tracks(const Cut& cut) const {
      return apply<ParticleFinder>(currentEvent(), "Tracks").particlesByPt(cut);
    }
 

    const std::vector<std::string> jetnames() const {
      std::vector<std::string> rtn;
      rtn.reserve(_jets.size());
      for (const auto& kv : _jets) rtn.push_back(kv.first);
      return rtn;
    }

    const Jets& jets(const std::string& jetname = "DEFAULT") const {
      if (_jets.find(jetname) == _jets.end()) { //< not cached
        _jets[jetname] = apply<JetFinder>(currentEvent(), "Jets" + jetname).jetsByPt(_jcuts.at(jetname));
      }
      return _jets.at(jetname);
    }

    Jets jets(const Cut& cut, const std::string& jetname = "DEFAULT") const {
      return select(jets(jetname), cut);
    }
 

    const Jets& bjets(const std::string& jetname = "DEFAULT") const {
      if (_bjets.find(jetname) == _bjets.end()) { //< not cached
        _bjets[jetname] = select(jets(jetname), hasBTag(_bcut, _bdeltaR));
      }
      return _bjets.at(jetname);
    }

    Jets bjets(const Cut& cut, const std::string& jetname = "DEFAULT") const {
      return select(bjets(jetname), cut);
    }
 

    Vector3 vmet() const {
      return apply<METFinder>(currentEvent(), "MET").vectorMissingPt();
    }

    double met() const {
      return vmet().mod();
    }

    double set() const {
      return apply<METFinder>(currentEvent(), "MET").scalarEt();
    }

    double metSignf() const {
      return apply<SmearedMET>(currentEvent(), "MET").missingEtSignf();
    }

 


    void doSimpleOverlapRemoval(IDClass eclass = IDClass::MEDIUM,
                                IDClass muclass = IDClass::MEDIUM,
                                IDClass phoclass = IDClass::MEDIUM, //IDClass tauclass=IDClass::MEDIUM,
                                bool include_indirect = false,
                                const std::string& jetname = "DEFAULT") {
      // Get non-const handles for filtering (and ensure objects exist)
      Jets& myjets = const_cast<Jets&>(jets(jetname));
      Jets& mybjets = const_cast<Jets&>(bjets(jetname));
      Particles& myes = const_cast<Particles&>(electrons(eclass, include_indirect));
      Particles& myms = const_cast<Particles&>(muons(muclass, include_indirect));
      //Particles& myts = const_cast<Particles&>(taus(tauclass, include_indirect));
      Particles& myys = const_cast<Particles&>(photons(phoclass, include_indirect));
      // Remove all jets within dR < 0.2 of a dressed lepton
      idiscardIfAnyDeltaRLess(myjets, myes, 0.2);
      idiscardIfAnyDeltaRLess(myjets, myms, 0.2);
      idiscardIfAnyDeltaRLess(mybjets, myes, 0.2);
      idiscardIfAnyDeltaRLess(mybjets, myms, 0.2);
      // Remove all isolated photons within dR < 0.2 of a dressed electron
      idiscardIfAnyDeltaRLess(myys, myes, 0.2);
      // Remove dressed leptons within the radius of a remaining jet
      idiscardIfAnyDeltaRLess(myes, myjets, _jdefs.at(jetname).R);
      idiscardIfAnyDeltaRLess(myms, myjets, _jdefs.at(jetname).R);
    }
 

    // void doJetRecluster(/* jet collection names (in and out), recluster-alg enum and dR */) {
    //   /// @todo ...
    // }

 


    template <typename T>
    void scaleToIntLumi(T& ao, double intlumi) {
      scale(ao, intlumi * crossSection() / sumOfWeights());
    }

 


    const ParticleFinder& electronsProj(IDClass eclass = IDClass::MEDIUM, bool include_indirect = false) {
      return getProjection<ParticleFinder>(_pckey(include_indirect, eclass, "Electrons"));
    }

    const ParticleFinder& muonsProj(IDClass muclass = IDClass::MEDIUM, bool include_indirect = false) {
      return getProjection<ParticleFinder>(_pckey(include_indirect, muclass, "Muons"));
    }

    const ParticleFinder& tausProj(IDClass tauclass = IDClass::MEDIUM, bool include_indirect = false) {
      return getProjection<ParticleFinder>(_pckey(include_indirect, tauclass, "Taus"));
    }

    const ParticleFinder& photonsProj(IDClass phoclass = IDClass::MEDIUM, bool include_indirect = false) {
      return getProjection<ParticleFinder>(_pckey(include_indirect, phoclass, "Photons"));
    }

    const ParticleFinder& tracksProj() {
      return getProjection<ParticleFinder>("Tracks");
    }

    const JetFinder& jetsProj(const std::string& jetname = "DEFAULT") {
      return getProjection<JetFinder>("Jets_" + jetname);
    }

    const METFinder& metProj() {
      return getProjection<METFinder>("MET");
    }

 
 
  private:

    std::map<IDClass, ParticleEffFn> _effs_electron{{IDClass::LOOSE, PARTICLE_EFF_ONE},
                                                    {IDClass::MEDIUM, PARTICLE_EFF_ONE},
                                                    {IDClass::TIGHT, PARTICLE_EFF_ONE}};
    std::map<IDClass, ParticleEffFn> _effs_muon{{IDClass::LOOSE, PARTICLE_EFF_ONE},
                                                {IDClass::MEDIUM, PARTICLE_EFF_ONE},
                                                {IDClass::TIGHT, PARTICLE_EFF_ONE}};
    std::map<IDClass, ParticleEffFn> _effs_photon{{IDClass::LOOSE, PARTICLE_EFF_ONE},
                                                  {IDClass::MEDIUM, PARTICLE_EFF_ONE},
                                                  {IDClass::TIGHT, PARTICLE_EFF_ONE}};
    std::map<IDClass, ParticleEffFn> _effs_tau{{IDClass::LOOSE, PARTICLE_EFF_ONE},
                                               {IDClass::MEDIUM, PARTICLE_EFF_ONE},
                                               {IDClass::TIGHT, PARTICLE_EFF_ONE}};
    ParticleSmearFn _smear_electron{PARTICLE_SMEAR_IDENTITY}, _smear_muon{PARTICLE_SMEAR_IDENTITY},
        _smear_photon{PARTICLE_SMEAR_IDENTITY}, _smear_tau{PARTICLE_SMEAR_IDENTITY};
    ParticleEffFn _eff_trk{PARTICLE_EFF_ONE};
    ParticleSmearFn _smear_trk{PARTICLE_SMEAR_IDENTITY};
    std::map<std::string, JetSmearFn> _smears_jet{{"DEFAULT", JET_SMEAR_IDENTITY}};
    std::map<std::string, JetEffFn> _effs_btag{{"DEFAULT", JET_BTAG_IDENTITY}};
    METSmearParamsFn _smearps_met{MET_SMEARPARAMS_IDENTITY};
    METSmearFn _smear_met{MET_SMEAR_IDENTITY};

    Cut _acut = Cuts::abseta < 5.0; //< base acceptance, applied to all
    Cut _ecut = Cuts::pT > 10 * GeV;
    Cut _mcut = Cuts::pT > 10 * GeV;
    Cut _tcut = Cuts::pT > 10 * GeV;
    Cut _trkcut = Cuts::abseta < 2.5;
    Cut _phocut = Cuts::pT > 10 * GeV;
    std::map<std::string, Cut> _jcuts = {{"DEFAULT", Cuts::pT > 20 * GeV&& Cuts::abseta < 5.0}};
    std::map<std::string, JetScheme> _jdefs = {{"DEFAULT", JetScheme(JetAlg::ANTIKT, 0.4)}};
 

    Cut _bcut = Cuts::pT > 5 * GeV&& Cuts::abseta < 2.5;
    double _bdeltaR = 0.3;
 

 
    // Internal typedefs
    // using ParticlesMap = std::map<std::string, Particles>;
    // using JetsMap = std::map<std::string, Jets>;

    mutable std::map<std::string, Particles> _electrons, _muons, _taus, _photons;

    mutable std::map<std::string, Jets> _jets, _bjets;

 

    //
    const std::map<IDClass, std::string>& _idclasses() const {
      // static const std::vector<IDClass> _idclasses{ {IDClass::LOOSE, IDClass::MEDIUM, IDClass::TIGHT} };
      static const auto* map = new std::map<IDClass, std::string>{{IDClass::LOOSE, "Loose"},
                                                                  {IDClass::MEDIUM, "Medium"},
                                                                  {IDClass::TIGHT, "Tight"}};
      return *map;
    }
    const std::string& _idclassname(const IDClass& idc) const {
      return _idclasses().at(idc);
    }

    std::string _pckey(bool include_indirect, IDClass idc, const std::string& prefix = "") const {
      return prefix + (include_indirect ? "All" : "Direct") + _idclassname(idc);
    }
  };
 
 
}
 

#define RIVET_DEFAULT_SIMPLEANALYSIS_CTOR(clsname) clsname() : SimpleAnalysis(# clsname) {}
 

#define setDetSmearing(DET, BTAG) do {                                  \
  setElectronReco(ELECTRON_EFF_ ## DET ## _LOOSE, ELECTRON_SMEAR_ ## DET, IDClass::LOOSE) ; \
  setElectronReco(ELECTRON_EFF_ ## DET ## _MEDIUM, ELECTRON_SMEAR_ ## DET, IDClass::MEDIUM) ; \
  setElectronReco(ELECTRON_EFF_ ## DET ## _TIGHT, ELECTRON_SMEAR_ ## DET, IDClass::TIGHT) ; \
  setMuonReco(MUON_EFF_ ## DET ## _LOOSE, MUON_SMEAR_ ## DET, IDClass::LOOSE) ; \
  setMuonReco(MUON_EFF_ ## DET ## _MEDIUM, MUON_SMEAR_ ## DET, IDClass::MEDIUM) ; \
  setMuonReco(MUON_EFF_ ## DET ## _TIGHT, MUON_SMEAR_ ## DET, IDClass::TIGHT) ; \
  setPhotonReco(PHOTON_EFF_ ## DET ## _LOOSE, ELECTRON_SMEAR_ ## DET, IDClass::LOOSE) ; \
  setPhotonReco(PHOTON_EFF_ ## DET ## _MEDIUM, ELECTRON_SMEAR_ ## DET, IDClass::MEDIUM) ; \
  setPhotonReco(PHOTON_EFF_ ## DET ## _TIGHT, ELECTRON_SMEAR_ ## DET, IDClass::TIGHT) ; \
  setTauReco(TAU_EFF_ ## DET ## _LOOSE, TAU_SMEAR_ ## DET, IDClass::LOOSE) ; \
  setTauReco(TAU_EFF_ ## DET ## _MEDIUM, TAU_SMEAR_ ## DET, IDClass::MEDIUM) ; \
  setTauReco(TAU_EFF_ ## DET ## _TIGHT, TAU_SMEAR_ ## DET, IDClass::TIGHT) ; \
  setTrkReco(TRK_EFF_ ## DET, TRK_SMEAR_ ## DET) ;                      \
  for (const std::string& jname : jetnames())                           \
    setJetReco(JET_SMEAR_ ## DET, JET_BTAG_ ## DET ## _ ## BTAG, jname) ; \
  setMETReco(MET_SMEARPARAMS_ ## DET, MET_SMEAR_ ## DET) ;                                \
  } while (0)