NusVars.cxx

Go to the documentation of this file.

#include "3FlavorAna/Vars/NueEnergy2020.h"
#include "NuXAna/Vars/NusVars.h"

#include "StandardRecord/Proxy/SRProxy.h"


namespace ana
{

  // Nus20 base variables for selection
  const Var kIsElastic       = SIMPLEVAR(vtx.elastic.IsValid);
  const Var kNFuzzyProng     = SIMPLEVAR(vtx.elastic.fuzzyk.npng);
  const Var kNContPlanes     = SIMPLEVAR(slc.ncontplanes);
  const Var kCVNnc_looseptp  = SIMPLEVAR(sel.cvnloosepreselptp.ncid);
  const Var kCVNnc_oldpresel = SIMPLEVAR(sel.cvnoldpresel.ncid);

  // Nus20 BDT variables
  const Var kNSliceHits      = SIMPLEVAR(slc.nhit);
  const Var kShwNHitX        = SIMPLEVAR(vtx.elastic.fuzzyk.png[0].shwlid.nhitx);
  const Var kShwNHitY        = SIMPLEVAR(vtx.elastic.fuzzyk.png[0].shwlid.nhity);
  const Var kShwDirY         = SIMPLEVAR(vtx.elastic.fuzzyk.png[0].shwlid.dir.y);
  const Var kNShwLid         = SIMPLEVAR(vtx.elastic.fuzzyk.nshwlid);
  const Var kNMipHit         = SIMPLEVAR(slc.nmiphit);
 
  const Var kShwNHitXPlusY([](const caf::SRProxy* sr)
    {
      double nHits = sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx +
                     sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;

      return nHits;
    });
 
  const Var kShwNHitXMinY([](const caf::SRProxy* sr)
    {
      double nHits = sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx -
                     sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;

      return nHits;
    });
  
  const Var kShwNHitXOverY([](const caf::SRProxy* sr)
    {
      double hitRatio = (double)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx /
                        (double)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;

      return hitRatio;
    });

  // Nus17 reconstructed (visible) energy - see docdb 17797 and 17985
  const double FDscaleCalE17 = 0.6377;
  const double NDscaleCalE17 = 0.6484;

  const Var kNus17Energy([](const caf::SRProxy* sr)
    {
      double cale  = sr->slc.calE;
      double recoE = 0.;
      if(sr->hdr.det == caf::kFARDET)  recoE = FDscaleCalE17*cale;
      if(sr->hdr.det == caf::kNEARDET) recoE = NDscaleCalE17*cale;
      return recoE;
    });


  const double NDscaleCalE18RHC = 1.15;
  const double NDscaleCalE18 = 1.11;
  const double FDscaleCalE18RHC = 1.18;
  const double FDscaleCalE18 = 1.20;

  const Var kNus18Energy([](const caf::SRProxy* sr) {
      bool h_FHC = sr->spill.isFHC;
      bool h_RHC = sr->spill.isRHC;
      double cale  = sr->slc.calE;
      double recoE = 0.;
      if (h_FHC && sr->hdr.det == caf::kFARDET) recoE = FDscaleCalE18*cale;
      else if (h_FHC && sr->hdr.det == caf::kNEARDET) recoE = NDscaleCalE18*cale;
      else if (h_RHC && sr->hdr.det == caf::kFARDET) recoE = FDscaleCalE18RHC*cale;
      else if (h_RHC && sr->hdr.det == caf::kNEARDET) recoE = NDscaleCalE18RHC*cale;
      return recoE ; // default case
    });


