Functions | Variables
make_nueFDprediction_kinematics_FHC_REW.C File Reference
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/SALoaders.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "3FlavorAna/Decomp/BENDecomp.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "OscLib/OscCalcDumb.h"
#include "StandardRecord/Proxy/SRProxy.h"
#include <iostream>
#include <iomanip>

Go to the source code of this file.

Functions

void CheckFileOverwrite (TString)
 
void make_nueFDprediction_kinematics_FHC_REW (const std::string sample, const std::string wfile, const std::string &outfilename="FDprediction_kinematics_REW.root")
 

Variables

osc::OscCalcDumb calc
 

Function Documentation

void CheckFileOverwrite ( TString  outfilename)

Definition at line 302 of file make_nueFDprediction_kinematics_FHC_REW.C.

References om::cout.

Referenced by make_nueFDprediction_kinematics_FHC_REW().

302  {
303  if(!gSystem->AccessPathName(outfilename.Data())){
304  std::cout << "\n\nThis will overwrite " << outfilename
305  << "\n\nAre you sure you want to continue? y/n ";
306  string input = ""; getline(cin,input);
307  if (input!="y") abort();
308  }
309 }
make_template_knob_config.outfilename
string outfilename
knobs that need extra care
Definition: make_template_knob_config.py:120
om::cout
OStream cout
Definition: OStream.cxx:6
void make_nueFDprediction_kinematics_FHC_REW ( const std::string  sample,
const std::string  wfile,
const std::string &  outfilename = "FDprediction_kinematics_REW.root" 
)

Definition at line 67 of file make_nueFDprediction_kinematics_FHC_REW.C.

References ana::HistDef::axis, axisNue, axisNumu, ana::bins, CheckFileOverwrite(), om::cout, ana::defs, allTimeWatchdog::endl, file, make_syst_table_plots::h, MECModelEnuComparisons::i, ana::kBeam, ana::kCosNumi, ana::kCosTheta, cmf::kData, caf::kFARDET, kMC, caf::kNEARDET, ana::kNoShift, ana::kNue2018Axis(), ana::kNue2018Binning, ana::kNue2018CVNCut, ana::kNue2018FDAllSamples, ana::kNue2018NDPresel, ana::kNumSels, ana::kNumuCCAxis, ana::kNumuCutND2018, ana::kNumVars, ana::kPPFXFluxCVWgt, ana::kPtP, ana::kRecoQ2, ana::kStandardSpillCuts, cmf::kTauSwap, ana::kTrueQ2, ana::kTrueW, ana::kXSecCVWgt2018, loaders, ana::HistDef::name, nbins, NueExtrap(), make_template_knob_config::outfilename, ana::PredictionExtrap::SaveTo(), ana::PredictionExtendToPeripheral::SaveTo(), ana::PredictionNoExtrap::SaveTo(), ana::selNames, ana::Binning::Simple(), sr, and string.

