make_nueFDprediction_kinematics_FHC_REW.C
Go to the documentation of this file.
1 
2 // Analysis
5 #include "CAFAna/Analysis/SALoaders.h"
7 
8 // Core
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Core/Cut.h"
11 #include "CAFAna/Core/HistAxis.h"
13 #include "CAFAna/Core/Spectrum.h"
16 
17 // Cuts
18 #include "CAFAna/Cuts/Cuts.h"
21 #include "CAFAna/Cuts/NueCutsSecondAna.h"
22 #include "CAFAna/Cuts/SpillCuts.h"
25 
26 
27 // Decomp
32 
33 // Extrap
35 
36 // Predict
40 
41 // Systs
42 #include "CAFAna/Systs/Systs.h"
43 
44 // Vars
48 #include "CAFAna/Vars/Vars.h"
51 #include "CAFAna/Vars/HistAxes.h"
52 
53 #include "OscLib/OscCalcDumb.h"
54 
56 
57 
58 #include <iostream>
59 #include <iomanip>
60 
61 using namespace ana;
62 
63 void CheckFileOverwrite(TString);
64 
66 
67 void make_nueFDprediction_kinematics_FHC_REW(const std::string sample, const std::string wfile, const std::string& outfilename = "FDprediction_kinematics_REW.root")
68 
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 ********
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] = {
105  };
106  const std::string selNames[kNumSels] = {
107  sample
108  };
109 
110  const Cut cutNDNumu = {
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 
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);
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];
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
300 
301 //****************************************
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 }
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
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Color_t kMC
std::string name
Definition: NuePlotLists.h:12
const Cut kNue2018CVNCut([](const caf::SRProxy *sr){if(kIsRHC(sr)) return kNue2018RHCCVNCut(sr);else return kNue2018FHCCVNCut(sr);})
Definition: NueCuts2018.h:49
const Cut kNue2018NDPresel
Definition: NueCuts2018.h:145
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 HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
Definition: HistAxes.h:8
osc::OscCalcDumb calc
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 Binning kNue2018Binning
Definition: NueCuts2018.cxx:93
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:22
OStream cout
Definition: OStream.cxx:6
const HistAxis axisNue
Definition: syst_header.h:378
void make_nueFDprediction_kinematics_FHC_REW(const std::string sample, const std::string wfile, const std::string &outfilename="FDprediction_kinematics_REW.root")
const Binning bins
Definition: NumuCC_CPiBin.h:8
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.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void CheckFileOverwrite(TString)
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
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
const std::string selNames[kNumSels]
Definition: vars.h:46
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
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 HistAxis kNue2018Axis("NuE Energy / Analysis Bin", kNue2018Binning, kNue2018AnaBin)
Use this Axis for Ana2018, official Axis.
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
enum BeamMode string