NumuCCIncVars.cxx

Go to the documentation of this file.

#include "NDAna/numucc_inc/NumuCCIncVars.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "sys/stat.h"
#include "TCanvas.h"
#include "TFile.h"

namespace ana{

  TH2F *hprob_dedx_vs_scat_sig;
  TH2F *hprob_dedx_vs_scat_bg;

  float MuonLLR(const float dedxll,
                const float scatll)
  {

    if(!hprob_dedx_vs_scat_sig || !hprob_dedx_vs_scat_bg){
      const char* path = getenv("SRT_PUBLIC_CONTEXT");
      std::string fname = std::string(path) + "/NDAna/numucc_inc/muon_pdf.root";

      TFile *muon_prob_file = new TFile( fname.c_str(), "r");

      hprob_dedx_vs_scat_sig = (TH2F*)muon_prob_file->Get("prob_sig");
      hprob_dedx_vs_scat_bg  = (TH2F*)muon_prob_file->Get("prob_bg");

    }

    const int dedx_vs_scat_nxbins = hprob_dedx_vs_scat_bg->GetXaxis()->GetNbins();
    const int dedx_vs_scat_nybins = hprob_dedx_vs_scat_bg->GetYaxis()->GetNbins();

    int xbin_sig = hprob_dedx_vs_scat_sig->GetXaxis()->FindBin(dedxll);
    int ybin_sig = hprob_dedx_vs_scat_sig->GetYaxis()->FindBin(scatll);

    int xbin_bg = hprob_dedx_vs_scat_bg->GetXaxis()->FindBin(dedxll);
    int ybin_bg = hprob_dedx_vs_scat_bg->GetYaxis()->FindBin(scatll);

    float llr = -1e12;
    // if the bins are underflow or overflow, it's beyond the range of our model
    // skip over to the next track.
    if (xbin_sig == 0 || xbin_sig > dedx_vs_scat_nxbins ||
        ybin_sig == 0 || ybin_sig > dedx_vs_scat_nybins ||
        xbin_bg == 0  || xbin_bg  > dedx_vs_scat_nxbins ||
        ybin_bg == 0  || ybin_bg  > dedx_vs_scat_nybins)
      return llr;

    float psig = std::log(hprob_dedx_vs_scat_sig->GetBinContent(xbin_sig, ybin_sig));
    float pbg  = std::log(hprob_dedx_vs_scat_bg->GetBinContent(xbin_bg, ybin_bg));

    llr = psig - pbg;
    return llr;
  }


  //---------------------------------
  // Muon Energy for Active region
  //---------------------------------

  float MuonEAct(double trklenact){

    double p0    =  1.67012e-01;
    double p1    =  1.79305e-01;
    double p2    =  3.74708e-03;
    double p3    = -1.54232e-04;

    float MuonE = 0.0;

    if (trklenact <= 0.0) return 0.0;

    MuonE = p0
      + p1 * trklenact
      + p2 * std::pow(trklenact, 2)
      + p3 * std::pow(trklenact, 3);

    return MuonE;
  }


  //---------------------------------------
  // Muon Energy for Active+Catcher region
  //(for Active region Track length)
  //---------------------------------------

  float MuonEActandCat(double trklenactandcat){

    double  p0   =  1.21130e-02;
    double  p1   =  1.97903e-01;
    double  p2   =  7.82459e-04;

    float MuonE = 0.0;

    if (trklenactandcat <= 0.0) return 0.0;

    MuonE = p0
      + p1 * trklenactandcat
      + p2 * std::pow(trklenactandcat, 2);

    return MuonE;
  }


  //---------------------------------------
  // Muon Energy for Active+Catcher region
  //(for catcher region Track length)
  //---------------------------------------

  float MuonECat(double trklencat){

    double  offset  =  1.31325e-01;
    double  slope  =  5.35146e-01;

    float MuonE = 0.0;

    if (trklencat <= 0.0) return 0.0;

    MuonE = slope*trklencat + offset;

    return MuonE;
  }

  //---------------------------------------
  // Haronic Energy
  //----------------------------------------

