NumuCCIncVars.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuEFxs.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Utilities.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/TruthVars.h"

#include "Utilities/func/MathUtil.h"

#include "NDAna/numucc_inc/NumuCCIncBins.h"
#include "NDAna/numucc_inc/NDXSecMuonPID.h"

#include "TH2.h"
#include "TMath.h"
#include "TH2.h"
#include "TH1.h"
#include "TAxis.h"
#include "TVector3.h"
#include "TParticlePDG.h"
#include "TDatabasePDG.h"
#include "TLorentzVector.h"

namespace ana
{
  //-------------------------------------------------------------------------------------------------
  // 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);

  extern const Var kMostMuonLikeTrackIdx;

  extern const Var kMuonLLRPID;

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

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

  extern const Var kMuStartX;
  extern const Var kMuStartY;
  extern const Var kMuStartZ;

  extern const Var kMuStopX;
  extern const Var kMuStopY;
  extern const Var kMuStopZ;

  extern const Var kMuonLength;

  extern const Var kVisHadE;

  extern const Var kHadEonTrack;

  extern const Var kAveragededx40;

//==============================================================================
// 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.

  extern const Var kTrkLenAct;

  extern const Var kTrkLenCat;

  const Var kTrkLenActandCat = kTrkLenAct + kTrkLenCat;


  // Muon Energy for Inclusive numuCC
  //------------------------------------
  extern const Var kIncXsecMuonE;

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

  extern const Var kIncXsecHadE;

  // Neutrino Energy for Inclusive numuCC
  //---------------------------------------

  const Var kIncXsecNuE = kIncXsecMuonE + kIncXsecHadE;

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

  const Var kRecoE = kIncXsecNuE;

  inline double MuonMass()
  {
    static TDatabasePDG* pdgdb = TDatabasePDG::Instance();
    static double mass = pdgdb->GetParticle(13)->Mass();
    return mass;
  }

  const Var kMuonMass([](const caf::SRProxy* sr)
  {
    return MuonMass();
  });

  const Var kMuonMass2([](const caf::SRProxy* sr)
  {
    return util::sqr(MuonMass());
  });

  inline double NeutronMass()
  {
    static TDatabasePDG* pdgdb = TDatabasePDG::Instance();
    static double mass = pdgdb->GetParticle(2112)->Mass();
    return mass;
  }

  const Var kNeutronMass([](const caf::SRProxy* sr)
  {
    return NeutronMass();
  });

  const Var kNeutronMass2([](const caf::SRProxy* sr)
  {
    return util::sqr(NeutronMass());
  });

  const Var kRecoMuKE =  kIncXsecMuonE - kMuonMass;

  // Reconstructed candidate muon angle wrt beam
  extern const Var kRecoMuCostheta;

  // Angle systematics (+/- 0.0025 mrad wrt beam)
  extern const Var kRecoMuTheta;

  const Var kThetaSystShift = Constant(0.0025);
  const Var kRecoMuThetaUp = kRecoMuTheta + kThetaSystShift;
  const Var kRecoMuThetaDw = kRecoMuTheta - kThetaSystShift;

  extern const Var kRecoMuCosthetaUp;
  extern const Var kRecoMuCosthetaDw;

  extern const NuTruthVar kTrueEST;
  const Var kTrueNuE = VarFromNuTruthVar(kTrueEST);

  extern const NuTruthVar kTrueMuEST;
  const Var kTrueMuE = VarFromNuTruthVar(kTrueMuEST);

  extern const NuTruthVar kTrueMuKEST;
  const Var kTrueMuKE = VarFromNuTruthVar(kTrueMuKEST);


  // True muon cos theta wrt beam direction

  extern const NuTruthVar kTrueMuCosthetaST;
  const Var kTrueMuCostheta = VarFromNuTruthVar(kTrueMuCosthetaST);

  // True hadronic y
  extern const NuTruthVar kTrueYST;

  // True hadronic energy
  extern const NuTruthVar kTrueHadEST;
  const Var kTrueHadE = VarFromNuTruthVar(kTrueHadEST);

  // Reconstructed hadronic energy
  const Var kRecoHadE = kIncXsecHadE;

  const NuTruthVar kTrueEavailST = kTrueEavail_NT;

