Nus18SystsLightLevelAna.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // Light Level systematics for Nus18
3 //
4 // Author: Sijith Edayath
5 // Modified by Mike Wallbank (wallbank@fnal.gov)
6 /////////////////////////////////////////////////////////////////////
7 
12 #include "CAFAna/Vars/HistAxes.h"
13 #include "NuXAna/Vars/HistAxes.h"
14 
17 
18 using namespace ana;
19 
21  TH1* hNC_ND;
22  TH1* hBG_ND;
23  TH1* hNC_FD;
24  TH1* hBG_FD;
25  TH1* hNC_EX;
26  TH1* hBG_EX;
28  if (type == "NC") {
29  if (det == "ND") return hNC_ND;
30  else if (det == "FD") return hNC_FD;
31  else if (det == "EX") return hNC_EX;
32  else std::cout << "Detector must be ND, FD or EX" << std::endl;
33  } else if (type == "BG") {
34  if (det == "ND") return hBG_ND;
35  else if (det == "FD") return hBG_FD;
36  else if (det == "EX") return hBG_EX;
37  else std::cout << "Detector must be ND, FD or EX" << std::endl;
38  } else
39  std::cout << "Type must be NC or BG" << std::endl;
40  return nullptr;
41  }
42 };
43 
45 
46  TH1::AddDirectory(0);
47 
48  // Set up path to file with filled CAFAna objects and for output of analysis
49  std::string fileName = inFile.substr(0, inFile.find(".root"));
50  std::string saveName = fileName + "_ana.root";
51  std::string textName = fileName + "_ana.txt";
52 
53  // Get options
54  std::string polarity = "";
55  bool extrap = false;
56  if (opt.Contains("rhc", TString::kIgnoreCase) and
57  !opt.Contains("fhc", TString::kIgnoreCase))
58  polarity = "rhc";
59  else if (opt.Contains("fhc", TString::kIgnoreCase) and
60  !opt.Contains("rhc", TString::kIgnoreCase))
61  polarity = "fhc";
62  else {
63  std::cout << "Flux required: fhc or rhc" << std::endl;
64  abort();
65  }
66  if (opt.Contains("extrap", TString::kIgnoreCase))
67  extrap = true;
68 
69  double kNDScale, kFDScale;
70  if (polarity == "fhc") {
71  kNDScale = kAna2018RHCNDPOT;
72  kFDScale = kAna2018RHCPOT;
73  } else if (polarity == "rhc") {
74  kNDScale = kAna2018FHCNDPOT;
75  kFDScale = kAna2018FHCPOT;
76  }
77 
78  // A map of systematic labels
79  std::map<std::string, std::string> shiftLabels;
80  shiftLabels["ULLevel"] = "Light level Up and Calibration Down";
81  shiftLabels["DLLevel"] = "Light level Down and Calibration Up";;
82  shiftLabels["Cherenkov"] = "Cherenkov variation";
83 
84  // Get the vector of samples to analyse
85  std::vector<std::string> cut_samples;
86  cut_samples.push_back("Nus18");
87 
88  // Set up maps to the decompositions/predictions (each flavor component)
89  // Nominal maps are indexed purely by sample label
90  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
91  std::map<std::string, NDPredictionSterile*> predsND_nominal;
92  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
93  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
94  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
95  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
96  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
97 
98  // Fill the maps from the data in the input file
99  TFile* rootL = new TFile(inFile.c_str(), "READ");
100  LoadMaps(rootL, &predsND_nominal, &predsND_shifted);
101  LoadMaps(rootL, &predsFD_nominal, &predsFD_shifted);
102  if (extrap)
103  LoadMaps(rootL, &predsExtrap_nominal, &predsExtrap_shifted);
104  rootL->Close();
105 
106  // Get output files
107  TFile* rootF = new TFile(saveName.c_str(), "RECREATE");
108  FILE* textF;
109  textF = fopen(textName.c_str(), "w+");
110 
111  // Shifted spectra
112  std::map<std::string, ShiftedLightLevel> shiftedSpectra;
113 
114  // Look through all samples
115  for (const auto& sample : cut_samples) {
116 
117  strings strs;
118  strs.fComponent = "NC";
119  strs.fDet = "ND";
120  strs.fPOT = Form("%f", kNDScale);
121  strs.fSample = sample;
122  strs.fSystType = "LightLevel";
123  strs.fXLabel = kNus18AxisE.GetLabels()[0];
124 
125  InitializeSystText(textF, strs);
126 
127  // Get nominal spectra
128  TH1* hNomNC_predND = GetNC(predsND_nominal[sample], kNDScale);
129  TH1* hNomBG_predND = GetBG(predsND_nominal[sample], kNDScale);
130  TH1* hNomNC_predFD = GetNC(predsFD_nominal[sample], kFDScale);
131  TH1* hNomBG_predFD = GetBG(predsFD_nominal[sample], kFDScale);
132  TH1 *hNomNC_predExtrap, *hNomBG_predExtrap;
133  if (extrap) {
134  hNomNC_predExtrap = GetNC(predsExtrap_nominal[sample], kFDScale);
135  hNomBG_predExtrap = GetBG(predsExtrap_nominal[sample], kFDScale);
136  }
137 
138  // Look through all systematic shifts
139  for (const auto& shifts : shiftLabels) {
140 
141  strs.fSystS = shifts.first;
142  strs.fSystL = shifts.second;
143 
144  strs.fDet = "ND";
145  strs.fComponent = "NC";
146  PlotSyst(hNomNC_predND,
147  predsND_shifted[sample][shifts.first][1],
148  rootF, textF, strs, kNDScale, true, false);
149  shiftedSpectra[shifts.first].hNC_ND = GetNC(predsND_shifted[sample][shifts.first][1], kNDScale);
150 
151  strs.fComponent = "BG";
152  PlotSyst(hNomBG_predND,
153  predsND_shifted[sample][shifts.first][1],
154  rootF, textF, strs, kNDScale, false, false);
155  shiftedSpectra[shifts.first].hBG_ND = GetBG(predsND_shifted[sample][shifts.first][1], kNDScale);
156 
157  strs.fDet = "FD";
158  strs.fComponent = "NC";
159  PlotSyst(hNomNC_predFD,
160  predsFD_shifted[sample][shifts.first][1],
161  rootF, textF, strs, kFDScale, true, false);
162  shiftedSpectra[shifts.first].hNC_FD = GetNC(predsFD_shifted[sample][shifts.first][1], kFDScale);
163 
164  strs.fComponent = "BG";
165  PlotSyst(hNomBG_predFD,
166  predsFD_shifted[sample][shifts.first][1],
167  rootF, textF, strs, kFDScale, false, false);
168  shiftedSpectra[shifts.first].hBG_FD = GetBG(predsFD_shifted[sample][shifts.first][1], kFDScale);
169 
170  if (extrap) {
171  strs.fDet = "EX";
172  strs.fComponent = "NC";
173  PlotSyst(hNomNC_predExtrap,
174  predsExtrap_shifted[sample][shifts.first][1],
175  rootF, textF, strs, kFDScale, true, false);
176  shiftedSpectra[shifts.first].hNC_EX = GetNC(predsExtrap_shifted[sample][shifts.first][1], kFDScale);
177 
178  strs.fComponent = "BG";
179  PlotSyst(hNomBG_predExtrap,
180  predsExtrap_shifted[sample][shifts.first][1],
181  rootF, textF, strs, kFDScale, false, false);
182  shiftedSpectra[shifts.first].hBG_EX = GetBG(predsExtrap_shifted[sample][shifts.first][1], kFDScale);
183 
184  }
185 
186  }
187 
188  // for (const auto& det : {"ND", "FD"})
189  // for (const auto& type : {"NC", "BG", "EX"})
190  // for (const auto& syst : shiftLabels) {
191 
192  strs.fDet = "ND";
193  strs.fComponent = "NC";
194  PlotSyst(hNomNC_predND,
195  shiftedSpectra["ULLevel"].hNC_ND,
196  shiftedSpectra["DLLevel"].hNC_ND,
197  shiftedSpectra["Cherenkov"].hNC_ND,
198  rootF, textF, strs);
199 
200  strs.fComponent = "BG";
201  PlotSyst(hNomBG_predND,
202  shiftedSpectra["ULLevel"].hBG_ND,
203  shiftedSpectra["DLLevel"].hBG_ND,
204  shiftedSpectra["Cherenkov"].hBG_ND,
205  rootF, textF, strs);
206 
207  strs.fDet = "FD";
208  strs.fComponent = "NC";
209  PlotSyst(hNomNC_predFD,
210  shiftedSpectra["ULLevel"].hNC_FD,
211  shiftedSpectra["DLLevel"].hNC_FD,
212  shiftedSpectra["Cherenkov"].hNC_FD,
213  rootF, textF, strs);
214 
215  strs.fComponent = "BG";
216  PlotSyst(hNomBG_predFD,
217  shiftedSpectra["ULLevel"].hBG_FD,
218  shiftedSpectra["DLLevel"].hBG_FD,
219  shiftedSpectra["Cherenkov"].hBG_FD,
220  rootF, textF, strs);
221 
222  if (extrap) {
223  strs.fDet = "EX";
224  strs.fComponent = "NC";
225  PlotSyst(hNomNC_predExtrap,
226  shiftedSpectra["ULLevel"].hNC_EX,
227  shiftedSpectra["DLLevel"].hNC_EX,
228  shiftedSpectra["Cherenkov"].hNC_EX,
229  rootF, textF, strs);
230 
231  strs.fComponent = "BG";
232  PlotSyst(hNomBG_predExtrap,
233  shiftedSpectra["ULLevel"].hBG_EX,
234  shiftedSpectra["DLLevel"].hBG_EX,
235  shiftedSpectra["Cherenkov"].hBG_EX,
236  rootF, textF, strs);
237  }
238 
239  }
240 
241  fclose(textF);
242  rootF->Close();
243 
244  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
245 
246 }
std::string fDet
Definition: PPFXHelper.h:104
fileName
Definition: plotROC.py:78
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string fComponent
Definition: PPFXHelper.h:103
void PlotSyst(TH1 *h, TH1 *hp1, TH1 *hm1, TH1 *hp2, TH1 *hm2, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:662
TH1 * GetSpectra(std::string type, std::string det)
const double kAna2018RHCPOT
Definition: Exposures.h:208
const HistAxis kNus18AxisE("Energy Deposited in Scintillator (GeV)", kNus18EnergyBinning, kNus18Energy)
Axes used in Nus18 analysis by nus group.
Definition: HistAxes.h:16
TH1 * GetNC(IDecomp *specs, double POT)
Definition: PPFXHelper.h:211
fclose(fg1)
ifstream inFile
Definition: AnaPlotMaker.h:34
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
std::string fXLabel
Definition: PPFXHelper.h:110
void Nus18SystsLightLevelAna(std::string inFile, TString opt)
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
const double kAna2018FHCPOT
Definition: Exposures.h:207
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
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