Nus18SystsCalibAbsAna.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // Absolute Calibration 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["CalFlatDown"] = "Flat Miscalibration Down";
81  shiftLabels["CalFlatUp"] = "Flat Miscalibration Up";
82  shiftLabels["CalShape"] = "Calibration shape";
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, ShiftedCalibAbs> 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 = "CalibAbs";
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,
149  shifts.first.find("Flat") != std::string::npos); // use max for flat systs
150  shiftedSpectra[shifts.first].hNC_ND = GetNC(predsND_shifted[sample][shifts.first][1], kNDScale);
151 
152  strs.fComponent = "BG";
153  PlotSyst(hNomBG_predND,
154  predsND_shifted[sample][shifts.first][1],
155  rootF, textF, strs, kNDScale, false,
156  shifts.first.find("Flat") != std::string::npos);
157  shiftedSpectra[shifts.first].hBG_ND = GetBG(predsND_shifted[sample][shifts.first][1], kNDScale);
158 
159  strs.fDet = "FD";
160  strs.fComponent = "NC";
161  PlotSyst(hNomNC_predFD,
162  predsFD_shifted[sample][shifts.first][1],
163  rootF, textF, strs, kFDScale, true,
164  shifts.first.find("Flat") != std::string::npos);
165  shiftedSpectra[shifts.first].hNC_FD = GetNC(predsFD_shifted[sample][shifts.first][1], kFDScale);
166 
167  strs.fComponent = "BG";
168  PlotSyst(hNomBG_predFD,
169  predsFD_shifted[sample][shifts.first][1],
170  rootF, textF, strs, kFDScale, false,
171  shifts.first.find("Flat") != std::string::npos);
172  shiftedSpectra[shifts.first].hBG_FD = GetBG(predsFD_shifted[sample][shifts.first][1], kFDScale);
173 
174  if (extrap) {
175  strs.fDet = "EX";
176  strs.fComponent = "NC";
177  PlotSyst(hNomNC_predExtrap,
178  predsExtrap_shifted[sample][shifts.first][1],
179  rootF, textF, strs, kFDScale, true,
180  shifts.first.find("Flat") != std::string::npos);
181  shiftedSpectra[shifts.first].hNC_EX = GetNC(predsExtrap_shifted[sample][shifts.first][1], kFDScale);
182 
183  strs.fComponent = "BG";
184  PlotSyst(hNomBG_predExtrap,
185  predsExtrap_shifted[sample][shifts.first][1],
186  rootF, textF, strs, kFDScale, false,
187  shifts.first.find("Flat") != std::string::npos);
188  shiftedSpectra[shifts.first].hBG_EX = GetBG(predsExtrap_shifted[sample][shifts.first][1], kFDScale);
189 
190  }
191 
192  }
193 
194  // for (const auto& det : {"ND", "FD"})
195  // for (const auto& type : {"NC", "BG", "EX"})
196  // for (const auto& syst : shiftLabels) {
197 
198  strs.fDet = "ND";
199  strs.fComponent = "NC";
200  PlotSyst(hNomNC_predND,
201  shiftedSpectra["CalFlatUp"].hNC_ND,
202  shiftedSpectra["CalFlatDown"].hNC_ND,
203  shiftedSpectra["CalShape"].hNC_ND,
204  rootF, textF, strs);
205 
206  strs.fComponent = "BG";
207  PlotSyst(hNomBG_predND,
208  shiftedSpectra["CalFlatUp"].hBG_ND,
209  shiftedSpectra["CalFlatDown"].hBG_ND,
210  shiftedSpectra["CalShape"].hBG_ND,
211  rootF, textF, strs);
212 
213  strs.fDet = "FD";
214  strs.fComponent = "NC";
215  PlotSyst(hNomNC_predFD,
216  shiftedSpectra["CalFlatUp"].hNC_FD,
217  shiftedSpectra["CalFlatDown"].hNC_FD,
218  shiftedSpectra["CalShape"].hNC_FD,
219  rootF, textF, strs);
220 
221  strs.fComponent = "BG";
222  PlotSyst(hNomBG_predFD,
223  shiftedSpectra["CalFlatUp"].hBG_FD,
224  shiftedSpectra["CalFlatDown"].hBG_FD,
225  shiftedSpectra["CalShape"].hBG_FD,
226  rootF, textF, strs);
227 
228  if (extrap) {
229  strs.fDet = "EX";
230  strs.fComponent = "NC";
231  PlotSyst(hNomNC_predExtrap,
232  shiftedSpectra["CalFlatUp"].hNC_EX,
233  shiftedSpectra["CalFlatDown"].hNC_EX,
234  shiftedSpectra["CalShape"].hNC_EX,
235  rootF, textF, strs);
236 
237  strs.fComponent = "BG";
238  PlotSyst(hNomBG_predExtrap,
239  shiftedSpectra["CalFlatUp"].hBG_EX,
240  shiftedSpectra["CalFlatDown"].hBG_EX,
241  shiftedSpectra["CalShape"].hBG_EX,
242  rootF, textF, strs);
243  }
244 
245  }
246 
247  fclose(textF);
248  rootF->Close();
249 
250  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
251 
252 }
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
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 * GetSpectra(std::string type, std::string det)
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
void Nus18SystsCalibAbsAna(std::string inFile, TString opt)
std::string fPOT
Definition: PPFXHelper.h:105
OStream cout
Definition: OStream.cxx:6
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
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