25 const std::vector<std::string>
vars3D{
"EAvail"};
26 const std::vector<std::string>
vars1D 39 {
"angleup",
"Angle Up"},
40 {
"angledw",
"Angle Down"},
42 {
"neutron_Up",
"Neutron Up"},
43 {
"neutron_Dw",
"Neutron Down"}
47 {
"neutron",
"Neutron"}
64 TAxis*
xaxis =
new TAxis(1, 0.5, 1.0);
65 TAxis*
yaxis =
new TAxis(1, 0.5, 2.5);
66 TAxis*
zaxis =
new TAxis(1, -0.25, 0.25);
67 TAxis*
eaxis =
new TAxis(1, 0.0, 2.5);
72 TH2F * rat =
new TH2F(*syst - *nominal);
78 TH1F * rat =
new TH1F(*syst - *nominal);
85 std::string outputFile = outputDir +
"/UnfoldingResults.root";
86 std::string outputPDF = outputDir +
"/all_uncertainties.pdf";
88 TFile *
fin =
new TFile(inputFile.c_str());
89 TFile *
fout =
new TFile(outputFile.c_str(),
"RECREATE");
90 TDirectory * systDir = fout->mkdir(
"systematics");
93 TCanvas *
c1 =
new TCanvas(
"c1",
"c1", 1600, 1200);
94 c1->SaveAs((outputPDF +
"[").c_str());
97 gStyle->SetPaintTextFormat(
".3f");
102 TH3F * nom3D = (TH3F*) fin->Get((
xs_mode +
"/hXS3Dsec_neutron_nom_" +
var).c_str());
103 TH2F * nom2D = (TH2F*) fin->Get((
xs_mode +
"/hXS2Dsec_neutron_nom_" +
var).c_str());
104 TH1F * nom1D = (TH1F*) fin->Get((
xs_mode +
"/hXS1Dsec_neutron_nom_" +
var).c_str());
105 TH2F * nomEff = (TH2F*) fin->Get((
"h2DEff_nom_" +
var).c_str());
107 for (std::pair<std::string, std::string> syst :
systNames){
110 TH2F * systHist2D = (TH2F*) fin->Get((
xs_mode +
"/hXS2Dsec_" + systName +
"_" +
var).c_str());
112 systUncert2D->SetName((systName +
"_fracuncert").c_str());
113 systUncert2D->SetTitle((
"(" + systTitle +
" - nominal) / nominal" +
yxaxistitle).c_str());
114 systUncert2D->GetXaxis()->SetRangeUser(
xaxis->GetXmin(),
xaxis->GetXmax());
115 systUncert2D->GetYaxis()->SetRangeUser(
yaxis->GetXmin(),
yaxis->GetXmax());
116 systUncert2D->GetZaxis()->SetRangeUser(
zaxis->GetXmin(),
zaxis->GetXmax());
117 systUncert2D->Draw(
"COLZ");
118 systUncert2D->Write();
119 c1->SaveAs((outputDir +
"/" + systName +
"_fracuncert" + outputFormat).c_str());
120 c1->SaveAs(outputPDF.c_str());
122 TH1F * systHist1D = (TH1F*) fin->Get((
xs_mode +
"/hXS1Dsec_" + systName +
"_" +
var).c_str());
124 systUncert1D->SetName((systName +
"_fracuncert_1d").c_str());
125 systUncert1D->SetTitle((
"(" + systTitle +
" - nominal) / nominal" +
zaxistitle).c_str());
126 systUncert1D->GetXaxis()->SetRangeUser(
eaxis->GetXmin(),
eaxis->GetXmax());
127 systUncert1D->Draw(
"hist");
128 systUncert1D->Write();
129 c1->SaveAs((outputDir +
"/" + systName +
"_fracuncert_1d" + outputFormat).c_str());
130 c1->SaveAs(outputPDF.c_str());
134 for (
unsigned int iupdown = 0; iupdown <
updownlabel.size(); iupdown++){
136 std::string systTitle = syst + updownlabel[iupdown];
137 TH2F * systHist2D = (TH2F*) fin->Get((
xs_mode +
"/hXS2Dsec_" + systName +
"_" +
var).c_str());
139 systUncert2D->SetName((systName +
"_fracuncert").c_str());
140 systUncert2D->SetTitle((
"(" + systTitle +
" - nominal) / nominal" +
yxaxistitle).c_str());
141 systUncert2D->GetXaxis()->SetRangeUser(
xaxis->GetXmin(),
xaxis->GetXmax());
142 systUncert2D->GetYaxis()->SetRangeUser(
yaxis->GetXmin(),
yaxis->GetXmax());
143 systUncert2D->GetZaxis()->SetRangeUser(
zaxis->GetXmin(),
zaxis->GetXmax());
144 systUncert2D->Draw(
"COLZ");
145 systUncert2D->Write();
146 c1->SaveAs((outputDir +
"/" + systName +
"_fracuncert" + outputFormat).c_str());
147 c1->SaveAs(outputPDF.c_str());
150 nom2D->SetName(
"nominal_2Dxsec");
151 nom2D->SetTitle(
"Nominal Cross-section;Cos(#theta_{#mu});Muon Kinetic Energy (GeV)");
152 nom2D->GetXaxis()->SetRangeUser(
xaxis->GetXmin(),
xaxis->GetXmax());
153 nom2D->GetYaxis()->SetRangeUser(
yaxis->GetXmin(),
yaxis->GetXmax());
154 gPad->SetRightMargin(0.15);
157 c1->SaveAs((outputDir +
"/nominal_xsec" + outputFormat).c_str());
159 nom1D->SetName(
"nominal_1Dxsec");
160 nom1D->SetTitle(
"Nominal Cross-section;Available Energy (GeV)");
161 nom1D->GetXaxis()->SetRangeUser(
eaxis->GetXmin(),
eaxis->GetXmax());
164 c1->SaveAs((outputDir +
"/nominal_xsec_1d" + outputFormat).c_str());
166 nomEff->SetName(
"nominal_eff");
167 nomEff->SetTitle((
"Nominal Efficiency" +
yxaxistitle).c_str());
168 nomEff->GetXaxis()->SetRangeUser(
xaxis->GetXmin(),
xaxis->GetXmax());
169 nomEff->GetYaxis()->SetRangeUser(
yaxis->GetXmin(),
yaxis->GetXmax());
170 nomEff->GetZaxis()->SetRangeUser(0.0, 1.0);
171 nomEff->Draw(
"COLZ");
173 c1->SaveAs((outputDir +
"/nominal_eff" + outputFormat).c_str());
175 for (std::pair<std::string, std::string> syst : systNames){
178 TH2F * systEff = (TH2F*) fin->Get((
"h2DEff_" + systName +
"_" +
var).c_str());
179 systEff->SetName((systName +
"_eff").c_str());
180 systEff->SetTitle((systTitle +
" Efficiency" +
yxaxistitle).c_str());
181 systEff->GetXaxis()->SetRangeUser(
xaxis->GetXmin(),
xaxis->GetXmax());
182 systEff->GetYaxis()->SetRangeUser(
yaxis->GetXmin(),
yaxis->GetXmax());
183 systEff->GetZaxis()->SetRangeUser(0.0, 1.0);
184 systEff->Draw(
"COLZ");
186 c1->SaveAs((outputDir +
"/" + systName +
"_eff" + outputFormat).c_str());
190 for (
unsigned int iupdown = 0; iupdown <
updownlabel.size(); iupdown++){
192 std::string systTitle = syst + updownlabel[iupdown];
193 TH2F * systEff = (TH2F*) fin->Get((
"h2DEff_" + systName +
"_" +
var).c_str());
194 systEff->SetName((systName +
"_eff").c_str());
195 systEff->SetTitle((systTitle +
" Efficiency" +
yxaxistitle).c_str());
196 systEff->GetXaxis()->SetRangeUser(
xaxis->GetXmin(),
xaxis->GetXmax());
197 systEff->GetYaxis()->SetRangeUser(
yaxis->GetXmin(),
yaxis->GetXmax());
198 systEff->GetZaxis()->SetRangeUser(0.0, 1.0);
199 systEff->Draw(
"COLZ");
201 c1->SaveAs((outputDir +
"/" + systName +
"_eff" + outputFormat).c_str());
205 c1->SaveAs((outputPDF +
"]").c_str());
const std::vector< std::string > vars1D
const std::string xs_mode0
std::vector< std::string > beamSysts
TH1F * makeSystematic1D(TH1F *syst, TH1F *nominal)
std::map< std::string, std::string > systUpDownPairs
const std::string xs_mode1
const std::string zaxistitle
const std::string xs_mode
const std::vector< std::string > vars3D
void PlotUnfolding(std::string inputFile, std::string outputDir)
std::vector< std::string > updownlabel
std::map< std::string, std::string > systNames
const std::string yxaxistitle
TH2F * makeSystematic2D(TH2F *syst, TH2F *nominal)