SystsBeamLoad.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Loaders.h"
4 #include "CAFAna/Cuts/Cuts.h"
5 #include "NuXAna/Cuts/NusCuts.h"
12 #include "CAFAna/Systs/Systs.h"
14 #include "NuXAna/Vars/HistAxes.h"
17 
18 using namespace ana;
19 
21 {
22  TH1::AddDirectory(0);
23 
25 
31 
33 
34  // Vector of sigmas, levels to set each systematic
35  std::vector<int> sigmas{-3, -2, -1, +1, +2, +3};
36 
37  // A map to set up systematics, starting with a label for the systematic,
38  // and finishing with a vector for items necessary to set up the Syst object
39  std::map<std::string, std::vector<std::string> > shiftlabels;
40  shiftlabels["BeamHornCurrent"] = std::vector<std::string>();
41  shiftlabels["BeamSpotSize"] = std::vector<std::string>();
42  shiftlabels["BeamPosX"] = std::vector<std::string>();
43  shiftlabels["BeamPosY"] = std::vector<std::string>();
44  shiftlabels["BeamH1PosXY"] = std::vector<std::string>();
45  shiftlabels["BeamH2PosXY"] = std::vector<std::string>();
46  shiftlabels["BeamTarget"] = std::vector<std::string>();
47  shiftlabels["BeamExp"] = std::vector<std::string>();
48  shiftlabels["BeamGeomWater"] = std::vector<std::string>();
49  shiftlabels["BeamFlukaVersion"] = std::vector<std::string>();
50  shiftlabels["BeamG4"] = std::vector<std::string>();
51  shiftlabels["BeamNA49"] = std::vector<std::string>();
52 
53  shiftlabels["BeamHornCurrent"] .push_back("HornCurrent");
54  shiftlabels["BeamHornCurrent"] .push_back("Horn Current");
55  shiftlabels["BeamSpotSize"] .push_back("BmSpotSize");
56  shiftlabels["BeamSpotSize"] .push_back("Spot Size");
57  shiftlabels["BeamPosX"] .push_back("BmPosX");
58  shiftlabels["BeamPosX"] .push_back("Beam Position X");
59  shiftlabels["BeamPosY"] .push_back("BmPosY");
60  shiftlabels["BeamPosY"] .push_back("Beam Position Y");
61  shiftlabels["BeamH1PosXY"] .push_back("Horn1PosXY");
62  shiftlabels["BeamH1PosXY"] .push_back("Horn 1 Position");
63  shiftlabels["BeamH2PosXY"] .push_back("Horn2PosXY");
64  shiftlabels["BeamH2PosXY"] .push_back("Horn 2 Position");
65  shiftlabels["BeamTarget"] .push_back("TargetPos");
66  shiftlabels["BeamTarget"] .push_back("Target Position");
67  shiftlabels["BeamExp"] .push_back("ExpBField");
68  shiftlabels["BeamExp"] .push_back("Exponential B Field");
69  shiftlabels["BeamGeomWater"] .push_back("GeomWater");
70  shiftlabels["BeamGeomWater"] .push_back("New Horn Geometry and 1mm water");
71  shiftlabels["BeamFlukaVersion"].push_back("FlukaVersion");
72  shiftlabels["BeamFlukaVersion"].push_back("Fluka version");
73  shiftlabels["BeamG4"] .push_back("G4NuMI");
74  shiftlabels["BeamG4"] .push_back("G4NuMI");
75  shiftlabels["BeamNA49"] .push_back("NA49");
76  shiftlabels["BeamNA49"] .push_back("NA49");
77 
78  // A map of systematic objects, indexed by the same label as above
79  std::map<std::string, std::map<int, SystShifts*> > systematics;
80 
81  std::string rootfile = FindCAFAnaDir() + "/data/beam/beamSyst_histos.root";
82  for(const auto& shiftlabel : shiftlabels) {
83  for(const auto& sigma : sigmas) {
84  systematics[shiftlabel.first][sigma] = new SystShifts(
85  new BeamSyst(rootfile, shiftlabel.second[0], shiftlabel.second[1]),
86  sigma
87  );
88  }
89  }
90 
91  // Define a map of samples and cuts/axes to set them up
92  // Each sample, essentially, is a way to set up the analysis
93  std::string labelRecoE = "Calorimetric Energy (GeV)";
94  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
95  std::pair<Cut*, Cut*> Ana01(new Cut(kNusND), new Cut(kNusFD));
96  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
98  Ana01
99  );
100 
101  // Set up maps to the decompositions/predictions (each flavor component)
102  // Nominal maps are indexed purely by sample label
103  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
104  std::map<std::string, IDecomp*> decomps_nominal;
105  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
106  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
107  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
108  std::map<std::string, PredictionSterile*> predsSt_nominal;
109  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
110 
111  // Set up the actual decompositions/predictions
112  // sample.first = the sample label
113  // shiftlabel.first = the systematic shift label
114  // sample.second.first = the sample HistAxis
115  // sample.second.second = the sample Cut
116  // syst.first = the integer number of sigmas for the shift
117  // syst.second = the systematic being shifted
118  for(const auto& sample : cut_samples) {
119  SterileGenerator gen(
120  *sample.second.first, kNCBinsNumuCCAxis,
121  *sample.second.second.second, *sample.second.second.first,
123  );
124 
125  // Set up the nominal decompositions
126  // Only one is needed per sample
127  decomps_nominal[sample.first] = new CheatDecomp(
129  *sample.second.first, *sample.second.second.first,
131  );
132  predsNE_nominal[sample.first] = new PredictionNoExtrap(
133  loaders,
134  *sample.second.first, *sample.second.second.second,
136  );
137  predsSt_nominal[sample.first] = (PredictionSterile*)gen.Generate(
138  loaders
139  ).release();
140 
141  // Set up the shifted decompositions
142  // One is needed for every sample, systematic, and sigma level
143  for(const auto& shiftlabel : shiftlabels) {
144  for(const auto& syst : systematics[shiftlabel.first]) {
145  decomps_shifted[sample.first][shiftlabel.first][syst.first] = new CheatDecomp(
147  *sample.second.first, *sample.second.second.first,
148  *syst.second, kTuftsWeightCC
149  );
150  predsNE_shifted[sample.first][shiftlabel.first][syst.first] = new PredictionNoExtrap(
151  loaders,
152  *sample.second.first, *sample.second.second.second,
153  *syst.second, kTuftsWeightCC
154  );
155  predsSt_shifted[sample.first][shiftlabel.first][syst.first] = (PredictionSterile*)gen.Generate(
156  loaders,
157  *syst.second
158  ).release();
159  }
160  }
161  }
162 
163  loaders.Go();
164 
165  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
166  std::string filenm = "SystsBeam";
167 
168  std::string fullLocation = folder + filenm + ".root";
169  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
170 
171  SaveMaps(rootF,
172  decomps_nominal, decomps_shifted,
173  predsNE_nominal, predsNE_shifted,
174  predsSt_nominal, predsSt_shifted
175  );
176 }
Near Detector underground.
Definition: SREnums.h:10
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const std::string fnametau_concat
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
Beam systematics:
Definition: BeamSysts.h:38
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kNusFD
Definition: NusCuts.h:46
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
void SystsBeamLoad()
Definition: SystsBeamLoad.C:20
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
Definition: PPFXHelper.h:1077
Generates extrapolated NC predictions using ProportionalDecomp.
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const SpillCut kOnly14DB([](const caf::SRSpillProxy *spill){if(spill->det!=caf::kFARDET) return true;std::bitset< 14 > binary(spill->dibmask);for(int i=0;i< 14;++i){if(!binary[i]) return false;}return true;})
const std::string fnamenear_concat
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
const std::string fnameneardata_concat
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
const Cut kNumuND
Definition: NumuCuts.h:55
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
double sigma(TH1F *hist, double percentile)
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const HistAxis kNCBinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNCDisappearanceEnergyBinning, kCCE)
Definition: HistAxes.h:9
const Var kTuftsWeightCC
Definition: XsecTunes.h:30
std::vector< Loaders * > loaders
Definition: syst_header.h:386
A prediction object compatible with sterile oscillations.
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Prediction that just uses FD MC, with no extrapolation.
const std::string fnameswap_concat
const std::string fnamefar_concat
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
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
const Cut kNusND
Definition: NusCuts.h:71
const Binning kNCDisappearanceEnergyBinning
Energy binnings used in Ana01 for nus extrapolation.
Definition: Binning.cxx:68
enum BeamMode string