ExptResolutionFunctions.hh

// -*- C++ -*-
#ifndef RIVET_ExptResolutionFunctions_HH
#define RIVET_ExptResolutionFunctions_HH
 
#include "Rivet/Tools/ResolutionFunctions.hh"
 
#include "YODA/IO.h"
#include "YODA/Estimate.h"
#include "YODA/BinnedEstimate.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