  // Reconstructed available energy (Determined by Travis Olson)
  extern const Var kRecoEavail;


  extern const Var kRecoMuKEVsCos;

  extern const NuTruthVar kTrueMuKEVsCosST;
  const Var kTrueMuKEVsCos = VarFromNuTruthVar(kTrueMuKEVsCosST);


  extern const NuTruthVar kTrueMuKEVsCosVsEnuST;

  extern const Var kRecoMuKEVsCosVsEnu;

  extern const NuTruthVar kTrueMuKEVsCosVsEavailST;

  extern const Var kRecoMuKEVsCosVsEavail;

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

  // Vars related to muon energy scale systematics:
  extern const Var kIncXsecMuonEUp;

  extern const Var kIncXsecMuonEDw;

  // Muon Energy Shift Vars
  const Var kIncXsecNuEUp = kIncXsecMuonEUp + kIncXsecHadE;
  const Var kIncXsecNuEDw = kIncXsecMuonEDw + kIncXsecHadE;

  const Var kRecoEUp      = kIncXsecNuEUp;
  const Var kRecoEDw      = kIncXsecNuEDw;

  const Var kRecoMuKEUp   = kIncXsecMuonEUp - kMuonMass;
  const Var kRecoMuKEDw   = kIncXsecMuonEDw - kMuonMass;

  // 3D vars for muon energy shift (Enu as 3rd axis)
  extern const Var kRecoMuKEVsCosVsEnu_onlyMuKEUp;
  extern const Var kRecoMuKEVsCosVsEnu_onlyMuKEDw;
  extern const Var kRecoMuKEVsCosVsEnuUp;
  extern const Var kRecoMuKEVsCosVsEnuDw;

  // 3D vars for muon energy shift (Eavail as 3rd axis)
  extern const Var kRecoMuKEVsCosVsEavailUp;
  extern const Var kRecoMuKEVsCosVsEavailDw;

  // 3D vars for muon angle shift (Eavail as 3rd axis)
  extern const Var kRecoMuKEVsCosVsEavailAngleUp;
  extern const Var kRecoMuKEVsCosVsEavailAngleDw;

  // Hist axes for p, costheta and energy
  const HistAxis kRecoMuKEVsCosStandardAxis(
    "Reconstructed T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}",
    angvsmukebins,
    kRecoMuKEVsCos);
  const NuTruthHistAxis kTrueMuKEVsCosStandardAxisST(
    "True T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}",
    angvsmukebins,
    kTrueMuKEVsCosST);

  const HistAxis kTrueMuKEVsCosStandardAxis = HistAxisFromNuTruthHistAxis(kTrueMuKEVsCosStandardAxisST);

  const HistAxis kRecoEStandardAxis(
    "Reconstructed Neutrino Energy (GeV)",
    enubins,
    kRecoE);
  const NuTruthHistAxis kTrueEStandardAxisST(
    "True Neutrino Energy (GeV)",
    enubins,
    kTrueEST);

  const HistAxis kTrueEStandardAxis = HistAxisFromNuTruthHistAxis(kTrueEStandardAxisST);


  // 3D axis, for 3D efficiency calculation
  const NuTruthHistAxis kTrueMuKEVsCosVsEnuStandardAxisST(
    "True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
    angvsmukevsebins,
    kTrueMuKEVsCosVsEnuST);

  const HistAxis kTrueMuKEVsCosVsEnuStandardAxis = HistAxisFromNuTruthHistAxis(kTrueMuKEVsCosVsEnuStandardAxisST);


  const NuTruthHistAxis kTrueMuKEVsCosVsEavailStandardAxisST(
    "True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
    angvsmukevseavailbins,
    kTrueMuKEVsCosVsEavailST);

  const HistAxis kTrueMuKEVsCosVsEavailStandardAxis = HistAxisFromNuTruthHistAxis(kTrueMuKEVsCosVsEavailStandardAxisST);

  // 3D reco axis, for 3D efficiency calculation
  const HistAxis kRecoMuKEVsCosVsEnuStandardAxis(
    "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
    angvsmukevsebins,
    kRecoMuKEVsCosVsEnu);

  const HistAxis kRecoMuKEVsCosVsEavailStandardAxis(
    "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
    angvsmukevseavailbins,
    kRecoMuKEVsCosVsEavail);

