7 #include "RooUnfoldResponse.h" 8 #include "RooUnfoldBayes.h" 15 const double pot = 8.09e20;
25 TH2F* hDD = (TH2F*)sig->Project3D(
"yx");
27 for(
int i = 1;
i <= sig->GetXaxis()->GetNbins();
i++){
28 for(
int j = 1;
j <= sig->GetYaxis()->GetNbins();
j++){
34 hBkgd1D->Add(hSig1D,-1);
35 int sigbin = hSig1D->FindBin(0.85);
37 float tot_sig = hSig1D->Integral(sigbin,-1);
38 float tot_bkgd = hBkgd1D->Integral(sigbin,-1);
40 if(tot_bkgd <= 0) hDD->SetBinContent(i,j,0);
41 else hDD->SetBinContent(i,j,tot_sig/tot_bkgd);
46 hDD->GetZaxis()->SetTitle(
"Signal/Background (CVNe > 0.85)");
47 hDD->GetXaxis()->SetTitle(
"cos #theta_{e}");
48 hDD->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
49 hDD->GetXaxis()->SetRangeUser(0.75,1);
50 hDD->GetYaxis()->SetRangeUser(1.0,10.0);
51 hDD->GetXaxis()->CenterTitle();
52 hDD->GetYaxis()->CenterTitle();
53 hDD->GetZaxis()->CenterTitle();
54 TCanvas* c88 =
new TCanvas(
"c88",
"c88");
55 c88->SetRightMargin(0.15);
56 c88->SetLeftMargin(0.15);
57 hDD->SetMarkerColor(kWhite);
58 gStyle->SetPaintTextFormat(
"4.2f");
59 hDD->Draw(
"COLZ TEXT");
61 c88->SaveAs(out.c_str());
71 bool shouldSave =
true)
75 TFile* fAnalysis =
new TFile((
sDir+
sAnalysis).c_str(),
"read");
84 std::vector<std::unique_ptr<ana::Spectrum>> ppfx_spects;
85 for(
int iuniv = 0; iuniv < 100; iuniv++){
87 sprintf(name,
"%s_%i_%s_%s",
"ppfx", iuniv,
88 pidName.c_str(), varName.c_str());
115 std::vector<std::unique_ptr<ana::Spectrum>>genie_spects;
116 for(
int iuniv = 0; iuniv < 100; iuniv++){
118 sprintf(name,
"%s_%i_%s_%s",
"genie", iuniv,
119 pidName.c_str(), varName.c_str());
143 std::vector<TH3F*> systs_hists;
144 std::vector<TH3F*> systs_hists_up;
145 std::vector<TH3F*> systs_hists_down;
148 sprintf(name,
"calibpos_%s_%s_%s", pidName.c_str(),
149 varName.c_str(),
chns[0].
name.c_str());
151 (fTemplates->GetDirectory(name)))
154 sprintf(name,
"lightup_%s_%s_%s", pidName.c_str(),
155 varName.c_str(),
chns[0].
name.c_str());
157 (fTemplates->GetDirectory(name)))
159 systs_hists_up.push_back(hGenieUpper);
160 systs_hists_up.push_back(hPPFXUpper);
162 sprintf(name,
"calibneg_%s_%s_%s", pidName.c_str(),
163 varName.c_str(),
chns[0].
name.c_str());
165 (fTemplates->GetDirectory(name)))
168 sprintf(name,
"lightdown_%s_%s_%s", pidName.c_str(),
169 varName.c_str(),
chns[0].
name.c_str());
171 (fTemplates->GetDirectory(name)))
173 systs_hists_down.push_back(hGenieLower);
174 systs_hists_down.push_back(hPPFXLower);
176 sprintf(name,
"calibshape_%s_%s_%s", pidName.c_str(),
177 varName.c_str(),
chns[0].
name.c_str());
179 (fTemplates->GetDirectory(name)))->ToTH3(
pot));
180 sprintf(name,
"ckv_%s_%s_%s", pidName.c_str(),
181 varName.c_str(),
chns[0].
name.c_str());
183 (fTemplates->GetDirectory(name)))->ToTH3(
pot));
185 std::vector<TH3F*> systs_calibshape = {systs_hists[0]};
186 std::vector<TH3F*> systs_ckv = {systs_hists[1]};
187 std::vector<TH3F*> systs_genie_up = {hGenieUpper};
188 std::vector<TH3F*> systs_genie_down = {hGenieLower};
189 std::vector<TH3F*> systs_ppfx_up = {hPPFXUpper};
190 std::vector<TH3F*> systs_ppfx_down = {hPPFXLower};
191 std::vector<TH3F*> systs_light_up = {systs_hists_up[1]};
192 std::vector<TH3F*> systs_light_down = {systs_hists_down[1]};
193 std::vector<TH3F*> systs_calib_up = {systs_hists_up[0]};
194 std::vector<TH3F*> systs_calib_down = {systs_hists_down[0]};
196 std::vector<TH3F*> systs_onesided;
197 std::vector<TH3F*> systs_upper;
198 std::vector<TH3F*> systs_lower;
200 if(systName ==
"tot"){
201 systs_onesided = {systs_hists[0],systs_hists[1]};
202 systs_upper = {systs_hists_up[0], systs_hists_up[1],
203 hGenieUpper,hPPFXUpper};
204 systs_lower = {systs_hists_down[0], systs_hists_down[1],
205 hGenieLower,hPPFXLower};
207 else if(systName ==
"calibshape"){
208 systs_onesided = {systs_hists[0]};
212 else if(systName ==
"ckv"){
213 systs_onesided = {systs_hists[1]};
217 else if(systName ==
"genie"){
219 systs_upper = {hGenieUpper};
220 systs_lower = {hGenieLower};
222 else if(systName ==
"ppfx"){
224 systs_upper = {hPPFXUpper};
225 systs_lower = {hPPFXLower};
227 else if(systName ==
"calib"){
229 systs_upper = {systs_hists_down[0]};
230 systs_lower = {systs_hists_down[0]};
232 else if(systName ==
"light"){
234 systs_upper = {systs_hists_up[1]};
235 systs_lower = {systs_hists_down[1]};
244 std::vector<TH3F*> nominal_hists;
246 sprintf(name,
"nominal_%s_%s_%s", pidName.c_str(),
249 (fTemplates->GetDirectory(name)))
253 TH2F* hPlot2D = (TH2F*)nominal_hists[0]->Project3D(
"yx");
255 hPlot2D->SetTitle(
"");
256 hPlot2D->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
257 hPlot2D->GetXaxis()->SetTitle(
"cos #theta_{e}");
258 hPlot2D->GetYaxis()->CenterTitle();
259 hPlot2D->GetYaxis()->SetTitleOffset(0.45);
260 hPlot2D->GetZaxis()->SetMaxDigits(3);
261 hPlot2D->GetXaxis()->CenterTitle();
262 hPlot2D->GetZaxis()->SetTitle(
"Events/8.09 #times 10^{20} POT");
263 hPlot2D->GetZaxis()->CenterTitle();
264 hPlot2D->GetXaxis()->SetRangeUser(0.75,1.00);
265 hPlot2D->GetYaxis()->SetRangeUser(1.0,6.0);
266 TH2F* hPlot2D_sig = (TH2F*)nominal_hists[1]->Project3D(
"yx");
267 TH2F* hPlot2D_sig_2 = (TH2F*)nominal_hists[7]->Project3D(
"yx");
269 hPlot2D_sig->GetXaxis()->SetRangeUser(0.75,1);
270 hPlot2D_sig->GetYaxis()->SetRangeUser(1.0,10.0);
272 hPlot2D_sig->SetTitle(
"");
273 hPlot2D_sig->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
274 hPlot2D_sig->GetXaxis()->SetTitle(
"cos #theta_{e}");
275 hPlot2D_sig->GetYaxis()->CenterTitle();
276 hPlot2D_sig->GetYaxis()->SetTitleOffset(0.45);
277 hPlot2D_sig->GetZaxis()->SetMaxDigits(3);
278 hPlot2D_sig->GetXaxis()->CenterTitle();
279 hPlot2D_sig->GetZaxis()->SetTitle(
"Signal Events/8.09 #times 10^{20} POT");
280 hPlot2D_sig->GetZaxis()->CenterTitle();
284 TCanvas* c29 =
new TCanvas(
"c28",
"c29");
285 c29->SetRightMargin(0.15);
286 hPlot2D_sig->SetMarkerColor(kWhite);
287 gStyle->SetPaintTextFormat(
"4.0f");
288 hPlot2D_sig->Draw(
"COLZ TEXT");
290 c29->SaveAs(
"signal_mc_2d.png");
294 TH3F* hNominalSignalLike =
296 hNominalSignalLike->Add(nominal_hists[2]);
297 hNominalSignalLike->Add(nominal_hists[7]);
300 "s_b_ratio_signal_region.png");
303 TCanvas *c88 =
new TCanvas(
"c88",
"c88");
305 c88->SetRightMargin(0.15);
307 hPlot2D->Draw(
"COLZ");
309 c88->SaveAs(
"total_mc_2d.png");
313 TH2F* hPlot2D_a = (TH2F*)nominal_hists[0]->Project3D(
"yz");
315 hPlot2D_a->SetTitle(
"");
316 hPlot2D_a->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
317 hPlot2D_a->GetXaxis()->SetTitle(
"Event CVNe");
318 hPlot2D_a->GetYaxis()->CenterTitle();
319 hPlot2D_a->GetYaxis()->SetTitleOffset(0.45);
320 hPlot2D_a->GetZaxis()->SetMaxDigits(3);
321 hPlot2D_a->GetXaxis()->CenterTitle();
322 hPlot2D_a->GetZaxis()->SetTitle(
"Events/8.09 #times 10^{20} POT");
323 hPlot2D_a->GetZaxis()->CenterTitle();
324 hPlot2D_a->GetXaxis()->SetRangeUser(0.,1.00);
325 hPlot2D_a->GetYaxis()->SetRangeUser(1.0,6.0);
327 TCanvas *c89 =
new TCanvas(
"c89",
"c89");
329 c89->SetRightMargin(0.15);
331 hPlot2D_a->Draw(
"COLZ");
333 c89->SaveAs(
"total_mc_2d_electron_pid.png");
337 std::vector<int> intvect = {};
344 std::vector<TH3F*> data_hists;
345 if (
dataName.find(
"ppfx") != std::string::npos) {
346 sprintf(name,
"mc_%s_%s_%s",
dataName.c_str(),
349 else if (
dataName.find(
"genie") != std::string::npos) {
350 sprintf(name,
"mc_%s_%s_%s",
dataName.c_str(),
354 sprintf(name,
"mc_%s_%s_%s",
dataName.c_str(),
358 GetDirectory(name)))->
362 std::vector<TH3F*> nominal_hists_3d;
364 sprintf(name,
"%s_%s_%s_%s",
"mc",
"nominal",
365 "analysis",
chns[
i].name.c_str());
377 std::vector<TH3F*> weighted_hists_3d =
379 systs_upper,systs_lower,nominal_hists_3d,
380 pidName,varName,out);
387 if(weighted_hists_3d.size() < 4){
396 sprintf(name,
"mc_%s_%s_signal_postfit",
dataName.c_str(), systName.c_str());
397 TH3F* hSignalOut = (TH3F*)weighted_hists_3d[1]->
Clone(name);
398 hSignalOut->SetName(name);
402 const int nx = hSignalOut->GetXaxis()->GetNbins();
403 const int ny = hSignalOut->GetYaxis()->GetNbins();
404 const int nz = hSignalOut->GetZaxis()->GetNbins();
407 TH1F* hHolder =
new TH1F(
"hHolder",
"",nx*ny*nz,0,nx*ny*nz);
409 for(
int i = 0;
i < hHolder->GetNbinsX(); ++
i){
410 const int nynz = ny*nz;
411 const int nmodnynz =
i%nynz;
412 const int ix =
i/nynz;
413 const int iy = nmodnynz/nz;
416 const double val = hSignalOut->GetBinContent(ix+1, iy+1, iz+1);
417 const double err = hSignalOut->GetBinError (ix+1, iy+1, iz+1);
419 hHolder->SetBinContent(
i+1,val);
420 hHolder->SetBinError(
i+1,err);
423 std::vector<TH1F*> vHolder;
424 for(
int j = 0;
j < (
int)weighted_hists_3d.size();
j++){
425 TH1F* hHold =
new TH1F(
ana::UniqueName().c_str(),
"",nx*ny*nz,0,nx*ny*nz);
427 for(
int i = 0;
i < hHold->GetNbinsX(); ++
i){
428 const int nynz = ny*nz;
429 const int nmodnynz =
i%nynz;
430 const int ix =
i/nynz;
431 const int iy = nmodnynz/nz;
434 const double val = weighted_hists_3d[
j]->GetBinContent(ix+1, iy+1, iz+1);
435 const double err = weighted_hists_3d[
j]->GetBinError (ix+1, iy+1, iz+1);
437 hHold->SetBinContent(
i+1,val);
438 hHold->SetBinError(
i+1,err);
440 vHolder.push_back(hHold);
441 sprintf(name,
"mc_%s_%s_%s_postfit",
dataName.c_str(),
443 vHolder[
j]->SetName(name);
446 if(shouldSave ==
true){
450 "CrossSection_TemplateFitResults_fakedata_development_stat.root";
451 TFile *
outf =
new TFile((
sDir+outfile).c_str(),
"update");
453 sprintf(name,
"mc_%s_%s_signal_postfit",
dataName.c_str(), systName.c_str());
454 sSignalOut.
SaveTo(outf->mkdir(name));
456 std::vector<Spectrum*> sHolder;
457 for(
int j = 0;
j < (
int)vHolder.size();
j++){
458 sprintf(name,
"mc_%s_%s_%s_postfit",
dataName.c_str(),
461 sHolder[
j]->SaveTo(outf->mkdir(name));
const std::string sTemplates
void PlotSignalBackground(TH3F *sig, TH3F *tot, std::string out)
Cuts and Vars for the 2020 FD DiF Study.
const ana::Binning eelecbins
void makeXSecPlots_TemplateFit(std::string pidName="prongCVN", std::string varName="double", std::string dataName="nominal", std::string systName="tot", bool shouldSave=true)
const SelDef chns[knumchns]
const ana::Binning costhetabins
Representation of a spectrum in any variable, with associated POT.
const ana::Binning pidbins
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
void SaveTo(TDirectory *dir, const std::string &name) const
std::vector< float > Spectrum
std::vector< TH3F * > BinByBinTemplateFit(TH3F *data, std::vector< TH3F * > templates, std::vector< TH3F * > systs_hists, std::vector< TH3F * > systs_hists_up, std::vector< TH3F * > systs_hists_down, std::vector< TH3F * > analysis_templates, std::string pidName, std::string varName, std::string dataName)
const std::string sAnalysis
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
std::string UniqueName()
Return a different string each time, for creating histograms.