69 {
70  //********** 0. Basic Setup *********************
71 
72  // ND loaders
73  std::string loaderNDData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
74  std::string loaderNDMC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
75 
76  // FD loaders
77  std::string loaderFDnonswapMC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1";
78  std::string loaderFDfluxswapMC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_fluxswap_fhc_nova_v08_full_v1";
79  std::string loaderFDtauswapMC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_tau_fhc_nova_v08_full_v1";
80 
81  //******* 1. Set up selection and variables ********
82  const HistAxis axisNumu = kNumuCCAxis;
83 
84  const int kNumSels = 1;
85  const int kNumVars = 5;
86 
87  const HistDef defs[kNumVars] = {
88  {"cvn2d", kNue2018Axis},
89  {"trueQ2", {"true Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kTrueQ2}},
90  {"trueW2", {"true W^{2} (GeV^{2})", Binning::Simple(100, 0, 2), kTrueW}},
91  {"PtP", {"p_{t}/p", Binning::Simple(500,0,1), kPtP}},
92  {"CosTheta", {"CosTheta", Binning::Simple(500,0,1), kCosTheta}}
93  };
94 
95  Cut nueFDsel = kNue2018FDAllSamples;
96 
97  /*if (sample == "CORE") nueFDsel = kNue2018FD;
98  else if (sample == "PERI") nueFDsel = kNue2018FDPeripheral;*/
99 
100  const Cut selFD[kNumSels] = {
101  nueFDsel
102  };
103  const Cut selND[kNumSels] = {
104  kNue2018NDPresel && kNue2018CVNCut
105  };
106  const std::string selNames[kNumSels] = {
107  sample
108  };
109 
110  const Cut cutNDNumu = {
111  kNumuCutND2018
112  };
113 
114  // open the weights file
115  TFile *weight_file = new TFile(wfile.c_str(),"READ");
116 
117  // construction of trueQ2 weight
118  TH1 *h;
119  gDirectory->GetObject(("trueQ2_weight_"+sample).c_str(),h);
120 
121  std::cout << h->GetNbinsX() << endl;
122  int nbins = h->GetNbinsX();
123  double bins[nbins], binslow[nbins], binshigh[nbins];
124  for(int i = 1; i <= nbins; ++i){
125  bins[i-1] = h->GetBinContent(i);
126  binslow[i-1] = h->GetBinLowEdge(i);
127  binshigh[i-1] = h->GetBinLowEdge(i) + h->GetBinWidth(i);
128  }
129 
130  const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](const caf::SRProxy *sr)
131  {
132  for(int i = 0; i < nbins; ++i){
133  if((kRecoQ2(sr) > binslow[i]) && (kRecoQ2(sr) < binshigh[i])){return bins[i];}
134  }
135  return 1.0;
136  }
137  );
138 
139  //construction of PtP weight
140  TH1 *hptp;
141  gDirectory->GetObject(("PtP_weight_"+sample).c_str(),hptp);
142 
143  std::cout << hptp->GetNbinsX() << endl;
144  int nbins_ptp = hptp->GetNbinsX();
145  double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp];
146  for(int i = 1; i <= nbins_ptp; ++i){
147  bins_ptp[i-1] = hptp->GetBinContent(i);
148  binslow_ptp[i-1] = hptp->GetBinLowEdge(i);
149  binshigh_ptp[i-1] = hptp->GetBinLowEdge(i) + hptp->GetBinWidth(i);
150  }
151 
152  const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](const caf::SRProxy *sr)
153  {
154  for(int i = 0; i < nbins_ptp; ++i){
155  if((kPtP(sr) > binslow_ptp[i]) && (kPtP(sr) < binshigh_ptp[i])){return bins_ptp[i];}
156  }
157  return 1.0;
158  }
159  );
160 
161  //construction of Cos weight
162  TH1 *hcos;
163  gDirectory->GetObject(("CosNumi_weight_"+sample).c_str(),hcos);
164 
165  std::cout << hcos->GetNbinsX() << endl;
166  int nbins_cos = hptp->GetNbinsX();
167  double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos];
168  for(int i = 1; i <= nbins_cos; ++i){
169  bins_cos[i-1] = hcos->GetBinContent(i);
170  binslow_cos[i-1] = hcos->GetBinLowEdge(i);
171  binshigh_cos[i-1] = hcos->GetBinLowEdge(i) + hcos->GetBinWidth(i);
172  }
173 
174  const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](const caf::SRProxy *sr)
175  {
176  for(int i = 0; i < nbins_cos; ++i){
177  if((kCosNumi(sr) > binslow_cos[i]) && (kCosNumi(sr) < binshigh_cos[i])){return bins_cos[i];}
178  }
179  return 1.0;
180  }
181  );
182 
183  // close weights file
184  weight_file->Close();
185 
186  CheckFileOverwrite(outfilename.c_str());
187 
188  //********** 2. Set up loaders *********************
189  SANueExtrapLoaders loaders;
190 
191  loaders.SetLoaderPath(loaderNDData, caf::kNEARDET, ana::Loaders::DataMC::kData);
192  loaders.SetLoaderPath(loaderNDMC, caf::kNEARDET, ana::Loaders::DataMC::kMC);
193 
194  loaders.SetLoaderPath(loaderFDnonswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam, ana::Loaders::SwappingConfig::kNonSwap);
195  loaders.SetLoaderPath(loaderFDfluxswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam,ana::Loaders::SwappingConfig::kFluxSwap);
196  loaders.SetLoaderPath(loaderFDtauswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam, ana::Loaders::SwappingConfig::kTauSwap);
197 
198  loaders.SetSpillCut(kStandardSpillCuts);
199 
200  //******* 3. Set up the extrapolations and go *******
201  SystShifts thisShift = kNoShift;
202 
203  // reweighted numudecomps
204  NumuDecomp * numuDecomp_trueQ2;
205  NumuDecomp * numuDecomp_PtP;
206  NumuDecomp * numuDecomp_Cos;
207 
208  // propdecomp
209  ProportionalDecomp * propDecomp[kNumSels][kNumVars];
210 
211  // extraps from different numudecomps
212  ModularExtrap * extrapProp_trueQ2[kNumSels][kNumVars];
213  ModularExtrap * extrapProp_PtP[kNumSels][kNumVars];
214  ModularExtrap * extrapProp_Cos[kNumSels][kNumVars];
215 
216  PredictionExtrap * predicProp_trueQ2[kNumSels][kNumVars];
217  PredictionExtrap * predicProp_PtP[kNumSels][kNumVars];
218  PredictionExtrap * predicProp_Cos[kNumSels][kNumVars];
219 
220  PredictionNoExtrap * predicNoXP[kNumSels][kNumVars];
221 
222 
223 
224  numuDecomp_trueQ2 = new NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kTrueQ2Weight);
225  numuDecomp_PtP = new NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kPtPWeight);
226  numuDecomp_Cos = new NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kCosWeight);
227 
228  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
229  for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
230  const HistAxis& axisNue = defs[varIdx].axis;
231 
232  propDecomp[selIdx][varIdx] = new ProportionalDecomp (loaders,
233  axisNue, selND[selIdx],
234  thisShift, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt );
235 
236  // trueQ2 reweighted extrap
237  extrapProp_trueQ2[selIdx][varIdx] = new NueExtrap (loaders,
238  *propDecomp[selIdx][varIdx], *numuDecomp_trueQ2,
239  axisNue, axisNumu,
240  selFD[selIdx], selND[selIdx], cutNDNumu,
241  thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
242 
243  // PtP reweighted extrap
244  extrapProp_PtP[selIdx][varIdx] = new NueExtrap (loaders,
245  *propDecomp[selIdx][varIdx], *numuDecomp_PtP,
246  axisNue, axisNumu,
247  selFD[selIdx], selND[selIdx], cutNDNumu,
248  thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt));
249  // Cos reweighted extrap
250  extrapProp_Cos[selIdx][varIdx] = new NueExtrap (loaders,
251  *propDecomp[selIdx][varIdx], *numuDecomp_Cos,
252  axisNue, axisNumu,
253  selFD[selIdx], selND[selIdx], cutNDNumu,
254  thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt));
255 
256 
257  predicProp_trueQ2[selIdx][varIdx] = new PredictionExtrap (extrapProp_trueQ2[selIdx][varIdx]);
258  predicProp_PtP[selIdx][varIdx] = new PredictionExtrap (extrapProp_PtP[selIdx][varIdx]);
259  predicProp_Cos[selIdx][varIdx] = new PredictionExtrap (extrapProp_Cos[selIdx][varIdx]);
260 
261  predicNoXP[selIdx][varIdx] = new PredictionNoExtrap (loaders,
262  axisNue,
263  selFD[selIdx], thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
264 
265  }//vars
266  }//sels
267 
268  loaders.Go();
269 
270  PredictionExtendToPeripheral * predicExte_trueQ2[kNumSels][kNumVars];
271  PredictionExtendToPeripheral * predicExte_PtP[kNumSels][kNumVars];
272  PredictionExtendToPeripheral * predicExte_Cos[kNumSels][kNumVars];
273 
274  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
275  for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
276  predicExte_trueQ2[selIdx][varIdx] = new PredictionExtendToPeripheral (predicProp_trueQ2[selIdx][varIdx], predicNoXP[selIdx][varIdx], kNue2018Binning, false);
277  predicExte_PtP[selIdx][varIdx] = new PredictionExtendToPeripheral (predicProp_PtP[selIdx][varIdx], predicNoXP[selIdx][varIdx], kNue2018Binning, false);
278  predicExte_Cos[selIdx][varIdx] = new PredictionExtendToPeripheral (predicProp_Cos[selIdx][varIdx], predicNoXP[selIdx][varIdx], kNue2018Binning, false);
279  }
280  }
281 
282  //********** 4. Save everything ********************
283  TFile * file = new TFile(outfilename.c_str(),"RECREATE");
284  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
285  TDirectory* dFDshi = file->mkdir((selNames[selIdx]+"_"+sample).c_str());
286  for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
287  TString varName = defs[varIdx].name.c_str();
288  predicExte_trueQ2[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_trueQ2_"+selNames[selIdx]+"_"+varName) ;
289  predicExte_PtP[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_PtP_"+selNames[selIdx]+"_"+varName) ;
290  predicExte_Cos[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_Cos_"+selNames[selIdx]+"_"+varName) ;
291  predicProp_trueQ2[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextend_trueQ2_"+selNames[selIdx]+"_"+varName) ;
292  predicProp_PtP[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextend_PtP_"+selNames[selIdx]+"_"+varName) ;
293  predicProp_Cos[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextend_Cos_"+selNames[selIdx]+"_"+varName) ;
294  predicNoXP[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextrap_"+selNames[selIdx]+"_"+varName) ;
295  }
296  }
297  file->Close();
298 
299 } // end
caf::kNEARDET
Near Detector underground.
Definition: SREnums.h:10
ana::kRecoQ2
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
ana::kNumVars
const int kNumVars
Definition: vars.h:14
caf::kFARDET
Far Detector at Ash River.
Definition: SREnums.h:11
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
ana::kNumuCutND2018
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
kMC
const Color_t kMC
Definition: make_extrap_figure.C:44
ana::HistDef::name
std::string name
Definition: NuePlotLists.h:12
ana::kNue2018CVNCut
const Cut kNue2018CVNCut([](const caf::SRProxy *sr){if(kIsRHC(sr)) return kNue2018RHCCVNCut(sr);else return kNue2018FHCCVNCut(sr);})
Definition: NueCuts2018.h:49
ana::kNue2018NDPresel
const Cut kNue2018NDPresel
Definition: NueCuts2018.h:145
ana::kTrueQ2
const Var kTrueQ2
Definition: TruthVars.h:27
make_template_knob_config.outfilename
string outfilename
knobs that need extra care
Definition: make_template_knob_config.py:120
nbins
const int nbins
Definition: cellShifts.C:15
ana::kPtP
const Var kPtP
Transverse momentum fraction in slice.
Definition: NueVars.cxx:90
ana::NumuDecomp
Uses MC for NC and CC components, assigns remainder of data to CC.
Definition: NumuDecomp.h:10
ana::kNumuCCAxis
const HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
Definition: HistAxes.h:8
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
bins
const Binning bins
Definition: nus17_fd_arrays.h:114
ana::defs
const HistDef defs[kNumVars]
Definition: vars.h:15
ana::kNue2018Binning
const Binning kNue2018Binning
Definition: NueCuts2018.cxx:93
ana::kNumSels
const int kNumSels
Definition: vars.h:43
ana::PredictionNoExtrap::SaveTo
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: PredictionNoExtrap.cxx:79
ana::HistDef::axis
HistAxis axis
Definition: NuePlotLists.h:13
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
om::cout
OStream cout
Definition: OStream.cxx:6
axisNue
const HistAxis axisNue
Definition: syst_header.h:378
NueExtrap
void NueExtrap(string beam="fhc", string cvntype="oldpresel")
Definition: NueExtrap.C:28
ana::PredictionExtrap::SaveTo
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: PredictionExtrap.cxx:139
ana::ProportionalDecomp
Splits Data proportionally according to MC.
Definition: ProportionalDecomp.h:20
ana::_HistAxis< Var >
ana::kXSecCVWgt2018
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
loaders
std::vector< Loaders * > loaders
Definition: syst_header.h:386
CheckFileOverwrite
void CheckFileOverwrite(TString)
Definition: make_nueFDprediction_kinematics_FHC_REW.C:302
ana::kNue2018FDAllSamples
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
file
TFile * file
Definition: cellShifts.C:17
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
ana::PredictionNoExtrap
Prediction that just uses FD MC, with no extrapolation.
Definition: PredictionNoExtrap.h:17
axisNumu
const HistAxis axisNumu
Definition: syst_variations.h:206
ana::PredictionExtrap
Take the output of an extrapolation and oscillate it as required.
Definition: PredictionExtrap.h:12
ana::_Cut< caf::SRProxy >
ana::kCosTheta
const Var kCosTheta
Definition: NueVarsExtra.cxx:240
ana::selNames
const std::string selNames[kNumSels]
Definition: vars.h:46
ana::_Var< caf::SRProxy >
ana::ModularExtrap
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
ana::PredictionExtendToPeripheral::SaveTo
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: PredictionExtendToPeripheral.cxx:420
ana::kTrueW
const Var kTrueW
Definition: TruthVars.h:22
ana::kNue2018Axis
const HistAxis kNue2018Axis("NuE Energy / Analysis Bin", kNue2018Binning, kNue2018AnaBin)
Use this Axis for Ana2018, official Axis.
ana::kCosNumi
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
string
enum BeamMode string

Variable Documentation

osc::OscCalcDumb calc

Definition at line 65 of file make_nueFDprediction_kinematics_FHC_REW.C.

Generated on Fri Apr 9 2021 14:13:56 for NOvA by   doxygen 1.8.11