14 #include "CAFAna/Analysis/SALoaders.h" 18 #include "CAFAna/Core/Binning.h" 19 #include "CAFAna/Core/Cut.h" 20 #include "CAFAna/Core/HistAxis.h" 29 #include "CAFAna/Cuts/NueCutsSecondAna.h" 73 const bool hastau =
false 80 std::string loaderNDData =
"prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_nue_or_numu_or_nus_contain_v1_goodruns";
81 std::string loaderNDMC =
"prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
84 std::string loaderFDnonswapMC =
"prod_decaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
85 std::string loaderFDfluxswapMC =
"prod_decaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
86 std::string loaderFDtauswapMC =
"prod_decaf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
94 TFile *weight_file =
new TFile(
"/nova/app/users/nostom66/nue/extrapsysts/combined/Weights.root",
"READ");
98 gDirectory->GetObject(
"trueQ2_weight_AllSamples",h);
101 int nbins = h->GetNbinsX();
105 bins[
i-1] = h->GetBinContent(
i);
106 binslow[
i-1] = h->GetBinLowEdge(
i);
107 binshigh[
i-1] = h->GetBinLowEdge(
i) + h->GetBinWidth(
i);
110 const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](
const caf::SRProxy *sr)
113 if((
kRecoQ2(sr) > binslow[
i]) && (
kRecoQ2(sr) < binshigh[i])){
return bins[
i];}
120 TH1 *hptp =
new TH1();
121 gDirectory->GetObject(
"PtP_weight_AllSamples",hptp);
124 int nbins_ptp = hptp->GetNbinsX();
125 double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp];
127 for(
int i = 1;
i <= nbins_ptp; ++
i){
128 bins_ptp[
i-1] = hptp->GetBinContent(
i);
129 binslow_ptp[
i-1] = hptp->GetBinLowEdge(
i);
130 binshigh_ptp[
i-1] = hptp->GetBinLowEdge(
i) + hptp->GetBinWidth(
i);
133 const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](
const caf::SRProxy *sr)
135 for(
int i = 0;
i < nbins_ptp; ++
i){
136 if((
kPtP(sr) > binslow_ptp[
i]) && (
kPtP(sr) < binshigh_ptp[i])){
return bins_ptp[
i];}
143 TH1 *hcos =
new TH1();
144 gDirectory->GetObject(
"CosNumi_weight_AllSamples",hcos);
147 int nbins_cos = hptp->GetNbinsX();
148 double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos];
150 for(
int i = 1;
i <= nbins_cos; ++
i){
151 bins_cos[
i-1] = hcos->GetBinContent(
i);
152 binslow_cos[
i-1] = hcos->GetBinLowEdge(
i);
153 binshigh_cos[
i-1] = hcos->GetBinLowEdge(
i) + hcos->GetBinWidth(
i);
156 const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](
const caf::SRProxy *sr)
158 for(
int i = 0;
i < nbins_cos; ++
i){
159 if((
kCosNumi(sr) > binslow_cos[
i]) && (
kCosNumi(sr) < binshigh_cos[i])){
return bins_cos[
i];}
166 weight_file->Close();
191 const Cut cutNDNumu = {
241 for(
unsigned int selIdx = 0; selIdx <
kNumSels; ++selIdx){
242 for(
unsigned int varIdx = 0; varIdx <
kNumVars; ++varIdx){
247 axisNue, selND[selIdx],
251 extrapProp_trueQ2[selIdx][varIdx] =
NueExtrap (loaders,
252 *propDecomp[selIdx][varIdx], *numuDecomp_trueQ2,
254 selFD[selIdx], selND[selIdx], cutNDNumu,
258 extrapProp_PtP[selIdx][varIdx] =
NueExtrap (loaders,
259 *propDecomp[selIdx][varIdx], *numuDecomp_PtP,
261 selFD[selIdx], selND[selIdx], cutNDNumu,
264 extrapProp_Cos[selIdx][varIdx] =
NueExtrap (loaders,
265 *propDecomp[selIdx][varIdx], *numuDecomp_Cos,
267 selFD[selIdx], selND[selIdx], cutNDNumu,
271 predicProp_trueQ2[selIdx][varIdx] =
new PredictionExtrap (extrapProp_trueQ2[selIdx][varIdx]);
272 predicProp_PtP[selIdx][varIdx] =
new PredictionExtrap (extrapProp_PtP[selIdx][varIdx]);
273 predicProp_Cos[selIdx][varIdx] =
new PredictionExtrap (extrapProp_Cos[selIdx][varIdx]);
285 TDirectory * dFDshi = file->mkdir(
"prediction");
286 for(
unsigned int selIdx = 0; selIdx <
kNumSels; ++selIdx){
287 for(
unsigned int varIdx = 0; varIdx <
kNumVars; ++varIdx){
288 TString varName = defs[varIdx].
name.c_str();
290 predicProp_trueQ2[selIdx][varIdx]->
SaveTo(dFDshi,
"nue_pred_trueQ2_"+selNames[selIdx]+
"_"+varName) ;
291 predicProp_PtP[selIdx][varIdx]->
SaveTo(dFDshi,
"nue_pred_PtP_"+selNames[selIdx]+
"_"+varName) ;
292 predicProp_Cos[selIdx][varIdx]->
SaveTo(dFDshi,
"nue_pred_Cos_"+selNames[selIdx]+
"_"+varName) ;
293 predicNoXP[selIdx][varIdx]->
SaveTo(dFDshi,
"nue_pred_noextrap_"+selNames[selIdx]+
"_"+varName) ;
302 if(!gSystem->AccessPathName(outfilename.Data())){
303 std::cout <<
"\n\nThis will overwrite " << outfilename
304 <<
"\n\nAre you sure you want to continue? y/n ";
305 string input =
""; getline(cin,input);
306 if (input!=
"y") abort();
Near Detector underground.
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 four-momentum transfer invariant (Q^2)
Far Detector at Ash River.
Cuts and Vars for the 2020 FD DiF Study.
Simple record of shifts applied to systematic parameters.
Proxy for caf::StandardRecord.
caf::Proxy< std::vector< caf::SRNeutrino > > nu
void make_nueFDprediction_kinematics_REW(const std::string &outfilename="FDprediction_kinematics_REW.root", const bool hastau=false)
void CheckFileOverwrite(TString)
string outfilename
knobs that need extra care
const Var kPtP
Transverse momentum fraction in slice.
Uses MC for NC and CC components, assigns remainder of data to CC.
const Var kNueEnergy2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2017FDFit(kCVNemE(sr), kCVNhadE(sr));})
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it's s...
const HistDef defs[kNumVars]
const SystShifts kNoShift
const Cut kNue2017NDPresel
caf::Proxy< caf::SRTruthBranch > mc
const Var kWOscDumb([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.f;return float(sr->mc.nu[0].woscdumb);})
Splits Data proportionally according to MC.
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
const std::string selNames[kNumSels]
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
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 Cut kNue2017FDAllSamples
Our FD selection including all samples, for making predictions, etc.