5 #include "CAFAna/Core/HistAxis.h" 10 #include "Utilities/func/MathUtil.h" 21 #include "TParticlePDG.h" 22 #include "TDatabasePDG.h" 23 #include "TLorentzVector.h" 101 static TDatabasePDG* pdgdb = TDatabasePDG::Instance();
102 static double mass = pdgdb->GetParticle(13)->Mass();
118 static TDatabasePDG* pdgdb = TDatabasePDG::Instance();
119 static double mass = pdgdb->GetParticle(2112)->Mass();
228 "Reconstructed T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}",
232 "True T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}",
239 "Reconstructed Neutrino Energy (GeV)",
243 "True Neutrino Energy (GeV)",
252 "True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
254 kTrueMuKEVsCosVsEnuST);
260 "True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
262 kTrueMuKEVsCosVsEavailST);
268 "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
270 kRecoMuKEVsCosVsEnu);
273 "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
275 kRecoMuKEVsCosVsEavail);
278 "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
280 kRecoMuKEVsCosVsEnu_onlyMuKEUp);
282 "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
284 kRecoMuKEVsCosVsEnu_onlyMuKEDw);
286 "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
288 kRecoMuKEVsCosVsEnuUp);
290 "Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)",
292 kRecoMuKEVsCosVsEnuDw);
295 "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
297 kRecoMuKEVsCosVsEavailUp);
299 "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
301 kRecoMuKEVsCosVsEavailDw);
305 "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
307 kRecoMuKEVsCosVsEavailAngleUp);
309 "Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)",
311 kRecoMuKEVsCosVsEavailAngleDw);
482 "Reconstructed Neutrino Energy (GeV)",
488 "Reconstructed Neutrino Energy (GeV)",
const HistAxis kRecoEStandardAxis_MuonEUp("Reconstructed Neutrino Energy (GeV)", enubins, kRecoEUp)
const Var kRecoMuKEVsCosVsEavailUp
const NuTruthHistAxis kTrueMuKEVsCosVsEnuStandardAxisST("True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kTrueMuKEVsCosVsEnuST)
const Binning angvsmukevseavailbins
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 HistAxis kTrueEStandardAxis
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;})
const HistAxis kRecoMuKEVsCosVsEnuStandardAxis_onlyMuKEDw("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnu_onlyMuKEDw)
const HistAxis kRecoQ2StandardAxis("Reco Q2 (GeV)", q2bins, kRecoq2)
const NuTruthHistAxis kTrueMuKEVsCosStandardAxisST("True T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}", angvsmukebins, kTrueMuKEVsCosST)
const Var kReMIdTrkLenAct([](const caf::SRProxy *sr){return(sr->energy.numu.ndtrklenact/100);})
Cuts and Vars for the 2020 FD DiF Study.
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 HistAxis kRecoQ2StandardAxis_MuonEDw("Reco Q2 (GeV)", q2bins, kRecoq2_MuonEDw)
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 kRecoq2_AngleDw
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);})
const Var kRecoq2_AngleUp
const Var kRecoMuKEVsCosVsEnu_onlyMuKEDw
const Var kNeutronMass([](const caf::SRProxy *sr){return NeutronMass();})
Proxy for caf::StandardRecord.
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;})
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
const HistAxis kTrueMuKEVsCosVsEnuStandardAxis
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)
const HistAxis kRecoMuKEVsCosStandardAxis("Reconstructed T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}", angvsmukebins, kRecoMuKEVsCos)
const Var kReMIdTrkLenCat([](const caf::SRProxy *sr){return(sr->energy.numu.ndtrklencat/100);})
const HistAxis kRecoQ2StandardAxis_AngleDw("Reco Q2 (GeV)", q2bins, kRecoq2_AngleDw)
Var Constant(double c)
Use to weight events up and down by some factor.
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 HistAxis kRecoQ2StandardAxis_AngleUp("Reco Q2 (GeV)", q2bins, kRecoq2_AngleUp)
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;})
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 HistAxis kRecoMuKEVsCosVsEavailStandardAxisAngleDw("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailAngleDw)
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 NuTruthHistAxis kTrueQ2StandardAxisST("True Q2 (GeV)", q2bins, kTrueQ2_NT)
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);})
float MuonEActandCat(double trklenactandcat)
const HistAxis kRecoMuKEVsCosVsEavailStandardAxisUp("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailUp)
const Var kRecoMuKEVsCosVsEnuDw
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 HistAxis kRecoMuKEVsCosVsEnuStandardAxisUp("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnuUp)
const Binning angvsmukebins
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 HistAxis kRecoEStandardAxis_MuonEDw("Reconstructed Neutrino Energy (GeV)", enubins, kRecoEDw)
const HistAxis kRecoMuKEVsCosVsEnuStandardAxisDw("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnuDw)
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;})
HistAxis HistAxisFromNuTruthHistAxis(NuTruthHistAxis ntha, double _default)
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
const HistAxis kTrueMuKEVsCosVsEavailStandardAxis
float MuonEAct(double trklenact)
const Var kRecoq2_MuonEUp
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 kTrkLenActandCat
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;})
const Var kNeutronMass2([](const caf::SRProxy *sr){return util::sqr(NeutronMass());})
const Var kRecoMuKEVsCosVsEavailAngleDw
const HistAxis kRecoMuKEVsCosVsEnuStandardAxis_onlyMuKEUp("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnu_onlyMuKEUp)
const Var kRecoMuKEVsCosVsEavailDw
const Var kThetaSystShift
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
const HistAxis kRecoMuKEVsCosVsEnuStandardAxis("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnu)
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 kMuonMass([](const caf::SRProxy *sr){return MuonMass();})
const HistAxis kRecoMuKEVsCosVsEavailStandardAxisAngleUp("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailAngleUp)
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 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 Binning angvsmukevsebins
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));})
const Var kDummy([](const caf::SRProxy *sr){return 1;})
const NuTruthHistAxis kTrueMuKEVsCosVsEavailStandardAxisST("True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevseavailbins, kTrueMuKEVsCosVsEavailST)
const HistAxis kTrueMuKEVsCosStandardAxis
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;})
const HistAxis kRecoMuKEVsCosVsEavailStandardAxis("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavail)
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();})
const HistAxis kRecoMuKEVsCosVsEavailStandardAxisDw("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailDw)
const NuTruthVar kTrueEavail_NT([](const caf::SRNeutrinoProxy *nu){double eAvail=0;for(const auto &particle:nu->prim){double particle_Eavail=0;if(particle.pdg==2212||abs(particle.pdg)==211){double gamma=particle.p.Gamma();particle_Eavail=(gamma-1)/gamma *particle.p.E;}else if(particle.pdg==111 or abs(particle.pdg)==11 or particle.pdg==22) particle_Eavail=particle.p.E;else if(particle.pdg >=2000000000){}else if(particle.pdg >=1000000000){}else if(particle.pdg >=2000 &&particle.pdg!=2212 &&particle.pdg!=2112){particle_Eavail=particle.p.E-0.9382;}else if(particle.pdg<=-2000){particle_Eavail=particle.p.E+0.9382;}else if(particle.pdg!=2112 &&(abs(particle.pdg)< 11||abs(particle.pdg) > 16)){particle_Eavail=particle.p.E;}eAvail+=particle_Eavail;}return eAvail;})
MINERvA "true Eavail" (used in 2p2h xsec)
const HistAxis kTrueQ2StandardAxis
const HistAxis kRecoEStandardAxis("Reconstructed Neutrino Energy (GeV)", enubins, kRecoE)
const NuTruthHistAxis kTrueEStandardAxisST("True Neutrino Energy (GeV)", enubins, kTrueEST)
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);})
const NuTruthVar kTrueQ2_NT([](const caf::SRNeutrinoProxy *nu){return nu->q2;})
True square of four-momentum transfer.
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 kRecoq2_MuonEDw
const Var kResVsTrueCostheta
const Var kRecoMuKEVsCosVsEnu
const Var kResVsRecoCostheta
const Var kMuonMass2([](const caf::SRProxy *sr){return util::sqr(MuonMass());})
const HistAxis kRecoQ2StandardAxis_MuonEUp("Reco Q2 (GeV)", q2bins, kRecoq2_MuonEUp)