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"
13 
14 #include <TFile.h>
15 #include <TString.h>
16 #include <vector>
17 #include <iostream>
18 #include <regex>
19 
20 
21 using namespace ana;
22 
23 
24 ana::Loaders* GetLoaders2020(const TString option,
25  const TString period="full")
26 {
27  Loaders * loaders = new Loaders();
28 
29  // Define CAF type
30  auto caftype = ana::ECAFType::kFullCAF;
31  if (option.Contains("nueconcat")) caftype = ana::ECAFType::kNueConcat;
32  if (option.Contains("numuconcat")) caftype = ana::ECAFType::kNumuConcat;
33 
34  // Get flux (FHC or RHC) from option
35  if (!(option.Contains("FHC") || option.Contains("RHC")))
36  {
37  std:: cerr << "Please mention FHC or RHC in option" << std::endl;
38  abort(); // unable to determine if FHC or RHC
39  }
40  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
41 
42  if(option.Contains("CalibrationUp"))
43  loaders = new Prod5AbsCalibLoaders(caftype, flux, +1,
44  period.Data(), period.Data());
45  else if(option.Contains("CalibrationDown"))
46  loaders = new Prod5AbsCalibLoaders(caftype, flux, -1,
47  period.Data(), period.Data());
48  else if(option.Contains("CalibShape"))
49  loaders = new Prod5CalibShapeLoaders(caftype, flux,
50  period.Data(), period.Data());
51  else if(option.Contains("CalibDrift"))
52  loaders = new Prod5CalibDriftLoaders(caftype, flux,
53  period.Data(), period.Data());
54  else if(option.Contains("Cherenkov"))
55  loaders = new Prod5CherenkovLoaders(caftype, flux,
56  period.Data(), period.Data());
57  else if(option.Contains("LightLevelUp"))
58  loaders = new Prod5LightLevelLoaders(caftype, flux, +1,
59  period.Data(), period.Data());
60  else if(option.Contains("LightLevelDown"))
61  loaders = new Prod5LightLevelLoaders(caftype, flux, -1,
62  period.Data(), period.Data());
63  else if(option.Contains("Light_Level_FD")){
64  int sign = (option.Contains("Up") ? +1:-1);
65  loaders = new Prod5LightLevelLoaders(caftype, flux, sign,
66  period.Data(), period.Data());
67  // Reset ND to nominal MC
68  auto temploaders = new Prod5NomLoaders(caftype, flux, period.Data(), period.Data());
69  auto NomNDMC = temploaders->GetLoaderPath(
71  loaders->SetLoaderPath(NomNDMC, caf::kNEARDET,
73  std::cout << "Swapped Light_Level_FD loader: " << NomNDMC << std::endl;
74  }
75  else if(option.Contains("Light_Level_ND")){
76  int sign = (option.Contains("Up") ? +1:-1);
77  loaders = new Prod5NomLoaders(caftype, flux, period.Data(), period.Data());
78  // Start from Nom and shift the ND since setting all
79  // the different FD loaders would be a pain
80  auto temploaders = new Prod5LightLevelLoaders(caftype, flux, sign,
81  period.Data(), period.Data());
82  auto ShiftNDMC = temploaders->GetLoaderPath(
84  loaders->SetLoaderPath(ShiftNDMC, caf::kNEARDET,
86  std::cout << "Swapped Light_Level_ND loader: " << ShiftNDMC << std::endl;
87  }
88  else if(option.Contains("RelativeCalib"))
89  {
90  int sign = (option.Contains("Up") ? +1:-1);
91  loaders = new Prod5AbsCalibLoaders(caftype, flux, sign,
92  period.Data(), period.Data());
93  auto temploaders = new Prod5AbsCalibLoaders(
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  else
103  loaders = new Prod5NomLoaders(caftype, flux, period.Data(), period.Data());
104 
105  return loaders;
106 }
107 
109  const TString option,
110  const TString period="full")
111 {
112  auto caftype = ana::ECAFType::kFullCAF;
113  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
114 
115  if (option.Contains("nueconcat")) caftype = ana::ECAFType::kNueConcat;
116  if (option.Contains("numuconcat")) caftype = ana::ECAFType::kNumuConcat;
117 
118  auto temploaders = new Prod5NomLoaders (caftype, flux, period.Data(), "full");
119  auto nominalNDMC = temploaders->GetLoaderPath(
121  delete temploaders;
122  std::cout << "\n\nSwapping ND Data for fake data \n\n";
123  std::cout << "Before " << loaders->GetLoaderPath(
125  loaders->SetLoaderPath(
127  std::cout << "After " << loaders->GetLoaderPath(
129  << std::endl << std::endl;
130 }
131 
132 //----------------------------------------------------------------------
133 // Parse option for Pt quantiles
134 
135 unsigned int GetPtNQuants( const TString option )
136 {
137  const std::string option_str( option.Data() );
138  const std::regex search_expr( "[a-zA-Z_0-9]*PtQuants([0-9])[a-zA-Z_0-9]*" );
139  std::smatch match;
140  if ( std::regex_match( option_str, match, search_expr ) )
141  {
142  const int nquants = std::stoi( match[1].str() );
143  if ( nquants <= 1 )
144  {
145  std::cout << "Wrong number of Pt quantiles: " << nquants << std::endl;
146  abort();
147  }
148  return nquants;
149  }
150  return 0;
151 }
152 
154 {
155  const unsigned int quant;
156  const unsigned int nquants;
157 };
158 
159 NumuPtQuant GetNumuPtQuant( const TString option )
160 {
161  const std::string option_str( option.Data() );
162  const std::regex search_expr("[a-zA-Z_0-9]*PtQuant([0-9])of([0-9])[a-zA-Z_0-9]*");
163  std::smatch match;
164  if ( std::regex_match( option_str, match, search_expr ) )
165  {
166  const int quant = std::stoi( match[1].str() );
167  const int nquants = std::stoi( match[2].str() );
168  if ( nquants <= 1 )
169  {
170  std::cout << "Wrong number of Numu Pt quantiles: " << nquants << std::endl;
171  abort();
172  }
173  if ( quant < 1 || quant > nquants )
174  {
175  std::cout << "Wrong Numu Pt quantile: " << quant << std::endl;
176  abort();
177  }
178  return NumuPtQuant({unsigned(quant), unsigned(nquants)});
179  }
180  return NumuPtQuant({0, 0});
181 }
182 
183 //////////////////////////////////////////////////////////////////////
184 
185  // Cuts
186  // Numu ND
188  // Nue ND
190  // Numu FD
192  // Nue FD
193  const Cut cutNueFDAll = kNue2020FDAllSamples; // core and peripheral
194  const Cut cutNueFDCore = kNue2020FD; // core only
195 
196  // Axes
197  // NuMu
199  // NuMu Reconstructed Muon Energy
200  const HistAxis axisMuEnergy( "Reconstructed Muon Energy (GeV)", kNumuCCEOptimisedBinning, kMuE);
201  // NuMu Hadronic energy fraction
202  const HistAxis axisHadEfrac( "Hadronic Energy Fraction", Binning::Simple(50, 0.0, 1.01), kHadEFrac);
204 
205 //----------------------------------------------------------------------
206 // Extrapolation defs
207 
208  struct ExtrapDef{
209  TString name;
210  const HistAxis axis;
211  const Cut fdcut;
212  const Cut nuendcut;
213  const Cut numundcut;
214  };
215  //all 2020 cuts, with 3 or 4 bins
216 
217  std::vector <ExtrapDef> GetExtrapolationDefs(
218  const TString analysis, const TString period, const TString option="")
219  {
220 
221  if ( !option.Contains("FHC") && !option.Contains("RHC") )
222  {
223  std:: cerr << "Please mention FHC or RHC in option" << std::endl;
224  abort();
225  }
226  const bool isRHC = option.Contains( "RHC" );
227 
228  std::vector <ExtrapDef> ret;
229  if ( analysis.Contains( "nue" ) )
230  {
231  std::string axisName = "Nue2020Axis";
232 
233  // Signal-only extrapolation with Pt quantiles
234  if ( analysis.Contains( "signal" ) )
235  {
236  const unsigned int PtNQuants = GetPtNQuants( option );
237  // Need 2 sets of quantile cuts, one for ND (numu) and one for FD (nue) due to nue/numu specific axis variables
238  std::vector<Cut> PtQuantCutsND = GetNueQuantCuts2020( isRHC, caf::kNEARDET, PtNQuants, ana::kExtrapPt );
239  std::vector<Cut> PtQuantCutsFD = GetNueQuantCuts2020( isRHC, caf::kFARDET , PtNQuants, ana::kExtrapPt );
240 
241  std::string axisNameSignal = axisName + "_Signal";
242  ret.push_back( { axisNameSignal, kNue2020Axis, cutNueFDAll, kNoCut, cutNumuND } ); // Signal-only extrap, no ND nue cut
243  for ( unsigned int quant = 0; quant < PtNQuants; ++quant )
244  {
245  std::cout << Form( "Adding Nue ExtrapDef for Signal Pt Quantile %d", quant+1 ) << std::endl;
246  std::string axisNameSignalPt = axisNameSignal + Form( "_PtQuant%d", quant+1 );
247  // Signal-only extrap, no ND nue cut
248  ret.push_back( { axisNameSignalPt, kNue2020Axis, cutNueFDAll && PtQuantCutsFD[quant], kNoCut, cutNumuND && PtQuantCutsND[quant] } );
249  }
250  }
251  // Regular and background-only extrapolation
252  else
253  {
254  ret.push_back( { axisName, kNue2020Axis, cutNueFDAll, cutNueND, cutNumuND } );
255  }
256  }
257  else if ( analysis.Contains( "numu" ) )
258  {
259  std::vector<Cut> HadEFracQuantCuts = GetNumuEhadFracQuantCuts2020( isRHC );
260  NumuPtQuant PtQuant = GetNumuPtQuant( option );
261 
262  std::string axisName = "NumuEnergy";
263  for ( unsigned int EhadQuant = 0; EhadQuant < HadEFracQuantCuts.size(); ++EhadQuant )
264  {
265  std::string axisNameEhad = axisName + Form( "_EhadFracQuant%d", EhadQuant+1 );
266  if ( PtQuant.quant > 0 )
267  {
268  // As of now we can only generate predictions for 1 Pt quantile for
269  // each of the EhadFrac quantiles at a time due to output file size and merging
270  std::cout << Form( "Adding Numu ExtrapDef for EhadFrac Quantile %d, Pt Quantile %d of %d",
271  EhadQuant+1, PtQuant.quant, PtQuant.nquants ) << std::endl;
272  std::vector<Cut> PtQuantCuts = GetNumuPtQuantCuts2020( isRHC, EhadQuant+1, PtQuant.nquants );
273  std::string axisNamePt = axisNameEhad + Form( "_PtQuant%dof%d", PtQuant.quant, PtQuant.nquants );
274  ret.push_back( { axisNamePt, axisNumu,
275  cutNumuFD && HadEFracQuantCuts[EhadQuant] && PtQuantCuts[PtQuant.quant-1], kNoCut,
276  cutNumuND && HadEFracQuantCuts[EhadQuant] && PtQuantCuts[PtQuant.quant-1] } );
277  }
278  else
279  {
280  ret.push_back( { axisNameEhad, axisNumu, cutNumuFD && HadEFracQuantCuts[EhadQuant],
281  kNoCut, cutNumuND && HadEFracQuantCuts[EhadQuant] } );
282  }
283  }
284  if ( PtQuant.nquants == 0 )
285  {
286  ret.push_back( { axisName, axisNumu, cutNumuFD, kNoCut, cutNumuND } );
287  }
288  }
289  else
290  std:: cerr << "Please mention nue or numu in analysis" << std::endl;
291 
292  std::cout << "ret has size " << ret.size() << std::endl;
293 
294  return ret;
295  }
unsigned int GetPtNQuants(const TString option)
Near Detector underground.
Definition: SREnums.h:10
const unsigned int quant
const HistAxis axisNumuForNueSig
const XML_Char * name
Definition: expat.h:151
const Cut cutNumuND
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut cutNumuFD
const Cut cutNueND
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
OStream cerr
Definition: OStream.cxx:7
const HistAxis kNumuCCOptimisedAxis2020("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kNumuE2020)
Definition: HistAxes.h:25
const Cut kNue2020ND
Definition: NueCuts2020.h:178
Loaders::FluxType flux
Loaders for absolute calibration paths/definitions.
Definition: Prod5Loaders.h:122
Loaders for Cherenkov paths/definitions.
Definition: Prod5Loaders.h:173
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)
std::vector< Cut > GetNumuPtQuantCuts2020(const bool isRHC, const unsigned int ehad_quant, const unsigned int nquants)
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod5Loaders.h:101
const Var kHadEFrac
Definition: NumuVars.h:24
const HistAxis kNue2020Axis("NuE Energy / Analysis Bin", kNue2020Binning, kNue2020AnaBin)
Use this Axis for Ana2020, official Axis.
Definition: NueCuts2020.h:195
std::vector< ExtrapDef > GetExtrapolationDefs(const TString analysis, const TString period)
Loaders for calibration drift (detector aging) paths/definitions.
Definition: Prod5Loaders.h:218
std::vector< Cut > GetNueQuantCuts2020(const bool isRHC, const caf::Det_t det, const unsigned int nquants, const ExtrapVar var)
const Cut kNumu2020FD
Definition: NumuCuts2020.h:59
ana::Loaders * GetLoaders2020(const TString option, const TString period="full")
const Cut kNue2020FD
Definition: NueCuts2020.h:65
OStream cout
Definition: OStream.cxx:6
Loaders for light level paths/definitions.
Definition: Prod5Loaders.h:147
const unsigned int nquants
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
int regex_match(const std::string &s, const std::string &p)
Definition: RegexMatch.cxx:7
NumuPtQuant GetNumuPtQuant(const TString option)
const HistAxis axisNumu
const Cut kNue2020FDAllSamples
Definition: NueCuts2020.h:84
const Var kMuE
Definition: NumuVars.h:22
Loaders for calibration shape paths/definitions.
Definition: Prod5Loaders.h:196
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)
void SwapNDDataLoader(Loaders *loaders, const TString option, const TString period="full")
std::vector< Cut > GetNumuEhadFracQuantCuts2020(const bool isRHC, const unsigned int nquants)
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
enum BeamMode string