make_nueFDprediction_kinematics_REW.C
Go to the documentation of this file.
1 // Macro to make nominal and shifted predictions
2 // Takes names and datasets as input
3 // Create files for nominal and shifted samples
4 // set systematicSample to false for beam or genie systematics
5 // USAGE
6 // cafe -bq make_nueFDprediction_shifted.C bool "ShortName" "Long Name"
7 // "filename" "nd data" "nd mc dataset" "FDMC nonswap" "swap" "tau"
8 // See full examples in run_nueFDprediction*.sh
9 //hastau,plus2,plus1,minu1,minu2,knoshift,compareto
10 
11 // Analysis
13 #include "CAFAna/Analysis/Plots.h"
14 #include "CAFAna/Analysis/SALoaders.h"
15 #include "CAFAna/Analysis/Style.h"
16 
17 // Core
18 #include "CAFAna/Core/Binning.h"
19 #include "CAFAna/Core/Cut.h"
20 #include "CAFAna/Core/HistAxis.h"
22 #include "CAFAna/Core/Spectrum.h"
25 
26 // Cuts
27 #include "CAFAna/Cuts/Cuts.h"
29 #include "CAFAna/Cuts/NueCutsSecondAna.h"
30 #include "CAFAna/Cuts/SpillCuts.h"
32 
33 
34 // Decomp
39 
40 // Extrap
42 
43 // Predict
46 
47 // Systs
48 #include "CAFAna/Systs/Systs.h"
49 
50 // Vars
52 #include "CAFAna/Vars/HistAxes.h"
54 #include "CAFAna/Vars/Vars.h"
57 
58 #include "OscLib/OscCalcDumb.h"
59 
61 
62 
63 #include <iostream>
64 #include <iomanip>
65 
66 using namespace ana;
67 
68 void CheckFileOverwrite(TString);
69 
71 
72 void make_nueFDprediction_kinematics_REW(const std::string& outfilename = "FDprediction_kinematics_REW.root",
73  const bool hastau = false
74  )
75 
76 {
77  //********** 0. Basic Setup *********************
78 
79  // ND loaders
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";
82 
83  // FD loaders
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";
87 
88  //******* 1. Set up selection and variables ********
89  const HistAxis axisNumu = kNumuNonQEAxisFirstAna;
90 
91  const Var kWOscDumb([](const caf::SRProxy* sr){return sr->mc.nu[0].woscdumb;});
92 
93  // open the weights file
94  TFile *weight_file = new TFile("/nova/app/users/nostom66/nue/extrapsysts/combined/Weights.root","READ");
95 
96  // construction of trueQ2 weight
97  TH1 *h = new TH1();
98  gDirectory->GetObject("trueQ2_weight_AllSamples",h);
99 
100  std::cout << h->GetNbinsX() << endl;
101  int nbins = h->GetNbinsX();
102  double bins[nbins], binslow[nbins], binshigh[nbins];
103 
104  for(int i = 1; i <= nbins; ++i){
105  bins[i-1] = h->GetBinContent(i);
106  binslow[i-1] = h->GetBinLowEdge(i);
107  binshigh[i-1] = h->GetBinLowEdge(i) + h->GetBinWidth(i);
108  }
109 
110  const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](const caf::SRProxy *sr)
111  {
112  for(int i = 0; i < nbins; ++i){
113  if((kRecoQ2(sr) > binslow[i]) && (kRecoQ2(sr) < binshigh[i])){return bins[i];}
114  }
115  return 1.0;
116  }
117  );
118 
119  //construction of PtP weight
120  TH1 *hptp = new TH1();
121  gDirectory->GetObject("PtP_weight_AllSamples",hptp);
122 
123  std::cout << hptp->GetNbinsX() << endl;
124  int nbins_ptp = hptp->GetNbinsX();
125  double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp];
126 
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);
131  }
132 
133  const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](const caf::SRProxy *sr)
134  {
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];}
137  }
138  return 1.0;
139  }
140  );
141 
142  //construction of Cos weight
143  TH1 *hcos = new TH1();
144  gDirectory->GetObject("CosNumi_weight_AllSamples",hcos);
145 
146  std::cout << hcos->GetNbinsX() << endl;
147  int nbins_cos = hptp->GetNbinsX();
148  double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos];
149 
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);
154  }
155 
156  const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](const caf::SRProxy *sr)
157  {
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];}
160  }
161  return 1.0;
162  }
163  );
164 
165  // close weights file
166  weight_file->Close();
167 
168 
169 
170  const int kNumVars = 5;
171 
172  const HistDef defs[kNumVars] = {
173  {"recoE", {"E_{reco} (GeV)", Binning::Simple(10, 0, 5), kNueEnergy2017}},
174  {"trueQ2", {"true Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kTrueQ2}},
175  {"trueW2", {"true W^{2} (GeV^{2})", Binning::Simple(100, 0, 2), kTrueW}},
176  {"PtP", {"p_{t}/p", Binning::Simple(500,0,1), kPtP}},
177  {"CosTheta", {"CosTheta", Binning::Simple(500,0,1), kCosTheta}}
178  };
179 
180  const int kNumSels = 1;
181 
182  const Cut selFD[kNumSels] = {
184  };
185  const Cut selND[kNumSels] = {
186  kNue2017NDPresel && kNueSecondAnaCVNeSsb
187  };
188  const std::string selNames[kNumSels] = {
189  "AllSamples"
190  };
191  const Cut cutNDNumu = {
192  kNumuNDCvn
193  };
194 
196 
197  //********** 2. Set up loaders *********************
198  SANueExtrapLoaders loaders;
199 
200  loaders.SetLoaderPath(loaderNDData, caf::kNEARDET, ana::Loaders::DataMC::kData);
201  loaders.SetLoaderPath(loaderNDMC, caf::kNEARDET, ana::Loaders::DataMC::kMC);
202 
203  loaders.SetLoaderPath(loaderFDnonswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam, ana::Loaders::SwappingConfig::kNonSwap);
204  loaders.SetLoaderPath(loaderFDfluxswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam,ana::Loaders::SwappingConfig::kFluxSwap);
205  if(hastau)
207  else
209 
210  loaders.SetSpillCut(kStandardSpillCuts);
211 
212  //******* 3. Set up the extrapolations and go *******
213  SystShifts thisShift = kNoShift;
214 
215 
216  // reweighted numudecomps
217  NumuDecomp * numuDecomp_trueQ2;
218  NumuDecomp * numuDecomp_PtP;
219  NumuDecomp * numuDecomp_Cos;
220 
221  // propdecomp
222  ProportionalDecomp * propDecomp[kNumSels][kNumVars];
223 
224  // extraps from different numudecomps
225  ModularExtrap * extrapProp_trueQ2[kNumSels][kNumVars];
226  ModularExtrap * extrapProp_PtP[kNumSels][kNumVars];
227  ModularExtrap * extrapProp_Cos[kNumSels][kNumVars];
228 
229  PredictionExtrap * predicProp_trueQ2[kNumSels][kNumVars];
230  PredictionExtrap * predicProp_PtP[kNumSels][kNumVars];
231  PredictionExtrap * predicProp_Cos[kNumSels][kNumVars];
232 
233  PredictionNoExtrap * predicNoXP[kNumSels][kNumVars];
234 
235 
236 
237  numuDecomp_trueQ2 = new NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt*kTrueQ2Weight);
238  numuDecomp_PtP = new NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt*kPtPWeight);
239  numuDecomp_Cos = new NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt*kCosWeight);
240 
241  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
242  for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
243 
244  const HistAxis& axisNue = defs[varIdx].axis;
245 
246  propDecomp[selIdx][varIdx] = new ProportionalDecomp (loaders,
247  axisNue, selND[selIdx],
248  thisShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt );
249 
250  // trueQ2 reweighted extrap
251  extrapProp_trueQ2[selIdx][varIdx] = NueExtrap (loaders,
252  *propDecomp[selIdx][varIdx], *numuDecomp_trueQ2,
253  axisNue, axisNumu,
254  selFD[selIdx], selND[selIdx], cutNDNumu,
255  thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
256 
257  // PtP reweighted extrap
258  extrapProp_PtP[selIdx][varIdx] = NueExtrap (loaders,
259  *propDecomp[selIdx][varIdx], *numuDecomp_PtP,
260  axisNue, axisNumu,
261  selFD[selIdx], selND[selIdx], cutNDNumu,
262  thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
263  // Cos reweighted extrap
264  extrapProp_Cos[selIdx][varIdx] = NueExtrap (loaders,
265  *propDecomp[selIdx][varIdx], *numuDecomp_Cos,
266  axisNue, axisNumu,
267  selFD[selIdx], selND[selIdx], cutNDNumu,
268  thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
269 
270 
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]);
274 
275  predicNoXP[selIdx][varIdx] = new PredictionNoExtrap (loaders,
276  axisNue,
277  selFD[selIdx], thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
278 
279  }//vars
280  }//sels
281  loaders.Go();
282 
283  //********** 4. Save everything ********************
284  TFile * file = new TFile(outfilename.c_str(),"RECREATE");
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();
289 
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) ;
294  }
295  }
296  file->Close();
297 
298 } // end
299 
300 //****************************************
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();
307  }
308 }
Near Detector underground.
Definition: SREnums.h:10
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)
Definition: NumuVars.h:146
const int kNumVars
Definition: vars.h:14
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const Color_t kMC
std::string name
Definition: NuePlotLists.h:12
void make_nueFDprediction_kinematics_REW(const std::string &outfilename="FDprediction_kinematics_REW.root", const bool hastau=false)
void CheckFileOverwrite(TString)
const Var kTrueQ2
Definition: TruthVars.h:27
string outfilename
knobs that need extra care
const int nbins
Definition: cellShifts.C:15
const Var kPtP
Transverse momentum fraction in slice.
Definition: NueVars.cxx:90
Uses MC for NC and CC components, assigns remainder of data to CC.
Definition: NumuDecomp.h:10
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));})
Definition: NueEnergy2017.h:11
const Cut kNumuNDCvn
Definition: NumuCuts.h:62
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
caf::StandardRecord * sr
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
HistAxis axis
Definition: NuePlotLists.h:13
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const HistAxis axisNue
Definition: syst_header.h:378
const Cut kNue2017NDPresel
Definition: NueCuts2017.h:285
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Var kWOscDumb([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.f;return float(sr->mc.nu[0].woscdumb);})
Definition: TruthVars.h:10
void NueExtrap(string beam="fhc", string cvntype="oldpresel")
Definition: NueExtrap.C:28
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Splits Data proportionally according to MC.
std::vector< Loaders * > loaders
Definition: syst_header.h:386
TFile * file
Definition: cellShifts.C:17
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Prediction that just uses FD MC, with no extrapolation.
const HistAxis axisNumu
Take the output of an extrapolation and oscillate it as required.
const Var kCosTheta
osc::OscCalcDumb calc
const std::string selNames[kNumSels]
Definition: vars.h:46
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kTrueW
Definition: TruthVars.h:22
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;})
Definition: NumuVars.h:43
const Cut kNue2017FDAllSamples
Our FD selection including all samples, for making predictions, etc.
Definition: NueCuts2017.h:155
enum BeamMode string