Nus18SystsBeamTranspAna.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // Beam Transport 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["BeamHornCurrent"] = "2kA";
57  shiftLabels["BeamHornCurrent"] = "Horn Current";
58  shiftLabels["BeamSpotSize"] = "Spot Size";
59  shiftLabels["BeamposX"] = "Beam Position X";
60  shiftLabels["BeamposY"] = "Beam Position Y";
61  shiftLabels["BeamH1PosX"] = "Horn 1 X Position";
62  shiftLabels["BeamH1PosY"] = "Horn 1 Y Position";
63  shiftLabels["BeamH2PosX"] = "Horn 2 X Position";
64  shiftLabels["BeamH2PosY"] = "Horn 2 Y Position";
65  shiftLabels["TargetPosZ"] = "Target Z position";
66  shiftLabels["BeamExp"] = "Magnet Field in Decay Pipe";
67  shiftLabels["HornWater"] = "New Horn Geometry and 1mm water";
68  shiftLabels["totErr"] = "Combined Beam Transport Systematics";
69 
70  // Get the vector of samples to analyse
71  std::vector<std::string> cut_samples;
72  cut_samples.push_back("Nus18");
73 
74  // Set up maps to the decompositions/predictions (each flavor component)
75  // Nominal maps are indexed purely by sample label
76  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
77  std::map<std::string, NDPredictionSterile*> predsND_nominal;
78  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
79  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
80  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
81  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
82  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
83 
84  // Fill the maps from the data in the input file
85  TFile* rootL = new TFile(inFile.c_str(), "READ");
86  LoadMaps(rootL, &predsND_nominal, &predsND_shifted);
87  LoadMaps(rootL, &predsFD_nominal, &predsFD_shifted);
88  if (extrap)
89  LoadMaps(rootL, &predsExtrap_nominal, &predsExtrap_shifted);
90  rootL->Close();
91 
92  // Get output files
93  TFile* rootF = new TFile(saveName.c_str(), "RECREATE");
94  FILE* textF;
95  textF = fopen(textName.c_str(), "w+");
96 
97  // Look through all samples
98  for (const auto& sample : cut_samples) {
99 
100  strings strs;
101  strs.fComponent = "NC";
102  strs.fDet = "ND";
103  strs.fPOT = Form("%f", kNDScale);
104  strs.fSample = sample;
105  strs.fSystType = "Beam";
106  strs.fXLabel = kNus18AxisE.GetLabels()[0];
107 
108  InitializeSystText(textF, strs);
109 
110  // Get nominal spectra
111  TH1* hNomNC_predND = GetNC(predsND_nominal[sample], kNDScale);
112  TH1* hNomBG_predND = GetBG(predsND_nominal[sample], kNDScale);
113  TH1* hNomNC_predFD = GetNC(predsFD_nominal[sample], kFDScale);
114  TH1* hNomBG_predFD = GetBG(predsFD_nominal[sample], kFDScale);
115  TH1 *hNomNC_predExtrap, *hNomBG_predExtrap;
116  if (extrap) {
117  hNomNC_predExtrap = GetNC(predsExtrap_nominal[sample], kFDScale);
118  hNomBG_predExtrap = GetBG(predsExtrap_nominal[sample], kFDScale);
119  }
120 
121  // Look through all systematic shifts
122  for (const auto& shifts : shiftLabels) {
123 
124  strs.fSystS = shifts.first;
125  strs.fSystL = shifts.second;
126 
127  strs.fDet = "ND";
128  strs.fComponent = "NC";
129  PlotSyst(hNomNC_predND,
130  predsND_shifted[sample][shifts.first],
131  rootF, textF, strs, kNDScale, true, false);
132 
133  strs.fComponent = "BG";
134  PlotSyst(hNomBG_predND,
135  predsND_shifted[sample][shifts.first],
136  rootF, textF, strs, kNDScale, false, false);
137 
138  strs.fDet = "FD";
139  strs.fComponent = "NC";
140  PlotSyst(hNomNC_predFD,
141  predsFD_shifted[sample][shifts.first],
142  rootF, textF, strs, kFDScale, true, false);
143 
144  strs.fComponent = "BG";
145  PlotSyst(hNomBG_predFD,
146  predsFD_shifted[sample][shifts.first],
147  rootF, textF, strs, kFDScale, false, false);
148 
149  if (extrap) {
150  strs.fDet = "EX";
151  strs.fComponent = "NC";
152  PlotSyst(hNomNC_predExtrap,
153  predsExtrap_shifted[sample][shifts.first],
154  rootF, textF, strs, kFDScale, true, false);
155 
156  strs.fComponent = "BG";
157  PlotSyst(hNomBG_predExtrap,
158  predsExtrap_shifted[sample][shifts.first],
159  rootF, textF, strs, kFDScale, false, false);
160  }
161 
162  }
163 
164  strs.fDet = "ND";
165  strs.fComponent = "NC";
166  PlotSystBand(hNomNC_predND, rootF, textF, strs);
167 
168  strs.fComponent = "BG";
169  PlotSystBand(hNomBG_predND, rootF, textF, strs);
170 
171  strs.fDet = "FD";
172  strs.fComponent = "NC";
173  PlotSystBand(hNomNC_predFD, rootF, textF, strs);
174 
175  strs.fComponent = "BG";
176  PlotSystBand(hNomBG_predFD, rootF, textF, strs);
177 
178  if (extrap) {
179  strs.fDet = "EX";
180  strs.fComponent = "NC";
181  PlotSystBand(hNomNC_predExtrap, rootF, textF, strs);
182 
183  strs.fComponent = "BG";
184  PlotSystBand(hNomBG_predExtrap, rootF, textF, strs);
185  }
186 
187  }
188 
189  fclose(textF);
190  rootF->Close();
191 
192  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
193 
194 }
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 Nus18SystsBeamTranspAna(std::string inFile, TString opt)
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
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