Nus18SystsXSecOnOffAna.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // Cross-section on/off 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 
22  TH1::AddDirectory(0);
23 
24  // Set up path to file with filled CAFAna objects and for output of analysis
25  std::string fileName = inFile.substr(0, inFile.find(".root"));
26  std::string saveName = fileName + "_ana.root";
27  std::string textName = fileName + "_ana.txt";
28 
29  // Get options
30  std::string polarity = "";
31  bool extrap = false;
32  if (opt.Contains("rhc", TString::kIgnoreCase) and
33  !opt.Contains("fhc", TString::kIgnoreCase))
34  polarity = "rhc";
35  else if (opt.Contains("fhc", TString::kIgnoreCase) and
36  !opt.Contains("rhc", TString::kIgnoreCase))
37  polarity = "fhc";
38  else {
39  std::cout << "Flux required: fhc or rhc" << std::endl;
40  abort();
41  }
42  if (opt.Contains("extrap", TString::kIgnoreCase))
43  extrap = true;
44 
45  double kNDScale, kFDScale;
46  if (polarity == "fhc") {
47  kNDScale = kAna2018RHCNDPOT;
48  kFDScale = kAna2018RHCPOT;
49  } else if (polarity == "rhc") {
50  kNDScale = kAna2018FHCNDPOT;
51  kFDScale = kAna2018FHCPOT;
52  }
53 
54  // A map of systematic labels
55  std::map<std::string, std::string> shiftLabels;
56  shiftLabels["NoXSec"] = "Cross-section off";
57 
58  // Get the vector of samples to analyse
59  std::vector<std::string> cut_samples;
60  cut_samples.push_back("Nus18");
61 
62  // Set up maps to the decompositions/predictions (each flavor component)
63  // Nominal maps are indexed purely by sample label
64  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
65  std::map<std::string, NDPredictionSterile*> predsND_nominal;
66  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
67  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
68  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
69  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
70  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
71 
72  // Fill the maps from the data in the input file
73  TFile* rootL = new TFile(inFile.c_str(), "READ");
74  LoadMaps(rootL, &predsND_nominal, &predsND_shifted);
75  LoadMaps(rootL, &predsFD_nominal, &predsFD_shifted);
76  if (extrap)
77  LoadMaps(rootL, &predsExtrap_nominal, &predsExtrap_shifted);
78  rootL->Close();
79 
80  // Get output files
81  TFile* rootF = new TFile(saveName.c_str(), "RECREATE");
82  FILE* textF;
83  textF = fopen(textName.c_str(), "w+");
84 
85  // Look through all samples
86  for (const auto& sample : cut_samples) {
87 
88  strings strs;
89  strs.fComponent = "NC";
90  strs.fDet = "ND";
91  strs.fPOT = Form("%f", kNDScale);
92  strs.fSample = sample;
93  strs.fSystType = "XSecOnOff";
94  strs.fXLabel = kNus18AxisE.GetLabels()[0];
95 
96  InitializeSystText(textF, strs);
97 
98  // Get nominal spectra
99  TH1* hNomNC_predND = GetNC(predsND_nominal[sample], kNDScale);
100  TH1* hNomBG_predND = GetBG(predsND_nominal[sample], kNDScale);
101  TH1* hNomNC_predFD = GetNC(predsFD_nominal[sample], kFDScale);
102  TH1* hNomBG_predFD = GetBG(predsFD_nominal[sample], kFDScale);
103  TH1 *hNomNC_predExtrap, *hNomBG_predExtrap;
104  if (extrap) {
105  hNomNC_predExtrap = GetNC(predsExtrap_nominal[sample], kFDScale);
106  hNomBG_predExtrap = GetBG(predsExtrap_nominal[sample], kFDScale);
107  }
108 
109  // Look through all systematic shifts
110  for (const auto& shifts : shiftLabels) {
111 
112  strs.fSystS = shifts.first;
113  strs.fSystL = shifts.second;
114 
115  strs.fDet = "ND";
116  strs.fComponent = "NC";
117  PlotSyst(hNomNC_predND,
118  predsND_shifted[sample][shifts.first][1],
119  rootF, textF, strs, kNDScale, true, false);
120 
121  strs.fComponent = "BG";
122  PlotSyst(hNomBG_predND,
123  predsND_shifted[sample][shifts.first][1],
124  rootF, textF, strs, kNDScale, false, false);
125 
126  strs.fDet = "FD";
127  strs.fComponent = "NC";
128  PlotSyst(hNomNC_predFD,
129  predsFD_shifted[sample][shifts.first][1],
130  rootF, textF, strs, kFDScale, true, false);
131 
132  strs.fComponent = "BG";
133  PlotSyst(hNomBG_predFD,
134  predsFD_shifted[sample][shifts.first][1],
135  rootF, textF, strs, kFDScale, false, false);
136 
137  if (extrap) {
138  strs.fDet = "EX";
139  strs.fComponent = "NC";
140  PlotSyst(hNomNC_predExtrap,
141  predsExtrap_shifted[sample][shifts.first][1],
142  rootF, textF, strs, kFDScale, true, false);
143 
144  strs.fComponent = "BG";
145  PlotSyst(hNomBG_predExtrap,
146  predsExtrap_shifted[sample][shifts.first][1],
147  rootF, textF, strs, kFDScale, false, false);
148  }
149 
150  }
151 
152  strs.fDet = "ND";
153  strs.fComponent = "NC";
154  PlotSystBand(hNomNC_predND, rootF, textF, strs);
155 
156  strs.fComponent = "BG";
157  PlotSystBand(hNomBG_predND, rootF, textF, strs);
158 
159  strs.fDet = "FD";
160  strs.fComponent = "NC";
161  PlotSystBand(hNomNC_predFD, rootF, textF, strs);
162 
163  strs.fComponent = "BG";
164  PlotSystBand(hNomBG_predFD, rootF, textF, strs);
165 
166  if (extrap) {
167  strs.fDet = "EX";
168  strs.fComponent = "NC";
169  PlotSystBand(hNomNC_predExtrap, rootF, textF, strs);
170 
171  strs.fComponent = "BG";
172  PlotSystBand(hNomBG_predExtrap, rootF, textF, strs);
173  }
174 
175  }
176 
177  fclose(textF);
178  rootF->Close();
179 
180  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
181 
182 }
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 * 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
void PlotSystBand(TH1 *h, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:968
std::string fXLabel
Definition: PPFXHelper.h:110
void Nus18SystsXSecOnOffAna(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