Nus18SystsPPFXAna.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // PPFX 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 = kAna2018FHCNDPOT;
48  kFDScale = kAna2018FHCPOT;
49  } else if (polarity == "rhc") {
50  kNDScale = kAna2018RHCNDPOT;
51  kFDScale = kAna2018RHCPOT;
52  }
53 
54  // Get the vector of samples to analyse
55  std::vector<std::string> cut_samples;
56  cut_samples.push_back("Nus18");
57 
58  // Set up maps to the decompositions/predictions (each flavor component)
59  // Nominal maps are indexed purely by sample label
60  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
61  std::map<std::string, NDPredictionSterile*> predsND_nominal;
62  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
63  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
64  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
65  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
66  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
67 
68  // Containers for each of the random PPFX universes
69  std::vector<NDPredictionSterile*> predsND_shifted_Idx;
70  std::vector<FDPredictionSterile*> predsFD_shifted_Idx;
71  std::vector<PredictionSterile*> predsExtrap_shifted_Idx;
72 
73  // Fill the maps from the data in the input file
74  TFile* rootL = new TFile(inFile.c_str(), "READ");
75  LoadMaps(rootL, &predsND_nominal, &predsND_shifted);
76  LoadMaps(rootL, &predsFD_nominal, &predsFD_shifted);
77  if (extrap)
78  LoadMaps(rootL, &predsExtrap_nominal, &predsExtrap_shifted);
79  rootL->Close();
80 
81  // Get output files
82  TFile* rootF = new TFile(saveName.c_str(), "RECREATE");
83  FILE* textF;
84  textF = fopen(textName.c_str(), "w+");
85 
86  // Look through all samples
87  for (const auto& sample : cut_samples) {
88 
89  strings strs;
90  strs.fComponent = "NC";
91  strs.fDet = "ND";
92  strs.fPOT = Form("%f", kNDScale);
93  strs.fSample = sample;
94  strs.fSystType = "PPFX";
95  strs.fXLabel = kNus18AxisE.GetLabels()[0];
96 
97  InitializeSystText(textF, strs);
98 
99  // Get nominal spectra
100  TH1* hNomNC_predND = GetNC(predsND_nominal[sample], kNDScale);
101  TH1* hNomBG_predND = GetBG(predsND_nominal[sample], kNDScale);
102  TH1* hNomNC_predFD = GetNC(predsFD_nominal[sample], kFDScale);
103  TH1* hNomBG_predFD = GetBG(predsFD_nominal[sample], kFDScale);
104  TH1 *hNomNC_predExtrap, *hNomBG_predExtrap;
105  if (extrap) {
106  hNomNC_predExtrap = GetNC(predsExtrap_nominal[sample], kFDScale);
107  hNomBG_predExtrap = GetBG(predsExtrap_nominal[sample], kFDScale);
108  }
109 
110  // Get the shifted decompositions
111  // One is needed for every sample and PPFX universe
112  for (unsigned int Idx = 0; Idx < 100; ++Idx) {
113  predsND_shifted_Idx.push_back(predsND_shifted[sample][Form("%d",Idx)][1]);
114  predsFD_shifted_Idx.push_back(predsFD_shifted[sample][Form("%d",Idx)][1]);
115  if (extrap)
116  predsExtrap_shifted_Idx.push_back(predsExtrap_shifted[sample][Form("%d",Idx)][1]);
117  }
118 
119  strs.fSystS = "PPFXWGT";
120  strs.fSystL = "PPFX";
121 
122  strs.fDet = "ND";
123  strs.fComponent = "NC";
124  PlotMultiSyst(hNomNC_predND,
125  predsND_shifted_Idx,
126  rootF, textF, strs, kNDScale, true, true);
127 
128  strs.fComponent = "BG";
129  PlotMultiSyst(hNomBG_predND,
130  predsND_shifted_Idx,
131  rootF, textF, strs, kNDScale, false, true);
132 
133  strs.fDet = "FD";
134  strs.fComponent = "NC";
135  PlotMultiSyst(hNomNC_predFD,
136  predsFD_shifted_Idx,
137  rootF, textF, strs, kFDScale, true, true);
138 
139  strs.fComponent = "BG";
140  PlotMultiSyst(hNomBG_predFD,
141  predsFD_shifted_Idx,
142  rootF, textF, strs, kFDScale, false, true);
143 
144  if (extrap) {
145  strs.fDet = "EX";
146  strs.fComponent = "NC";
147  PlotMultiSyst(hNomNC_predExtrap,
148  predsExtrap_shifted_Idx,
149  rootF, textF, strs, kFDScale, true, true);
150 
151  strs.fComponent = "BG";
152  PlotMultiSyst(hNomBG_predExtrap,
153  predsExtrap_shifted_Idx,
154  rootF, textF, strs, kFDScale, false, true);
155  }
156 
157  }
158 
159  fclose(textF);
160  rootF->Close();
161 
162  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
163 
164 }
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
void PlotMultiSyst(TH1 *h, std::vector< TH1 * > hp, TDirectory *out, FILE *text, strings strs)
Plot multiple systematics.
Definition: PPFXHelper.h:464
std::string fComponent
Definition: PPFXHelper.h:103
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
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
void Nus18SystsPPFXAna(std::string inFile, TString opt)
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