ExptResolutionFunctions.hh

// -*- C++ -*-
#ifndef RIVET_ExptResolutionFunctions_HH
#define RIVET_ExptResolutionFunctions_HH
 
#include "Rivet/Tools/ResolutionFunctions.hh"
 
#include "YODA/BinnedEstimate.h"
#include "YODA/Estimate.h"
#include "YODA/IO.h"
 
namespace Rivet {
 
  class ATLAS_RUN2_EMTOPO_PT_RES : public ResolutionFunctor<Jet> {
  public:

    ATLAS_RUN2_EMTOPO_PT_RES() {
      _refFileName = "ExptResolutionFunctions/ATLAS_RUN2_RES_JET";
      _h["hist_JERMC16EMTopo"] = refData<YODA::Estimate2D>("hist_JERMC16EMTopo");
      _h["hist_JERData16EMTopo"] = refData<YODA::Estimate2D>("hist_JERData16EMTopo");
    }
 
    // Get resolution based on jet et, abseta.
    // Return max of MC/Data results.
    double resolution(const double et, const double abseta) const {
      const double et_adj = max(min(et, 3000), 17);
      const double abseta_adj = min(abseta, 4.5);
 
      return max(_h.at("hist_JERMC16EMTopo").binAt(et_adj, abseta_adj).val(),
                 _h.at("hist_JERData16EMTopo").binAt(et_adj, abseta_adj).val());
    }

    double resolution(const Jet& j) const override {
      return resolution(j.Et(), j.abseta());
    }
 
  private:
 
    map<string, YODA::Estimate2D> _h;
  };
 
  class ATLAS_RUN2_JET_PHI_RES : public ResolutionFunctor<Jet> {
  public:

    ATLAS_RUN2_JET_PHI_RES() {
      _refFileName = "ExptResolutionFunctions/ATLAS_RUN2_RES_JET";
      for (const string s : {"20", "50", "100"}) {
        const string histname = "hist_jetphi" + s;
        _h[histname] = refData<YODA::Estimate2D>(histname);
      }
    }

    double getRes(const double pt, const double eta, const double phi) const {
      if (pt < 20 * GeV)
        return 0.;
      else if (pt < 50 * GeV)
        return _h.at("hist_jetphi20").binAt(eta, phi).val();
      else if (pt < 100 * GeV)
        return _h.at("hist_jetphi50").binAt(eta, phi).val();
      return _h.at("hist_jetphi100").binAt(eta, phi).val();
    }

    double resolution(const Jet& j) const override {
      return getRes(j.pt(), j.eta(), j.phi(MINUSPI_PLUSPI));
    }
 
  private:
 
    map<string, YODA::Estimate2D> _h;
  };

  class ATLAS_RUN2_ELECTRON_ET_RES : public ResolutionFunctor<Particle> {
  public:
 
    ATLAS_RUN2_ELECTRON_ET_RES() {
      _refFileName = "ExptResolutionFunctions/ATLAS_RUN2_RES_ELECTRON";
      // TODO: These histograms are small enough (17 bins IIRC) that we could hardcode them if preferred.
      _h["hist_electronNoise90"] = refData("hist_electronNoise90");
      _h["hist_electronConst90"] = refData("hist_electronConst90");
      _h["hist_electronSampling90"] = refData("hist_electronSampling90");
    }

    double resolution(const double et, const double e, const double abseta) const {
      // This et range seems almost pointlessly small but hey ho.
      const double et_adj = max(min(abseta, 50), 5);
 
      const double rsampling = _h.at("hist_electronSampling90").binAt(abseta).val();
      const double rnoise = _h.at("hist_electronNoise90").binAt(abseta).val();
      const double rconst = _h.at("hist_electronConst90").binAt(abseta).val();
      const double sigma2 = rsampling * rsampling / e + rnoise * rnoise / e / e + rconst * rconst;
      const double pileupNoiseMeV = sqrt(32.) * (60. + 40. * log(et_adj / 10.) / log(5.));
      const double pileupSigma2 = (pileupNoiseMeV / 1000. / et_adj) * (pileupNoiseMeV / 1000. / et_adj);
      return sqrt(sigma2 + pileupSigma2);
    }

    double resolution(const Particle& e) const override {
      // TODO: In typical ATLAS fashion, there's a good chance they mean pT
      return resolution(e.Et(), e.E(), e.abseta());
    }
 
  private:
 
    map<string, YODA::Estimate1D> _h;
  };

  class ATLAS_RUN2_MUON_PT_RES : public ResolutionFunctor<Particle> {
  public:
 
    ATLAS_RUN2_MUON_PT_RES() {
      _refFileName = "ExptResolutionFunctions/ATLAS_RUN2_RES_MUON";
      for (const string& s : _histosbarrel) {
        _h[s] = refData<YODA::Estimate2D>("hist_" + s);
      }
      for (const string& s : _histosendcap) {
        _h[s] = refData<YODA::Estimate2D>("hist_" + s);
      }
    }

    double resolution(const double pT2, const double eta, const double phi) const {
      vector<string> whichhistos = abs(eta) < 1.05 ? _histosbarrel : _histosendcap;
      vector<double> pars = {0., 0., 0., 0., 0.};
      for (size_t i = 0; i < 5; ++i) {
        pars[i] = _h.at(whichhistos[i]).binAt(phi, eta).val();
      }
      double IDResSq = pars[0] * pars[0] + pars[1] * pars[1] * pT2;
      double MSResSq = pars[2] * pars[2] / pT2 + pars[3] * pars[3] + pars[4] * pars[4] * pT2;
      return sqrt(IDResSq * MSResSq / (IDResSq + MSResSq));
    }

