xsec_extrap_plots.C
Go to the documentation of this file.
3 #include "CAFAna/Core/Loaders.h"
4 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Cuts/NueCutsSecondAna.h"
12 #include "CAFAna/Systs/Systs.h"
14 #include "CAFAna/Vars/HistAxes.h"
16 #include "temp_draw_plots_util.h"
17 
18 #include <TFile.h>
19 #include <TString.h>
20 #include <vector>
21 #include <iostream>
22 
23 using namespace ana;
24 
25 void SetLoaderPaths2017(Loaders& loaders, const bool usedecafs = false)
26 {
27  std::string ldrNDData = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_nue2017";
28  std::string ldrNDMC = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_nue2017";
29  std::string ldrNonSwap = "prod_sumdecaf_R17-03-01-prod3reco.g_fd_genie_nonswap_fhc_nova_v08_full_v1_nue2017";
30  std::string ldrFlxSwap = "prod_sumdecaf_R17-03-01-prod3reco.g_fd_genie_fluxswap_fhc_nova_v08_full_v1_nue2017";
31  std::string ldrTauSwap = "prod_sumdecaf_R17-03-01-prod3reco.g_fd_genie_tau_fhc_nova_v08_full_v1_nue2017";
32 
33  if (usedecafs) {
34  ldrNDData = "prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_nue_or_numu_or_nus_contain_v1_goodruns";
35  ldrNDMC = "prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
36  ldrNonSwap = "prod_decaf_R17-03-01-prod3reco.g_fd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
37  ldrFlxSwap = "prod_decaf_R17-03-01-prod3reco.g_fd_genie_fluxswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
38  ldrTauSwap = "prod_decaf_R17-03-01-prod3reco.g_fd_genie_tau_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
39  }
40 
41 
42  loaders.SetLoaderPath(ldrNDData, caf::kNEARDET, Loaders::kData,
44  loaders.SetLoaderPath(ldrNDMC, caf::kNEARDET, Loaders::kMC,
46  loaders.SetLoaderPath(ldrNonSwap, caf::kFARDET, Loaders::kMC,
48  loaders.SetLoaderPath(ldrFlxSwap, caf::kFARDET, Loaders::kMC,
50  loaders.SetLoaderPath(ldrTauSwap, caf::kFARDET, Loaders::kMC,
52 
53 }
54 
55 
56 void xsec_extrap_plots(bool createFile = false,
57  TString suffix = "test")
58 {
59  TString filename = "pred_xsec_fhc_" + suffix + ".root";
60 
61  struct WeiDef {
62  TString name;
63  Var weight;
64  };
65  struct ShiDef {
66  TString name;
67  SystShifts shift;
68  };
69 
70  std::vector <WeiDef> weights;
71  weights.push_back({"ppfx_xsec", kPPFXFluxCVWgt * kXSecCVWgt2017});
72  weights.push_back({"ppfx_only", kPPFXFluxCVWgt});
73 
74  std::vector <ShiDef> shifts;
75  shifts.push_back({"no_shift", kNoShift });
76 
77  for (auto sigma:{+1, -1}) {
78  TString sig_str = (sigma == 1 ? "plus" : "minus");
79 
80  shifts.push_back({"MEC_enushape_" + sig_str,
82  shifts.push_back({"MEC_q0shape_" + sig_str,
84  shifts.push_back({"MEC_NPfrac_" + sig_str,
86  shifts.push_back({"MEC_scaleSA_" + sig_str,
88  shifts.push_back({"MEC_all3_" + sig_str,
89  SystShifts({
91  { &kMECq0ShapeSyst2017 , sigma},
93  })
94  });
95  } //end sigma
96 
97  if (createFile) {
98 
100  bool use_decafs = (suffix.Contains("test") ||
101  suffix.Contains("decaf"));
102  SetLoaderPaths2017(loaders, use_decafs);
104 
105  const HistAxis axisNumu = kNumuNonQEAxisFirstAna;
106  const HistAxis axisNue = kNueAxisSecondAnaCVN;
107 
108  std::vector <IPrediction*> v_pred_prop;
109  std::vector <IPrediction*> v_pred_nxp;
110  std::vector <TString> names;
111 
112  // Loop over weights and shifts
113  // Don't do all variations: e.g. only nominal weight for shifts
114 
115  for (auto weight:weights) {
116 
117  auto generator_nxp = NoExtrapGenerator(axisNue,
118  kNueSAFDCVNSsb,
119  weight.weight);
120 
121  for (auto shift:shifts) {
122 
123  if (shift.name == "no_shift" || weight.name == "ppfx_xsec") {
124 
125  auto generator_prop = NuePropExtrapGenerator(axisNue,
126  axisNumu,
127  kNueSAFDCVNSsb,
128  kNueNDCVNSsb,
129  kNumuND,
130  kNoShift,
131  weight.weight);
132 
133  auto pred_prop = generator_prop.Generate(loaders,
134  shift.shift);
135  auto pred_nxp = generator_nxp .Generate(loaders,
136  shift.shift);
137 
138  v_pred_prop.emplace_back(pred_prop.release());
139  v_pred_nxp. emplace_back(pred_nxp. release());
140 
141  names.push_back (shift.name + "_" + weight.name);
142  }//end if
143  }//end shifts
144  }//end weights
145 
146  loaders.Go();
147 
148  TFile* file = new TFile(filename, "recreate");
149 
150  for (int i = 0 ; i < (int) names.size() ; ++i) {
151  v_pred_nxp [i]->SaveTo(file, "pred_nxp_" + names[i]);
152  v_pred_prop[i]->SaveTo(file, "pred_xp_prop_" + names[i]);
153  }
154 
155  delete file;
156  std::cout << "Saved predictions to " << filename << std::endl;
157 
158  } //end create file
159  else {
160  // Make some plots
161  TFile* file = new TFile (filename, "READ");
162  TFile* outfile_xpplots = new TFile("plots/xp/modExtrapNuePlots.root",
163  "RECREATE") ;
164  bool printtable_numu = 1;
165  bool printtable_nue = 1;
166  bool pidaxes = true;
167  bool printtable = 1;
168 
169  // Weights are on/off. Compare 1 to nominal
170 
171  TDirectory* dir_pred_no, * dir_pred_sh;
172  TDirectory* dir_decomp_no, * dir_decomp_sh;
173  TString dirname_decomp, plot_title, plot_name;
174 
175  for (auto weight:weights) {
176  if (weight.name != "ppfx_xsec") // Skip nominal vs nominal
177  for (TString xp: {"nxp", "xp_prop"}) {
178 
179  dir_pred_no = file->GetDirectory ("pred_" + xp +
180  "_no_shift_ppfx_xsec");
181  dir_pred_sh = file->GetDirectory ("pred_" + xp +
182  "_no_shift_" +
183  weight.name);
184  //ND plots
185  if (xp == "xp_prop") {
186  //Numu
187  dirname_decomp = "extrap/MEextrap/Decomp";
188  dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
189  dir_decomp_sh = dir_pred_sh->GetDirectory (dirname_decomp);
190 
191  plot_title = weight.name + " (NearDet MC);" +
192  " Non-Quasielastic Energy Estimator (GeV); Events / Bin";
193  plot_name = "compare_numu_ND_" + weight.name;
194 
195  CompareNDDataMC (dir_decomp_no, dir_decomp_sh,
196  plot_title, plot_name, "Numu",
197  printtable_numu);
198 
199  //nue
200  dirname_decomp = "extrap/EEextrap/Decomp/";
201  plot_title = weight.name + " (NearDet MC)"+
202  "; Neutrino Energy Estimator (GeV) X CVN; Events / Bin";
203  plot_name = "compare_nue_ND_" + weight.name;
204 
205  dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
206  dir_decomp_sh = dir_pred_sh->GetDirectory (dirname_decomp);
207  CompareNDDataMC (dir_decomp_no, dir_decomp_sh,
208  plot_title, plot_name, "CVN",
209  printtable_nue, pidaxes);
210  } //end nd plots
211 
212  //FD Plots
213  plot_title = weight.name +
214  (xp != "nxp" ? " FarDet Extrap":" (FarDet MC)") +
215  ";Neutrino Energy Estimator (GeV) X CVN;"+ "Events / Bin";
216  plot_name = "compare_" + xp + "_" + weight.name;
217 
218  auto pred_no = ana::LoadFrom<IPrediction> (dir_pred_no);
219  auto pred_sh = ana::LoadFrom<IPrediction> (dir_pred_sh);
220  ComparePredictions (pred_no.get(), pred_sh.get(),
221  plot_title, plot_name, "CVN",
222  printtable, pidaxes);
223 
224  //Additional XP plots
225  if (xp == "xp_prop") {
226  auto xp = (ModularExtrap*)
227  ((PredictionExtrap*)pred_sh.release())->GetExtrap();
228  xp->SavePlotsNue (outfile_xpplots->mkdir (
229  "plots_xp_" + weight.name), 9E20);
230 
231  if ("weight.name" == "no_weight") {
232  xp = (ModularExtrap*)
233  ((PredictionExtrap*)pred_no.release())->GetExtrap();
234  xp->SavePlotsNue (outfile_xpplots->mkdir (
235  "plots_xp_ppfx_xsec"), 9E20);
236  } //end no weight
237  } //end XP plots
238 
239  } //end xp
240  } //end weights
241 
242  TDirectory* dir_pred_pl, * dir_pred_mi;
243  TDirectory* dir_decomp_pl, * dir_decomp_mi;
244 
245  // Compare plus and minus shifts to nominal at once
246  for (TString shiftname: {"MEC_enushape", "MEC_q0shape",
247  "MEC_NPfrac", "MEC_scaleSA",
248  "MEC_all3"}) {
249 
250  for (TString xp: {"nxp","xp_prop"}) {
251 
252  auto dir_pred_no = file->GetDirectory ("pred_" + xp +
253  "_no_shift_ppfx_xsec");
254  auto dir_pred_pl = file->GetDirectory ("pred_" + xp +
255  "_" + shiftname +
256  "_plus_ppfx_xsec");
257  auto dir_pred_mi = file->GetDirectory ("pred_" + xp +
258  "_" + shiftname +
259  "_minus_ppfx_xsec");
260 
261  // ND plots
262  if (xp == "xp_prop") {
263  dirname_decomp = "extrap/MEextrap/Decomp/";
264  plot_title = shiftname + " (NearDet MC)" +
265  "; Non-Quasielastic Energy Estimator (GeV); Events / Bin";
266  plot_name = "compare_numu_ND_" + shiftname;
267 
268  dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
269  dir_decomp_pl = dir_pred_pl->GetDirectory (dirname_decomp);
270  dir_decomp_mi = dir_pred_mi->GetDirectory (dirname_decomp);
271 
272  CompareNDDataMC (dir_decomp_no, dir_decomp_pl, dir_decomp_mi,
273  plot_title, plot_name, "Numu",
274  printtable_numu);
275 
276 
277  dirname_decomp = "extrap/EEextrap/Decomp/";
278  plot_title = shiftname + " (NearDet MC)" +
279  "; Neutrino Energy Estimator (GeV) X CVN; Events / Bin";
280  plot_name = "compare_nue_ND_" + shiftname;
281 
282  dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
283  dir_decomp_pl = dir_pred_pl->GetDirectory (dirname_decomp);
284  dir_decomp_mi = dir_pred_mi->GetDirectory (dirname_decomp);
285 
286  CompareNDDataMC (dir_decomp_no, dir_decomp_pl, dir_decomp_mi,
287  plot_title, plot_name, "CVN",
288  printtable_nue, pidaxes);
289  }
290 
291  // FD Plots
292  plot_title = shiftname +
293  (xp != "nxp" ? " FarDet Extrap" : " (FarDet MC)") +
294  ";Neutrino Energy Estimator (GeV) X CVN;"+ "Events / Bin";
295  plot_name = "compare_" + xp + "_" + shiftname;
296 
297  auto pred_no = ana::LoadFrom <IPrediction> (dir_pred_no);
298  auto pred_pl = ana::LoadFrom <IPrediction> (dir_pred_pl);
299  auto pred_mi = ana::LoadFrom <IPrediction> (dir_pred_mi);
300  ComparePredictions (pred_no.get(), pred_pl.get(), pred_mi.get(),
301  plot_title, plot_name, "CVN",
302  printtable, pidaxes);
303 
304  } // end xp
305  } // end weights
306  } // end plots
307 }
Near Detector underground.
Definition: SREnums.h:10
void xsec_extrap_plots(bool createFile=false, TString suffix="test")
const XML_Char * name
Definition: expat.h:151
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
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
string filename
Definition: shutoffs.py:106
const NOvARwgtSyst kMECScaleSystSA("mecScale","MEC Scale", novarwgt::kMECScaleSystSA)
Definition: MECSysts.h:48
Generates FD-only predictions (no extrapolation)
const NOvARwgtSyst kMECEnuShapeSyst2017("MECEnuShape","MEC E_{#nu} shape", novarwgt::kMECEnuShapeSyst2017)
Definition: MECSysts.h:43
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
Generates extrapolated Nue predictions using ProportionalDecomp.
void CompareNDDataMC(TDirectory *d_no, TDirectory *d_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=0)
const Cut kNumuND
Definition: NumuCuts.h:55
Var weights
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:18
double sigma(TH1F *hist, double percentile)
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const HistAxis axisNue
Definition: syst_header.h:378
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
Definition: MECSysts.h:41
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void GetExtrap(PredictionNoExtrap *all, PredictionNoExtrap *sel, osc::IOscCalc *calc, TString title)
TFile * file
Definition: cellShifts.C:17
void SetLoaderPaths2017(Loaders &loaders, const bool usedecafs=false)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const HistAxis axisNumu
Take the output of an extrapolation and oscillate it as required.
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
void ComparePredictions(IPrediction *pred_no, IPrediction *pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=false)
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
const NOvARwgtSyst kMECInitStateNPFracSyst2017("MECInitStateNPFrac","MEC initial state np fraction", novarwgt::kMECInitStateNPFracSyst2017)
Definition: MECSysts.h:42
enum BeamMode string