SystsGENIEAna.C
Go to the documentation of this file.
2 #include "CAFAna/Systs/Systs.h"
5 
6 using namespace ana;
7 
9 {
10  TH1::AddDirectory(0);
11 
12  // A map of systematic labels
13  std::map<std::string, rwgt::ReweightLabel_t> shiftlabels; // Easier to get label information like this...
14  // NCEL parameters:
15  shiftlabels["ReweightMaNCEL"] = rwgt::fReweightMaNCEL;
16  shiftlabels["ReweightEtaNCEL"] = rwgt::fReweightEtaNCEL;
17  //shiftlabels["ReweightMaCCQE"] = rwgt::fReweightMaCCQE;
18  shiftlabels["ReweightMaCCQEshape"] = rwgt::fReweightMaCCQEshape;
19  shiftlabels["ReweightNormCCQE"] = rwgt::fReweightNormCCQE;
20  //shiftlabels["ReweightNormCCQEenu"] = rwgt::fReweightNormCCQEenu;
21  shiftlabels["ReweightCCQEPauliSupViaKF"] = rwgt::fReweightCCQEPauliSupViaKF;
22  shiftlabels["ReweightVecCCQEshape"] = rwgt::fReweightVecCCQEshape;
23  shiftlabels["ReweightCCQEMomDistroFGtoSF"] = rwgt::fReweightCCQEMomDistroFGtoSF;
24  //shiftlabels["ReweightNormCCRES"] = rwgt::fReweightNormCCRES;
25  //shiftlabels["ReweightNormNCRES"] = rwgt::fReweightNormNCRES;
26  //shiftlabels["ReweightMaCCRESshape"] = rwgt::fReweightMaCCRESshape;
27  //shiftlabels["ReweightMvCCRESshape"] = rwgt::fReweightMvCCRESshape;
28  //shiftlabels["ReweightMaNCRESshape"] = rwgt::fReweightMaNCRESshape;
29  //shiftlabels["ReweightMvNCRESshape"] = rwgt::fReweightMvNCRESshape;
30  shiftlabels["ReweightMaCCRES"] = rwgt::fReweightMaCCRES;
31  shiftlabels["ReweightMvCCRES"] = rwgt::fReweightMvCCRES;
32  shiftlabels["ReweightMaNCRES"] = rwgt::fReweightMaNCRES;
33  shiftlabels["ReweightMvNCRES"] = rwgt::fReweightMvNCRES;
34  shiftlabels["ReweightMaCOHpi"] = rwgt::fReweightMaCOHpi;
35  shiftlabels["ReweightR0COHpi"] = rwgt::fReweightR0COHpi;
36  shiftlabels["ReweightRvpCC1pi"] = rwgt::fReweightRvpCC1pi;
37  shiftlabels["ReweightRvpCC2pi"] = rwgt::fReweightRvpCC2pi;
38  shiftlabels["ReweightRvnCC1pi"] = rwgt::fReweightRvnCC1pi;
39  shiftlabels["ReweightRvnCC2pi"] = rwgt::fReweightRvnCC2pi;
40  shiftlabels["ReweightRvpNC1pi"] = rwgt::fReweightRvpNC1pi;
41  shiftlabels["ReweightRvpNC2pi"] = rwgt::fReweightRvpNC2pi;
42  shiftlabels["ReweightRvnNC1pi"] = rwgt::fReweightRvnNC1pi;
43  shiftlabels["ReweightRvnNC2pi"] = rwgt::fReweightRvnNC2pi;
44  shiftlabels["ReweightRvbarpCC1pi"] = rwgt::fReweightRvbarpCC1pi;
45  shiftlabels["ReweightRvbarpCC2pi"] = rwgt::fReweightRvbarpCC2pi;
46  shiftlabels["ReweightRvbarnCC1pi"] = rwgt::fReweightRvbarnCC1pi;
47  shiftlabels["ReweightRvbarnCC2pi"] = rwgt::fReweightRvbarnCC2pi;
48  shiftlabels["ReweightRvbarpNC1pi"] = rwgt::fReweightRvbarpNC1pi;
49  shiftlabels["ReweightRvbarpNC2pi"] = rwgt::fReweightRvbarpNC2pi;
50  shiftlabels["ReweightRvbarnNC1pi"] = rwgt::fReweightRvbarnNC1pi;
51  shiftlabels["ReweightRvbarnNC2pi"] = rwgt::fReweightRvbarnNC2pi;
52  shiftlabels["ReweightAhtBY"] = rwgt::fReweightAhtBY;
53  shiftlabels["ReweightBhtBY"] = rwgt::fReweightBhtBY;
54  shiftlabels["ReweightCV1uBY"] = rwgt::fReweightCV1uBY;
55  shiftlabels["ReweightCV2uBY"] = rwgt::fReweightCV2uBY;
56  //shiftlabels["ReweightAhtBYshape"] = rwgt::fReweightAhtBYshape;
57  //shiftlabels["ReweightBhtBYshape"] = rwgt::fReweightBhtBYshape;
58  //shiftlabels["ReweightCV1uBYshape"] = rwgt::fReweightCV1uBYshape;
59  //shiftlabels["ReweightCV2uBYshape"] = rwgt::fReweightCV2uBYshape;
60  shiftlabels["ReweightNormDISCC"] = rwgt::fReweightNormDISCC;
61  shiftlabels["ReweightRnubarnuCC"] = rwgt::fReweightRnubarnuCC;
62  shiftlabels["ReweightDISNuclMod"] = rwgt::fReweightDISNuclMod;
63  shiftlabels["ReweightAGKY_xF1pi"] = rwgt::fReweightAGKY_xF1pi;
64  shiftlabels["ReweightAGKY_pT1pi"] = rwgt::fReweightAGKY_pT1pi;
65  shiftlabels["ReweightFormZone"] = rwgt::fReweightFormZone;
66  shiftlabels["ReweightTheta_Delta2Npi"] = rwgt::fReweightTheta_Delta2Npi;
67  shiftlabels["ReweightBR1gamma"] = rwgt::fReweightBR1gamma;
68  shiftlabels["ReweightBR1eta"] = rwgt::fReweightBR1eta;
69  shiftlabels["ReweightMFP_N"] = rwgt::fReweightMFP_N;
70  shiftlabels["ReweightFrCEx_N"] = rwgt::fReweightFrCEx_N;
71  shiftlabels["ReweightFrElas_N"] = rwgt::fReweightFrElas_N;
72  shiftlabels["ReweightFrInel_N"] = rwgt::fReweightFrInel_N;
73  shiftlabels["ReweightFrAbs_N"] = rwgt::fReweightFrAbs_N;
74  shiftlabels["ReweightFrPiProd_N"] = rwgt::fReweightFrPiProd_N;
75  shiftlabels["ReweightMFP_pi"] = rwgt::fReweightMFP_pi;
76  shiftlabels["ReweightFrCEx_pi"] = rwgt::fReweightFrCEx_pi;
77  shiftlabels["ReweightFrElas_pi"] = rwgt::fReweightFrElas_pi;
78  shiftlabels["ReweightFrInel_pi"] = rwgt::fReweightFrInel_pi;
79  shiftlabels["ReweightFrAbs_pi"] = rwgt::fReweightFrAbs_pi;
80  shiftlabels["ReweightFrPiProd_pi"] = rwgt::fReweightFrPiProd_pi;
81  //shiftlabels["ReweightNC"] = rwgt::fReweightNC;
82  // These two shifts are GENIE systs, but not in the form of reweights as above
83  // Unfortunately they need a rwgt label, so give them the null one.
84  shiftlabels["mecScale"] = rwgt::kReweightNull;
85  shiftlabels["RPA"] = rwgt::kReweightNull;
86 
87  std::string labelRecoE = "Calorimetric Energy (GeV)";
88  std::vector<std::string> cut_samples;
89  cut_samples.push_back("Ana01");
90 
91  // Set up maps to the decompositions/predictions (each flavor component)
92  // Nominal maps are indexed purely by sample label
93  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
94  std::map<std::string, IDecomp*> decomps_nominal;
95  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
96  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
97  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
98  std::map<std::string, PredictionSterile*> predsSt_nominal;
99  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
100 
101  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
102  std::string filenm = "SystsGENIE";
103 
104  std::string loadLocation = folder + filenm + ".root";
105  std::string saveLocation = folder + filenm + "Ana.root";
106  std::string textLocation = folder + filenm + "Ana.txt";
107 
108  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
109  LoadMaps(rootL, &decomps_nominal, &decomps_shifted);
110  LoadMaps(rootL, &predsNE_nominal, &predsNE_shifted);
111  LoadMaps(rootL, &predsSt_nominal, &predsSt_shifted);
112  rootL->Close();
113 
114  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
115  FILE* textF;
116  textF = fopen(textLocation.c_str(), "w");
117 
118  for(const auto& sample : cut_samples) {
119  strings strs;
120  strs.fComponent = "NC";
121  strs.fDet = "ND";
122  strs.fPOT = StringFromDouble(NCSCALE);
123  strs.fSample = sample;
124  strs.fSystType = "GENIE";
125  strs.fXLabel = labelRecoE;
126 
127  InitializeSystText(textF, strs);
128 
129  TH1* hNomNC_decomp = GetNC(decomps_nominal[sample], NCSCALE);
130  TH1* hNomBG_decomp = GetBG(decomps_nominal[sample], NCSCALE);
131 
132  TH1* hNomNC_predNE = GetNC(predsNE_nominal[sample], NCSCALE);
133  TH1* hNomBG_predNE = GetBG(predsNE_nominal[sample], NCSCALE);
134 
135  TH1* hNomNC_predSt = GetNC(predsSt_nominal[sample], NCSCALE);
136  TH1* hNomBG_predSt = GetBG(predsSt_nominal[sample], NCSCALE);
137 
138  for(const auto& shifts : shiftlabels) {
139  if(shifts.first.compare("RPA") != 0 &&
140  shifts.first.compare("mecScale") != 0) {
141  ISyst* temp = new GenieRwgtTableSyst(shifts.second);
142  strs.fSystS = temp->ShortName();
143  strs.fSystL = temp->LatexName();
144  }
145  else {
146  if(shifts.first.compare("RPA") == 0) {
147  strs.fSystS = "RPA";
148  strs.fSystL = "RPA Weight";
149  }
150  else {
151  strs.fSystS = "mecScale";
152  strs.fSystL = "MEC Scale";
153  }
154  }
155 
156  strs.fComponent = "NC";
157  strs.fDet = "ND";
158  PlotSyst(hNomNC_decomp,
159  decomps_shifted[sample][shifts.first],
160  rootF, textF, strs, NCSCALE, true, false);
161 
162  strs.fComponent = "BG";
163  PlotSyst(hNomBG_decomp,
164  decomps_shifted[sample][shifts.first],
165  rootF, textF, strs, NCSCALE, false, false);
166 
167  strs.fComponent = "NC";
168  strs.fDet = "FD";
169  PlotSyst(hNomNC_predNE,
170  predsNE_shifted[sample][shifts.first],
171  rootF, textF, strs, NCSCALE, true, false);
172 
173  strs.fComponent = "BG";
174  PlotSyst(hNomBG_predNE,
175  predsNE_shifted[sample][shifts.first],
176  rootF, textF, strs, NCSCALE, false, false);
177 
178  strs.fComponent = "NC";
179  strs.fDet = "EX";
180  PlotSyst(hNomNC_predSt,
181  predsSt_shifted[sample][shifts.first],
182  rootF, textF, strs, NCSCALE, true, false);
183 
184  strs.fComponent = "BG";
185  PlotSyst(hNomBG_predSt,
186  predsSt_shifted[sample][shifts.first],
187  rootF, textF, strs, NCSCALE, false, false);
188  }
189 
190  strs.fComponent = "NC";
191  strs.fDet = "ND";
192  PlotSystBand(hNomNC_decomp, rootF, textF, strs);
193 
194  strs.fComponent = "BG";
195  PlotSystBand(hNomBG_decomp, rootF, textF, strs);
196 
197  strs.fComponent = "NC";
198  strs.fDet = "FD";
199  PlotSystBand(hNomNC_predNE, rootF, textF, strs);
200 
201  strs.fComponent = "BG";
202  PlotSystBand(hNomBG_predNE, rootF, textF, strs);
203 
204  strs.fComponent = "NC";
205  strs.fDet = "EX";
206  PlotSystBand(hNomNC_predSt, rootF, textF, strs);
207 
208  strs.fComponent = "BG";
209  PlotSystBand(hNomBG_predSt, rootF, textF, strs);
210  }
211 
212  fclose(textF);
213  rootF->Close();
214 
215  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
216 }
std::string fDet
Definition: PPFXHelper.h:104
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string fComponent
Definition: PPFXHelper.h:103
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
void PlotSyst(TH1 *h, TH1 *hp1, TH1 *hm1, TH1 *hp2, TH1 *hm2, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:662
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
TH1 * GetNC(IDecomp *specs, double POT)
Definition: PPFXHelper.h:211
fclose(fg1)
void InitializeSystText(FILE *text, strings strs)
Print some initial text about a systematic – the systematic type and sample.
Definition: PPFXHelper.h:263
std::string fSystL
Definition: PPFXHelper.h:109
std::string fPOT
Definition: PPFXHelper.h:105
OStream cout
Definition: OStream.cxx:6
void PlotSystBand(TH1 *h, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:968
virtual const std::string & LatexName() const final
The name used on plots (ROOT&#39;s TLatex syntax)
Definition: ISyst.h:30
std::string fXLabel
Definition: PPFXHelper.h:110
std::string fSystType
Definition: PPFXHelper.h:107
std::string fSystS
Definition: PPFXHelper.h:108
void LoadMaps(TDirectory *dir, std::map< std::string, IDecomp * > *nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > *shifted)
Definition: PPFXHelper.h:273
static double NCSCALE
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
void SystsGENIEAna()
Definition: SystsGENIEAna.C:8
std::string StringFromDouble(double pot)
Definition: PPFXHelper.h:128
A helper structure to contain a group of string for plotting.
Definition: PPFXHelper.h:101
std::string fSample
Definition: PPFXHelper.h:106
enum BeamMode string