SystsMCStatsLoad.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"
8 #include "CAFAna/Extrap/ExtrapSterile.h"
11 #include "CAFAna/Systs/Systs.h"
13 #include "NuXAna/Vars/HistAxes.h"
16 
17 using namespace ana;
18 
20 {
21  TH1::AddDirectory(0);
22 
23  std::string fnear = "defname: prod_decaf_S16-04-08_nd_genie_nonswap_genierw_fhc_nova_v08_full_nus_or_numu_contain_v1 with stride 5";
24  std::string ffar = "defname: prod_decaf_S16-04-08_fd_genie_nonswap_fhc_nova_v08_full_nus_contain_v1 with stride 5";
25  std::string fswap = "defname: prod_decaf_S16-04-08_fd_genie_fluxswap_fhc_nova_v08_full_nus_contain_v1 with stride 5";
26  std::string ftau = "defname: prod_decaf_S16-04-08_fd_genie_tau_fhc_nova_v08_full_nus_contain_v1 with stride 5";
27 
29 
35 
37 
38  // A map to set up systematics, starting with a label for the systematic,
39  // and finishing with a vector for items necessary to set up the Syst object
40  std::map<std::string, std::string> shiftlabels;
41  shiftlabels["Set1"] = "Set 1";
42  shiftlabels["Set2"] = "Set 2";
43  shiftlabels["Set3"] = "Set 3";
44  shiftlabels["Set4"] = "Set 4";
45 
46  const Cut kMod0(
47  [](const caf::SRProxy* sr)
48  {
49  if(sr->hdr.ismc == false) { return true; }
50  int evtno = sr->hdr.evt;
51  if((evtno % 5) == 0) { return true; }
52  return false;
53  });
54  const Cut kMod1(
55  [](const caf::SRProxy* sr)
56  {
57  if(sr->hdr.ismc == false) { return true; }
58  int evtno = sr->hdr.evt;
59  if((evtno % 5) == 1) { return true; }
60  return false;
61  });
62  const Cut kMod2(
63  [](const caf::SRProxy* sr)
64  {
65  if(sr->hdr.ismc == false) { return true; }
66  int evtno = sr->hdr.evt;
67  if((evtno % 5) == 2) { return true; }
68  return false;
69  });
70  const Cut kMod3(
71  [](const caf::SRProxy* sr)
72  {
73  if(sr->hdr.ismc == false) { return true; }
74  int evtno = sr->hdr.evt;
75  if((evtno % 5) == 3) { return true; }
76  return false;
77  });
78  const Cut kMod4(
79  [](const caf::SRProxy* sr)
80  {
81  if(sr->hdr.ismc == false) { return true; }
82  int evtno = sr->hdr.evt;
83  if((evtno % 5) == 4) { return true; }
84  return false;
85  });
86 
87  // A map of systematic objects, indexed by the same label as above
88  std::map<std::string, Cut*> systematics;
89  systematics["Set1"] = new Cut(kMod1);
90  systematics["Set2"] = new Cut(kMod2);
91  systematics["Set3"] = new Cut(kMod3);
92  systematics["Set4"] = new Cut(kMod4);
93 
94  // Define a map of samples and cuts/axes to set them up
95  // Each sample, essentially, is a way to set up the analysis
96  std::string labelRecoE = "Calorimetric Energy (GeV)";
97  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
98  std::pair<Cut*, Cut*> Ana01(new Cut(kNusND), new Cut(kNusFD));
99  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
101  Ana01
102  );
103 
104  // Set up maps to the decompositions/predictions (each flavor component)
105  // Nominal maps are indexed purely by sample label
106  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
107  std::map<std::string, IDecomp*> decomps_nominal;
108  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
109  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
110  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
111  std::map<std::string, PredictionSterile*> predsSt_nominal;
112  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
113 
114  // Set up the actual decompositions/predictions
115  // sample.first = the sample label
116  // shiftlabel.first = the systematic shift label
117  // sample.second.first = the sample HistAxis
118  // sample.second.second = the sample Cut
119  for(const auto& sample : cut_samples) {
120  // Set up the nominal decompositions
121  // Only one is needed per sample
122  decomps_nominal[sample.first] = new ProportionalDecomp(
123  loaders,
124  *sample.second.first, *sample.second.second.first && kMod0,
126  );
127  predsNE_nominal[sample.first] = new PredictionNoExtrap(
128  loaders,
129  *sample.second.first, *sample.second.second.second && kMod0,
131  );
132 
133  ProportionalDecomp* decompnc_n = new ProportionalDecomp(
134  loaders,
135  *sample.second.first, *sample.second.second.first && kMod0,
137  );
138  ProportionalDecomp* decompnumu_n = new ProportionalDecomp(
139  loaders,
140  kNCBinsNumuCCAxis, kNumuND && kMod0,
142  );
144  loaders, *decompnc_n, *decompnumu_n,
145  *sample.second.first, kNCBinsNumuCCAxis,
146  *sample.second.second.second && kMod0,
147  *sample.second.second.first && kMod0,
148  kNumuND && kMod0,
150  ));
151  predsSt_nominal[sample.first] = new PredictionSterile(extrap_n);
152 
153  // Set up the shifted decompositions
154  // One is needed for every sample, systematic, and sigma level
155  for(const auto& systlabel : systematics) {
156  int ndpos = 0;
157  int fdpos = 1;
158  int fspos = 2;
159  int ftpos = 3;
160 
161  decomps_shifted[sample.first][systlabel.first][1] = new ProportionalDecomp(
162  loaders,
163  *sample.second.first, *sample.second.second.first && *systlabel.second,
165  );
166  predsNE_shifted[sample.first][systlabel.first][1] = new PredictionNoExtrap(
167  loaders,
168  *sample.second.first, *sample.second.second.second && *systlabel.second,
170  );
171 
172  ProportionalDecomp* decompnc_s = new ProportionalDecomp(
173  loaders,
174  *sample.second.first, *sample.second.second.first && *systlabel.second,
176  );
177  ProportionalDecomp* decompnumu_s = new ProportionalDecomp(
178  loaders,
179  kNCBinsNumuCCAxis, kNumuND && *systlabel.second,
181  );
183  loaders, *decompnc_s, *decompnumu_s,
184  *sample.second.first, kNCBinsNumuCCAxis,
185  *sample.second.second.second && *systlabel.second,
186  *sample.second.second.first && *systlabel.second,
187  kNumuND && *systlabel.second,
189  ));
190  predsSt_shifted[sample.first][systlabel.first][1] = new PredictionSterile(
191  extrap_s
192  );
193  }
194  }
195 
196  loaders.Go();
197 
198  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
199  std::string filenm = "SystsMCStats";
200 
201  std::string fullLocation = folder + filenm + ".root";
202  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
203 
204  SaveMaps(rootF,
205  decomps_nominal, decomps_shifted,
206  predsNE_nominal, predsNE_shifted,
207  predsSt_nominal, predsSt_shifted
208  );
209 }
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
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
const Cut kNusFD
Definition: NusCuts.h:46
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
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
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
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
caf::StandardRecord * sr
void SystsMCStatsLoad()
caf::Proxy< unsigned int > evt
Definition: SRProxy.h:237
const SystShifts kNoShift
Definition: SystShifts.cxx:21
const HistAxis kNCBinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNCDisappearanceEnergyBinning, kCCE)
Definition: HistAxes.h:9
Splits Data proportionally according to MC.
const Var kTuftsWeightCC
Definition: XsecTunes.h:31
std::vector< Loaders * > loaders
Definition: syst_header.h:386
caf::Proxy< bool > ismc
Definition: SRProxy.h:242
bool systematics
Definition: fnexvscaf.C:31
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
static ModularExtrapSterile NCDisappearance(Loaders &loaders, const IDecomp &NCSurvDecomp, const IDecomp &NumuOscDecomp, const HistAxis &axis, const HistAxis &axisNumuND, const Cut &fdcut, const Cut &NCNDcut, const Cut &NumuNDcut, const SystShifts &shiftMC=kNoShift, const Var &weight=kUnweighted)
Creates a NC disappearance extrapolation.
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