18 if(!hprob_dedx_vs_scat_sig || !hprob_dedx_vs_scat_bg){
19 const char*
path =
getenv(
"SRT_PUBLIC_CONTEXT");
22 TFile *muon_prob_file =
new TFile( fname.c_str(),
"r");
24 hprob_dedx_vs_scat_sig = (TH2F*)muon_prob_file->Get(
"prob_sig");
25 hprob_dedx_vs_scat_bg = (TH2F*)muon_prob_file->Get(
"prob_bg");
29 const int dedx_vs_scat_nxbins = hprob_dedx_vs_scat_bg->GetXaxis()->GetNbins();
30 const int dedx_vs_scat_nybins = hprob_dedx_vs_scat_bg->GetYaxis()->GetNbins();
32 int xbin_sig = hprob_dedx_vs_scat_sig->GetXaxis()->FindBin(dedxll);
33 int ybin_sig = hprob_dedx_vs_scat_sig->GetYaxis()->FindBin(scatll);
35 int xbin_bg = hprob_dedx_vs_scat_bg->GetXaxis()->FindBin(dedxll);
36 int ybin_bg = hprob_dedx_vs_scat_bg->GetYaxis()->FindBin(scatll);
41 if (xbin_sig == 0 || xbin_sig > dedx_vs_scat_nxbins ||
42 ybin_sig == 0 || ybin_sig > dedx_vs_scat_nybins ||
43 xbin_bg == 0 || xbin_bg > dedx_vs_scat_nxbins ||
44 ybin_bg == 0 || ybin_bg > dedx_vs_scat_nybins)
47 float psig =
std::log(hprob_dedx_vs_scat_sig->GetBinContent(xbin_sig, ybin_sig));
48 float pbg =
std::log(hprob_dedx_vs_scat_bg->GetBinContent(xbin_bg, ybin_bg));
61 double p0 = 1.67012e-01;
62 double p1 = 1.79305e-01;
63 double p2 = 3.74708e-03;
64 double p3 = -1.54232e-04;
68 if (trklenact <= 0.0)
return 0.0;
86 double p0 = 1.21130e-02;
87 double p1 = 1.97903e-01;
88 double p2 = 7.82459e-04;
92 if (trklenactandcat <= 0.0)
return 0.0;
95 + p1 * trklenactandcat
109 double offset = 1.31325e-01;
110 double slope = 5.35146e-01;
114 if (trklencat <= 0.0)
return 0.0;
116 MuonE = slope*trklencat +
offset;
127 double p0 = 5.85254e-02;
128 double p1 = 1.27796e+00;
129 double p2 = 3.75457e-01;
130 double p3 = -5.45618e-01;
131 double p4 = 1.65975e-01;
154 float maxllr = -1e12;
304 float vishadE = CalhadE + ExtraHadE;
418 float muonEact = 0.0;
419 float muonEcat = 0.0;
420 float muonEactandcat = 0.0;
422 float trkLenAct = 0.f;
423 float trkLenCat = 0.f;
442 muonE = muonEactandcat + muonEcat;
445 return muonE + muonEact;
469 float vishadE = CalhadE + ExtraHadE;
493 return (
float)dir.Dot(beamdir);
509 return (
float)dir.Angle(beamdir);
534 int nprims = nu->
prim.size();
535 for(
int iprim = 0; iprim < nprims; iprim++){
536 if(
abs(nu->
prim[iprim].pdg) == 13){
537 MuE = nu->
prim[iprim].p.T();
552 int nprims = nu->
prim.size();
553 for(
int iprim = 0; iprim < nprims; iprim++){
554 if(
abs(nu->
prim[iprim].pdg) == 13){
555 double E= nu->
prim[iprim].p.T();
570 int nprims = nu->
prim.size();
571 for(
int iprim = 0; iprim < nprims; iprim++){
572 if(
abs(nu->
prim[iprim].pdg) == 13){
573 TVector3 mudir = nu->
prim[iprim].p.Vect();
576 return mudir.Unit().Dot(beamdir.Unit());
655 float mue_act_up = 1.0079;
656 float mue_cat_up = 1.0120;
659 float muonEact = 0.0;
660 float muonEcat = 0.0;
661 float muonEactandcat = 0.0;
663 float trkLenAct = 0.f;
664 float trkLenCat = 0.f;
678 muonEact = mue_act_up *
MuonEAct(trkLenAct);
683 muonE = mue_act_up*muonEactandcat + mue_cat_up*muonEcat;
686 return muonE + muonEact;
692 float mue_act_dw = 0.9921;
693 float mue_cat_dw = 0.9880;
696 float muonEact = 0.0;
697 float muonEcat = 0.0;
698 float muonEactandcat = 0.0;
700 float trkLenAct = 0.f;
701 float trkLenCat = 0.f;
715 muonEact = mue_act_dw *
MuonEAct(trkLenAct);
720 muonE = mue_act_dw*muonEactandcat + mue_cat_dw*muonEcat;
723 return muonE + muonEact;
849 return (sr->
mc.
nnu == 0) ? 0. :
float(sr->
mc.
nu[0].q2);
856 return (sr->
mc.
nnu == 0) ? 0. :
float(sr->
mc.
nu[0].W2);
868 if(
abs(sr->
mc.
nu[0].pdg) != 14 || !sr->
mc.
nu[0].iscc)
871 int nprims = sr->
mc.
nu[0].prim.size();
872 for(
int iprim = 0; iprim < nprims; iprim++){
873 if(
abs(sr->
mc.
nu[0].prim[iprim].pdg) == 13){
875 return sr->
mc.
nu[0].prim[iprim].p.Vect().Mag();
889 if(
abs(sr->
mc.
nu[0].pdg) != 14 || !sr->
mc.
nu[0].iscc)
893 int nprims = sr->
mc.
nu[0].prim.size();
895 for(
int iprim = 0; iprim < nprims; iprim++){
896 if(
abs(sr->
mc.
nu[0].prim[iprim].pdg) == 13){
897 TVector3
p = sr->
mc.
nu[0].prim[iprim].p.Vect();
898 return (p-p.Dot(beamdir)*beamdir).
Mag();
911 if(
abs(sr->
mc.
nu[0].pdg) != 14 || !sr->
mc.
nu[0].iscc)
915 int nprims = sr->
mc.
nu[0].prim.size();
917 for(
int iprim = 0; iprim < nprims; iprim++){
918 if(
abs(sr->
mc.
nu[0].prim[iprim].pdg) == 13){
919 TVector3 p = sr->
mc.
nu[0].prim[iprim].p.Vect();
920 return p.Dot(beamdir);
929 return (
std::max(muonpsqr, (
float)0.0));
934 return (
std::max(muonpsqr, (
float)0.0));
939 return (
std::max(muonpsqr, (
float)0.0));
961 return float((pvec - pvec.Dot(beamdir)*beamdir).
Mag());
976 return pvec.Dot(beamdir);
985 if(
abs(sr->
mc.
nu[0].pdg) != 14 ||
997 float mup = -5, muRecoE;
998 int nprims = sr->
mc.
nu[0].prim.size();
999 for(
int iprim = 0; iprim < nprims; iprim++){
1000 if(
abs(sr->
mc.
nu[0].prim[iprim].pdg) == 13){
1001 return float(sr->
mc.
nu[0].prim[iprim].p.Vect().Mag());
1014 return float((recop-mup)/mup);
1024 if(
abs(sr->
mc.
nu[0].pdg) != 14 ||
1038 float mupt = -5, muRecoE;
1039 int nprims = sr->
mc.
nu[0].prim.size();
1040 for(
int iprim = 0; iprim < nprims; iprim++){
1041 if(
abs(sr->
mc.
nu[0].prim[iprim].pdg) == 13){
1042 TVector3 p = sr->
mc.
nu[0].prim[iprim].p.Vect();
1043 mupt = (p-p.Dot(beamdir)*beamdir).
Mag();
1055 double recopt = (recop-recop.Dot(beamdir)*beamdir).
Mag();
1056 return float((recopt-mupt)/mupt);
1066 if(
abs(sr->
mc.
nu[0].pdg) != 14 ||
1072 float mupz = -5, muRecoE;
1073 int nprims = sr->
mc.
nu[0].prim.size();
1074 for(
int iprim = 0; iprim < nprims; iprim++){
1075 if(
abs(sr->
mc.
nu[0].prim[iprim].pdg) == 13){
1076 TVector3 p = sr->
mc.
nu[0].prim[iprim].p.Vect();
1077 mupz = p.Dot(beamdir);
1090 double recopz = recop.Dot(beamdir);
1091 return float((recopz-mupz)/mupz);
1129 if(
abs(sr->
mc.
nu[0].pdg) != 14 ||
1140 double mu = 0, muReco = 0;
1142 int nprims = sr->
mc.
nu[0].prim.size();
1143 for(
int iprim = 0; iprim < nprims; iprim++){
1144 if(
abs(sr->
mc.
nu[0].prim[iprim].pdg) == 13){
1145 TVector3 mudir = sr->
mc.
nu[0].prim[iprim].p.Vect();
1146 mu = beamdir.Dot(mudir);
1147 mu = beamdir.Dot(sr->
mc.
nu[0].prim[iprim].p.Vect().Unit());
1156 muReco = beamdir.Dot(dir);
1158 return (muReco-mu)/mu;
1171 return float(sr->
mc.
nu[0].visE/sr->
mc.
nu[0].E);
const Var kRecoMuKEVsCosVsEavailUp
caf::Proxy< size_t > npng
Near Detector underground.
T max(const caf::Proxy< T > &a, T b)
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();}}return MuE;})
const Var kRecoMuKEVsCosVsEavail
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();}}return-5.0;})
caf::Proxy< caf::SRFuzzyK > fuzzyk
const Var kReMIdTrkLenAct([](const caf::SRProxy *sr){return(sr->energy.numu.ndtrklenact/100);})
Cuts and Vars for the 2020 FD DiF Study.
TH2F * hprob_dedx_vs_scat_bg
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());}}if(mu==0) return-5.0;TVector3 dir=sr->trk.kalman.tracks[0].dir;muReco=beamdir.Dot(dir);return(muReco-mu)/mu;})
const Var kIncXsecMuonEUp([](const caf::SRProxy *sr){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;})
Systematic Shift Vars and HistAxes.
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());}}return-5.0;})
const Var kTrueMuCostheta
const NuTruthVar kTrueHadEST([](const caf::SRNeutrinoProxy *sr){return float(sr->y *sr->E);})
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 kQsqr([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].q2);})
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;float ExtraHadE=sr->trk.kalman.tracks[ibesttrk].overlapE;float CalhadE=sr->slc.calE-sr->trk.kalman.tracks[ibesttrk].calE;float vishadE=CalhadE+ExtraHadE;return vishadE;})
const NuTruthVar kTrueMuKEVsCosST
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; 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;}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);})
caf::Proxy< unsigned int > idxremid
const Var kRecoMuKEVsCosVsEnu_onlyMuKEDw
TH2F * hprob_dedx_vs_scat_sig
caf::Proxy< size_t > ntracks
Proxy for caf::SRNeutrino.
caf::Proxy< caf::SRHeader > hdr
caf::Proxy< caf::SRNumuEnergy > numu
Proxy for caf::StandardRecord.
_Var< T > Var2D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb)
Variable formed from two input variables.
caf::Proxy< std::vector< caf::SRNeutrino > > nu
const Var kRecoEavail([](const caf::SRProxy *sr){return 1.68 *kNumuHadVisE(sr)+0.02 *pow(kNumuHadVisE(sr), 2);})
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; 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());} }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);})
const Var kRecoMuPsqrUp([](const caf::SRProxy *sr){float muonpsqr=std::pow(kRecoMuEUp(sr), 2)-kMuonMass2(sr);return(std::max(muonpsqr,(float) 0.0));})
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);})
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;}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);})
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();}}return ke;})
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 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;})
caf::Proxy< short int > nnu
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;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 sr->mc.nu[0].prim[iprim].p.Vect().Mag();}}return-5.0;})
const NuTruthVar kTrueMuKEVsCosVsEnuST
caf::Proxy< caf::SREnergyBranch > energy
const Var kRecoMuCosthetaUp([](const caf::SRProxy *sr){return(float) std::cos(kRecoMuThetaUp(sr));})
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;llr=MuonLLR(sr->trk.kalman.tracks[bestidx].dedxllh, sr->trk.kalman.tracks[bestidx].scatllh);return llr;})
T sqr(T x)
More efficient square function than pow(x,2)
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
const Var kReMIdTrkLenCat([](const caf::SRProxy *sr){return(sr->energy.numu.ndtrklencat/100);})
_Var< T > Var3D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb, const _Var< T > &c, const Binning &binsc)
This is just like a Var2D, but useful for 3D Spectra.
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
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);})
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);}}return-5.0;})
const NuTruthVar kTrueEST([](const caf::SRNeutrinoProxy *nu){return nu->E;})
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(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat< 0) return 0.f; 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(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;})
caf::Proxy< caf::SRElastic > elastic
const Var kIncXsecMuonEDw([](const caf::SRProxy *sr){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;})
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 kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
float VisibleHadE(float vishadE)
const Var kRecoMuKEVsCosVsEnu_onlyMuKEUp
const Var kRecoMuKEVsCosVsEavailAngleUp
const NuTruthVar kTrueMuKEVsCosVsEavailST
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 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);})
const Var kRemid([](const caf::SRProxy *sr){return(sr->sel.remid.pid);})
caf::Proxy< caf::SRTrackBranch > trk
float MuonEActandCat(double trklenactandcat)
const Var kRecoMuKEVsCosVsEnuDw
std::string getenv(std::string const &name)
const Var kweightedten([](const caf::SRProxy *){return 1.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.;})
const Var kMostMuonLikeTrackIdx([](const caf::SRProxy *sr){float maxllr=-1e12;int bestidx=-5;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;}}return bestidx;})
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(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(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float(sr->trk.kalman.tracks[ibesttrk].leninact/100.); if(sr->trk.kalman.tracks[ibesttrk].leninact< 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return 0.f;return-1000.f;})
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;return(float) sr->trk.kalman.tracks[ibesttrk].overlapE;})
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;float ExtraHadE=sr->trk.kalman.tracks[ibesttrk].overlapE;float CalhadE=sr->slc.calE-sr->trk.kalman.tracks[ibesttrk].calE;float vishadE=CalhadE+ExtraHadE;float hadE=0.0;hadE=VisibleHadE(vishadE);return hadE;})
const Var kRecoMuKEVsCosVsEnuUp
caf::Proxy< caf::SRRemid > remid
float TAMuE(double trklen)
float MuonEAct(double trklenact)
const NuTruthVar kTrueYST([](const caf::SRNeutrinoProxy *sr){return float(sr->y);})
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());})
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;})
const Var kNProngs([](const caf::SRProxy *sr){return(sr->vtx.elastic[0].fuzzyk.npng);})
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;})
caf::Proxy< caf::SRTruthBranch > mc
const Var kRecoMuKEVsCosVsEavailAngleDw
const Var kRecoMuKEVsCosVsEavailDw
caf::Proxy< caf::SRKalman > kalman
caf::Proxy< caf::SRSlice > slc
const Var kRecoMuCosthetaDw([](const caf::SRProxy *sr){return(float) std::cos(kRecoMuThetaDw(sr));})
Var Sqrt(const Var &v)
Use to take sqrt of a var.
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();})
caf::Proxy< float > ndtrklenact
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);})
float MuonECat(double trklencat)
const Var kRecoMuPsqr([](const caf::SRProxy *sr){float muonpsqr=std::pow(kRecoMuE(sr), 2)-kMuonMass2(sr);return(std::max(muonpsqr,(float) 0.0));})
float MuonLLR(const float dedxll, const float scatll)
const Var kRecoMuPsqrDw([](const caf::SRProxy *sr){float muonpsqr=std::pow(kRecoMuEDw(sr), 2)-kMuonMass2(sr);return(std::max(muonpsqr,(float) 0.0));})
assert(nhit_max >=nhit_nbins)
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
const Var kDummy([](const caf::SRProxy *sr){return 1;})
caf::Proxy< short int > pdg
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 kMuLength([](const caf::SRProxy *sr){return sr->trk.kalman.tracks[0].len/100;})
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
caf::Proxy< caf::SRIDBranch > sel
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 NuTruthVar kTrueEavailST
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();})
caf::Proxy< caf::SRVertexBranch > vtx
const Binning angbinsCustom
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 kWsqr([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].W2);})
caf::Proxy< float > ndtrklencat
caf::Proxy< std::vector< caf::SRTrueParticle > > prim
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);})
const Var kResVsTrueCostheta
const Var kRecoMuKEVsCosVsEnu
const Var kResVsRecoCostheta
const Var kMuonMass2([](const caf::SRProxy *sr){return util::sqr(MuonMass());})