datamc_ND_numu_kinematics_RHC_REW.C
Go to the documentation of this file.
1 
5 
6 // Core
9 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Core/HistAxis.h"
12 #include "CAFAna/Core/Binning.h"
13 #include "CAFAna/Core/Cut.h"
14 
15 // Cuts
16 #include "CAFAna/Cuts/Cuts.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
22 #include "CAFAna/Cuts/NumuCuts.h"
24 #include "CAFAna/Cuts/TruthCuts.h"
25 
26 // Vars
27 #include "CAFAna/Vars/Vars.h"
30 #include "CAFAna/Vars/NueVars.h"
31 #include "CAFAna/Vars/NumuVars.h"
33 #include "CAFAna/Vars/HistAxes.h"
34 
35 // Systs
36 #include "CAFAna/Systs/Systs.h"
37 #include "CAFAna/Analysis/Style.h"
38 #include "CAFAna/Decomp/IDecomp.h"
40 
42 
43 
44 #include "CAFAna/Core/SystShifts.h"
45 #include "CAFAna/Cuts/NumuCuts.h"
51 #include "TFile.h"
54 #include "CAFAna/Core/Spectrum.h"
55 #include "TH1.h"
57 #include "TCanvas.h"
59 #include "CAFAna/Analysis/Plots.h"
60 #include "CAFAna/Vars/FitVars.h"
61 #include "CAFAna/Cuts/Cuts.h"
62 #include "CAFAna/Cuts/NumuCuts.h"
63 
64 #include "TStopwatch.h"
65 #include "CAFAna/Core/Utilities.h"
66 
67 
68 
69 #include <fstream>
70 #include <iostream>
71 #include <cmath>
72 #include <string>
73 
74 using namespace ana;
75 
77 
78 // ROOT
79 #include "TCanvas.h"
80 #include "TFile.h"
81 #include "TH1.h"
82 #include "TLegend.h"
83 #include "TPad.h"
84 #include "TMath.h"
85 #include "TH2.h"
86 #include "TLine.h"
87 #include "TText.h"
88 #include "TString.h"
89 
90 #include "Utilities/func/MathUtil.h"
91 
92 
93 //define some useful variables:
94 
96 {
97 
98  std::string nameData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns";
99  std::string nameMC = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
100 
101  //choose your samweb datasets or concat files:
102  SpectrumLoader loaderData(nameData);
103  SpectrumLoader loaderMC(nameMC);
104 
105  loaderData.SetSpillCut(kStandardSpillCuts);
106  loaderMC.SetSpillCut(kStandardDQCuts);
107 
108  //which variables do you want to look at?
109  const int kNumVars = 8;
110 
111  const HistDef defs[kNumVars] = {
112  {"numu", kNumuCCAxis},
113  {"recoE", {"E_{reco} (GeV)", Binning::Simple(100, 0, 5), kCCE}},
114  {"trueQ2", {"true Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kTrueQ2}},
115  {"recoQ2", {"reco Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kRecoQ2}},
116  {"trueW2", {"true W^{2} (GeV^{2})", Binning::Simple(100, 0, 2), kTrueW}},
117  {"trueE", {"true E (GeV)", Binning::Simple(100, 0, 5), kTrueE}},
118  {"PtP", {"p_{t}/p", Binning::Simple(500,0,1), kPtP}},
119  {"CosNumi", {"CosNumi", Binning::Simple(500,0,1), kCosNumi}}
120  };
121 
122  const int kNumSels = 1;
123  const Cut sels[kNumSels] = {
125  };
126  const std::string selNames[kNumSels] = {
127  sample
128  };
129 
130  TFile *weight_file = new TFile(("/nova/app/users/nostom66/nue/ExtrapSyst/NDFDkinematics/RHC/Weights_"+sample+"_RHC.root").c_str(),"READ");
131 
132  //trueQ2 weight
133  TH1 *h;
134  gDirectory->GetObject(("trueQ2_weight_"+sample).c_str(),h);
135 
136  std::cout << h->GetNbinsX() << endl;
137  int nbins = h->GetNbinsX();
138  double bins[nbins], binslow[nbins], binshigh[nbins];
139 
140  for(int i = 1; i <= nbins; ++i){
141  bins[i-1] = h->GetBinContent(i);
142  binslow[i-1] = h->GetBinLowEdge(i);
143  binshigh[i-1] = h->GetBinLowEdge(i) + h->GetBinWidth(i);
144  }
145 
146  const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](const caf::SRProxy *sr)
147  {
148  for(int i = 0; i < nbins; ++i){
149  if((kRecoQ2(sr) > binslow[i]) && (kRecoQ2(sr) < binshigh[i])){return bins[i];}
150  }
151  return 1.0;
152  }
153  );
154 
155  //PtP weight
156  TH1 *hptp;;
157  gDirectory->GetObject(("PtP_weight_"+sample).c_str(),hptp);
158 
159  std::cout << hptp->GetNbinsX() << endl;
160  int nbins_ptp = hptp->GetNbinsX();
161  double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp];
162 
163  for(int i = 1; i <= nbins_ptp; ++i){
164  bins_ptp[i-1] = hptp->GetBinContent(i);
165  binslow_ptp[i-1] = hptp->GetBinLowEdge(i);
166  binshigh_ptp[i-1] = hptp->GetBinLowEdge(i) + hptp->GetBinWidth(i);
167  }
168 
169  const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](const caf::SRProxy *sr)
170  {
171  for(int i = 0; i < nbins_ptp; ++i){
172  if((kPtP(sr) > binslow_ptp[i]) && (kPtP(sr) < binshigh_ptp[i])){return bins_ptp[i];}
173  }
174  return 1.0;
175  }
176  );
177 
178  //Cos weight
179  TH1 *hcos;
180  gDirectory->GetObject(("CosNumi_weight_"+sample).c_str(),hcos);
181 
182  std::cout << hcos->GetNbinsX() << endl;
183  int nbins_cos = hptp->GetNbinsX();
184  double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos];
185 
186  for(int i = 1; i <= nbins_cos; ++i){
187  bins_cos[i-1] = hcos->GetBinContent(i);
188  binslow_cos[i-1] = hcos->GetBinLowEdge(i);
189  binshigh_cos[i-1] = hcos->GetBinLowEdge(i) + hcos->GetBinWidth(i);
190  }
191 
192  const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](const caf::SRProxy *sr)
193  {
194  for(int i = 0; i < nbins_cos; ++i){
195  if((kCosNumi(sr) > binslow_cos[i]) && (kCosNumi(sr) < binshigh_cos[i])){return bins_cos[i];}
196  }
197  return 1.0;
198  }
199  );
200 
201  weight_file->Close();
202 
204  Spectrum* spects_trueQ2[kNumSels][kNumVars];
205  Spectrum* spects_PtP[kNumSels][kNumVars];
206  Spectrum* spects_Cos[kNumSels][kNumVars];
207  IPrediction* preds[kNumSels][kNumVars];
208  IPrediction* preds_trueQ2[kNumSels][kNumVars];
209  IPrediction* preds_PtP[kNumSels][kNumVars];
210  IPrediction* preds_Cos[kNumSels][kNumVars];
211 
212  //loop over the selectors and variables to create a full selection
213  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
214  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
215  const HistAxis& axis = defs[varIdx].axis;
216 
217  // non-reweighted
218  spects[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift);
219  preds[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
220  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
221  sels[selIdx], kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
222 
223  // trueQ2 reweighted
224  spects_trueQ2[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift, kTrueQ2Weight);
225  preds_trueQ2[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
226  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
227  sels[selIdx], kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kTrueQ2Weight);
228  // PtP reweighted
229  spects_PtP[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift, kPtPWeight);
230  preds_PtP[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
231  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
232  sels[selIdx], kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kPtPWeight);
233  // CosTheta reweighted
234  spects_Cos[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift, kCosWeight);
235  preds_Cos[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
236  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
237  sels[selIdx], kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kCosWeight);
238 
239  }
240  }
241 
242  loaderMC.Go();
243  loaderData.Go();
244 
245  TString fname = "datamc_ND_numu_kinematics_REW.root";
246  TFile* fout = new TFile(fname, "RECREATE");
247  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
248  TDirectory* d = fout->mkdir(selNames[selIdx].c_str());
249  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
250  const char* name = defs[varIdx].name.c_str();
251  spects[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_noweight_%s", name)));
252  preds[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("pred_noweight_%s", name)));
253 
254  spects_trueQ2[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_trueQ2_%s", name)));
255  preds_trueQ2[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("pred_trueQ2_%s", name)));
256 
257  spects_PtP[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_PtP_%s", name)));
258  preds_PtP[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("pred_PtP_%s", name)));
259 
260  spects_Cos[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_Cos_%s", name)));
261  preds_Cos[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("pred_Cos_%s", name)));
262  }
263  }
264 }
265 
const XML_Char * name
Definition: expat.h:151
const int kNumVars
Definition: vars.h:14
Oscillation analysis framework, runs over CAF files outside of ART.
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
virtual void SaveTo(TDirectory *dir) const
std::vector< T > GetVars() const
Definition: HistAxis.h:79
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:40
HistAxis axis
Definition: NuePlotLists.h:13
Proxy for StandardRecord.
Definition: SRProxy.h:2237
void SetSpillCut(const SpillCut &cut)
const Var kTrueQ2
Definition: TruthVars.h:27
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const int nbins
Definition: cellShifts.C:15
const Var kPtP
Transverse momentum fraction in slice.
Definition: NueVars.cxx:90
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:87
const Cut sels[kNumSels]
Definition: vars.h:44
const HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
void datamc_ND_numu_kinematics_RHC_REW(const std::string sample)
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:12
const Var kCCE
Definition: Vars.h:90
Float_t d
Definition: plot.C:236
virtual void Go() override
Load all the registered spectra.
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
std::string name
Definition: NuePlotLists.h:12
const SystShifts kNoShift
Definition: SystShifts.h:112
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
const ana::Var kRecoQ2([](const caf::SRProxy *sr){const double M_mu_sqrd=util::sqr(0.1056);double E_mu=kNumuTrackE(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:128
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:50
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:28
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void spects(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC)
Definition: sensitivity.C:100
Prediction that just uses FD MC, with no extrapolation.
const std::string selNames[kNumSels]
Definition: vars.h:46
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Var kTrueW
Definition: TruthVars.h:22
std::vector< std::string > GetLabels() const
Definition: HistAxis.h:77
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:27
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
std::vector< Binning > GetBinnings() const
Definition: HistAxis.h:78
static constexpr Double_t sr
Definition: Munits.h:164
h
Definition: demo3.py:41
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1029