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 <vector>
17 
18 namespace ana
19 {
23  kNumuSummedSmallGENIESyst(FindCAFAnaDir()+"/data/numu/summed_genie_systs_nd.root",
24  FindCAFAnaDir()+"/data/numu/summed_genie_systs_fd.root");
25 
32 
33 
34  //----------------------------------------------------------------------
36  Shift(double sigma, caf::SRProxy* /*sr*/, double& weight) const
37  {
38  weight *= 1+.01*sigma; // 1% shift
39  }
40 
41  //----------------------------------------------------------------------
43  Shift(double sigma, caf::SRProxy* sr, double& weight) const
44  {
45  const double scale = 1 + 0.01*sigma; //shifts 1% across boundary
46 
47  sr->energy.numu.trkccE *= scale;
48  }
49 
50  //----------------------------------------------------------------------
52  Shift(double sigma, caf::SRProxy* sr, double& weight) const
53  {
54  if(sr->mc.nnu == 0) return;
55  if(!sr->mc.nu[0].iscc) weight *= 1+sigma;
56  if(weight < 0) weight = 0;
57  }
58 
59  //----------------------------------------------------------------------
61  Shift(double sigma, caf::SRProxy* sr, double& weight) const
62  {
63  if(sr->mc.nnu == 0) return;
64  if(sr->mc.nu[0].iscc && abs(sr->mc.nu[0].pdg) == 16) weight *= 1+sigma;
65  }
66 
67  //----------------------------------------------------------------------
68  void NumuNormSyst::
69  Shift(double sigma, caf::SRProxy* /*sr*/, double& weight) const
70  {
71  weight *= 1+.01*sigma; // 1% absolute error
72  }
73 
74  //----------------------------------------------------------------------
76  Shift(double sigma, caf::SRProxy* sr, double& weight) const
77  {
78  // 1% difference between detectors + 0.7% from intensity + 3% ND norm (rock) + 3.7% (cosmic overlay) all in quadrature
79  if (sr->hdr.det == caf::kFARDET) return;
80  weight *= 1 + .049 * sigma;
81  }
82 
83  //----------------------------------------------------------------------
85  Shift(double sigma, caf::SRProxy* sr, double& weight) const
86  {
87  const double scale = 1 + .05*sigma;
88 
89  sr->slc.calE *= scale;
90  sr->energy.numu.trknonqeE *= scale;
91  sr->energy.numu.trkqeE *= scale;
92  sr->energy.numu.trkccE *= scale;
93  }
94 
95  //----------------------------------------------------------------------
97  Shift(double sigma, caf::SRProxy* sr, double& weight) const
98  {
99  if (sr->hdr.det != caf::kFARDET) return;
100 
101  const double scale = 1 + .025*sigma /// 5% difference between detectors
102  *( (sr->hdr.det==caf::kFARDET) ? 1 : -1 );
103 
104  sr->slc.calE *= scale;
105  sr->energy.numu.trknonqeE *= scale;
106  sr->energy.numu.trkqeE *= scale;
107  sr->energy.numu.trkccE *= scale;
108  }
109 
110 
111  //----------------------------------------------------------------------
112  // 2018 systematics list -- still a work in progress
113  std::vector<const ISyst*> getAllDirectNumuSysts2018() {
114  // This is the list of (mostly) *not-tiny* systematics for numu,
115  // for making things like the ND MC uncertainty bands and various cross checks.
116  // The latest greatest xsec systs are taken from getAllXsecSysts_2018()
117  // Beam stuff will need updating (not sure RHC is working???)
118  // The direct absolute hadE is no longer commented out
119  // (we will use the special calib files but that's not what this function is for)
120  // The really tiny GENIE syst that were included in the 2017 version have been removed
121 
122  std::vector<const ISyst*> systs =
123  {
129 
135 
139 
141 
148 
149  &kMuEScaleSyst2017, // correlated muon E scale applied to muon track length
150  &kRelMuEScaleSyst2017, // uncorrelated muon E scale applied to muon track length
151  &kDirectRelHadEScaleSyst2017, // relative near/far hadE syst
152  &kDirectHadEScaleSyst2017, // absolute hadE syst, correlated between detectors
153  };
154 
155  std::set<rwgt::ReweightLabel_t> largeGENIEParameters;
156  largeGENIEParameters.insert(rwgt::fReweightMaCCRES);
157  largeGENIEParameters.insert(rwgt::fReweightMvCCRES);
158  largeGENIEParameters.insert(rwgt::fReweightMaNCRES);
159 
160  for (const auto & parameter:largeGENIEParameters)
161  systs.push_back(GetGenieKnobSyst(parameter));
162 
163  return systs;
164  }
165 
166  //----------------------------------------------------------------------
167  // 2017 systematics list;
168  std::vector<const ISyst*> getAllDirectNumuSysts2017()
169  {
170  std::vector<const ISyst*> systs =
171  {
180  &kRPACCQEEnhSyst2017, // new RPA error (1-sided more enhancement)
181  &kRPACCQESuppSyst2017, // new RPA error (1-sided more suppression)
183  &kMuEScaleSyst2017, // correlated muon E scale applied to muon track length
184  &kRelMuEScaleSyst2017, // uncorrelated muon E scale applied to muon track length
185  &kDirectRelHadEScaleSyst2017, // directly applied 5% hadronic E uncorrelated syst
186  // kDirectHadEScaleSyst2017 will be provided elsewhere soon.
187  &kMAQEGenieReducedSyst2017, // replacing fReweightNormCCQE and fReweightMaCCQEshape from SA. Reduced M_A uncertainty (as compared to GENIE M_A knob)
188  };
189 
190  systs.push_back(&kNumuSummedSmallGENIESyst); // horrible and needs to be re-looked at for TA
191  std::set<rwgt::ReweightLabel_t> largeGENIEParameters;
192  largeGENIEParameters.insert(rwgt::fReweightMaNCEL);
193  largeGENIEParameters.insert(rwgt::fReweightMaCCRES);
194  largeGENIEParameters.insert(rwgt::fReweightMvCCRES);
195  largeGENIEParameters.insert(rwgt::fReweightMaNCRES);
196  largeGENIEParameters.insert(rwgt::fReweightMvNCRES);
197  largeGENIEParameters.insert(rwgt::fReweightCCQEPauliSupViaKF);
198 
199  for (const auto & parameter:largeGENIEParameters)
200  systs.push_back(GetGenieKnobSyst(parameter));
201 
202  return systs;
203  }
204 
205  //----------------------------------------------------------------------
208  const std::string& fnameFD)
209  : ISyst("numuSumSmallGENIE", "Summed Small GENIE")
210  {
211  fHistoFnames[0] = fnameND;
212  fHistoFnames[1] = fnameFD;
213 
214  // Signal that these aren't loaded yet
215  fHists[kFD][kSig] = 0;
216  }
217 
218 
219  //----------------------------------------------------------------------
221  {
222 
223  // Bail if we've already loaded histograms
224  if(fHists[kFD][kSig]) return;
225 
226 
227  for(int iDet = kND; iDet < kNumDets; ++iDet){
228  TFile fin(fHistoFnames[iDet].c_str(), "read");
229  if(fin.IsZombie()){
230  std::cerr << "Warning: couldn't open " <<
231  fHistoFnames[iDet] << std::endl;
232  return;
233  }
234 
235  fHists[iDet][kSig] = (TH1D*)fin.Get("allcc_sig")->Clone();
236  fHists[iDet][kSig]->SetDirectory(0);
237  fHists[iDet][kBkg] = (TH1D*)fin.Get("allcc_bkg")->Clone();
238  fHists[iDet][kBkg]->SetDirectory(0);
239  }
240  }
241 
242 
243  //----------------------------------------------------------------------
245  caf::SRProxy* sr, double& weight) const
246  {
248  /// Initialize the scale factor for the weight
249  double shift = 0;
250  if(sr->mc.nnu == 0) return;
251 
252  int det;
253  switch (sr->hdr.det){
254  case 1 : det = kND; break;
255  case 2 : det = kFD; break;
256  default: std::cerr <<"Wrong detector; ignore" << std::endl; return;
257  }
258 
259  TH1D* hist = NULL;
260  if(abs(sr->mc.nu[0].pdg) == 14 && sr->mc.nu[0].iscc)
261  hist = fHists[det][kSig];
262  else
263  hist = fHists[det][kBkg];
264 
265  double energy = sr->energy.numu.trkccE;
266 
267  if (energy > hist->GetXaxis()->GetXmin() &&
268  energy < hist->GetXaxis()->GetXmax() )
269  shift = hist->Interpolate(energy) - 1;
270  weight *= 1 + shift * sigma;
271 
272  }
273 
274  //----------------------------------------------------------------------
276  {
277  for(int i = 0; i < kNumDets; ++i){
278  for(int j = 0; j < kNumComponents; ++j){
279  delete fHists[i][j];
280  fHists[i][j] = 0;
281  }
282  }
283  }
284 
285 }
BeamSyst * GetPPFXPrincipals(int PCIdx)
Definition: BeamSysts.cxx:63
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:26
caf::Proxy< float > trkccE
Definition: SRProxy.h:194
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:97
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:30
const NumuTauContaminationSyst kNumuTauContaminationSyst
Definition: NumuSysts.cxx:27
const NumuNormSyst kNumuNormSyst
Definition: NumuSysts.cxx:28
caf::Proxy< float > trknonqeE
Definition: SRProxy.h:195
const Var weight
GeniePCASyst * GetGeniePrincipals2018(int PCIdx)
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2119
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:213
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:597
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:29
caf::Proxy< short int > nnu
Definition: SRProxy.h:596
const NumuGEANTNormSyst kNumuGEANTNormSyst
Definition: NumuSysts.cxx:20
std::string fHistoFnames[kNumDets]
Definition: NumuSysts.h:56
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2118
const NOvARwgtSyst kMECEnuShapeSyst2018Nu("MECEnuShapeNu","MEC E_{#nu} shape, neutrinos", novarwgt::kMECEnuShapeSyst2017_NuOnly)
Definition: MECSysts.h:31
std::string FindCAFAnaDir()
Definition: Utilities.cxx:208
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:80
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
std::vector< const ISyst * > getAllDirectNumuSysts2018()
Definition: NumuSysts.cxx:113
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:69
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:21
caf::Proxy< float > trkqeE
Definition: SRProxy.h:196
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:76
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:125
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:85
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:244
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:31
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2120
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const DirectRelHadEScaleSyst2017 kDirectRelHadEScaleSyst2017(0.05)
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2124
caf::Proxy< float > calE
Definition: SRProxy.h:1271
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:207
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:36
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:52
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:61
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:168
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: NumuSysts.cxx:43
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:220
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231
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