make_nueFDprediction_kinematics_RHC_REW.C
Go to the documentation of this file.
1 // Analysis
6 
7 // Core
8 #include "CAFAna/Core/Binning.h"
9 #include "CAFAna/Core/Cut.h"
10 #include "CAFAna/Core/HistAxis.h"
12 #include "CAFAna/Core/Spectrum.h"
15 
16 // Cuts
17 #include "CAFAna/Cuts/Cuts.h"
20 #include "CAFAna/Cuts/SpillCuts.h"
22 
23 // Decomp
28 
29 // Extrap
31 
32 // Predict
37 
38 // Systs
39 #include "CAFAna/Systs/Systs.h"
40 
41 // Vars
42 #include "CAFAna/Vars/XsecTunes.h"
43 #include "CAFAna/Vars/TruthVars.h"
46 #include "CAFAna/Vars/Vars.h"
49 
50 #include "OscLib/OscCalcDumb.h"
52 
53 #include <iostream>
54 #include <iomanip>
55 
56 using namespace ana;
57 
59 
60 //sample : ALL/CORE/PERI and the corresponding Weight file name.
62 
63 {
64  const std::string outfilename = ("FDprediction_kinematics_RHC_REW_"+sample+".root").c_str();
65 
66  //For grid
67  // std::string wfile_name = ("/pnfs/nova/scratch/users/zvallari/AcceptSyst/"+wfile).c_str();
68  // TFile *weight_file = TFile::Open(pnfs2xrootd(wfile_name).c_str() );
69 
70  //For non-grid
71  std::string wfile_name = ("/nova/ana/users/zvallari/AcceptSyst2020/"+wfile).c_str();
72  TFile *weight_file = TFile::Open(wfile_name.c_str());
73 
74  // construction of trueQ2 weight
75  TH1 *h;
76  gDirectory->GetObject(("trueQ2_weight_"+sample).c_str(),h);
77 
78  std::cout << h->GetNbinsX() << endl;
79  int nbins = h->GetNbinsX();
80  double bins[nbins], binslow[nbins], binshigh[nbins];
81  for(int i = 1; i <= nbins; ++i){
82  bins[i-1] = h->GetBinContent(i);
83  binslow[i-1] = h->GetBinLowEdge(i);
84  binshigh[i-1] = h->GetBinLowEdge(i) + h->GetBinWidth(i);
85  }
86 
87  const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](const caf::SRProxy *sr){
88  if(sr->hdr.det == caf::kFARDET) return 1.0;
89  if(!kNumu2020ND(sr)) {
90  std::cout << "THERE IS AN EVENT WHICH DOES NOT PASS NUMU ND CUT IN NUE SIGNAL EXTRAP!!!!" << std::endl;
91  exit(0);
92  return 1.0;
93  }
94  for(int i = 0; i < nbins; ++i){
95  if((kRecoQ2(sr) > binslow[i]) && (kRecoQ2(sr) < binshigh[i])){
96  return bins[i];
97  }
98  }
99  return 1.0;
100  }
101  );
102 
103  //construction of PtP weight
104  TH1 *hptp;
105  gDirectory->GetObject(("PtP_weight_"+sample).c_str(),hptp);
106 
107  std::cout << hptp->GetNbinsX() << endl;
108  int nbins_ptp = hptp->GetNbinsX();
109  double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp];
110  for(int i = 1; i <= nbins_ptp; ++i){
111  bins_ptp[i-1] = hptp->GetBinContent(i);
112  binslow_ptp[i-1] = hptp->GetBinLowEdge(i);
113  binshigh_ptp[i-1] = hptp->GetBinLowEdge(i) + hptp->GetBinWidth(i);
114  }
115 
116  const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](const caf::SRProxy *sr){
117  if(sr->hdr.det == caf::kFARDET) return 1.0;
118  for(int i = 0; i < nbins_ptp; ++i){
119  if((kPtP(sr) > binslow_ptp[i]) && (kPtP(sr) < binshigh_ptp[i])){return bins_ptp[i];}
120  }
121  return 1.0;
122  }
123  );
124 
125  //construction of Cos weight
126  TH1 *hcos;
127  gDirectory->GetObject(("CosNumi_weight_"+sample).c_str(),hcos);
128 
129  std::cout << hcos->GetNbinsX() << endl;
130  int nbins_cos = hptp->GetNbinsX();
131  double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos];
132  for(int i = 1; i <= nbins_cos; ++i){
133  bins_cos[i-1] = hcos->GetBinContent(i);
134  binslow_cos[i-1] = hcos->GetBinLowEdge(i);
135  binshigh_cos[i-1] = hcos->GetBinLowEdge(i) + hcos->GetBinWidth(i);
136  }
137 
138  const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](const caf::SRProxy *sr){
139  if(sr->hdr.det == caf::kFARDET) return 1.0;
140  for(int i = 0; i < nbins_cos; ++i){
141  if((kCosNumi(sr) > binslow_cos[i]) && (kCosNumi(sr) < binshigh_cos[i])){return bins_cos[i];}
142  }
143  return 1.0;
144  }
145  );
146 
147  // close weights file
148  weight_file->Close();
149 
152 
153  struct GenDef{
154  const IPredictionGenerator *gen;
155  const TString cutname;
156  const TString varname;
157  };
158 
159  std::vector <GenDef> gens;
160  std::vector <IPrediction * > predictions;
161 
163  sample.c_str(), "nueAxis_NoExtrap"});
164 
165 
169  kNumu2020ND,
170  kNoShift,
171  kXSecCVWgt2020*kPPFXFluxCVWgt*kTrueQ2Weight),
172  sample.c_str(), "nueAxis_NueSignalExtrap_Q2Weight"});
173 
177  kNumu2020ND,
178  kNoShift,
179  kXSecCVWgt2020*kPPFXFluxCVWgt*kCosWeight),
180  sample.c_str(), "nueAxis_NueSignalExtrap_CosWeight"});
181 
185  kNumu2020ND,
186  kNoShift,
187  kXSecCVWgt2020*kPPFXFluxCVWgt*kPtPWeight),
188  sample.c_str(), "nueAxis_NueSignalExtrap_PtPWeight"});
189 
190 
191  for(auto & gen:gens){
192  predictions.push_back(gen.gen->Generate(loaders).release());
193  }
194 
195  loaders.Go();
196 
197  TFile* file = new TFile(outfilename.c_str(),"RECREATE");
198  for (int i = 0 ; i < (int)gens.size(); ++i){
199  auto dir = file->GetDirectory(gens[i].cutname);
200  if(!dir) dir = file->mkdir(gens[i].cutname);
201  predictions[i]->SaveTo(dir, gens[i].varname);
202  }
203 
204  file->Close();
205 
206 }
207 
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
Far Detector at Ash River.
Definition: SREnums.h:11
std::map< TString, IPredictionGenerator * > gens
Definition: syst_header.h:387
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
const HistAxis kNumuCCOptimisedAxis2020("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kNumuE2020)
Definition: HistAxes.h:25
Generates FD-only predictions (no extrapolation)
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
string outfilename
knobs that need extra care
const int nbins
Definition: cellShifts.C:15
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod5Loaders.h:101
const HistAxis kNue2020Axis("NuE Energy / Analysis Bin", kNue2020Binning, kNue2020AnaBin)
Use this Axis for Ana2020, official Axis.
Definition: NueCuts2020.h:195
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
Generates extrapolated Nue signal-only predictions.
caf::StandardRecord * sr
osc::OscCalcDumb calc
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
TDirectory * dir
Definition: macro.C:5
std::vector< Loaders * > loaders
Definition: syst_header.h:386
exit(0)
TFile * file
Definition: cellShifts.C:17
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Given loaders and an MC shift, Generate() generates an IPrediction.
const Cut kNue2020FDAllSamples
Definition: NueCuts2020.h:84
void make_nueFDprediction_kinematics_RHC_REW(const std::string sample, const std::string wfile, const std::string &outfilename="FDprediction_kinematics_REW.root")
const Var kXSecCVWgt2020
Definition: XsecTunes.h:106
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
enum BeamMode string