  float VisibleHadE(float vishadE) {

    double  p0    =   5.85254e-02;
    double  p1    =   1.27796e+00;
    double  p2    =   3.75457e-01;
    double  p3    =  -5.45618e-01;
    double  p4    =   1.65975e-01;

    float HadE = 0.0;

    HadE = p0
      + p1 * vishadE
      + p2 * std::pow(vishadE, 2)
      + p3 * std::pow(vishadE, 3)
      + p4 * std::pow(vishadE, 4);

    return HadE;
  }


  //-------------------------------------------------------------------------------------------------
  // PID related Vars

  /// Return the muon to non-muon likelihod ration, given
  /// the dedx and scattering likelihood ratios for the track
  float MuonLLR(const float dedxll, const float scatll);

  const Var kMostMuonLikeTrackIdx([](const caf::SRProxy* sr)
  {
    float maxllr = -1e12;
    int bestidx = -5;
    // find the track with max LLR for being a muon
    for(unsigned int itrk = 0; itrk < sr->trk.kalman.ntracks; itrk++){
      float llr = MuonLLR(sr->trk.kalman.tracks[itrk].dedxllh,
                          sr->trk.kalman.tracks[itrk].scatllh);
      if(llr > maxllr){
        maxllr  = llr;
        bestidx =  itrk;
      }
    } // end loop over tracks
    return bestidx;
  });

  const Var kMuonLLRPID([](const caf::SRProxy* sr)
  {
    int bestidx = kMostMuonLikeTrackIdx(sr);
    float llr = -1e12;

    if(bestidx < 0 || bestidx >= int(sr->trk.kalman.ntracks))
      return llr;

    // if we get here, the bestidx is a valid index into the tracks vector
    llr =  MuonLLR(sr->trk.kalman.tracks[bestidx].dedxllh,
    sr->trk.kalman.tracks[bestidx].scatllh);
    return llr;
  });

  //----------------------------
  // Energy Estimation
  //----------------------------

  float MuonEActandCat(double trklenactandcat);
  float MuonEAct(double trklenact);
  float MuonECat(double trklencat);
  float VisibleHadE(float vishadE);

  const Var kMuStartX([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    return sr->trk.kalman.tracks[ibesttrk].start.X();
  });


  const Var kMuStartY([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    return sr->trk.kalman.tracks[ibesttrk].start.Y();
  });


  const Var kMuStartZ([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    return sr->trk.kalman.tracks[ibesttrk].start.Z();
  });

  const Var kMuStopX([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    return sr->trk.kalman.tracks[ibesttrk].stop.X();
  });


  const Var kMuStopY([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    return sr->trk.kalman.tracks[ibesttrk].stop.Y();
  });


  const Var kMuStopZ([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    return sr->trk.kalman.tracks[ibesttrk].stop.Z();
  });


  const Var kMuonLength([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    return float(sr->trk.kalman.tracks[ibesttrk].len/100.);
  });

  const Var kVisHadE([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    // Extra energy (hadronic contamination) associated with muon track
    float ExtraHadE = sr->trk.kalman.tracks[ibesttrk].overlapE;

    // calorimetric slice  energy - Calorimetric Energy of Muon Track
    float CalhadE = sr->slc.calE - sr->trk.kalman.tracks[ibesttrk].calE;

    // Add calorimetric hadE and Overlap energy in the muon track
    float vishadE = CalhadE + ExtraHadE;

    return vishadE;
  });

  const Var kHadEonTrack([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    // Extra energy (hadronic contamination) associated with muon track
    return (float)sr->trk.kalman.tracks[ibesttrk].overlapE;
  });


  const Var kAveragededx40([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    if (sr->trk.kalman.tracks[ibesttrk].avedEdxlast40cm == -5)
      return -1000.f;

    return float(sr->trk.kalman.tracks[ibesttrk].avedEdxlast40cm);
  });

//==============================================================================
// Redefine track lengths in active and catcher volumes with the new PID.
// Note that at this moment what is called leninact is actually distance from
// the start of the track to the transition plane. If the transition plane is
// outside of the track, that number is the track length plus the length of the
// projection from the last hit along the direction determined by the last two
// hits to the transition plane. The lenincat behaves similarly.
// In the next production the two variables should be changed to what they are
// supposed to mean literally.

