SystsGENIELoad.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 
24 
29  //loaders.SetLoaderPath(fnameneardata_concat, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);
30 
32 
33  // Vector of sigmas, levels to set each systematic
34  std::vector<int> sigmas{-2, -1, +1, +2};
35 
36 /*
37 redundant
38 "ReweightMaCCQE"
39 "ReweightMaCCRESshape"
40 "ReweightMvCCRESshape"
41 "ReweightMaNCRESshape"
42 "ReweightMvNCRESshape"
43 "ReweightAhtBYshape"
44 "ReweightBhtBYshape"
45 "ReweightCV1uBYshape"
46 "ReweightCV2uBYshape"
47 "ReweightNormCCQEenu"
48 "ReweightNormCCRES"
49 "ReweightNormNCRES"
50 "ReweightNC"
51 */
52 
53  // A map to set up systematics, starting with a label for the systematic,
54  // and finishing with a vector for items necessary to set up the Syst object
55  std::map<std::string, rwgt::ReweightLabel_t> shiftlabels;
56  // NCEL parameters:
57  shiftlabels["ReweightMaNCEL"] = rwgt::fReweightMaNCEL; // Ma NCEL, affects dsigma(NCEL)/dQ2 both in shape and normalization
58  shiftlabels["ReweightEtaNCEL"] = rwgt::fReweightEtaNCEL; // NCEL strange axial form factor eta, affects dsigma(NCEL)/dQ2 both in shape and normalization
59  // CCQE parameters:
60  //shiftlabels["ReweightMaCCQE"] = rwgt::fReweightMaCCQE; // Ma CCQE, affects dsigma(CCQE)/dQ2 both in shape and normalization
61  shiftlabels["ReweightMaCCQEshape"] = rwgt::fReweightMaCCQEshape; // Ma CCQE, affects dsigma(CCQE)/dQ2 in shape only (normalized to constant integral)
62  shiftlabels["ReweightNormCCQE"] = rwgt::fReweightNormCCQE; // CCQE normalization (energy independent)
63  //shiftlabels["ReweightNormCCQEenu"] = rwgt::fReweightNormCCQEenu; // CCQE normalization (maintains dependence on neutrino energy)
64  shiftlabels["ReweightCCQEPauliSupViaKF"] = rwgt::fReweightCCQEPauliSupViaKF; //
65  shiftlabels["ReweightVecCCQEshape"] = rwgt::fReweightVecCCQEshape; // elastic nucleon form factors (BBA/default -> dipole) - shape only effect of dsigma(CCQE)/dQ2
66  shiftlabels["ReweightCCQEMomDistroFGtoSF"] = rwgt::fReweightCCQEMomDistroFGtoSF; //
67  // Resonance neutrino-production parameters:
68  //shiftlabels["ReweightNormCCRES"] = rwgt::fReweightNormCCRES; // CCRES normalization
69  //shiftlabels["ReweightNormNCRES"] = rwgt::fReweightNormNCRES; // NCRES normalization
70  //shiftlabels["ReweightMaCCRESshape"] = rwgt::fReweightMaCCRESshape; // Ma CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
71  //shiftlabels["ReweightMvCCRESshape"] = rwgt::fReweightMvCCRESshape; // Mv CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
72  //shiftlabels["ReweightMaNCRESshape"] = rwgt::fReweightMaNCRESshape; // Ma NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
73  //shiftlabels["ReweightMvNCRESshape"] = rwgt::fReweightMvNCRESshape; // Mv NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
74  shiftlabels["ReweightMaCCRES"] = rwgt::fReweightMaCCRES; // Ma CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
75  shiftlabels["ReweightMvCCRES"] = rwgt::fReweightMvCCRES; // Mv CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
76  shiftlabels["ReweightMaNCRES"] = rwgt::fReweightMaNCRES; // Ma NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
77  shiftlabels["ReweightMvNCRES"] = rwgt::fReweightMvNCRES; // Mv NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
78  // Coherent pion production parameters:
79  shiftlabels["ReweightMaCOHpi"] = rwgt::fReweightMaCOHpi; // Ma for COH pion production
80  shiftlabels["ReweightR0COHpi"] = rwgt::fReweightR0COHpi; // R0 for COH pion production
81  // Non-resonance background parameters:
82  shiftlabels["ReweightRvpCC1pi"] = rwgt::fReweightRvpCC1pi; // the 1pi non-RES bkg in the RES region, for v+p CC
83  shiftlabels["ReweightRvpCC2pi"] = rwgt::fReweightRvpCC2pi; // the 2pi non-RES bkg in the RES region, for v+p CC
84  shiftlabels["ReweightRvnCC1pi"] = rwgt::fReweightRvnCC1pi; // the 1pi non-RES bkg in the RES region, for v+n CC
85  shiftlabels["ReweightRvnCC2pi"] = rwgt::fReweightRvnCC2pi; // the 2pi non-RES bkg in the RES region, for v+n CC
86  shiftlabels["ReweightRvpNC1pi"] = rwgt::fReweightRvpNC1pi; // the 1pi non-RES bkg in the RES region, for v+p NC
87  shiftlabels["ReweightRvpNC2pi"] = rwgt::fReweightRvpNC2pi; // the 2pi non-RES bkg in the RES region, for v+p NC
88  shiftlabels["ReweightRvnNC1pi"] = rwgt::fReweightRvnNC1pi; // the 1pi non-RES bkg in the RES region, for v+n NC
89  shiftlabels["ReweightRvnNC2pi"] = rwgt::fReweightRvnNC2pi; // the 2pi non-RES bkg in the RES region, for v+n NC
90  shiftlabels["ReweightRvbarpCC1pi"] = rwgt::fReweightRvbarpCC1pi; // the 1pi non-RES bkg in the RES region, for vbar+p CC
91  shiftlabels["ReweightRvbarpCC2pi"] = rwgt::fReweightRvbarpCC2pi; // the 2pi non-RES bkg in the RES region, for vbar+p CC
92  shiftlabels["ReweightRvbarnCC1pi"] = rwgt::fReweightRvbarnCC1pi; // the 1pi non-RES bkg in the RES region, for vbar+n CC
93  shiftlabels["ReweightRvbarnCC2pi"] = rwgt::fReweightRvbarnCC2pi; // the 2pi non-RES bkg in the RES region, for vbar+n CC
94  shiftlabels["ReweightRvbarpNC1pi"] = rwgt::fReweightRvbarpNC1pi; // the 1pi non-RES bkg in the RES region, for vbar+p NC
95  shiftlabels["ReweightRvbarpNC2pi"] = rwgt::fReweightRvbarpNC2pi; // the 2pi non-RES bkg in the RES region, for vbar+p NC
96  shiftlabels["ReweightRvbarnNC1pi"] = rwgt::fReweightRvbarnNC1pi; // the 1pi non-RES bkg in the RES region, for vbar+n NC
97  shiftlabels["ReweightRvbarnNC2pi"] = rwgt::fReweightRvbarnNC2pi; // the 2pi non-RES bkg in the RES region, for vbar+n NC
98  // DIS parameters - applied for DIS events with (Q2>Q2o, W>Wo),
99  // typically Q2okReweight =1GeV^2, WokReweight =1.7-2.0GeV
100  shiftlabels["ReweightAhtBY"] = rwgt::fReweightAhtBY; // the Bodek-Yang model parameter A_{ht} - incl. both shape and normalization effect
101  shiftlabels["ReweightBhtBY"] = rwgt::fReweightBhtBY; // the Bodek-Yang model parameter B_{ht} - incl. both shape and normalization effect
102  shiftlabels["ReweightCV1uBY"] = rwgt::fReweightCV1uBY; // the Bodek-Yang model parameter CV1u - incl. both shape and normalization effect
103  shiftlabels["ReweightCV2uBY"] = rwgt::fReweightCV2uBY; // the Bodek-Yang model parameter CV2u - incl. both shape and normalization effect
104  //shiftlabels["ReweightAhtBYshape"] = rwgt::fReweightAhtBYshape; // the Bodek-Yang model parameter A_{ht} - shape only effect to d2sigma(DIS)/dxdy
105  //shiftlabels["ReweightBhtBYshape"] = rwgt::fReweightBhtBYshape; // the Bodek-Yang model parameter B_{ht} - shape only effect to d2sigma(DIS)/dxdy
106  //shiftlabels["ReweightCV1uBYshape"] = rwgt::fReweightCV1uBYshape; // the Bodek-Yang model parameter CV1u - shape only effect to d2sigma(DIS)/dxdy
107  //shiftlabels["ReweightCV2uBYshape"] = rwgt::fReweightCV2uBYshape; // the Bodek-Yang model parameter CV2u - shape only effect to d2sigma(DIS)/dxdy
108  shiftlabels["ReweightNormDISCC"] = rwgt::fReweightNormDISCC; // the inclusive DIS CC normalization (not currently working in genie)
109  shiftlabels["ReweightRnubarnuCC"] = rwgt::fReweightRnubarnuCC; // the ratio of \sigma(\bar\nu CC) / \sigma(\nu CC) (not currently working in genie)
110  shiftlabels["ReweightDISNuclMod"] = rwgt::fReweightDISNuclMod; // DIS nuclear modification (shadowing, anti-shadowing, EMC). Not working in GENIE at the moment
111 
112  // Hadronization (free nucleon target)
113  shiftlabels["ReweightAGKY_xF1pi"] = rwgt::fReweightAGKY_xF1pi; // xF distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
114  shiftlabels["ReweightAGKY_pT1pi"] = rwgt::fReweightAGKY_pT1pi; // pT distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
115  shiftlabels["ReweightFormZone"] = rwgt::fReweightFormZone; // Hadron formation zone
116  // Resonance decays
117  shiftlabels["ReweightTheta_Delta2Npi"] = rwgt::fReweightTheta_Delta2Npi; // distort pi angular distribution in Delta -> N + pi
118  shiftlabels["ReweightBR1gamma"] = rwgt::fReweightBR1gamma; // Resonance -> X + gamma branching ratio, eg Delta+(1232) -> p gamma
119  shiftlabels["ReweightBR1eta"] = rwgt::fReweightBR1eta; // Resonance -> X + eta branching ratio, eg N+(1440) -> p eta
120 
121  // Intranuclear rescattering systematics. There are 2 sets of parameters:
122  // - parameters that control the total rescattering probability, P(total)
123  // - parameters that control the fraction of each process (`fate'), given a total rescat. prob., P(fate|total)
124  // These parameters are considered separately for pions and nucleons.
125  shiftlabels["ReweightMFP_N"] = rwgt::fReweightMFP_N; // mean free path for nucleons
126  shiftlabels["ReweightFrCEx_N"] = rwgt::fReweightFrCEx_N; // charge exchange probability for nucleons, for given total rescattering probability
127  shiftlabels["ReweightFrElas_N"] = rwgt::fReweightFrElas_N; // elastic probability for nucleons, for given total rescattering probability
128  shiftlabels["ReweightFrInel_N"] = rwgt::fReweightFrInel_N; // inelastic probability for nucleons, for given total rescattering probability
129  shiftlabels["ReweightFrAbs_N"] = rwgt::fReweightFrAbs_N; // absorption probability for nucleons, for given total rescattering probability
130  shiftlabels["ReweightFrPiProd_N"] = rwgt::fReweightFrPiProd_N; // pion production probability for nucleons, for given total rescattering probability
131  shiftlabels["ReweightMFP_pi"] = rwgt::fReweightMFP_pi; // mean free path for pions
132  shiftlabels["ReweightFrCEx_pi"] = rwgt::fReweightFrCEx_pi; // charge exchange probability for pions, for given total rescattering probability
133  shiftlabels["ReweightFrElas_pi"] = rwgt::fReweightFrElas_pi; // elastic probability for pions, for given total rescattering probability
134  shiftlabels["ReweightFrInel_pi"] = rwgt::fReweightFrInel_pi; // inelastic probability for pions, for given total rescattering probability
135  shiftlabels["ReweightFrAbs_pi"] = rwgt::fReweightFrAbs_pi; // absorption probability for pions, for given total rescattering probability
136  shiftlabels["ReweightFrPiProd_pi"] = rwgt::fReweightFrPiProd_pi; // pion production probability for pions, for given total rescattering probability
137  // Nuclear model
138 
139  //shiftlabels["ReweightNC"] = rwgt::fReweightNC; //
140 
141  // A map of systematic objects, indexed by the same label as above
142  std::map<std::string, std::map<int, SystShifts*> > systematics;
143 
144  for(const auto& shiftlabel : shiftlabels){
145  for(const auto& sigma : sigmas) {
146  systematics[shiftlabel.first][sigma] = new SystShifts(
147  new GenieRwgtTableSyst(shiftlabel.second), sigma
148  );
149  }
150  }
151 
152  // These two shifts are GENIE systs, but not in the form of reweights as above
153  // Unfortunately they need a rwgt label, so give them the null one.
154  shiftlabels["mecScale"] = rwgt::kReweightNull;
155  shiftlabels["RPA"] = rwgt::kReweightNull;
156 
157  for(const auto& sigma : sigmas) {
158  systematics["mecScale"][sigma] = new SystShifts(new MECScaleSystSA(), sigma);
159  systematics["RPA"][sigma] = new SystShifts(new RPACCQESystSA(), sigma);
160  }
161 
162  // Define a map of samples and cuts/axes to set them up
163  // Each sample, essentially, is a way to set up the analysis
164  std::string labelRecoE = "Calorimetric Energy (GeV)";
165  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
166  std::pair<Cut*, Cut*> Ana01(new Cut(kNusND), new Cut(kNusFD));
167  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
169  Ana01
170  );
171 
172  // Cuts to split MC into two groups
173  const Cut kMod0(
174  [](const caf::SRProxy* sr)
175  {
176  if(sr->hdr.ismc == false) { return true; }
177  int evtno = sr->hdr.evt;
178  if((evtno % 2) == 0) { return true; }
179  return false;
180  });
181  const Cut kMod1(
182  [](const caf::SRProxy* sr)
183  {
184  if(sr->hdr.ismc == false) { return true; }
185  int evtno = sr->hdr.evt;
186  if((evtno % 2) == 1) { return true; }
187  return false;
188  });
189 
190  // Set up maps to the decompositions/predictions (each flavor component)
191  // Nominal maps are indexed purely by sample label
192  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
193  std::map<std::string, IDecomp*> decomps_nominal;
194  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
195  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
196  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
197  std::map<std::string, PredictionSterile*> predsSt_nominal;
198  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
199 
200  // Set up the actual decompositions/predictions
201  // sample.first = the sample label
202  // shiftlabel.first = the systematic shift label
203  // sample.second.first = the sample HistAxis
204  // sample.second.second = the sample Cut
205  // syst.first = the integer number of sigmas for the shift
206  // syst.second = the systematic being shifted
207  for(const auto& sample : cut_samples) {
208  // Set up the nominal decompositions
209  // Only one is needed per sample
210  decomps_nominal[sample.first] = new CheatDecomp(
212  *sample.second.first, *sample.second.second.first && kMod0,
214  );
215  predsNE_nominal[sample.first] = new PredictionNoExtrap(
216  loaders,
217  *sample.second.first, *sample.second.second.second,
218  kNoShift, kTuftsWeightCC
219  );
220 
221  CheatDecomp* decompnc_n = new CheatDecomp(
223  *sample.second.first, *sample.second.second.first && kMod1,
225  );
226  CheatDecomp* decompnumu_n = new CheatDecomp(
228  kNCBinsNumuCCAxis, kNumuND && kMod1,
230  );
232  loaders, *decompnc_n, *decompnumu_n,
233  *sample.second.first, kNCBinsNumuCCAxis,
234  *sample.second.second.second, *sample.second.second.first && kMod0, kNumuND && kMod0,
235  kNoShift, kTuftsWeightCC
236  ));
237  predsSt_nominal[sample.first] = new PredictionSterile(extrap_n);
238 
239  // Set up the shifted decompositions
240  // One is needed for every sample, systematic, and sigma level
241  for(const auto& shiftlabel : shiftlabels) {
242  for(const auto& syst : systematics[shiftlabel.first]) {
243  decomps_shifted[sample.first][shiftlabel.first][syst.first] = new CheatDecomp(
245  *sample.second.first, *sample.second.second.first && kMod0,
246  *syst.second, kTuftsWeightCC
247  );
248  predsNE_shifted[sample.first][shiftlabel.first][syst.first] = new PredictionNoExtrap(
249  loaders,
250  *sample.second.first, *sample.second.second.second,
251  *syst.second, kTuftsWeightCC
252  );
253 
254  CheatDecomp* decompnc_s = new CheatDecomp(
256  *sample.second.first, *sample.second.second.first && kMod1,
258  );
259  CheatDecomp* decompnumu_s = new CheatDecomp(
261  kNCBinsNumuCCAxis, kNumuND && kMod1,
263  );
265  loaders, *decompnc_s, *decompnumu_s,
266  *sample.second.first, kNCBinsNumuCCAxis,
267  *sample.second.second.second, *sample.second.second.first && kMod0, kNumuND && kMod0,
268  *syst.second, kTuftsWeightCC
269  ));
270  predsSt_shifted[sample.first][shiftlabel.first][syst.first] = new PredictionSterile(
271  extrap_s
272  );
273  }
274  }
275  }
276 
277  loaders.Go();
278 
279  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
280  std::string filenm = "SystsGENIE";
281 
282  std::string fullLocation = folder + filenm + ".root";
283  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
284 
285  SaveMaps(rootF,
286  decomps_nominal, decomps_shifted,
287  predsNE_nominal, predsNE_shifted,
288  predsSt_nominal, predsSt_shifted
289  );
290 }
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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 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
caf::StandardRecord * sr
double sigma(TH1F *hist, double percentile)
caf::Proxy< unsigned int > evt
Definition: SRProxy.h:237
const SystShifts kNoShift
Definition: SystShifts.cxx:21
void SystsGENIELoad()
const HistAxis kNCBinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNCDisappearanceEnergyBinning, kCCE)
Definition: HistAxes.h:9
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.
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