NumuVars.cxx

Go to the documentation of this file.

#include <vector>
#include <iostream>
#include "TH1.h"
#include "TFile.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuVars.h"

#include "CAFAna/Core/Utilities.h"

#include "Utilities/func/MathUtil.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <cassert>

#include "TVector3.h"

namespace
{
  const TVector3 beamDirFD = ana::NuMIBeamDirection(caf::kFARDET);
  const TVector3 beamDirND = ana::NuMIBeamDirection(caf::kNEARDET);
}

namespace ana
{
  const Var kNonQEE = SIMPLEVAR(energy.numu.trknonqeE);

  const Var kQEE = SIMPLEVAR(energy.numu.trkqeE);

  const Var kSliceTime(
           [](const caf::SRProxy *sr)
           {
             return sr->slc.meantime/1000;
           });

  const Var kSliceDuration(
         [](const caf::SRProxy *sr)
         {
           return (sr->slc.endtime - sr->slc.starttime);
         });

  const Var kDirX([](const caf::SRProxy *sr)
                  {
                    if(sr->trk.kalman.ntracks < 1) return -5.f;
                    return sr->trk.kalman.tracks[0].dir.X();
                  });

  const Var kDirY([](const caf::SRProxy *sr)
                  {
                    if(sr->trk.kalman.ntracks < 1) return -5.f;
                    return sr->trk.kalman.tracks[0].dir.Y();
                  });

  const Var kcosDirY([](const caf::SRProxy *sr)
                  {
                    if(sr->trk.kalman.ntracks < 1) return -5.f;
                    return (float)cos( kDirY(sr) );
                  });

  const Var kDirZ([](const caf::SRProxy *sr)
                  {
                    if(sr->trk.kalman.ntracks < 1) return -5.f;
                    return sr->trk.kalman.tracks[0].dir.Z();
                  });

