SystsFidContLoad.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"
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  // A map to set up systematics, starting with a label for the systematic,
35  // and finishing with a vector for items necessary to set up the Syst object
36  std::map<std::string, std::string> shiftlabels;
37  shiftlabels["FidXPos"] = "Fiducial Positive X";
38  shiftlabels["FidXNeg"] = "Fiducial Negative X";
39  shiftlabels["FidYPos"] = "Fiducial Positive Y";
40  shiftlabels["FidYNeg"] = "Fiducial Negative Y";
41  shiftlabels["FidZLow"] = "Fiducial Low Z";
42  shiftlabels["FidZHgh"] = "Fiducial High Z";
43  shiftlabels["FidInner"] = "Fiducial Inner Region";
44  shiftlabels["FidOuter"] = "Fiducial Outer Region";
45 
46 /* const Cut kXPos(
47  [](const caf::SRProxy* sr)
48  {
49  if(sr->mc.nnu < 1) { return false; }
50  double xDiv = NDL + (NDR - NDL)/2.;
51  if(sr->mc.nu[0].vtx.X() < xDiv) { return false; }
52  return true;
53  });
54 
55  const Cut kXNeg(
56  [](const caf::SRProxy* sr)
57  {
58  if(sr->mc.nnu < 1) { return false; }
59  double xDiv = NDL + (NDR - NDL)/2.;
60  if(sr->mc.nu[0].vtx.X() > xDiv) { return false; }
61  return true;
62  });
63 
64  const Cut kYPos(
65  [](const caf::SRProxy* sr)
66  {
67  if(sr->mc.nnu < 1) { return false; }
68  double yDiv = NDB + (NDT - NDB)/2.;
69  if(sr->mc.nu[0].vtx.Y() < yDiv) { return false; }
70  return true;
71  });
72 
73  const Cut kYNeg(
74  [](const caf::SRProxy* sr)
75  {
76  if(sr->mc.nnu < 1) { return false; }
77  double yDiv = NDB + (NDT - NDB)/2.;
78  if(sr->mc.nu[0].vtx.Y() > yDiv) { return false; }
79  return true;
80  });
81 
82  const Cut kZLow(
83  [](const caf::SRProxy* sr)
84  {
85  if(sr->mc.nnu < 1) { return false; }
86  double zDiv = NDF + (NDE - NDF)/2.;
87  if(sr->mc.nu[0].vtx.Z() > zDiv) { return false; }
88  return true;
89  });
90 
91  const Cut kZHgh(
92  [](const caf::SRProxy* sr)
93  {
94  if(sr->mc.nnu < 1) { return false; }
95  double zDiv = NDF + (NDE - NDF)/2.;
96  if(sr->mc.nu[0].vtx.Z() < zDiv) { return false; }
97  return true;
98  });
99 
100  const Cut kInner(
101  [](const caf::SRProxy* sr)
102  {
103  if(sr->mc.nnu < 1) { return false; }
104  double x = NDR - NDL;
105  double y = NDT - NDB;
106  double eps = (x + y)/4. - sqrt(x*x + y*y)/4.;
107  double xL = NDL + eps;
108  double xR = NDR - eps;
109  double yB = NDB + eps;
110  double yT = NDT - eps;
111  if(sr->mc.nu[0].vtx.X() < xL) { return false; }
112  if(sr->mc.nu[0].vtx.X() > xR) { return false; }
113  if(sr->mc.nu[0].vtx.Y() < yB) { return false; }
114  if(sr->mc.nu[0].vtx.Y() > yT) { return false; }
115  return true;
116  });
117 
118  const Cut kOuter(
119  [](const caf::SRProxy* sr)
120  {
121  if(sr->mc.nnu < 1) { return false; }
122  double x = NDR - NDL;
123  double y = NDT - NDB;
124  double eps = (x + y)/4. - sqrt(x*x + y*y)/4.;
125  double xL = NDL + eps;
126  double xR = NDR - eps;
127  double yB = NDB + eps;
128  double yT = NDT - eps;
129  if(sr->mc.nu[0].vtx.X() > xL &&
130  sr->mc.nu[0].vtx.X() < xR) { return false; }
131  if(sr->mc.nu[0].vtx.Y() > yB &&
132  sr->mc.nu[0].vtx.Y() < yT) { return false; }
133  return true;
134  });
135 */
136  const Cut kXPos(
137  [](const caf::SRProxy* sr)
138  {
139  if(!sr->vtx.elastic.IsValid ) { return false; }
140  if(sr->vtx.elastic.vtx.X() < 0.) { return false; }
141  return true;
142  });
143 
144  const Cut kXNeg(
145  [](const caf::SRProxy* sr)
146  {
147  if(!sr->vtx.elastic.IsValid) { return false; }
148  if(sr->vtx.elastic.vtx.X() > 0.) { return false; }
149  return true;
150  });
151 
152  const Cut kYPos(
153  [](const caf::SRProxy* sr)
154  {
155  if( !sr->vtx.elastic.IsValid ) { return false; }
156  if(sr->vtx.elastic.vtx.Y() < 0.) { return false; }
157  return true;
158  });
159 
160  const Cut kYNeg(
161  [](const caf::SRProxy* sr)
162  {
163  if(!sr->vtx.elastic.IsValid ) { return false; }
164  if(sr->vtx.elastic.vtx.Y() > 0.) { return false; }
165  return true;
166  });
167 
168  const Cut kZLow(
169  [](const caf::SRProxy* sr)
170  {
171  if(!sr->vtx.elastic.IsValid ) { return false; }
172  if(sr->vtx.elastic.vtx.Z() > 600.) { return false; }
173  return true;
174  });
175 
176  const Cut kZHgh(
177  [](const caf::SRProxy* sr)
178  {
179  if(!sr->vtx.elastic.IsValid ) { return false; }
180  if(sr->vtx.elastic.vtx.Z() < 600.) { return false; }
181  return true;
182  });
183 
184  const Cut kInner(
185  [](const caf::SRProxy* sr)
186  {
187  if(!sr->vtx.elastic.IsValid ) { return false; }
188  double rDiv = 100./sqrt(2.);
189  if(fabs(sr->vtx.elastic.vtx.X()) > rDiv) { return false; }
190  if(fabs(sr->vtx.elastic.vtx.Y()) > rDiv) { return false; }
191  return true;
192  });
193 
194  const Cut kOuter(
195  [](const caf::SRProxy* sr)
196  {
197  if( !sr->vtx.elastic.IsValid ) { return false; }
198  double rDiv = 100./sqrt(2.);
199  if(fabs(sr->vtx.elastic.vtx.X()) < rDiv &&
200  fabs(sr->vtx.elastic.vtx.Y()) < rDiv) { return false; }
201  return true;
202  });
203 
204  // A map of systematic objects, indexed by the same label as above
205  std::map<std::string, Cut*> systematics;
206  systematics["FidXPos"] = new Cut(kXPos);
207  systematics["FidXNeg"] = new Cut(kXNeg);
208  systematics["FidYPos"] = new Cut(kYPos);
209  systematics["FidYNeg"] = new Cut(kYNeg);
210  systematics["FidZLow"] = new Cut(kZLow);
211  systematics["FidZHgh"] = new Cut(kZHgh);
212  systematics["FidInner"] = new Cut(kInner);
213  systematics["FidOuter"] = new Cut(kOuter);
214 
215  // Define a map of samples and cuts/axes to set them up
216  // Each sample, essentially, is a way to set up the analysis
217  std::string labelRecoE = "Calorimetric Energy (GeV)";
218  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
219  std::pair<Cut*, Cut*> Ana01(new Cut(kNusND), new Cut(kNusFD));
220  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
222  Ana01
223  );
224 
225  // Set up maps to the decompositions/predictions (each flavor component)
226  // Nominal maps are indexed purely by sample label
227  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
228  std::map<std::string, IDecomp*> decomps_nominal;
229  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
230  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
231  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
232  std::map<std::string, PredictionSterile*> predsSt_nominal;
233  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
234 
235  // Set up the actual decompositions/predictions
236  // sample.first = the sample label
237  // shiftlabel.first = the systematic shift label
238  // sample.second.first = the sample HistAxis
239  // sample.second.second = the sample Cut
240  // syst.first = the integer number of sigmas for the shift
241  // syst.second = the systematic being shifted
242  for(const auto& sample : cut_samples) {
243  SterileGenerator gen_n(
244  *sample.second.first, kNCBinsNumuCCAxis,
245  *sample.second.second.second, *sample.second.second.first,
247  );
248 
249  // Set up the nominal decompositions
250  // Only one is needed per sample
251  decomps_nominal[sample.first] = new ProportionalDecomp(
252  loaders,
253  *sample.second.first, *sample.second.second.first,
255  );
256  predsNE_nominal[sample.first] = new PredictionNoExtrap(
257  loaders,
258  *sample.second.first, *sample.second.second.second,
260  );
261  predsSt_nominal[sample.first] =
262  (PredictionSterile*)gen_n.Generate(loaders).release();
263 
264  // Set up the shifted decompositions
265  // One is needed for every sample, systematic, and sigma level
266  for(const auto& shiftlabel : shiftlabels) {
267  decomps_shifted[sample.first][shiftlabel.first][1] = new ProportionalDecomp(
268  loaders,
269  *sample.second.first, *sample.second.second.first && *systematics[shiftlabel.first],
271  );
272  predsNE_shifted[sample.first][shiftlabel.first][1] = new PredictionNoExtrap(
273  loaders,
274  *sample.second.first, *sample.second.second.second,
276  );
277 
278  SterileGenerator gen_s(
279  *sample.second.first, kNCBinsNumuCCAxis,
280  *sample.second.second.second,
281  *sample.second.second.first && *systematics[shiftlabel.first],
282  kNumuND && *systematics[shiftlabel.first],
284  );
285  predsSt_shifted[sample.first][shiftlabel.first][1] =
286  (PredictionSterile*)gen_s.Generate(loaders).release();
287  }
288  }
289 
290  loaders.Go();
291 
292  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
293  std::string filenm = "SystsFidCont";
294 
295  std::string fullLocation = folder + filenm + ".root";
296  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
297 
298  SaveMaps(rootF,
299  decomps_nominal, decomps_shifted,
300  predsNE_nominal, predsNE_shifted,
301  predsSt_nominal, predsSt_shifted
302  );
303 }
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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
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.
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
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
const SystShifts kNoShift
Definition: SystShifts.cxx:21
const HistAxis kNCBinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNCDisappearanceEnergyBinning, kCCE)
Definition: HistAxes.h:9
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
void SystsFidContLoad()
Splits Data proportionally according to MC.
const Var kTuftsWeightCC
Definition: XsecTunes.h:31
std::vector< Loaders * > loaders
Definition: syst_header.h:386
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
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.
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
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
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