datamc_ND_numu_kinematics_RHC_REW_pTBins.C
Go to the documentation of this file.
1 #include "CAFAna/Analysis/Plots.h"
2 #include "CAFAna/Prediction/PredictionNoExtrap.h"
3 #include "StandardRecord/Proxy/SRProxy.h"
4 
5 // Core
6 #include "CAFAna/Core/EventList.h"
7 #include "CAFAna/Core/SpectrumLoader.h"
8 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Core/SpectrumLoaderBase.h"
10 #include "CAFAna/Core/HistAxis.h"
11 #include "CAFAna/Core/Binning.h"
12 #include "CAFAna/Core/Cut.h"
13 
14 // Cuts
15 #include "CAFAna/Cuts/Cuts.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
17 #include "3FlavorAna/Cuts/NueCuts2020.h"
18 #include "3FlavorAna/Cuts/NumuCuts.h"
19 #include "3FlavorAna/Cuts/NumuCuts2020.h"
20 #include "CAFAna/Cuts/TruthCuts.h"
21 
22 // Vars
23 #include "CAFAna/Vars/Vars.h"
24 #include "CAFAna/Vars/XsecTunes.h"
25 #include "CAFAna/Vars/GenieWeights.h"
26 #include "3FlavorAna/Vars/NueVarsExtra.h"
27 #include "3FlavorAna/Vars/NueVars.h"
28 #include "3FlavorAna/Vars/NumuVars.h"
29 #include "CAFAna/Vars/PPFXWeights.h"
30 #include "CAFAna/Vars/HistAxes.h"
31 
32 // Systs
33 #include "CAFAna/Systs/Systs.h"
34 #include "CAFAna/Analysis/Style.h"
35 #include "CAFAna/Decomp/IDecomp.h"
36 #include "CAFAna/Decomp/ProportionalDecomp.h"
37 
38 #include "OscLib/OscCalcPMNSOpt.h"
39 
40 
41 #include "CAFAna/Core/SystShifts.h"
42 #include "3FlavorAna/Cuts/NumuCuts.h"
43 #include "3FlavorAna/Cuts/QuantileCuts2020.h"
44 #include "3FlavorAna/Decomp/NumuDecomp.h"
45 #include "CAFAna/Decomp/CheatDecomp.h"
46 #include "3FlavorAna/Decomp/MichelDecomp.h"
47 #include "CAFAna/Extrap/ModularExtrap.h"
48 #include "CAFAna/Prediction/PredictionExtrap.h"
49 #include "TFile.h"
50 #include "OscLib/OscCalcPMNSOpt.h"
51 #include "OscLib/OscCalcDumb.h"
52 #include "CAFAna/Core/Spectrum.h"
53 #include "TH1.h"
54 #include "CAFAna/Experiment/SingleSampleExperiment.h"
55 #include "TCanvas.h"
56 #include "CAFAna/Fit/FrequentistSurface.h"
57 #include "CAFAna/Analysis/Plots.h"
58 #include "3FlavorAna/Vars/HistAxes.h"
59 #include "CAFAna/Vars/FitVars.h"
60 #include "CAFAna/Cuts/Cuts.h"
61 #include "3FlavorAna/Cuts/NumuCuts.h"
62 
63 #include "TStopwatch.h"
64 #include "CAFAna/Core/Utilities.h"
65 
66 
67 
68 #include <fstream>
69 #include <iostream>
70 #include <cmath>
71 #include <string>
72 
73 using namespace ana;
74 
75 #include "OscLib/IOscCalc.h"
76 
77 // ROOT
78 #include "TCanvas.h"
79 #include "TFile.h"
80 #include "TH1.h"
81 #include "TLegend.h"
82 #include "TPad.h"
83 #include "TMath.h"
84 #include "TH2.h"
85 #include "TLine.h"
86 #include "TText.h"
87 #include "TString.h"
88 
89 #include "Utilities/func/MathUtil.h"
90 
91 //sample : ALL/CORE/PERI and the corresponding Weight file name.
92 void datamc_ND_numu_kinematics_RHC_REW_pTBins(const std::string sample, const std::string wfile)
93 {
94 
95  std::string ndData = "prod_caf_R19-11-18-prod5reco.g_nd_numi_fhc_full_v1_goodruns";
96  std::string ndMC = "prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1";
97 
98  SpectrumLoader loaderData(ndData);
99  SpectrumLoader loaderMC(ndMC);
100 
101  loaderData.SetSpillCut(kStandardSpillCuts);
102  loaderMC.SetSpillCut(kStandardDQCuts);
103 
104  //which variables do you want to look at?
105  const int kNumVars = 2;
106 
107  struct HistDef{
108  std::string name;
109  const HistAxis axis;
110  };
111 
112  const HistDef defs[kNumVars] = {
113  {"numu", kNumuCCAxis},
114  {"recoE", {"E_{reco} (GeV)", Binning::Simple(100, 0, 5), kCCE}},
115  {"trueQ2", {"true Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kTrueQ2}},
116  {"recoQ2", {"reco Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kRecoQ2}},
117  {"PtP", {"p_{t}/p", Binning::Simple(500,0,1), kPtP}},
118  {"CosNumi", {"CosNumi", Binning::Simple(500,0,1), kCosNumi}}
119  };
120 
121  const int kNumSels = 4;
122  const unsigned int nQuants = 3;
123  bool isRHC = true;
124  std::vector<Cut> PtQuantCutsND = GetNueQuantCuts2020( isRHC, caf::kNEARDET, nQuants, ana::kExtrapPt );
125  const Cut sels[kNumSels] = {
126  kNumu2020ND,
127  kNumu2020ND && PtQuantCutsND[0],
128  kNumu2020ND && PtQuantCutsND[1],
129  kNumu2020ND && PtQuantCutsND[2]
130  };
131  const std::string selNames[kNumSels] = {
132  (sample+"_AllBins").c_str(),
133  (sample+"_Quant1").c_str(),
134  (sample+"_Quant2").c_str(),
135  (sample+"_Quant3").c_str(),
136  };
137 
138  Spectrum* spects[kNumSels][kNumVars];
139  Spectrum* spects_trueQ2[kNumSels][kNumVars];
140  Spectrum* spects_PtP[kNumSels][kNumVars];
141  Spectrum* spects_Cos[kNumSels][kNumVars];
142  IPrediction* preds[kNumSels][kNumVars];
143  IPrediction* preds_trueQ2[kNumSels][kNumVars];
144  IPrediction* preds_PtP[kNumSels][kNumVars];
145  IPrediction* preds_Cos[kNumSels][kNumVars];
146 
147 
148  //For grid
149  // std::string wfile_name = ("/pnfs/nova/scratch/users/zvallari/AcceptSyst/"+wfile).c_str();
150  // TFile *weight_file = TFile::Open(pnfs2xrootd(wfile_name).c_str() );
151  //For non-grid
152  std::string wfile_name = ("/nova/ana/users/zvallari/AcceptSyst2020/"+wfile).c_str();
153  TFile *weight_file = TFile::Open(wfile_name.c_str());
154 
155  TH1 *h;
156  gDirectory->GetObject(("trueQ2_weight_"+sample).c_str(),h);
157 
158  std::cout << h->GetNbinsX() << endl;
159  int nbins = h->GetNbinsX();
160  double bins[nbins], binslow[nbins], binshigh[nbins];
161 
162  for(int i = 1; i <= nbins; ++i){
163  bins[i-1] = h->GetBinContent(i);
164  binslow[i-1] = h->GetBinLowEdge(i);
165  binshigh[i-1] = h->GetBinLowEdge(i) + h->GetBinWidth(i);
166  }
167 
168  const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](const caf::SRProxy *sr)
169  {
170  for(int i = 0; i < nbins; ++i){
171  if((kRecoQ2(sr) > binslow[i]) && (kRecoQ2(sr) < binshigh[i])){return bins[i];}
172  }
173  return 1.0;
174  }
175  );
176 
177  //PtP weight
178  TH1 *hptp;;
179  gDirectory->GetObject(("PtP_weight_"+sample).c_str(),hptp);
180 
181  std::cout << hptp->GetNbinsX() << endl;
182  int nbins_ptp = hptp->GetNbinsX();
183  double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp];
184 
185  for(int i = 1; i <= nbins_ptp; ++i){
186  bins_ptp[i-1] = hptp->GetBinContent(i);
187  binslow_ptp[i-1] = hptp->GetBinLowEdge(i);
188  binshigh_ptp[i-1] = hptp->GetBinLowEdge(i) + hptp->GetBinWidth(i);
189  }
190 
191  const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](const caf::SRProxy *sr)
192  {
193  for(int i = 0; i < nbins_ptp; ++i){
194  if((kPtP(sr) > binslow_ptp[i]) && (kPtP(sr) < binshigh_ptp[i])){return bins_ptp[i];}
195  }
196  return 1.0;
197  }
198  );
199 
200 
201  //Cos weight
202  TH1 *hcos;
203  gDirectory->GetObject(("CosNumi_weight_"+sample).c_str(),hcos);
204 
205  std::cout << hcos->GetNbinsX() << endl;
206  int nbins_cos = hptp->GetNbinsX();
207  double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos];
208 
209  for(int i = 1; i <= nbins_cos; ++i){
210  bins_cos[i-1] = hcos->GetBinContent(i);
211  binslow_cos[i-1] = hcos->GetBinLowEdge(i);
212  binshigh_cos[i-1] = hcos->GetBinLowEdge(i) + hcos->GetBinWidth(i);
213  }
214 
215  const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](const caf::SRProxy *sr)
216  {
217  for(int i = 0; i < nbins_cos; ++i){
218  if((kCosNumi(sr) > binslow_cos[i]) && (kCosNumi(sr) < binshigh_cos[i])){return bins_cos[i];}
219  }
220  return 1.0;
221  }
222  );
223 
224  weight_file->Close();
225 
226 
227 
228  //loop over the selectors and variables to create a full selection
229  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
230  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
231  const HistAxis& axis = defs[varIdx].axis;
232 
233  // non-reweighted
234  spects[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift);
235  preds[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
236  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
237  sels[selIdx], kNoShift, kXSecCVWgt2020*kPPFXFluxCVWgt);
238 
239  // trueQ2 reweighted
240  spects_trueQ2[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift, kTrueQ2Weight);
241  preds_trueQ2[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
242  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
243  sels[selIdx], kNoShift, kXSecCVWgt2020*kPPFXFluxCVWgt*kTrueQ2Weight);
244  // PtP reweighted
245  spects_PtP[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift, kPtPWeight);
246  preds_PtP[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
247  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
248  sels[selIdx], kNoShift, kXSecCVWgt2020*kPPFXFluxCVWgt*kPtPWeight);
249  // CosTheta reweighted
250  spects_Cos[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift, kCosWeight);
251  preds_Cos[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
252  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
253  sels[selIdx], kNoShift, kXSecCVWgt2020*kPPFXFluxCVWgt*kCosWeight);
254 
255  }
256  }
257 
258  loaderMC.Go();
259  loaderData.Go();
260 
261  TString fname = "datamc_ND_numu_kinematics_RHC_REW_pTBins.root";
262  TFile* fout = new TFile(fname, "RECREATE");
263  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
264  TDirectory* d = fout->mkdir(selNames[selIdx].c_str());
265  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
266  const char* name = defs[varIdx].name.c_str();
267  spects[selIdx][varIdx]->SaveTo(d, TString::Format("spect_noweight_%s", name));
268  preds[selIdx][varIdx]->SaveTo(d, TString::Format("pred_noweight_%s", name));
269 
270  spects_trueQ2[selIdx][varIdx]->SaveTo(d, TString::Format("spect_trueQ2_%s", name));
271  preds_trueQ2[selIdx][varIdx]->SaveTo(d, TString::Format("pred_trueQ2_%s", name));
272 
273  spects_PtP[selIdx][varIdx]->SaveTo(d, TString::Format("spect_PtP_%s", name));
274  preds_PtP[selIdx][varIdx]->SaveTo(d, TString::Format("pred_PtP_%s", name));
275 
276  spects_Cos[selIdx][varIdx]->SaveTo(d, TString::Format("spect_Cos_%s", name));
277  preds_Cos[selIdx][varIdx]->SaveTo(d, TString::Format("pred_Cos_%s", name));
278  }
279  }
280 }
281 
caf::kNEARDET
Near Detector underground.
Definition: SREnums.h:10
name
const XML_Char * name
Definition: expat.h:151
plot_validation_datamc.fname
fname
Definition: plot_validation_datamc.py:180
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
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
datamc_ND_numu_kinematics_RHC_REW_pTBins
void datamc_ND_numu_kinematics_RHC_REW_pTBins(const std::string sample, const std::string wfile)
Definition: datamc_ND_numu_kinematics_RHC_REW_pTBins.C:92
ana::_HistAxis::GetVars
const std::vector< T > & GetVars() const
Definition: HistAxis.h:92
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
make_mec_shifts_plots.isRHC
string isRHC
Definition: make_mec_shifts_plots.py:100
submit_syst.fout
fout
Definition: submit_syst.py:244
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
ana::kNumu2020ND
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
ana::HistDef::name
std::string name
Definition: NuePlotLists.h:12
ana::IPrediction::SaveTo
virtual void SaveTo(TDirectory *dir, const std::string &name) const
Definition: IPrediction.cxx:126
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
ana::kTrueQ2
const Var kTrueQ2
Definition: TruthVars.h:27
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
nbins
const int nbins
Definition: cellShifts.C:15
ana::kPtP
const Var kPtP
Transverse momentum fraction in slice.
Definition: NueVars.cxx:90
ana::sels
const Cut sels[kNumSels]
Definition: vars.h:44
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
ana::kStandardDQCuts
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:16
ana::kCCE
const Var kCCE
Definition: NumuVars.h:21
d
Float_t d
Definition: plot.C:236
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
ana::Spectrum::SaveTo
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
ana::defs
const HistDef defs[kNumVars]
Definition: vars.h:15
ana::GetNueQuantCuts2020
std::vector< Cut > GetNueQuantCuts2020(const bool isRHC, const caf::Det_t det, const unsigned int nquants, const ExtrapVar var)
Definition: QuantileCuts2020.cxx:90
ana::kNumSels
const int kNumSels
Definition: vars.h:43
cmf::Spectrum
std::vector< float > Spectrum
Definition: Constants.h:759
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
ana::LabelsAndBins::GetBinnings
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
ana::bins
const Binning bins
Definition: NumuCC_CPiBin.h:8
ana::_HistAxis< Var >
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::IPrediction
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
ana::PredictionNoExtrap
Prediction that just uses FD MC, with no extrapolation.
Definition: PredictionNoExtrap.h:17
ana::_Cut< caf::SRProxy >
ana::selNames
const std::string selNames[kNumSels]
Definition: vars.h:46
ana::_Var< caf::SRProxy >
ana::Binning::Simple
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
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
ana::kNullLoader
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
Definition: SpectrumLoaderBase.h:278
ana::LabelsAndBins::GetLabels
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
ana::kXSecCVWgt2020
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105
string
enum BeamMode string
Generated on Mon Jan 25 2021 14:10:52 for NOvA by   doxygen 1.8.11