syst_variations.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/Loaders.h"
6 #include "CAFAna/Cuts/Cuts.h"
7 #include "3FlavorAna/Cuts/NueCuts2018.h" // kNue2018NDCVNSsb, kNue2018FDAllSamples,
8  // kNue2018FD
9 #include "3FlavorAna/Cuts/NumuCuts2018.h" // kNumuCutND2018, kNumuCutFD2018,
10  // kNue2018Axis,
11  // kNue2018AxisMergedPeripheral
14 #include "CAFAna/Systs/XSecSysts.h"
16 #include "CAFAna/Systs/BeamSysts.h"
24 #include "3FlavorAna/Systs/GeniePCASyst.h" //Adding GENIEPCA systematic prediction function (GetGeniePrincipals2018).
25 
26 #include <TFile.h>
27 #include <TString.h>
28 #include <vector>
29 #include <iostream>
30 
31 
32 using namespace ana;
33 
34 
35 ana::Loaders* GetLoaders2018(const TString option,
36  const TString period="full")
37 {
38  Loaders * loaders = new Loaders();
39 
40  // Define CAF type
41  auto caftype = ana::ECAFType::kFullCAF;
42  if (option.Contains("nueconcat")) caftype = ana::ECAFType::kNueConcat;
43  if (option.Contains("numuconcat")) caftype = ana::ECAFType::kNumuConcat;
44 
45  // Get flux (FHC or RHC) from option
46  if (!(option.Contains("FHC") || option.Contains("RHC")))
47  {
48  std:: cerr << "Please mention FHC or RHC in option" << std::endl;
49  abort(); // unable to determine if FHC or RHC
50  }
51  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
52 
53  if(option.Contains("CalibrationUp"))
54  loaders = new Prod4AbsCalibLoaders(caftype, flux, +1,
55  period.Data(), period.Data());
56  else if(option.Contains("CalibrationDown"))
57  loaders = new Prod4AbsCalibLoaders(caftype, flux, -1,
58  period.Data(), period.Data());
59  else if(option.Contains("CalibShape"))
60  loaders = new Prod4CalibShapeLoaders(caftype, flux,
61  period.Data(), period.Data());
62  else if(option.Contains("Cherenkov"))
63  loaders = new Prod4CherenkovLoaders(caftype, flux,
64  period.Data(), period.Data());
65  else if(option.Contains("LightLevelUp"))
66  loaders = new Prod4LightLevelLoaders(caftype, flux, +1,
67  period.Data(), period.Data());
68  else if(option.Contains("LightLevelDown"))
69  loaders = new Prod4LightLevelLoaders(caftype, flux, -1,
70  period.Data(), period.Data());
71  else if(option.Contains("LightLevelNom")) // FD ONLY
72  {
73  loaders = new Prod4LightLevelLoaders(caftype, flux, 0,
74  period.Data(), period.Data());
75  auto temploaders = new Prod4NomLoaders(caftype, flux,
76  period.Data(), "full");
77  auto nominalNDMC = temploaders->GetLoaderPath(
79  delete temploaders;
80  std::cout << "LightLevelNom ONLY for FD" << std::endl
81  << " --> adding loader path for standard ND nominal\n\n";
82  loaders->SetLoaderPath(
84  std::cout << "Now: " << loaders->GetLoaderPath(
86  << std::endl << std::endl;
87  }
88  else if(option.Contains("RelativeCalib"))
89  {
90  int sign = (option.Contains("Up") ? +1:-1);
91  loaders = new Prod4AbsCalibLoaders(caftype, flux, sign,
92  period.Data(), period.Data());
93  auto temploaders = new Prod4AbsCalibLoaders(
94  caftype, flux, - sign, period.Data(), period.Data());
95  auto oppositeNDMC = temploaders->GetLoaderPath(
97  loaders->SetLoaderPath(oppositeNDMC, caf::kNEARDET,
99  std::cout << "Swapped calibration loader" << oppositeNDMC << std::endl;
100  delete temploaders;
101  }
102  // TODO: Do we still need this?
103  else if (option.Contains("NoTau"))
104  {
105  std::cout << "WARNING: Option NoTau currently does nothing" << std::endl;
106  //loaders->DisableLoader(caf::kFARDET, Loaders::kMC,
107  // ana::kBeam, Loaders::kTauSwap);
108  /*
109  if (option.Contains("2sigmaNoTau"))
110  {
111  std::string posneg = (option.Contains("Up")) ? "pos" : "neg";
112  std::string NonSwa = "prod_decaf_R17-03-01-prod3reco.m_fd_genie_nonswap_fhc_nova_v08_full_calib-shift-fd-xyview-"+posneg+"-offset-2sigma_nue_or_numu_or_nus_contain_v1";
113  std::string FluSwa = "prod_decaf_R17-03-01-prod3reco.m_fd_genie_fluxswap_fhc_nova_v08_full_calib-shift-fd-xyview-"+posneg+"-offset-2sigma_nue_or_numu_or_nus_contain_v1";
114  loaders->SetLoaderPath( NonSwa, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap );
115  loaders->SetLoaderPath( FluSwa, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap );
116  }*/
117  }
118  else loaders = new Prod4NomLoaders(caftype, flux,
119  period.Data(), period.Data());
120 
121  return loaders;
122 }
123 
124 // Use the shifted ND MC as fake data
125 // Use Swap** below to use nominal NDMC as fake data instead
126 /*
127 ana::Loaders * GetLoaders2017FakeData (const TString option, const TString period="full"){
128 
129  Loaders * loaders, *temploaders;
130 
131  auto caftype = ana::ECAFType::kDecaf ;
132 
133  if (option.Contains("nueconcat")) caftype = ana::ECAFType::kNueConcat;
134  if (option.Contains("numuconcat")) caftype = ana::ECAFType::kNumuConcat;
135 
136  // Get flux (FHC or RHC) from option
137  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
138 
139  if(option.Contains("LightLevel")) loaders = new Prod4LightLevelLoaders(caftype, flux, 0,period.Data(), period.Data());
140  else loaders = new Prod4NomLoaders (caftype, flux, period.Data(), period.Data());
141 
142  if(option.Contains("Nom")) {
143  std::cerr << "You picked fake data = nominal MC. "
144  <<"This is a silly in and out test, and all your ratios will be one. "
145  <<"\nEdit the code if this is truly what you want" << std::endl;
146  exit(1);
147  }
148 
149  if(option.Contains("CalibrationUp")) temploaders = new Prod4AbsCalibLoaders(caftype, flux, +1, period.Data(), period.Data());
150  else if(option.Contains("CalibrationDown")) temploaders = new Prod4AbsCalibLoaders(caftype, flux, -1, period.Data(), period.Data());
151  else if(option.Contains("CalibShape")) temploaders = new Prod4CalibShapeLoaders(caftype, flux, period.Data(), period.Data());
152  else if(option.Contains("Cherenkov")) temploaders = new Prod4CherenkovLoaders(caftype, flux, period.Data(), period.Data());
153  else if(option.Contains("LightLevelUp")) temploaders = new Prod4LightLevelLoaders(caftype, flux,+1, period.Data(), period.Data());
154  else if(option.Contains("LightLevelDown")) temploaders = new Prod4LightLevelLoaders(caftype, flux,-1,period.Data(), period.Data());
155  else temploaders = new Prod4NomLoaders (caftype, flux, period.Data(), period.Data());
156 
157  if(option.Contains("RHC")) temploaders = new Prod4NomLoaders(caftype, flux, period.Data(), "full", Loaders::kRHC);
158 
159  auto shiftedMCpath = temploaders->GetLoaderPath(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
160  loaders->SetLoaderPath(shiftedMCpath, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);
161  std::cout << "Swapped real data loader for " << shiftedMCpath << std::endl;
162  return loaders;
163 }
164 
165 */
166 
168  const TString option,
169  const TString period="full")
170 {
171  auto caftype = ana::ECAFType::kFullCAF;
172  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
173 
174  if (option.Contains("nueconcat")) caftype = ana::ECAFType::kNueConcat;
175  if (option.Contains("numuconcat")) caftype = ana::ECAFType::kNumuConcat;
176 
177  auto temploaders = new Prod4NomLoaders (caftype, flux, period.Data(), "full");
178  auto nominalNDMC = temploaders->GetLoaderPath(
180  delete temploaders;
181  std::cout << "\n\nSwapping ND Data for fake data \n\n";
182  std::cout << "Before " << loaders->GetLoaderPath(
184  loaders->SetLoaderPath(
186  std::cout << "After " << loaders->GetLoaderPath(
188  << std::endl << std::endl;
189 }
190 
191 //////////////////////////////////////////////////////////////////////
192 
193  // Cuts
194  // Numu ND
196  // Nue ND
198  // Numu FD
200  // Nue FD
201  const Cut cutNueFDAll = kNue2018FDAllSamples; // core and peripheral
202  const Cut cutNueFDCore = kNue2018FD; // core only
203 
204  // Axes
205  // NuMu
207  // NuMu Reconstructed Muon Energy
208  const HistAxis axisMuEnergy( "Reconstructed Muon Energy (GeV)", kNumuCCEOptimisedBinning, kMuE);
209  // NuMu Hadronic energy fraction
210  const HistAxis axisHadEfrac( "Hadronic Energy Fraction", Binning::Simple(50, 0.0, 1.01), kHadEFrac);
212 
213  // Nue
214  // Mainly using Axes defined in NueCuts2018.h
215  // Only extra axis we need is for reconstructed energy
216  const HistAxis axisNueEnergy("Reconstructed neutrino energy (GeV)",
218 
219  struct ExtrapDef{
220  TString name;
221  const HistAxis axis;
222  const Cut fdcut;
223  const Cut nuendcut;
224  const Cut numundcut;
225  };
226  //all 2017 cuts, with 3 or 4 bins
227 
228  std::vector <ExtrapDef> GetExtrapolationDefs(
229  const TString analysis, const TString period, const TString option="")
230  {
231  std::vector <ExtrapDef> ret;
232  if (analysis.Contains("nue"))
233  {
234  // Not making core only predictions for 2018 analysis
235  // Nue 2018 analysis axis/binning
236  // This is the only axis we need for Nue systematics 2018
237  ret.push_back(
238  {"Nue2018Axis", kNue2018Axis,
239  cutNueFDAll, cutNueND, cutNumuND});
240  // Nue 2018 analysis axis/binning with single peripheral bin
241  //ret.push_back(
242  // {"Nue2018AxisMergedPeripheral", kNue2018AxisMergedPeripheral,
243  // cutNueFDAll, cutNueND, cutNumuND});
244  //ret.push_back(
245  // {"NueEnergy", axisNueEnergy, cutNueFDCore, cutNueND, cutNumuND});
246  }
247  else if(analysis.Contains("numu"))
248  {
249  // Want to load the quantile files from disk...
250  std::string sFHC_RHC = (analysis.Contains("FHC") == true ? "fhc" :"rhc");
251  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/official/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
252  std::cout << "\nI'm going to try and get my quantile boundaries from " << fdspecfile << ".\n\n" << std::endl;
253  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile).c_str() );
254  // Check that the file is valid.
255  if (inFile->IsZombie()) {
256  std::cout << "Problem with file " << fdspecfile << std::endl;
257  abort();
258  }
259  // Get the hadronic cut boundaries
260  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
261  const unsigned int nquantiles = 4;
262  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2( FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, nquantiles );
263 
264  ret.push_back({"NumuEnergy", axisNumu, cutNumuFD, kNoCut, cutNumuND});
265  for ( unsigned int quant=0; quant<HadEFracQuantCuts.size(); ++quant ) {
266  ret.push_back({"Energy_Quant" +std::to_string(quant+1), axisNumu , cutNumuFD && HadEFracQuantCuts[quant], kNoCut, cutNumuND && HadEFracQuantCuts[quant]});
267  //ret.push_back({"MuonEn_Quant" +std::to_string(quant+1), axisMuEnergy, cutNumuFD && HadEFracQuantCuts[quant], kNoCut, cutNumuND && HadEFracQuantCuts[quant]});
268  //ret.push_back({"HadEFrac_Quant"+std::to_string(quant+1), axisHadEfrac, cutNumuFD && HadEFracQuantCuts[quant], kNoCut, cutNumuND && HadEFracQuantCuts[quant]});
269  }
270  }
271  else
272  std:: cerr << "Please mention nue or numu in analysis" << std::endl;
273 
274  std::cout << "ret has size " << ret.size() << std::endl;
275 
276  return ret;
277  }
Near Detector underground.
Definition: SREnums.h:10
const HistAxis axisNumuForNueSig
const XML_Char * name
Definition: expat.h:151
const Cut cutNumuND
ana::Loaders * GetLoaders2018(const TString option, const TString period="full")
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut cutNumuFD
const Cut cutNueND
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
const Cut kNue2018NDCVNSsb
Definition: NueCuts2018.h:153
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
OStream cerr
Definition: OStream.cxx:7
Loaders for Cherenkov paths/definitions.
Definition: Prod4Loaders.h:172
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
Loaders::FluxType flux
const Binning kNueSAEnergyBinning
The energy part of the SA 2D binning.
Definition: Binnings.cxx:17
Loaders for calibration shape paths/definitions.
Definition: Prod4Loaders.h:196
std::string GetLoaderPath(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap) const
Definition: Loaders.cxx:89
const HistAxis axisHadEfrac("Hadronic Energy Fraction", Binning::Simple(50, 0.0, 1.01), kHadEFrac)
ifstream inFile
Definition: AnaPlotMaker.h:34
Loaders for absolute calibration paths/definitions.
Definition: Prod4Loaders.h:119
Loaders for light level paths/definitions.
Definition: Prod4Loaders.h:145
const Var kNueEnergy2018([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC(sr);else return kNueEnergyFHC(sr);})
Definition: NueEnergy2018.h:25
const Var kHadEFrac
Definition: NumuVars.h:24
std::vector< ExtrapDef > GetExtrapolationDefs(const TString analysis, const TString period)
const Cut kNue2018FD
Use this cut for the full FD Core selection, apply for both RHC and FHC.
Definition: NueCuts2018.h:54
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
OStream cout
Definition: OStream.cxx:6
const Cut cutNueFDAll
std::vector< Loaders * > loaders
Definition: syst_header.h:386
const Binning kNumuCCEOptimisedBinning
Optimised binning for numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39; in that talk...
Definition: Binnings.cxx:28
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const HistAxis axisNumu
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
const Var kMuE
Definition: NumuVars.h:22
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const HistAxis axisMuEnergy("Reconstructed Muon Energy (GeV)", kNumuCCEOptimisedBinning, kMuE)
const HistAxis kNue2018Axis("NuE Energy / Analysis Bin", kNue2018Binning, kNue2018AnaBin)
Use this Axis for Ana2018, official Axis.
void SwapNDDataLoader(Loaders *loaders, const TString option, const TString period="full")
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
def sign(x)
Definition: canMan.py:197
const Cut cutNueFDCore
const HistAxis axisNueEnergy("Reconstructed neutrino energy (GeV)", kNueSAEnergyBinning, kNueEnergy2018)
enum BeamMode string