  const Var kCosNumi(
                    [](const caf::SRProxy *sr)
                    {
                      if (sr->trk.kalman.ntracks > 0 && sr->trk.kalman.idxremid != 999){
                        if (sr->hdr.det == 1){
                          return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);
                        }
                        if (sr->hdr.det == 2){
                          return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);
                        }
                      }
                      return -5.f;
                    });

  const Var kRecoThetaMu( 
                        []( const caf::SRProxy* sr )
                        {
                          return acos( kCosNumi( sr ) ) * 180.0 / 3.1416;
                        });

  const Var kNumuMuonPtP([](const caf::SRProxy *sr) {
      if (sr->trk.kalman.ntracks > 0 && sr->trk.kalman.idxremid != 999){
        if (sr->hdr.det == 1){
          double Zbeam = sr->trk.kalman.tracks[0].dir.Dot(beamDirND);
          double ptp = sqrt(1 - Zbeam*Zbeam);
          return (float)ptp;
        }
        if (sr->hdr.det == 2){
          double Zbeam =  sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);
          double ptp = sqrt(1 - Zbeam*Zbeam);
          return (float)ptp;
        }
      }
      return -5.f;
    });

  const Var kNumuMuonPt( 
            []( const caf::SRProxy* sr )
            {
              if ( sr->trk.kalman.ntracks > 0 && sr->trk.kalman.idxremid != 999 )
              {
                // Important: Use kMuE to pick up muon energy scale systematics
                double pmu = sqrt( util::sqr( kMuE( sr ) ) - util::sqr( 0.1056 ) );
                double Zbeam = kCosNumi( sr );
                double pt = pmu * sqrt( 1 - Zbeam * Zbeam );
                return (float)pt;
              }
              return -5.f;
            });

  const Var kTrkStartX(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1)
               return -10.0f;
             return sr->trk.kalman.tracks[0].start.X()/100;
           });
  const Var kTrkStartY(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1)
               return -10.0f;
             return sr->trk.kalman.tracks[0].start.Y()/100;
           });
  const Var kTrkStartZ(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1)
               return -10.0f;
             return sr->trk.kalman.tracks[0].start.Z()/100;
           });

  const Var kTrkEndX(
         [](const caf::SRProxy *sr)
         {
           if(sr->trk.kalman.ntracks < 1)
             return -10.0f;
           return sr->trk.kalman.tracks[0].stop.X()/100;
         });
  const Var kTrkEndY(
         [](const caf::SRProxy *sr)
         {
           if(sr->trk.kalman.ntracks < 1)
             return -10.0f;
           return sr->trk.kalman.tracks[0].stop.Y()/100;
         });
  const Var kTrkEndZ(
         [](const caf::SRProxy *sr)
         {
           if(sr->trk.kalman.ntracks < 1)
             return -10.0f;
           return sr->trk.kalman.tracks[0].stop.Z()/100;
         });

  const Var kTrkNhits(
          [](const caf::SRProxy *sr)
          {
            if(sr->trk.kalman.ntracks < 1)
              return 65535;
            return int(sr->trk.kalman.tracks[0].nhit);
          });

  const Var kHadNHit(
          [](const caf::SRProxy *sr)
          {
            unsigned int nought = 0;
            if(sr->trk.kalman.ntracks < 1) return nought;
            return sr->slc.nhit - sr->trk.kalman.tracks[0].nhit;

          });

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

            if(sr->trk.kalman.ntracks < 1) return  0.0f;
            int nHit= sr->slc.nhit - sr->trk.kalman.tracks[0].nhit;
            if(nHit<=0) return 0.0f;
            float hadE = sr->energy.numu.hadcalE;
            return hadE/nHit;

          });

  const Var kTrkLength(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return  -5.f;
             return sr->trk.kalman.tracks[0].len / 100;
           });

  const Var kTrkCalE(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return 0.f;
             return float(sr->trk.kalman.tracks[0].calE);
           });

  const Var kTrkCalEPerNHit(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return 0.f;
             return sr->trk.kalman.tracks[0].calE /
                          sr->trk.kalman.tracks[0].nhit;
           });

  const Var kMuEPerNHit(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return 0.f;
             return (float)kMuE(sr) / sr->trk.kalman.tracks[0].nhit;
           });

  const Var kCosBackDist(
            [](const caf::SRProxy *sr)
            {
              return sr->sel.contain.cosbakdist / 100.; // Nicer in m
            });
  const Var kCosFwdDist(
            [](const caf::SRProxy *sr)
            {
              return sr->sel.contain.cosfwddist / 100.; // Nicer in m
            });

  const Var kKalmanFwdDist(
            [](const caf::SRProxy *sr)
            {
              return sr->sel.contain.kalfwddist / 100.; // Nicer in m
            });
  const Var kKalmanBackDist(
            [](const caf::SRProxy *sr)
            {
              return sr->sel.contain.kalbakdist / 100.; // Nicer in m
            });

  const Var kNonQeHadE(
            [](const caf::SRProxy *sr)
            {
              if (sr->trk.kalman.ntracks > 0 && sr->trk.kalman.idxremid != 999)
                return float(sr->energy.numu.recotrkcchadE);
              return -5.f;
            });

  const Var kQeHadE(
            [](const caf::SRProxy *sr)
            {
              if (sr->trk.kalman.ntracks > 0 && sr->trk.kalman.idxremid != 999)
                return sr->energy.numu.trkqeE - sr->energy.numu.recomuonE;
              return -5.f;
            });

  const Var kSlcUnCalibNHit(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.nhit - sr->slc.ncalhit;
            });

  const Var kSlcMeanTime(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.meantime / 1000.; // Nicer in \mu s
            });

  const Var kSlcStartTime(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.starttime / 1000.; // Nicer in \mu s
            });

  const Var kSlcEndTime(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.endtime / 1000.; // Nicer in \mu s
            });

  const Var kSlcMinX(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.boxmin.X() / 100.; // Nicer in m
            });

  const Var kSlcMinY(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.boxmin.Y() / 100.; // Nicer in m
            });

  const Var kSlcMinZ(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.boxmin.Z() / 100.; // Nicer in m
            });

  const Var kSlcMaxX(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.boxmax.X() / 100.; // Nicer in m
            });

  const Var kSlcMaxY(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.boxmax.Y() / 100.; // Nicer in m
            });

  const Var kSlcMaxZ(
            [](const caf::SRProxy *sr)
            {
              return sr->slc.boxmax.Z() / 100.; // Nicer in m
            });

  const Var kSlcExtentX(
            [](const caf::SRProxy *sr)
            {
              return (sr->slc.boxmax.X() - sr->slc.boxmin.X())/ 100.; // Nicer in m
            });

  const Var kSlcExtentY(
            [](const caf::SRProxy *sr)
            {
              return (sr->slc.boxmax.Y() - sr->slc.boxmin.Y())/ 100.; // Nicer in m
            });

  const Var kSlcExtentZ(
            [](const caf::SRProxy *sr)
            {
              return (sr->slc.boxmax.Z() - sr->slc.boxmin.Z())/ 100.; // Nicer in m
            });

  const Var kTrkNPlaneGap(
            [](const caf::SRProxy *sr)
            {
              if (sr->trk.kalman.ntracks > 0 && sr->trk.kalman.idxremid != 999)
                return int(sr->trk.kalman.tracks[0].nplanegap);
              return 500;
            });

  const Var kSlcCalEPerNHit(
            [](const caf::SRProxy *sr)
            {
              if (sr->slc.nhit > 0) return sr->slc.calE / (1.78 * sr->slc.nhit);
              return -5.;
            });

  //  const Var kNumuHadVisE(
  //            [](const caf::SRProxy *sr)
  //            {
  //             return sr->energy.numu.hadcalE + sr->energy.numu.hadtrkE;
  //            });

  //**********Variables useful for making muon energy resolution plots

  const Var kVisTrkCalE(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return 0.;
             if(sr->energy.numu.trkccE <= 0) return 0.;
             if(sr->energy.numu.hadtrkE < 0) return 0.;
             return ((sr->trk.kalman.tracks[0].calE)/1.78)-sr->energy.numu.hadtrkE;
           });

  const Var kTrueMuonE(
           [](const caf::SRProxy *sr)
           {
             if(sr->mc.nnu == 0) return 0.f;
             if(sr->mc.nu[0].prim.empty()) return 0.f;
             if(std::abs(sr->mc.nu[0].prim[0].pdg) != 13) return 0.f;
             return float(sr->mc.nu[0].prim[0].p.E);
           });

  const Var kRecoMuonEFromTrkLen(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return 0.f;
             if(sr->energy.numu.trkccE <= 0) return 0.f;
             if(sr->energy.numu.hadtrkE < 0) return 0.f;
             return float(sr->energy.numu.recomuonE);
           });

  const Var kRecoMuonEFromVisEFD(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return 0.;
             if(sr->energy.numu.trkccE <= 0) return 0.;
             if(sr->energy.numu.hadtrkE < 0) return 0.;
             double visE = ((sr->trk.kalman.tracks[0].calE)/1.78)-sr->energy.numu.hadtrkE;
             double offset = 0.32;
             double slope1 = 1.01;
             double slope2 = 1.504;
             double stitch = 0.315;
             if (visE<stitch){
               return(visE*slope1 + offset);
             }
             else{
               return(visE*slope2 + ((slope1-slope2)*stitch + offset));
             }
           });

  const Var kRecoMuonEFromVisEND(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return 0.;
             if(sr->energy.numu.trkccE <= 0) return 0.;
             if(sr->energy.numu.hadtrkE < 0) return 0.;
             double visE = ((sr->trk.kalman.tracks[0].calE)/1.78)-sr->energy.numu.hadtrkE;
             double offset = 0.222;
             double slope1 = 1.14;
             double slope2 = 1.475;
             double stitch = 0.312;
             if (visE<stitch){
               return(visE*slope1 + offset);
             }
             else{
               return(visE*slope2 + ((slope1-slope2)*stitch + offset));
             }
           });

  const Var kMuonRecoRelTrueVisEFD(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return -1.;
             if(sr->energy.numu.trkccE <= 0) return -1.;
             if(sr->energy.numu.hadtrkE < 0) return -1.;
             if(sr->mc.nnu == 0) return -1.;
             if(sr->mc.nu[0].prim.empty()) return -1.;
             if(std::abs(sr->mc.nu[0].prim[0].pdg) != 13) return -1.;
             double visE = ((sr->trk.kalman.tracks[0].calE)/1.78)-sr->energy.numu.hadtrkE;
             double offset = 0.32;
             double slope1 = 1.01;
             double slope2 = 1.504;
             double stitch = 0.315;
             double recoE = -5.0;
             if (visE<stitch){
               recoE = visE*slope1 + offset;
             }
             else{
               recoE = visE*slope2 + ((slope1-slope2)*stitch + offset);
             }
             double trueE = sr->mc.nu[0].prim[0].p.E;
             if (trueE <= 0.0) return -1.;
             return((recoE - trueE)/trueE);
           });

  const Var kMuonRecoRelTrueVisEND(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return -1.;
             if(sr->energy.numu.trkccE <= 0) return -1.;
             if(sr->energy.numu.hadtrkE < 0) return -1.;
             if(sr->mc.nnu == 0) return -1.;
             if(sr->mc.nu[0].prim.empty()) return -1.;
             if(std::abs(sr->mc.nu[0].prim[0].pdg) != 13) return -1.;
             double visE = ((sr->trk.kalman.tracks[0].calE)/1.78)-sr->energy.numu.hadtrkE;
             double offset = 0.222;
             double slope1 = 1.14;
             double slope2 = 1.475;
             double stitch = 0.312;
             double recoE = -5.0;
             if (visE<stitch){
               recoE = visE*slope1 + offset;
             }
             else{
               recoE = visE*slope2 + ((slope1-slope2)*stitch + offset);
             }
             double trueE = sr->mc.nu[0].prim[0].p.E;
             if (trueE <= 0.0) return -1.;
             return((recoE - trueE)/trueE);
           });

  const Var kMuonRecoRelTrueTrkLen(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return -1.;
             if(sr->energy.numu.trkccE <= 0) return -1.;
             if(sr->energy.numu.hadtrkE < 0) return -1.;
             if(sr->mc.nnu == 0) return -1.;
             if(sr->mc.nu[0].prim.empty()) return -1.;
             if(std::abs(sr->mc.nu[0].prim[0].pdg) != 13) return -1.;
             double recoE = sr->energy.numu.recomuonE;
             double trueE = sr->mc.nu[0].prim[0].p.E;
             if (trueE <= 0.0) return -1.;
             return((recoE - trueE)/trueE);
           });

  const Var kMuonRecoRelRecoFD(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return -1.;
             if(sr->energy.numu.trkccE <= 0) return -1.;
             if(sr->energy.numu.hadtrkE < 0) return -1.;
             if(sr->mc.nnu == 0) return -1.;
             if(sr->mc.nu[0].prim.empty()) return -1.;
             if(std::abs(sr->mc.nu[0].prim[0].pdg) != 13) return -1.;
             double recoETrkLen = sr->energy.numu.recomuonE;
             double visE = ((sr->trk.kalman.tracks[0].calE)/1.78)-sr->energy.numu.hadtrkE;
             double offset = 0.32;
             double slope1 = 1.01;
             double slope2 = 1.504;
             double stitch = 0.315;
             double recoE = -5.0;
             if (visE<stitch){
               recoE = visE*slope1 + offset;
             }
             else{
               recoE = visE*slope2 + ((slope1-slope2)*stitch + offset);
             }
             if (recoETrkLen <= 0.0) return -1.;
             return((recoE - recoETrkLen)/recoETrkLen);
           });

  const Var kMuonRecoRelRecoND(
           [](const caf::SRProxy *sr)
           {
             if(sr->trk.kalman.ntracks < 1) return -1.;
             if(sr->energy.numu.trkccE <= 0) return -1.;
             if(sr->energy.numu.hadtrkE < 0) return -1.;
             double recoETrkLen = sr->energy.numu.recomuonE;
             double visE = ((sr->trk.kalman.tracks[0].calE)/1.78)-sr->energy.numu.hadtrkE;
             double offset = 0.222;
             double slope1 = 1.14;
             double slope2 = 1.475;
             double stitch = 0.312;
             double recoE = -5.0;
             if (visE<stitch){
               recoE = visE*slope1 + offset;
             }
             else{
               recoE = visE*slope2 + ((slope1-slope2)*stitch + offset);
             }
             if (recoETrkLen <= 0.0) return -1.;
             return((recoE - recoETrkLen)/recoETrkLen);
           });

  const Var kNumuHadCalE        = SIMPLEVAR(energy.numu.hadcalE);
  const Var kNumuHadTrkE        = SIMPLEVAR(energy.numu.hadtrkE);
  const Var kNKalman            = SIMPLEVAR(trk.kalman.ntracks);
  const Var kNumuHadVisE(
            [](const caf::SRProxy *sr)
            {
              return kNumuHadCalE(sr) + kNumuHadTrkE(sr);
            });

  const Var kNumuCosRejAngleKal = SIMPLEVAR(sel.cosrej.anglekal);
  const Var kNumuCosRejAngleCos = SIMPLEVAR(sel.cosrej.anglecos);
  const Var kCosmicForwardCell  = SIMPLEVAR(sel.contain.cosfwdcell);
  const Var kCosmicBackwardCell = SIMPLEVAR(sel.contain.cosbakcell);
  const Var kKalmanForwardCell  = SIMPLEVAR(sel.contain.kalfwdcell);
  const Var kKalmanBackwardCell = SIMPLEVAR(sel.contain.kalbakcell);
  const Var kNumuContPID        = SIMPLEVAR(sel.cosrej.numucontpid2020);
  const Var kNumuContPID2019    = SIMPLEVAR(sel.cosrej.numucontpid2019);
  const Var kReMIdScatLLH       = SIMPLEVAR(sel.remid.scatllh);
  const Var kReMIdDEDxLLH       = SIMPLEVAR(sel.remid.dedxllh);
  const Var kReMIdMeasFrac      = SIMPLEVAR(sel.remid.measfrac);
  const Var kNplanesToFront     = SIMPLEVAR(sel.contain.nplanestofront);
  const Var kNplanesToBack      = SIMPLEVAR(sel.contain.nplanestoback);
  const Var kNcellsFromEdge     = SIMPLEVAR(slc.ncellsfromedge);
  const Var kNonQeE             = SIMPLEVAR(energy.numu.trknonqeE);
  const Var kQeE                = SIMPLEVAR(energy.numu.trkqeE);
  const Var kQeAngleE           = SIMPLEVAR(energy.numu.angleE);
  const Var kDirScore           = SIMPLEVAR(sel.cosrej.kdirscore);

  /// Reconstructed four-momentum transfer invariant (Q^2)
  const ana::Var kRecoQ2
                 ([] (const caf::SRProxy * sr)
                  {
                    const double M_mu_sqrd = util::sqr(0.1056);
                    double E_mu = kMuE(sr);
                    double p_mu = sqrt(util::sqr(E_mu) - M_mu_sqrd);
                    return 2 * kCCE(sr) * ( E_mu - p_mu * kCosNumi(sr) ) - M_mu_sqrd;
                  }
  );

  /// Reconstructed three-momentum transfer
  const Var kRecoQmag
                 ([] (const caf::SRProxy * sr)
                  {
                    // because q0 is just the energy transferred to the hadronic system
                    return sqrt(kRecoQ2(sr) + util::sqr(kHadE(sr)));
                  }
  );

  /// Reconstructed invariant mass (W)
  const Var kRecoW(
                 [](const caf::SRProxy *sr)
                 {
                   const double M_p = 0.938272;
                   double WSq = M_p * M_p + 2 * M_p * (kCCE(sr) - kMuE(sr)) - kRecoQ2(sr);
                   if (WSq > 0)
                     return sqrt(WSq);
                   else
                     return -5.;
                 });
}