EnergySysts2019.cxx

Go to the documentation of this file.

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

#include "StandardRecord/Proxy/SRProxy.h"

namespace ana
{
  // Error fully correlated between FD, main (non-MuCat) ND, MuCat
  const CorrMuEScaleSyst2019 kCorrMuEScaleSyst2019(0.0074, 0.0074, 0.0012);

  // Error for FD uncorrelated with main ND or muon catcher
  const UnCorrFDMuEScaleSyst2019 kUnCorrFDMuEScaleSyst2019(0.0015);

  // Error for main ND uncorrelated with FD or muon catcher
  const UnCorrNDMuEScaleSyst2019 kUnCorrNDMuEScaleSyst2019(0.0015);

  // Error for muon catcher uncorrelated with FD or main ND
  const UnCorrMuCatMuESyst2019 kUnCorrMuCatMuESyst2019(0.010);

  // One-sided error for neutron pile-up in the main ND and the muon catcher
  const PileupMuESyst2019 kPileupMuESyst2019(0.46, 1.3); // in cm

  // One-sided error for Bodek and Coulomb corrections
  const BolombMuESyst2019 kBolombMuESyst2019(-9.); // In MeV


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

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

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

    if(sr->hdr.det == caf::kFARDET)
      sr->trk.kalman.tracks[0].len *= 1+fScale*sigma;
  }

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

    if(sr->hdr.det == caf::kNEARDET)
      //                      act, not cat!
      sr->energy.numu.ndtrklenact *= 1+fScale*sigma;
  }

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

    if(sr->hdr.det == caf::kNEARDET && sr->energy.numu.ndtrklencat > 0){
      //                      cat, not act!
      sr->energy.numu.ndtrklencat *= 1+fScale*sigma;
    }
  }

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

    // This is a one-sided systematic, it can never *shorten* a track.
    // Implementing a one-sided error like this works fine, as it leaves
    // the objective function smooth and continuous.
    if(sigma < 0) return;

    if(sr->hdr.det == caf::kNEARDET){
      // Note subtle difference between "ndtrklenact" and "ndtrklencat"...
      if(sr->energy.numu.ndtrklenact > 0 && sr->energy.numu.ndtrklencat == 0){
        // Add for both ends of the track in the main ND
        sr->energy.numu.ndtrklenact += 2*fShiftND*sigma;
      }
      else if(sr->energy.numu.ndtrklenact == 0 && sr->energy.numu.ndtrklencat > 0){
        // Add for both ends of the track in the Muon Catcher
        sr->energy.numu.ndtrklencat += 2*fShiftMuCat*sigma;
      }
      else if(sr->energy.numu.ndtrklenact > 0 && sr->energy.numu.ndtrklencat > 0){
        // Add for one end in the main detector and one in the Muon Catcher
        sr->energy.numu.ndtrklenact += fShiftND   *sigma;
        sr->energy.numu.ndtrklencat += fShiftMuCat*sigma;
      }
      // and of couse if both are zero because there's no track, do nothing
    }
  }

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

    // This is a one-sided systematic, it can never *lengthen* a track.
    if(sigma < 0) return;

    // The systematic has to be applied as a length, not as an energy, so
    // here's the conversion constant, which is precise enough for these
    // purposes.  It is as per plots in doc-26649, discarding subtlety.  This
    // is the full physical length in the detector, averaged over air and
    // whatnot.  Since this correction pertains to the interaction vertex, we
    // don't need to worry about the Muon Catcher.
    const double cm_per_mev = 0.44;

    const double cmshift = fShift * cm_per_mev;

    // In principle, this error applies equally to nu_e and nu_tau, but no
    // attempt is made here to adjust anything but the muon.
    if(sr->hdr.det == caf::kNEARDET){
      sr->energy.numu.ndtrklenact =
        std::max(sr->energy.numu.ndtrklenact + cmshift*sigma, 0.);
    }
    else if(sr->hdr.det == caf::kFARDET){
      sr->trk.kalman.tracks[0].len =
        std::max(sr->trk.kalman.tracks[0].len + cmshift*sigma, 0.);
    }
  }
}