  const HistAxis kRecoMuKEVsCosVsEnuStandardAxis_onlyMuKEUp(
    "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
    angvsmukevsebins,
    kRecoMuKEVsCosVsEnu_onlyMuKEUp);
  const HistAxis kRecoMuKEVsCosVsEnuStandardAxis_onlyMuKEDw(
    "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
    angvsmukevsebins,
    kRecoMuKEVsCosVsEnu_onlyMuKEDw);
  const HistAxis kRecoMuKEVsCosVsEnuStandardAxisUp(
    "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
    angvsmukevsebins,
    kRecoMuKEVsCosVsEnuUp);
  const HistAxis kRecoMuKEVsCosVsEnuStandardAxisDw(
    "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
    angvsmukevsebins,
    kRecoMuKEVsCosVsEnuDw);

  const HistAxis kRecoMuKEVsCosVsEavailStandardAxisUp(
    "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
    angvsmukevseavailbins,
    kRecoMuKEVsCosVsEavailUp);
  const HistAxis kRecoMuKEVsCosVsEavailStandardAxisDw(
    "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
    angvsmukevseavailbins,
    kRecoMuKEVsCosVsEavailDw);

  // Angle systematic HistAxes
  const HistAxis kRecoMuKEVsCosVsEavailStandardAxisAngleUp(
    "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
    angvsmukevseavailbins,
    kRecoMuKEVsCosVsEavailAngleUp);
  const HistAxis kRecoMuKEVsCosVsEavailStandardAxisAngleDw(
    "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
    angvsmukevseavailbins,
    kRecoMuKEVsCosVsEavailAngleDw);

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


  extern const Var kDummy;

  extern const Var ksigResdown;
  extern const Var ksigDISup;

  extern const Var kbkgDISdown;

  //Number of prongs
  extern const Var kNProngs;
  // True Q2

  extern const Var kQsqr;

  // True W2

  extern const Var kWsqr;

  const Var kW = Sqrt(kWsqr);
  const Var kQ = Sqrt(kQsqr);


  extern const Var kweightedten;

 // Primary muon momentum

  extern const Var kTrueMuP;


 // Primary muon transverse momentum wrt beam direction

  extern const Var kTrueMuPt;

 // Primary muon momentum along the beam direction

  extern const Var kTrueMuPz;

  // Reconstructed muon momentum
  extern const Var kRecoMuPsqr;
  extern const Var kRecoMuPsqrUp;
  extern const Var kRecoMuPsqrDw;
  const Var kRecoMuP   = Sqrt(kRecoMuPsqr);
  const Var kRecoMuPUp = Sqrt(kRecoMuPsqrUp);
  const Var kRecoMuPDw = Sqrt(kRecoMuPsqrDw);


  // Reconstructed muon tranvserse momentum wrt beam direction

  extern const Var kRecoMuPt;

  // Reconstructed muon momentum along beam direction

  extern const Var kRecoMuPz;

  // Fractional neutrino energy resolution
  const Var kResE = (kNumuCCXsecNuE - kTrueNuE)/kTrueNuE;

  // Fractional muon momentum resolution
  extern const Var kResMuP;

  // Fractional muon transverse momentum resolution
  extern const Var kResMuPt;

  // Fractional resolution of muon momentum along beam direction
  extern const Var kResMuPz;

  // Remid
  extern const Var kRemid;

  // Remid variables
  extern const Var kScatLL;
  extern const Var kDedxLL;
  extern const Var kMeasfrac;

  // Length of candidate muon track
  extern const Var kMuLength;

  // Energy per length of the candidate muon track
  extern const Var kEPerLength;

  // Resolution of candidate muon angle wrt beam
  extern const Var kResMuCostheta;

  // True visible energy fraction of the muon
  extern const Var kTrueEFrac;

  // 2D var : energy resolution vs true energy
  extern const Var kResVsTrueE;

  // 2D var : energy resolution vs reco energy
  extern const Var kResVsRecoE;

  // 2D var : muon momentum resolution vs true muon momentum
  extern const Var kResVsTrueMuP;

  // 2D var : muon momentum resolution vs reco muon momentum
  extern const Var kResVsRecoMuP;

