EnergySysts2017.cxx

Go to the documentation of this file.

#include "3FlavorAna/Systs/EnergySysts2017.h"

#include "StandardRecord/Proxy/SRProxy.h"

namespace ana
{
  const DirectHadEScaleSyst2017 kDirectHadEScaleSyst2017(0.05);
  const DirectRelHadEScaleSyst2017 kDirectRelHadEScaleSyst2017(0.05);

  // These are 2019 values, shoehorned into 2017 2-systematic style
  const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012);
  const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5 /* cm */);


  //----------------------------------------------------------------------
  void MuEScaleSyst2017::
  Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    if (sr->trk.kalman.tracks.size()==0) return;

    if (sr->hdr.det == caf::kFARDET){
      sr->trk.kalman.tracks[0].len = (1-fScale*sigma)*sr->trk.kalman.tracks[0].len;
    }
    else if (sr->hdr.det == caf::kNEARDET){
      sr->energy.numu.ndtrklenact = (1-fScale*sigma)*sr->energy.numu.ndtrklenact;
      sr->energy.numu.ndtrklencat = (1-fScaleMC*sigma)*sr->energy.numu.ndtrklencat;
    }
  }

  //----------------------------------------------------------------------
  void RelMuEScaleSyst2017::
  Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    if (sr->trk.kalman.tracks.size()==0) return;

    // shift Near Det *down*, shift Far Det *up*
    if (sr->hdr.det == caf::kNEARDET){
      sr->energy.numu.ndtrklenact *= 1-fScale/2*sigma;
      if(sr->energy.numu.ndtrklencat > 0)
        sr->energy.numu.ndtrklencat =
          std::max(sr->energy.numu.ndtrklencat - fShiftMuCat/2*sigma, 0.);
    }
    else if (sr->hdr.det == caf::kFARDET){
      sr->trk.kalman.tracks[0].len *= 1+fScale/2*sigma;
    }
  }

  //----------------------------------------------------------------------
  void DirectHadEScaleSyst2017::
  Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    sr->energy.numu.hadtrkE = (1+fScale*sigma)*sr->energy.numu.hadtrkE;
    sr->energy.numu.hadcalE = (1+fScale*sigma)*sr->energy.numu.hadcalE;
  }

  //----------------------------------------------------------------------
  void DirectRelHadEScaleSyst2017::
  Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    if (sr->hdr.det == caf::kNEARDET){
      sr->energy.numu.hadtrkE = (1+(fScale/2)*sigma)*sr->energy.numu.hadtrkE;
      sr->energy.numu.hadcalE = (1+(fScale/2)*sigma)*sr->energy.numu.hadcalE;
    }

    if (sr->hdr.det == caf::kFARDET){
      sr->energy.numu.hadtrkE = (1-(fScale/2)*sigma)*sr->energy.numu.hadtrkE;
      sr->energy.numu.hadcalE = (1-(fScale/2)*sigma)*sr->energy.numu.hadcalE;
    }
  }
}