  const Var kTrkLenAct([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);
    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    // If leninact is positive and lenincat is negative,
    // the transition plane is to the right of the track.
    if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
       sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
      return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
                  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));

    // If leninact is positive and lenincat is positive,
    // the transition plane is in the middle of the track.
    if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
       sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
      return float(sr->trk.kalman.tracks[ibesttrk].leninact / 100.);

    // If leninact is negative and lenincat is positive,
    // the transition plane is to the left of the track.
    if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
       sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
      return 0.f;

    return -1000.f;
  });

  const Var kTrkLenCat([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);
    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    // If leninact is positive and lenincat is negative,
    // the transition plane is to the right of the track.
    if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
       sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
      return 0.f;

    // If leninact is positive and lenincat is positive,
    // the transition plane is in the middle of the track.
    if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
       sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
      return float(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.);

    // If leninact is negative and lenincat is positive,
    // the transition plane is to the left of the track.
    if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
       sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
      return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
                  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));

    return -1000.f;
  });


  // Muon Energy for Inclusive numuCC
  //------------------------------------
  const Var kIncXsecMuonE([](const caf::SRProxy *sr)
  {
    float muonE          = 0.0;
    float muonEact       = 0.0;
    float muonEcat       = 0.0;
    float muonEactandcat = 0.0;

    float trkLenAct = 0.f;
    float trkLenCat = 0.f;

    trkLenAct = kTrkLenAct(sr);
    trkLenCat = kTrkLenCat(sr);

    int ibesttrk = kBestTrack(sr);
    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
       sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
      muonEact = MuonEAct(trkLenAct);
    else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
            sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
      muonEcat = MuonECat(trkLenCat);
      muonEactandcat = MuonEActandCat(trkLenAct);
      muonE = muonEactandcat + muonEcat;
    }

    return muonE + muonEact;
  });


  // Hadronic Energy for Inclusive numuCC
  //---------------------------------------

  const Var kIncXsecHadE([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    // Extra energy (hadronic contamination) associated with muon track
    float ExtraHadE = sr->trk.kalman.tracks[ibesttrk].overlapE;

    // calorimetric slice  energy - Calorimetric Energy of Muon Track
    float CalhadE = sr->slc.calE - sr->trk.kalman.tracks[ibesttrk].calE;

    // Add calorimetric hadE and Overlap energy in the muon track
    float vishadE = CalhadE + ExtraHadE;

    float hadE = 0.0;
    hadE = VisibleHadE(vishadE);

    return hadE;
  });

  //-------------------------------------------------------------------------------------------------
  // Standard analysis vars and axes

  // Reconstructed candidate muon angle wrt beam
  const Var kRecoMuCostheta([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    TVector3 dir = sr->trk.kalman.tracks[ibesttrk].dir;
    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
    return (float)dir.Dot(beamdir);
  });

  // Angle systematics (+/- 0.0025 mrad wrt beam)
  const Var kRecoMuTheta([](const caf::SRProxy* sr)
  {
    int ibesttrk = kBestTrack(sr);

    if (sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if (ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    TVector3 dir = sr->trk.kalman.tracks[ibesttrk].dir;
    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
    return (float)dir.Angle(beamdir);
  });

  const Var kRecoMuCosthetaUp([](const caf::SRProxy* sr)
  {
    return (float)std::cos(kRecoMuThetaUp(sr));
  });

  const Var kRecoMuCosthetaDw([](const caf::SRProxy* sr)
  {
    return (float)std::cos(kRecoMuThetaDw(sr));
  });

  const NuTruthVar kTrueEST([](const caf::SRNeutrinoProxy* nu)
  {
    return nu->E;
  });

  const NuTruthVar kTrueMuEST([](const caf::SRNeutrinoProxy* nu)
  {
    float  MuE = -5.;

    if(abs(nu->pdg) != 14 || !nu->iscc)
      return MuE;

    int nprims = nu->prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(nu->prim[iprim].pdg) == 13){
        MuE = nu->prim[iprim].p.T();
      }
    }//end loop over primaries

    return MuE;
  });

  const NuTruthVar kTrueMuKEST([](const caf::SRNeutrinoProxy* nu)
  {

    float  ke = -5;

    if(abs(nu->pdg) != 14 || !nu->iscc)
      return ke;

    int nprims = nu->prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(nu->prim[iprim].pdg) == 13){
        double E= nu->prim[iprim].p.T();
        ke = E-MuonMass();
      }
    }//end loop over primaries

    return ke;
  });

  // True muon cos theta wrt beam direction

  const NuTruthVar kTrueMuCosthetaST([](const caf::SRNeutrinoProxy* nu)
  {
    if(abs(nu->pdg) != 14 || !nu->iscc)
      return -5.0;

    int nprims = nu->prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(nu->prim[iprim].pdg) == 13){
        TVector3 mudir = nu->prim[iprim].p.Vect();
        TVector3 beamdir = NuMIBeamDirection(caf::kNEARDET);

        return mudir.Unit().Dot(beamdir.Unit());
      }
    }//end loop over primaries
    return -5.0;
  });

  // True hadronic y
  const NuTruthVar kTrueYST([](const caf::SRNeutrinoProxy* sr)
  {
    return float(sr->y);
  });


  // True hadronic energy
  const NuTruthVar kTrueHadEST([](const caf::SRNeutrinoProxy* sr)
  {
    return float(sr->y * sr->E);
  });

  // Reconstructed available energy (Determined by Travis Olson)
  const Var kRecoEavail([] (const caf::SRProxy *sr)
  {
    return 1.68*kNumuHadVisE(sr)+0.02*pow(kNumuHadVisE(sr),2);
  });



  const Var kRecoMuKEVsCos = Var2D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKE,
    mukebins);

  const NuTruthVar kTrueMuKEVsCosST = Var2D(
    kTrueMuCosthetaST,
    angbinsCustom,
    kTrueMuKEST,
    mukebins);

  const NuTruthVar kTrueMuKEVsCosVsEnuST = Var3D(
    kTrueMuCosthetaST,
    angbinsCustom,
    kTrueMuKEST,
    mukebins,
    kTrueEST,
    enubins);

  const Var kRecoMuKEVsCosVsEnu = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKE,
    mukebins,
    kRecoE,
    enubins);

  const NuTruthVar kTrueMuKEVsCosVsEavailST = Var3D(
    kTrueMuCosthetaST,
    angbinsCustom,
    kTrueMuKEST,
    mukebins,
    kTrueEavailST,
    eavailbins);

  const Var kRecoMuKEVsCosVsEavail = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKE,
    mukebins,
    kRecoEavail,
    eavailbins);

