NumuSysts.cxx
Go to the documentation of this file.
2 
7 
8 // Names for the GENIE systs
9 #include "nugen/NuReweight/ReweightLabels.h"
10 
12 
13 #include "TFile.h"
14 #include "TH1.h"
15 
16 #include <iostream>
17 #include <vector>
18 
19 namespace ana
20 {
24  kNumuSummedSmallGENIESyst(FindCAFAnaDir()+"/data/numu/summed_genie_systs_nd.root",
25  FindCAFAnaDir()+"/data/numu/summed_genie_systs_fd.root");
26 
33 
34 
35  //----------------------------------------------------------------------
37  Shift(double sigma, caf::SRProxy* /*sr*/, double& weight) const
38  {
39  weight *= 1+.01*sigma; // 1% shift
40  }
41 
42  //----------------------------------------------------------------------
44  Shift(double sigma, caf::SRProxy* sr, double& weight) const
45  {
46  const double scale = 1 + 0.01*sigma; //shifts 1% across boundary
47 
48  sr->energy.numu.trkccE *= scale;
49  }
50 
51  //----------------------------------------------------------------------
53  Shift(double sigma, caf::SRProxy* sr, double& weight) const
54  {
55  if(sr->mc.nnu == 0) return;
56  if(!sr->mc.nu[0].iscc) weight *= 1+sigma;
57  if(weight < 0) weight = 0;
58  }
59 
60  //----------------------------------------------------------------------
62  Shift(double sigma, caf::SRProxy* sr, double& weight) const
63  {
64  if(sr->mc.nnu == 0) return;
65  if(sr->mc.nu[0].iscc && abs(sr->mc.nu[0].pdg) == 16) weight *= 1+sigma;
66  }
67 
68  //----------------------------------------------------------------------
69  void NumuNormSyst::
70  Shift(double sigma, caf::SRProxy* /*sr*/, double& weight) const
71  {
72  weight *= 1+.01*sigma; // 1% absolute error
73  }
74 
75  //----------------------------------------------------------------------
77  Shift(double sigma, caf::SRProxy* sr, double& weight) const
78  {
79  // 1% difference between detectors + 0.7% from intensity + 3% ND norm (rock) + 3.7% (cosmic overlay) all in quadrature
80  if (sr->hdr.det == caf::kFARDET) return;
81  weight *= 1 + .049 * sigma;
82  }
83 
84  //----------------------------------------------------------------------
86  Shift(double sigma, caf::SRProxy* sr, double& weight) const
87  {
88  const double scale = 1 + .05*sigma;
89 
90  sr->slc.calE *= scale;
91  sr->energy.numu.trknonqeE *= scale;
92  sr->energy.numu.trkqeE *= scale;
93  sr->energy.numu.trkccE *= scale;
94  }
95 
96  //----------------------------------------------------------------------
98  Shift(double sigma, caf::SRProxy* sr, double& weight) const
99  {
100  if (sr->hdr.det != caf::kFARDET) return;
101 
102  const double scale = 1 + .025*sigma /// 5% difference between detectors
103  *( (sr->hdr.det==caf::kFARDET) ? 1 : -1 );
104 
105  sr->slc.calE *= scale;
106  sr->energy.numu.trknonqeE *= scale;
107  sr->energy.numu.trkqeE *= scale;
108  sr->energy.numu.trkccE *= scale;
109  }
110 
111 
112  //----------------------------------------------------------------------
113  // 2018 systematics list -- still a work in progress
114  std::vector<const ISyst*> getAllDirectNumuSysts2018() {
115  // This is the list of (mostly) *not-tiny* systematics for numu,
116  // for making things like the ND MC uncertainty bands and various cross checks.
117  // The latest greatest xsec systs are taken from getAllXsecSysts_2018()
118  // Beam stuff will need updating (not sure RHC is working???)
119  // The direct absolute hadE is no longer commented out
120  // (we will use the special calib files but that's not what this function is for)
121  // The really tiny GENIE syst that were included in the 2017 version have been removed
122 
123  std::vector<const ISyst*> systs =
124  {
130 
136 
140 
142 
149 
150  &kMuEScaleSyst2017, // correlated muon E scale applied to muon track length
151  &kRelMuEScaleSyst2017, // uncorrelated muon E scale applied to muon track length
152  &kDirectRelHadEScaleSyst2017, // relative near/far hadE syst
153  &kDirectHadEScaleSyst2017, // absolute hadE syst, correlated between detectors
154  };
155 
156  std::set<rwgt::ReweightLabel_t> largeGENIEParameters;
157  largeGENIEParameters.insert(rwgt::fReweightMaCCRES);
158  largeGENIEParameters.insert(rwgt::fReweightMvCCRES);
159  largeGENIEParameters.insert(rwgt::fReweightMaNCRES);
160 
161  for (const auto & parameter:largeGENIEParameters)
162  systs.push_back(GetGenieKnobSyst(parameter));
163 
164  return systs;
165  }
166 
167  //----------------------------------------------------------------------
168  // 2017 systematics list;
169  std::vector<const ISyst*> getAllDirectNumuSysts2017()
170  {
171  std::vector<const ISyst*> systs =
172  {
181  &kRPACCQEEnhSyst2017, // new RPA error (1-sided more enhancement)
182  &kRPACCQESuppSyst2017, // new RPA error (1-sided more suppression)
184  &kMuEScaleSyst2017, // correlated muon E scale applied to muon track length
185  &kRelMuEScaleSyst2017, // uncorrelated muon E scale applied to muon track length
186  &kDirectRelHadEScaleSyst2017, // directly applied 5% hadronic E uncorrelated syst
187  // kDirectHadEScaleSyst2017 will be provided elsewhere soon.
188  &kMAQEGenieReducedSyst2017, // replacing fReweightNormCCQE and fReweightMaCCQEshape from SA. Reduced M_A uncertainty (as compared to GENIE M_A knob)
189  };
190 
191  systs.push_back(&kNumuSummedSmallGENIESyst); // horrible and needs to be re-looked at for TA
192  std::set<rwgt::ReweightLabel_t> largeGENIEParameters;
193  largeGENIEParameters.insert(rwgt::fReweightMaNCEL);
194  largeGENIEParameters.insert(rwgt::fReweightMaCCRES);
195  largeGENIEParameters.insert(rwgt::fReweightMvCCRES);
196  largeGENIEParameters.insert(rwgt::fReweightMaNCRES);
197  largeGENIEParameters.insert(rwgt::fReweightMvNCRES);
198  largeGENIEParameters.insert(rwgt::fReweightCCQEPauliSupViaKF);
199 
200  for (const auto & parameter:largeGENIEParameters)
201  systs.push_back(GetGenieKnobSyst(parameter));
202 
203  return systs;
204  }
205 
206  //----------------------------------------------------------------------
209  const std::string& fnameFD)
210  : ISyst("numuSumSmallGENIE", "Summed Small GENIE")
211  {
212  fHistoFnames[0] = fnameND;
213  fHistoFnames[1] = fnameFD;
214 
215  // Signal that these aren't loaded yet
216  fHists[kFD][kSig] = 0;
217  }
218 
219 
220  //----------------------------------------------------------------------
222  {
223 
224  // Bail if we've already loaded histograms
225  if(fHists[kFD][kSig]) return;
226 
227 
228  for(int iDet = kND; iDet < kNumDets; ++iDet){
229  TFile fin(fHistoFnames[iDet].c_str(), "read");
230  if(fin.IsZombie()){
231  std::cerr << "Warning: couldn't open " <<
232  fHistoFnames[iDet] << std::endl;
233  return;
234  }
235 
236  fHists[iDet][kSig] = (TH1D*)fin.Get("allcc_sig")->Clone();
237  fHists[iDet][kSig]->SetDirectory(0);
238  fHists[iDet][kBkg] = (TH1D*)fin.Get("allcc_bkg")->Clone();
239  fHists[iDet][kBkg]->SetDirectory(0);
240  }
241  }
242 
243 
244  //----------------------------------------------------------------------
246  caf::SRProxy* sr, double& weight) const
247  {
249  /// Initialize the scale factor for the weight
250  double shift = 0;
251  if(sr->mc.nnu == 0) return;
252 
253  int det;
254  switch (sr->hdr.det){
255  case 1 : det = kND; break;
256  case 2 : det = kFD; break;
257  default: std::cerr <<"Wrong detector; ignore" << std::endl; return;
258  }
259 
260  TH1D* hist = NULL;
261  if(abs(sr->mc.nu[0].pdg) == 14 && sr->mc.nu[0].iscc)
262  hist = fHists[det][kSig];
263  else
264  hist = fHists[det][kBkg];
265 
266  double energy = sr->energy.numu.trkccE;
267 
268  if (energy > hist->GetXaxis()->GetXmin() &&
269  energy < hist->GetXaxis()->GetXmax() )
270  shift = hist->Interpolate(energy) - 1;
271  weight *= 1 + shift * sigma;
272 
273  }
274 
275  //----------------------------------------------------------------------
277  {
278  for(int i = 0; i < kNumDets; ++i){
279  for(int j = 0; j < kNumComponents; ++j){
280  delete fHists[i][j];
281  fHists[i][j] = 0;
282  }
283  }
284  }
285 
286 }
BeamSyst * GetPPFXPrincipals(int PCIdx)
Definition: BeamSysts.cxx:72
const NOvARwgtSyst kMECShapeSyst2018AntiNu("MECShape2018AntiNu","MEC 2018 (q_{0}, |#vec{q}|) response, antineutrinos", novarwgt::kMECQ0Q3RespSystNubar2018)
Definition: MECSysts.h:26
TString fin
Definition: Style.C:24
SummedSmallGENIE systematics, from histograms:
Definition: NumuSysts.h:40
const NumuNCScaleSyst kNumuNCScaleSyst
Definition: NumuSysts.cxx:27
caf::Proxy< float > trkccE
Definition: SRProxy.h:195
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:98
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const NOvARwgtSyst kRPACCQEEnhSyst2018("RPAShapeenh2018","RPA shape: higher-Q^{2} enhancement (2018)", novarwgt::kRPACCQEEnhSyst2018)
Definition: RPASysts.h:13
const NumuEnergyScaleSyst kNumuEnergyScaleSyst
Definition: NumuSysts.cxx:31
const NumuTauContaminationSyst kNumuTauContaminationSyst
Definition: NumuSysts.cxx:28
const NumuNormSyst kNumuNormSyst
Definition: NumuSysts.cxx:29
caf::Proxy< float > trknonqeE
Definition: SRProxy.h:196
const Var weight
GeniePCASyst * GetGeniePrincipals2018(int PCIdx)
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
const NOvARwgtSyst kMAQEGenieReducedSyst2018(genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced_2018", genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced_2018", novarwgt::kMAQEGenieReducedSyst2018)
2018 &#39;reduced&#39; M_A^QE shift. See documentation in NOvARwgt
Definition: XSecSysts.h:51
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:214
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const NOvARwgtSyst kMECEnuShapeSyst2018AntiNu("MECEnuShapeAntiNu","MEC E_{#nu} shape, antineutrinos", novarwgt::kMECEnuShapeSyst2017_NubarOnly)
Definition: MECSysts.h:32
const NOvARwgtSyst kMAQEGenieReducedSyst2017(genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced", genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced", novarwgt::kMAQEGenieReducedSyst2017)
2017 &#39;reduced&#39; M_A^QE shift. See documentation in NOvARwgt.
Definition: XSecSysts.h:47
const NOvARwgtSyst kMECShapeSyst2018Nu("MECShape2018Nu","MEC 2018 (q_{0}, |#vec{q}|) response, neutrinos", novarwgt::kMECQ0Q3RespSystNu2018)
Definition: MECSysts.h:25
OStream cerr
Definition: OStream.cxx:7
const NOvARwgtSyst kMECInitStateNPFracSyst2018AntiNu("MECInitStateNPFracAntiNu","MEC initial state np fraction, antineutrinos", novarwgt::kMECInitStateNPFracSyst2017_NubarOnly)
Definition: MECSysts.h:35
void abs(TH1 *hist)
const NumuRelNormSyst kNumuRelNormSyst
Definition: NumuSysts.cxx:30
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const NumuGEANTNormSyst kNumuGEANTNormSyst
Definition: NumuSysts.cxx:21
std::string fHistoFnames[kNumDets]
Definition: NumuSysts.h:56
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2136
const NOvARwgtSyst kMECEnuShapeSyst2018Nu("MECEnuShapeNu","MEC E_{#nu} shape, neutrinos", novarwgt::kMECEnuShapeSyst2017_NuOnly)
Definition: MECSysts.h:31
std::string FindCAFAnaDir()
Definition: Utilities.cxx:203
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const NOvARwgtSyst kMECEnuShapeSyst2017("MECEnuShape","MEC E_{#nu} shape", novarwgt::kMECEnuShapeSyst2017)
Definition: MECSysts.h:43
BeamSyst * GetFluxPrincipals2018(int PCIdx)
Definition: BeamSysts.cxx:89
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
std::vector< const ISyst * > getAllDirectNumuSysts2018()
Definition: NumuSysts.cxx:114
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:70
const NOvARwgtSyst kRPACCQEEnhSyst2017("RPAShapeenh2017","RPA shape: higher-Q^{2} enhancement (2017)", novarwgt::kRPACCQEEnhSyst2017)
Definition: RPASysts.h:11
Double_t scale
Definition: plot.C:25
const NumuGEANTScaleSyst kNumuGEANTScaleSyst
Definition: NumuSysts.cxx:22
caf::Proxy< float > trkqeE
Definition: SRProxy.h:197
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:77
const NOvARwgtSyst kMECInitStateNPFracSyst2018Nu("MECInitStateNPFracNu","MEC initial state np fraction, neutrinos", novarwgt::kMECInitStateNPFracSyst2017_NuOnly)
Definition: MECSysts.h:34
double energy
Definition: plottest35.C:25
const BeamSyst kBeamAllTransport((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"totErr","beamTransportComb","Combined Beam Transport Systematics")
All Beam Transport systematics combined in quadratures.
Definition: BeamSysts.h:135
const NOvARwgtSyst kRPACCQESuppSyst2017("RPAShapesupp2017","RPA shape: low-Q^{2} suppression (2017)", novarwgt::kRPACCQESuppSyst2017)
Definition: RPASysts.h:12
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:86
caf::StandardRecord * sr
const DirectHadEScaleSyst2017 kDirectHadEScaleSyst2017(0.05)
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:245
const double j
Definition: BetheBloch.cxx:29
double sigma(TH1F *hist, double percentile)
Energy scale component of the numu GEANT systematic, see docdb 13539.
Definition: NumuSysts.h:27
const NumuRelEnergyScaleSyst kNumukRelEnergyScaleSyst
Definition: NumuSysts.cxx:32
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const DirectRelHadEScaleSyst2017 kDirectRelHadEScaleSyst2017(0.05)
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
caf::Proxy< float > calE
Definition: SRProxy.h:1292
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
Definition: MECSysts.h:41
NumuSummedSmallGENIESyst(const std::string &fnameND, const std::string &fnameFD)
Definition: NumuSysts.cxx:208
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:37
const NumuSummedSmallGENIESyst kNumuSummedSmallGENIESyst(FindCAFAnaDir()+"/data/numu/summed_genie_systs_nd.root", FindCAFAnaDir()+"/data/numu/summed_genie_systs_fd.root")
Definition: NumuSysts.h:62
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:53
const NOvARwgtSyst kRPARESSyst2017("RPAShapeRES2017","RPA shape: resonance events", novarwgt::kRPARESSyst2017)
Definition: RPASysts.h:20
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:62
const NOvARwgtSyst kRPACCQESuppSyst2018("RPAShapesupp2018","RPA shape: low-Q^{2} suppression (2018)", novarwgt::kRPACCQESuppSyst2018)
Definition: RPASysts.h:14
std::vector< const ISyst * > getAllDirectNumuSysts2017()
Get a vector of all the numu group systs.
Definition: NumuSysts.cxx:169
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:44
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
TH1D * fHists[kNumDets][kNumComponents]
Definition: NumuSysts.h:58
const NOvARwgtSyst kMECInitStateNPFracSyst2017("MECInitStateNPFrac","MEC initial state np fraction", novarwgt::kMECInitStateNPFracSyst2017)
Definition: MECSysts.h:42
void InitializeHistograms() const
Definition: NumuSysts.cxx:221
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
const NOvARwgtSyst kRPARESSyst2018("RPAShapeRES2018","RPA shape: resonance events (2018)", novarwgt::kRPARESSyst2018)
Definition: RPASysts.h:21
Normalization component of the numu GEANT systematic, see docdb 13539.
Definition: NumuSysts.h:14
enum BeamMode string