    double resolution(const Particle& mu) const override {
      return resolution(mu.pt2(), mu.eta(), mu.phi());
    }
 
  private:
 
    map<string, YODA::Estimate2D> _h;
    std::vector<string> _histosbarrel = {
        "r1_ID_MC_BARREL", "r2_ID_MC_BARREL", "r0_MS_MC_BARREL", "r1_MS_MC_BARREL", "r2_MS_MC_BARREL",
    };
    std::vector<string> _histosendcap = {"r1_ID_MC_ENDCAP", "r2_ID_MC_ENDCAP", "r0_MS_MC_ENDCAP",
                                         "r1_MS_MC_ENDCAP", "r2_MS_MC_ENDCAP"};
  };

  class ATLAS_RUN2_PHOTON_ET_RES : public ResolutionFunctor<Particle> {
  public:
 
    ATLAS_RUN2_PHOTON_ET_RES() {
      _refFileName = "ExptResolutionFunctions/ATLAS_RUN2_RES_PHOTON";
      _h["hist_photonConst90"] = refData<YODA::Estimate1D>("hist_photonConst90");
      _h["hist_photonNoise90"] = refData<YODA::Estimate1D>("hist_photonNoise90");
      _h["hist_photonSampling90"] = refData<YODA::Estimate1D>("hist_photonSampling90");
    }

    double resolution(const double E, const double pt, const double abseta) const {
      const double rsampling = abseta < 2.4 ? _h.at("hist_photonSampling90").binAt(abseta).val() : 0;
      const double rnoise = abseta < 2.4 ? _h.at("hist_photonNoise90").binAt(abseta).val() : 0;
      const double rconst = abseta < 2.4 ? _h.at("hist_photonConst90").binAt(abseta).val() : 0;
 
      double sigma2 = rsampling * rsampling / E + rnoise * rnoise / (E * E) + rconst * rconst;
 
      double et = std::min(50., std::max(5., pt)); // constraint 5-50 GeV
 
      double pileupNoiseMeV = sqrt(32.) * (60. + 40. * log(et / 10.) / log(5.));
      double pileupSigma2 = (pileupNoiseMeV / 1000. / pt) * (pileupNoiseMeV / 1000. / pt);
      // not clear why Egamma group uses Et and not E here? (comment left in from SA code)
      return sqrt(sigma2 + pileupSigma2);
    }
 
    double resolution(const Particle& gamma) const override {
      return resolution(gamma.E(), gamma.pt(), gamma.abseta());
    }
 
  private:
 
    map<string, YODA::Estimate1D> _h;
  };
 

  class ATLAS_RUN2_TAU_1p0n_PT_RES : public ResolutionFunctor<Jet> {
  public:
 
    ATLAS_RUN2_TAU_1p0n_PT_RES() {
      _refFileName = "ExptResolutionFunctions/ATLAS_RUN2_RES_TAU_1p0n";
      for (size_t i = 0; i < 5; ++i) {
        _h[i] = refData<YODA::Estimate1D>("hist_tauRes1p0nBin"s + to_string(i));
      }
    }

    double resolution(const double pt, const double abseta) const {
      const double ptMeV = 1000. * std::min(499., std::max(15., pt)); // only defined for 15-499 GeV
      const vector<double> etaEdges = {0.0, 0.3, 0.8, 1.3, 1.6};
      const double res = _h[binIndex(abseta, etaEdges, true)].binAt(ptMeV).val();
      return res;
    }
 
    double resolution(const Jet& hadTau) const override {
      // TODO: Should we consider an extra check to make sure it is a 1p0n decay?
      // And maybe a tau-tag check?
      return resolution(hadTau.pt(), hadTau.abseta());
    }
 
  private:
 
    vector<YODA::Estimate1D> _h;
  };
 


  class ATLAS_RUN2_MUON_PHI_RES : public ResolutionFunctor<Particle> {
  public:
 
    ATLAS_RUN2_MUON_PHI_RES() { }
    double resolution() const {
      return 0.001;
    }
    double resolution(const Particle& e) const override {
      return resolution();
    }
  };

  class ATLAS_RUN2_ELECTRON_PHI_RES : public ResolutionFunctor<Particle> {
  public:
 
    ATLAS_RUN2_ELECTRON_PHI_RES() { }
    double resolution() const {
      return 0.004;
    }
    double resolution(const Particle& e) const override {
      return resolution();
    }
  };
 

  class ATLAS_RUN2_PHOTON_PHI_RES : public ResolutionFunctor<Particle> {
  public:
 
    ATLAS_RUN2_PHOTON_PHI_RES() { }
    double resolution() const {
      return 0.004;
    }
    double resolution(const Particle& e) const override {
      return resolution();
    }
  };

  class ATLAS_RUN2_TAU_PHI_RES : public ResolutionFunctor<Jet> {
  public:
 
    ATLAS_RUN2_TAU_PHI_RES() { }
    double resolution() const {
      return 0.01;
    }
    double resolution(const Jet& j) const override {
      return resolution();
    }
  };

 
}
 
#endif