  const Var kNus20Energy([](const caf::SRProxy* sr) {
      if (sr->hdr.det != caf::kNEARDET && sr->hdr.det != caf::kFARDET) return (double) sr->slc.calE; // Only defined for ND & FD (eg, not test beam)
      double pars[4][6] = { {1.049, 0.795, 0.8409, 0.17,  0.82, -1.00},  // ND FHC (or 0HC)
                            {1.025, 0.797, 0.9162, 0.53, -0.26, -1.48},  // FD FHC (or 0HC)
                            {1.000, 1.000, 1.0000, 0.00,  0.00,  0.00},  // ND RHC (just return calE = e_EM + e_Had)
                            {1.000, 1.000, 1.0000, 0.00,  0.00,  0.00} };// FD RHC (just return calE = e_EM + e_Had)
      int detCur = (sr->hdr.det == caf::kFARDET) + ((sr->spill.isRHC)<<1); // bitshifting is magic
      double e_EM  = ana:: kEME_2020(sr);
      double e_Had = ana::kHADE_2020(sr);
      double e_Cal = sr->slc.calE;
      return (e_EM / pars[detCur][0] + e_Had / pars[detCur][1]) / (pars[detCur][2] + pars[detCur][3] * std::pow(e_Cal+pars[detCur][4],2) * std::exp(pars[detCur][5] * e_Cal));
    });

  /// NC Cosmic Rejection (G)BDT variables
  const Var kNCCosRej  = SIMPLEVAR(sel.nccosrej.cospidbdt);
  const Var kNCCosRejG = SIMPLEVAR(sel.nccosrej.cospidbdtg);

  const Var kPartPtp([](const caf::SRProxy* sr) {
      if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
         std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
      return float(sr->sel.nuecosrej.partptp);
    });

  const Var kDistTop([](const caf::SRProxy* sr) {
      return std::min(sr->sel.nuecosrej.starttop,
                      sr->sel.nuecosrej.stoptop);
    });

  const Var kMinDistFace([](const caf::SRProxy* sr) {
      const caf::SRNueCosRejProxy& cr = sr->sel.nuecosrej;
      float mindist = 10000.;
      mindist = std::min(mindist, std::min(cr.starteast,   cr.stopeast));
      mindist = std::min(mindist, std::min(cr.startwest,   cr.stopwest));
      mindist = std::min(mindist, std::min(cr.starttop,    cr.stoptop));
      mindist = std::min(mindist, std::min(cr.startbottom, cr.stopbottom));
      mindist = std::min(mindist, std::min(cr.startfront,  cr.stopfront));
      mindist = std::min(mindist, std::min(cr.startback,   cr.stopback));
      return (mindist > 9999. ? (float)-5. : mindist);
    });

  const Var kClosestSlcTime([](const caf::SRProxy* sr) {
      return (float)sr->slc.closestslicetime;
    });

  const Var kClosestSlcMinDist([](const caf::SRProxy* sr) {
      return (float)sr->slc.closestslicemindist;
    });

  const Var kClosestSlcMinTop([](const caf::SRProxy* sr) {
      return (float)sr->slc.closestsliceminfromtop;
    });

  //---------------------------------------------------------------------------------
  // Nus18 kaon
  const NuTruthVar kNus18BeamWeightP_NT([](const caf::SRNeutrinoProxy* nu) {
      if (abs(nu->beam.ptype) == 321 || nu->beam.ptype == 310 || nu->beam.ptype == 311 || nu->beam.ptype == 130)
        return 1.3;
      return 1.;
    });

  const NuTruthVar kNus18BeamWeightM_NT([](const caf::SRNeutrinoProxy* nu) {
      if (abs(nu->beam.ptype) == 321 || nu->beam.ptype == 310 || nu->beam.ptype == 311 || nu->beam.ptype == 130)
        return 0.7;
      return 1.;
    });

  //---------------------------------------------------------------------------------
  // Nus18 tau
  const NuTruthVar kFDTauWgtM_NT([](const caf::SRNeutrinoProxy* nu) {
      if (abs(nu->pdg) == 16 && nu->iscc)
        return 0.4;
      return 1.;
  });

  const NuTruthVar kFDTauWgtP_NT([](const caf::SRNeutrinoProxy* nu) {
      if (abs(nu->pdg) == 16 && nu->iscc)
        return 1.6;
      return 1.;
  });
}