  // 2D var : muon transverse momentum resolution vs true muon transverse momentum
  extern const Var kResVsTrueMuPt;

  // 2D var : muon transverse momentum resolution vs reco muon transverse momentum
  extern const Var kResVsRecoMuPt;

  // 2D var : muon longitudinal momentum resolution vs true muon longitudinal momentum
  extern const Var kResVsTrueMuPz;

  // 2D var : muon longitudinal momentum resolution vs reco muon longitudinal momentum
  extern const Var kResVsRecoMuPz;

  extern const Var kResVsTrueCostheta;

  extern const Var kResVsRecoCostheta;

  extern const Var kRecoPVsCos;

  extern const Var kRecoPtVsPz;

  extern const Var kTruePVsCos;

  extern const Var kTruePtVsPz;

  // background weights
  extern const Var kUncontainedBGWeight;

  //Reconstructed Muon energy
  const Var kRecoMuE   = kIncXsecMuonE;
  const Var kRecoMuEUp = kIncXsecMuonEUp;
  const Var kRecoMuEDw = kIncXsecMuonEDw;

  const Var kTwo = Constant(2);

  const Var kRecoq2 =         kTwo * kRecoE  *(kRecoMuE-  (kRecoMuP  *kRecoMuCostheta  ))-kMuonMass2;
  const Var kRecoq2_MuonEUp = kTwo * kRecoEUp*(kRecoMuEUp-(kRecoMuPUp*kRecoMuCostheta  ))-kMuonMass2;
  const Var kRecoq2_MuonEDw = kTwo * kRecoEDw*(kRecoMuEDw-(kRecoMuPDw*kRecoMuCostheta  ))-kMuonMass2;
  const Var kRecoq2_AngleUp = kTwo * kRecoE  *(kRecoMuE-  (kRecoMuP  *kRecoMuCosthetaUp))-kMuonMass2;
  const Var kRecoq2_AngleDw = kTwo * kRecoE  *(kRecoMuE-  (kRecoMuP  *kRecoMuCosthetaDw))-kMuonMass2;
  const Var kRecow2 =  kNeutronMass2+(kTwo*kNeutronMass*(kRecoE-kRecoMuE))-kRecoq2;
  const Var kRecoq  =  Sqrt(kRecoq2);
  const Var kRecow  =  Sqrt(kRecow2);


  // Var For Energy estimation
  extern const Var kTrueCatcherE;

  extern const Var kReMIdTrkLenAct;
  extern const Var kReMIdTrkLenCat;


  const Var kDeltatheta     = kRecoMuCostheta - kTrueMuCostheta;
  const Var kDeltamuKE      = kRecoMuKE - kTrueMuKE;

  // Q2 Variables
  const NuTruthHistAxis kTrueQ2StandardAxisST(
    "True Q2 (GeV)",
    q2bins,
    kTrueQ2_NT);

  const HistAxis kTrueQ2StandardAxis = HistAxisFromNuTruthHistAxis(kTrueQ2StandardAxisST);

  const HistAxis kRecoQ2StandardAxis(
    "Reco Q2 (GeV)",
    q2bins,
    kRecoq2);

  // Systematically shifted 1D axes
  const HistAxis kRecoEStandardAxis_MuonEUp(
    "Reconstructed Neutrino Energy (GeV)",
    enubins,
    kRecoEUp);

    // Systematically shifted 1D axes
  const HistAxis kRecoEStandardAxis_MuonEDw(
    "Reconstructed Neutrino Energy (GeV)",
    enubins,
    kRecoEDw);

  const HistAxis kRecoQ2StandardAxis_MuonEUp(
    "Reco Q2 (GeV)",
    q2bins,
    kRecoq2_MuonEUp);

  const HistAxis kRecoQ2StandardAxis_MuonEDw(
    "Reco Q2 (GeV)",
    q2bins,
    kRecoq2_MuonEDw);

  const HistAxis kRecoQ2StandardAxis_AngleUp(
    "Reco Q2 (GeV)",
    q2bins,
    kRecoq2_AngleUp);

  const HistAxis kRecoQ2StandardAxis_AngleDw(
    "Reco Q2 (GeV)",
    q2bins,
    kRecoq2_AngleDw);

}//end of namespace