PlotUnfolding.C
Go to the documentation of this file.
1 // PlotUnfolding.C
2 
3 #include "TH1.h"
4 #include "TH2.h"
5 #include "TH3.h"
6 #include "TLegend.h"
7 #include "TCanvas.h"
8 #include "TFile.h"
9 #include "TStyle.h"
10 
11 #include <string>
12 #include <vector>
13 
14 unsigned int nuniv = 0;
15 unsigned int nppfx = 0;
16 
17 const std::string xs_mode0 = "xs_nominal";
18 const std::string xs_mode1 = "xs_bkgcons";
19 
21 
22 const std::string yxaxistitle = ";Cos(#theta_{#mu});Muon Kinetic Energy (GeV)";
23 const std::string zaxistitle = "Available Energy";
24 
25 const std::vector<std::string> vars3D{"EAvail"};
26 const std::vector<std::string> vars1D
27 {
28  "ENu",
29  "Q2"
30 };
31 
32 std::map<std::string, std::string> systNames = {
33  // {"lightup", "syst_specs_lightup.root"},
34  // {"lightdown", "syst_specs_lightdw.root"},
35  // {"calibpos", "syst_specs_calibpos.root"},
36  // {"calibneg", "syst_specs_calibneg.root"},
37  // {"MuESDw", "Muon Energy Down"},
38  // {"MuESUp", "Muon Energy Up"},
39  {"angleup", "Angle Up"},
40  {"angledw", "Angle Down"},
41  {"ckv", "Cherenkov"},
42  {"neutron_Up", "Neutron Up"},
43  {"neutron_Dw", "Neutron Down"}
44 };
45 
46 std::map<std::string, std::string> systUpDownPairs = {
47  {"neutron", "Neutron"}
48 };
49 
50 std::vector<std::string> beamSysts = {
51  // "2kA",
52  // "02mmBeamSpotSize",
53  // "1mmBeamShiftX",
54  // "1mmBeamShiftY",
55  // "3mmHorn1X",
56  // "3mmHorn1Y",
57  // "3mmHorn2X",
58  // "3mmHorn2Y",
59  // "7mmTargetZ",
60  // "MagneticFieldinDecayPipe",
61  // "1mmHornWater"
62 };
63 
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);
68 
69 std::vector<std::string> updownlabel = {"_Up", "_Dw"};
70 
71 TH2F * makeSystematic2D(TH2F* syst, TH2F* nominal){
72  TH2F * rat = new TH2F(*syst - *nominal);
73  rat->Divide(nominal);
74  return rat;
75 }
76 
77 TH1F * makeSystematic1D(TH1F* syst, TH1F* nominal){
78  TH1F * rat = new TH1F(*syst - *nominal);
79  rat->Divide(nominal);
80  return rat;
81 }
82 
84  // Load the inputs
85  std::string outputFile = outputDir + "/UnfoldingResults.root";
86  std::string outputPDF = outputDir + "/all_uncertainties.pdf";
87  std::string outputFormat = ".png";
88  TFile * fin = new TFile(inputFile.c_str());
89  TFile * fout = new TFile(outputFile.c_str(), "RECREATE");
90  TDirectory * systDir = fout->mkdir("systematics");
91  fin->cd();
92 
93  TCanvas * c1 = new TCanvas("c1", "c1", 1600, 1200);
94  c1->SaveAs((outputPDF + "[").c_str());
95  fout->cd();
96 
97  gStyle->SetPaintTextFormat(".3f");
98 
99  for (std::string var : vars3D){
100 
101  // Load nominal
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());
106 
107  for (std::pair<std::string, std::string> syst : systNames){
108  std::string systName = syst.first;
109  std::string systTitle = syst.second;
110  TH2F * systHist2D = (TH2F*) fin->Get((xs_mode + "/hXS2Dsec_" + systName + "_" + var).c_str());
111  TH2F * systUncert2D = makeSystematic2D(systHist2D, nom2D);
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());
121 
122  TH1F * systHist1D = (TH1F*) fin->Get((xs_mode + "/hXS1Dsec_" + systName + "_" + var).c_str());
123  TH1F * systUncert1D = makeSystematic1D(systHist1D, nom1D);
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());
131  }
132 
133  for (std::string syst : beamSysts){
134  for (unsigned int iupdown = 0; iupdown < updownlabel.size(); iupdown++){
135  std::string systName = syst + updownlabel[iupdown];
136  std::string systTitle = syst + updownlabel[iupdown];
137  TH2F * systHist2D = (TH2F*) fin->Get((xs_mode + "/hXS2Dsec_" + systName + "_" + var).c_str());
138  TH2F * systUncert2D = makeSystematic2D(systHist2D, nom2D);
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());
148  }
149  }
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);
155  nom2D->Draw("COLZ");
156  nom2D->Write();
157  c1->SaveAs((outputDir + "/nominal_xsec" + outputFormat).c_str());
158 
159  nom1D->SetName("nominal_1Dxsec");
160  nom1D->SetTitle("Nominal Cross-section;Available Energy (GeV)");
161  nom1D->GetXaxis()->SetRangeUser(eaxis->GetXmin(), eaxis->GetXmax());
162  nom1D->Draw("hist");
163  nom1D->Write();
164  c1->SaveAs((outputDir + "/nominal_xsec_1d" + outputFormat).c_str());
165 
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");
172  nomEff->Write();
173  c1->SaveAs((outputDir + "/nominal_eff" + outputFormat).c_str());
174 
175  for (std::pair<std::string, std::string> syst : systNames){
176  std::string systName = syst.first;
177  std::string systTitle = syst.second;
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");
185  systEff->Write();
186  c1->SaveAs((outputDir + "/" + systName + "_eff" + outputFormat).c_str());
187  }
188 
189  for (std::string syst : beamSysts)
190  for (unsigned int iupdown = 0; iupdown < updownlabel.size(); iupdown++){
191  std::string systName = syst + updownlabel[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");
200  systEff->Write();
201  c1->SaveAs((outputDir + "/" + systName + "_eff" + outputFormat).c_str());
202  }
203  }
204 
205  c1->SaveAs((outputPDF + "]").c_str());
206  fout->Close();
207 }
TAxis * xaxis
Definition: PlotUnfolding.C:64
TString fin
Definition: Style.C:24
TFile * inputFile
Definition: PlotXSec.C:134
const std::vector< std::string > vars1D
Definition: PlotUnfolding.C:27
unsigned int nppfx
Definition: PlotUnfolding.C:15
const std::string xs_mode0
Definition: PlotUnfolding.C:17
std::vector< std::string > beamSysts
Definition: PlotUnfolding.C:50
TH1F * makeSystematic1D(TH1F *syst, TH1F *nominal)
Definition: PlotUnfolding.C:77
std::map< std::string, std::string > systUpDownPairs
Definition: PlotUnfolding.C:46
const std::string xs_mode1
Definition: PlotUnfolding.C:18
TAxis * yaxis
Definition: PlotUnfolding.C:65
const std::string zaxistitle
Definition: PlotUnfolding.C:23
TAxis * zaxis
Definition: PlotUnfolding.C:66
TAxis * eaxis
Definition: PlotUnfolding.C:67
const std::string xs_mode
Definition: PlotUnfolding.C:20
const std::vector< std::string > vars3D
Definition: PlotUnfolding.C:25
unsigned int nuniv
Definition: PlotUnfolding.C:14
void PlotUnfolding(std::string inputFile, std::string outputDir)
Definition: PlotUnfolding.C:83
std::vector< std::string > updownlabel
Definition: PlotUnfolding.C:69
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
c1
Definition: demo5.py:24
const std::string yxaxistitle
Definition: PlotUnfolding.C:22
TH2F * makeSystematic2D(TH2F *syst, TH2F *nominal)
Definition: PlotUnfolding.C:71
enum BeamMode string