/////////////////////////////////////////
/// \brief Systematic Shift Vars and HistAxes
/////////////////////////////////////////

  // Vars related to muon energy scale systematics:
  const Var kIncXsecMuonEUp([](const caf::SRProxy *sr)
  {
    //factors to make the muon energy up
    float mue_act_up = 1.0079;
    float mue_cat_up = 1.0120;

    float muonE          = 0.0;
    float muonEact       = 0.0;
    float muonEcat       = 0.0;
    float muonEactandcat = 0.0;

    float trkLenAct = 0.f;
    float trkLenCat = 0.f;

    trkLenAct = kTrkLenAct(sr);
    trkLenCat = kTrkLenCat(sr);

    int ibesttrk = kBestTrack(sr);
    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    if( sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
        sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
      muonEact = mue_act_up * MuonEAct(trkLenAct);
    else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
            sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
      muonEcat = MuonECat(trkLenCat);
      muonEactandcat = MuonEActandCat(trkLenAct);
      muonE = mue_act_up*muonEactandcat + mue_cat_up*muonEcat;
    }

    return muonE + muonEact;
  });

  const Var kIncXsecMuonEDw([](const caf::SRProxy *sr)
  {
    //factors to make the muon energy down
    float mue_act_dw = 0.9921;
    float mue_cat_dw = 0.9880;

    float muonE          = 0.0;
    float muonEact       = 0.0;
    float muonEcat       = 0.0;
    float muonEactandcat = 0.0;

    float trkLenAct = 0.f;
    float trkLenCat = 0.f;

    trkLenAct = kTrkLenAct(sr);
    trkLenCat = kTrkLenCat(sr);

    int ibesttrk = kBestTrack(sr);
    if(sr->trk.kalman.ntracks < 1)
      return -1000.f;

    if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
      return -1000.f;

    if( sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
        sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
      muonEact = mue_act_dw * MuonEAct(trkLenAct);
    else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
      sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
      muonEcat = MuonECat(trkLenCat);
      muonEactandcat = MuonEActandCat(trkLenAct);
      muonE = mue_act_dw*muonEactandcat + mue_cat_dw*muonEcat;
    }

    return muonE + muonEact;
  });

  // 3D vars for muon energy shift (Enu as 3rd axis)
  const Var kRecoMuKEVsCosVsEnu_onlyMuKEUp = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKEUp,
    mukebins,
    kRecoE,
    enubins);
  const Var kRecoMuKEVsCosVsEnu_onlyMuKEDw = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKEDw,
    mukebins,
    kRecoE,
    enubins);
  const Var kRecoMuKEVsCosVsEnuUp = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKEUp,
    mukebins,
    kRecoEUp,
    enubins);
  const Var kRecoMuKEVsCosVsEnuDw = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKEDw,
    mukebins,
    kRecoEDw,
    enubins);

  // 3D vars for muon energy shift (Eavail as 3rd axis)
  const Var kRecoMuKEVsCosVsEavailUp = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKEUp,
    mukebins,
    kRecoEavail,
    eavailbins);
  const Var kRecoMuKEVsCosVsEavailDw = Var3D(
    kRecoMuCostheta,
    angbinsCustom,
    kRecoMuKEDw,
    mukebins,
    kRecoEavail,
    eavailbins);

  // 3D vars for muon angle shift (Eavail as 3rd axis)
  const Var kRecoMuKEVsCosVsEavailAngleUp = Var3D(
    kRecoMuCosthetaUp,
    angbinsCustom,
    kRecoMuKE,
    mukebins,
    kRecoEavail,
    eavailbins);

  const Var kRecoMuKEVsCosVsEavailAngleDw = Var3D(
    kRecoMuCosthetaDw,
    angbinsCustom,
    kRecoMuKE,
    mukebins,
    kRecoEavail,
    eavailbins);

  //-------------------------------------------------------------------------------------------------


  const Var kDummy([](const caf::SRProxy* sr)
  {
    return 1;
  });



  const Var ksigResdown([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu == 0)
      return 0.;

    assert(sr->mc.nnu == 1);

    if( (sr->mc.nu[0].iscc) && sr->mc.nu[0].mode==caf::kRes)
      return 0.909090;

    else return 1.;
  });

  const Var ksigDISup([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu == 0)
      return 0.;

    assert(sr->mc.nnu == 1);

    if( (sr->mc.nu[0].iscc) && sr->mc.nu[0].mode==caf::kDIS)
      return 1.1;

    else return 1.;
  });

  const Var kbkgDISdown([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu == 0)
      return 0.;

    assert(sr->mc.nnu == 1);

    if( (!sr->mc.nu[0].iscc) && sr->mc.nu[0].mode==caf::kDIS)
      return 0.909090;

    else return 1.;
  });


  //Number of prongs
  const Var kNProngs([](const caf::SRProxy* sr)
  {
    return (sr->vtx.elastic[0].fuzzyk.npng);
  });

  // True Q2

  const Var kQsqr([](const caf::SRProxy* sr)
  {
    return (sr->mc.nnu == 0) ? 0. : float(sr->mc.nu[0].q2);
  });

  // True W2

  const Var kWsqr([](const caf::SRProxy* sr)
    {
      return (sr->mc.nnu == 0) ? 0. : float(sr->mc.nu[0].W2);
    });

  const Var kweightedten([](const caf::SRProxy*){return 1.1;});

 // Primary muon momentum

  const Var kTrueMuP([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0;

    if(abs(sr->mc.nu[0].pdg) != 14 || !sr->mc.nu[0].iscc)
      return -5.0;
    //     unsigned int ipartmu = 0;
    int nprims = sr->mc.nu[0].prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
      //  ipartmu = iprim;
        return sr->mc.nu[0].prim[iprim].p.Vect().Mag();
      }
    }//end loop over primaries
    return -5.0;
  });


 // Primary muon transverse momentum wrt beam direction

  const Var kTrueMuPt([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0;

    if(abs(sr->mc.nu[0].pdg) != 14 || !sr->mc.nu[0].iscc)
      return -5.0;

    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
    int nprims = sr->mc.nu[0].prim.size();

    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
        TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
        return (p-p.Dot(beamdir)*beamdir).Mag();
      }
    }//end loop over primaries
    return -5.0;
  });

 // Primary muon momentum along the beam direction

  const Var kTrueMuPz([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0;

    if(abs(sr->mc.nu[0].pdg) != 14 || !sr->mc.nu[0].iscc)
      return -5.0;

    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
    int nprims = sr->mc.nu[0].prim.size();

    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
        TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
        return p.Dot(beamdir);
      }
    }//end loop over primaries
    return -5.0;
  });

  // Reconstructed muon momentum
  const Var kRecoMuPsqr([](const caf::SRProxy * sr){
    float muonpsqr = std::pow(kRecoMuE(sr), 2) - kMuonMass2(sr);
    return (std::max(muonpsqr, (float)0.0)); // Return muonpsqr if > 0, otherwise 0.
  });

  const Var kRecoMuPsqrUp([](const caf::SRProxy * sr){
    float muonpsqr = std::pow(kRecoMuEUp(sr), 2) - kMuonMass2(sr);
    return (std::max(muonpsqr, (float)0.0)); // Return muonpsqr if > 0, otherwise 0.
  });

  const Var kRecoMuPsqrDw([](const caf::SRProxy * sr){
    float muonpsqr = std::pow(kRecoMuEDw(sr), 2) - kMuonMass2(sr);
    return (std::max(muonpsqr, (float)0.0)); // Return muonpsqr if > 0, otherwise 0.
  });


  // Reconstructed muon tranvserse momentum wrt beam direction

  const Var kRecoMuPt([](const caf::SRProxy* sr)
  {
    int bestidx = sr->trk.kalman.idxremid;
    int nkals = sr->trk.kalman.ntracks;

    if( nkals < 1 ||
        nkals < bestidx ||
        bestidx < 0 )
      return -5.0f;

    double E = TAMuE(sr->trk.kalman.tracks[0].len);

    TVector3 dir = sr->trk.kalman.tracks[0].dir;
    TVector3 pvec = TMath::Sqrt( E*E - util::sqr(MuonMass()))*dir;
    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);

    return float((pvec - pvec.Dot(beamdir)*beamdir).Mag());
  });


  // Reconstructed muon momentum along beam direction

  const Var kRecoMuPz([](const caf::SRProxy* sr)
  {

    double E = TAMuE(sr->trk.kalman.tracks[0].len);

    TVector3 dir = sr->trk.kalman.tracks[0].dir;
    TVector3 pvec = TMath::Sqrt( E*E - util::sqr(MuonMass()))*dir;
    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);

    return pvec.Dot(beamdir);
  });

  // Fractional muon momentum resolution
  const Var kResMuP([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0f;

    if( abs(sr->mc.nu[0].pdg) != 14 ||
        !sr->mc.nu[0].iscc)
      return -5.0f;
    int bestidx = sr->trk.kalman.idxremid;
    int nkals = sr->trk.kalman.ntracks;
    //int bestidx = sr->sel.remid.bestidx;
    //int nmups = sr->energy.mutrkE.size();
    if( nkals < 1 ||
        nkals < bestidx ||
        bestidx < 0 )
      return -5.0f;

    float mup = -5, muRecoE;
    int nprims = sr->mc.nu[0].prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
      return float(sr->mc.nu[0].prim[iprim].p.Vect().Mag());
      }
    // mup = sr->mc.nu[0].prim[iprim].p.P();
    // break;
    }//end loop over primaries

    if(mup == -5)
      return -5.0f;

    muRecoE = TAMuE(sr->trk.kalman.tracks[0].len);

    double recop = TMath::Sqrt( muRecoE*muRecoE - util::sqr(MuonMass()));

    return float((recop-mup)/mup);
  });


  // Fractional muon transverse momentum resolution
  const Var kResMuPt([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0f;

    if( abs(sr->mc.nu[0].pdg) != 14 ||
        !sr->mc.nu[0].iscc)
      return -5.0f;
    int bestidx = sr->trk.kalman.idxremid;
    int nkals = sr->trk.kalman.ntracks;
    //int bestidx = sr->sel.remid.bestidx;
    //int nmups = sr->energy.mutrkE.size();
    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);

    if( nkals < 1 ||
        nkals < bestidx ||
        bestidx < 0 )
    return -5.0f;

    float mupt = -5, muRecoE;
    int nprims = sr->mc.nu[0].prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
        TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
        mupt = (p-p.Dot(beamdir)*beamdir).Mag();
      }
      break;
    }//end loop over primaries

    if(mupt == -5)
      return -5.0f;

    muRecoE = TAMuE(sr->trk.kalman.tracks[0].len);

    TVector3 dir = sr->trk.kalman.tracks[0].dir;
    TVector3 recop = TMath::Sqrt( muRecoE*muRecoE - util::sqr(MuonMass()))*dir;
    double recopt = (recop-recop.Dot(beamdir)*beamdir).Mag();
    return float((recopt-mupt)/mupt);
  });


  // Fractional resolution of muon momentum along beam direction
  const Var kResMuPz([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0f;

    if( abs(sr->mc.nu[0].pdg) != 14 ||
        !sr->mc.nu[0].iscc)
      return -5.0f;

    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);

    float mupz = -5, muRecoE;
    int nprims = sr->mc.nu[0].prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
        TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
        mupz = p.Dot(beamdir);
      }
      break;
    }//end loop over primaries

    if(mupz == -5)
      return -5.0f;

    muRecoE = TAMuE(sr->trk.kalman.tracks[0].len);

    TVector3 dir = sr->trk.kalman.tracks[0].dir;
    TVector3 recop = TMath::Sqrt( muRecoE*muRecoE - util::sqr(MuonMass()))*dir;

    double recopz = recop.Dot(beamdir);
    return float((recopz-mupz)/mupz);
  });


  // Remid
  const Var kRemid([](const caf::SRProxy* sr){
    return (sr->sel.remid.pid);
  });

  // Remid variables
  const Var kScatLL = SIMPLEVAR(sel.remid.scatllh);
  const Var kDedxLL = SIMPLEVAR(sel.remid.dedxllh);
  const Var kMeasfrac = SIMPLEVAR(sel.remid.measfrac);

  // Length of candidate muon track
  const Var kMuLength([](const caf::SRProxy* sr)
  {
    //      return sr->trk.kalman[bestidx].len;
    return sr->trk.kalman.tracks[0].len/100;
  });

  // Energy per length of the candidate muon track
  const Var kEPerLength([](const caf::SRProxy* sr)
  {
    if (sr->trk.kalman.ntracks > 0 &&
        sr->trk.kalman.idxremid != 999)
      return TAMuE(sr->trk.kalman.tracks[0].len)/sr->trk.kalman.tracks[0].len;

    return -5.f;
  });


  // Resolution of candidate muon angle wrt beam
  const Var kResMuCostheta([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0;

    if( abs(sr->mc.nu[0].pdg) != 14 ||
        !sr->mc.nu[0].iscc)
      return -5.0;

    int bestidx = sr->trk.kalman.idxremid;
    int nkals = sr->trk.kalman.ntracks;
    if( nkals < 1 ||
        nkals < bestidx ||
        bestidx < 0 )
      return -5.0;

    double mu = 0, muReco = 0;
    TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
    int nprims = sr->mc.nu[0].prim.size();
    for(int iprim = 0; iprim < nprims; iprim++){
      if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
        TVector3 mudir = sr->mc.nu[0].prim[iprim].p.Vect();
        mu = beamdir.Dot(mudir);
        mu = beamdir.Dot(sr->mc.nu[0].prim[iprim].p.Vect().Unit());
      }
    }//end loop over primaries

    if(mu == 0)
      return -5.0;

    TVector3 dir = sr->trk.kalman.tracks[0].dir;

    muReco = beamdir.Dot(dir);

    return (muReco-mu)/mu;
  });





  // True visible energy fraction of the muon
  const Var kTrueEFrac([](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu < 1)
      return -5.0f;

    return float(sr->mc.nu[0].visE/sr->mc.nu[0].E);
  });

  // 2D var : energy resolution vs true energy
  const Var kResVsTrueE = Var2D(
    kResE,
    resbins,
    kTrueE,
    ebins);

  // 2D var : energy resolution vs reco energy
  const Var kResVsRecoE = Var2D(
    kResE,
    resbins,
    kRecoE,
    ebins);

  // 2D var : muon momentum resolution vs true muon momentum
  const Var kResVsTrueMuP = Var2D(
    kResMuP,
    resbins,
    kTrueMuP,
    ebins);

  // 2D var : muon momentum resolution vs reco muon momentum
  const Var kResVsRecoMuP = Var2D(
    kResMuP,
    resbins,
    kRecoMuP,
    ebins);

  // 2D var : muon transverse momentum resolution vs true muon transverse momentum
  const Var kResVsTrueMuPt = Var2D(
    kResMuPt,
    resbins,
    kTrueMuPt,
    ebins);

  // 2D var : muon transverse momentum resolution vs reco muon transverse momentum
  const Var kResVsRecoMuPt = Var2D(
    kResMuPt,
    resbins,
    kRecoMuPt,
    ebins);

  // 2D var : muon longitudinal momentum resolution vs true muon longitudinal momentum
  const Var kResVsTrueMuPz = Var2D(
    kResMuPz,
    resbins,
    kTrueMuPz,
    ebins);

  // 2D var : muon longitudinal momentum resolution vs reco muon longitudinal momentum
  const Var kResVsRecoMuPz = Var2D(
    kResMuPz,
    resbins,
    kRecoMuPz,
    ebins);


  const Var kResVsTrueCostheta = Var2D(
    kResMuCostheta,
    resbins,
    kTrueMuCostheta,
    angbins);

  const Var kResVsRecoCostheta = Var2D(
    kResMuCostheta,
    resbins,
    kRecoMuCostheta,
    angbins);

  const Var kRecoPVsCos = Var2D(
    kRecoMuP,
    ebins,
    kRecoMuCostheta,
    angbins);

  const Var kRecoPtVsPz = Var2D(
    kRecoMuPz,
    ebins,
    kRecoMuPt,
    ebins);

  const Var kTruePVsCos = Var2D(
    kTrueMuP,
    ebins,
    kTrueMuCostheta,
    angbins);

  const Var kTruePtVsPz = Var2D(
    kTrueMuPz,
    ebins,
    kTrueMuPt,
    ebins);

  // background weights
  const Var kUncontainedBGWeight([](const caf::SRProxy* sr)
  {
    if(sr->trk.kalman.ntracks < 1)
      return -1000.0;

    return -0.00032*sr->trk.kalman.tracks[0].start.Y() + 1.14;
  });


  // Var For Energy estimation
  const Var kTrueCatcherE = SIMPLEVAR(energy.numu.mc.truemuoncatcherE);

  const Var kReMIdTrkLenAct ([](const caf::SRProxy* sr)
  {
    return (sr->energy.numu.ndtrklenact / 100); // in m
  });

  const Var kReMIdTrkLenCat ([](const caf::SRProxy* sr)
  {
    return (sr->energy.numu.ndtrklencat / 100); // in m
  });

}// end of namespace