Functions | Variables
PlotXSec.C File Reference
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Analysis/Plots.h"
#include "NDAna/numucc_inc/Multiverse.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TFile.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TLegend.h"
#include <string>
#include <vector>
#include <map>
#include <iostream>

Go to the source code of this file.

Functions

const TAxis angleAxis (1, 0.5, 1.0)
 
const TAxis muoneAxis (1, 0.5, 2.5)
 
const TAxis uncertAxis (1,-0.2, 0.2)
 
const TAxis enuAxis (1, 0.0, 5.0)
 
const TAxis q2Axis (1, 0.0, 2.8)
 
void FixAxisRangeQ2 (TH1 *hist)
 
void FixAxisRangeENu (TH1 *hist)
 
void FixAxisRange2D (TH2 *hist)
 
void AddUpDwUncerts (TH1 *destHist, TH1 *upHist, TH1 *dwHist)
 
void AddUpDwUncerts2D (TH2 *destHist, TH2 *upHist, TH2 *dwHist)
 
TH3D * Load3DXSec (string syst)
 
TH2D * XSec2D (string syst)
 
TH1D * XSecENu (string syst)
 
TH1D * XSecQ2 (string syst)
 
TH2D * XSec2DUncert (string syst)
 
TH1D * XSecENuUncert (string syst)
 
TH1D * XSecQ2Uncert (string syst)
 
void LoadNominal (string syst)
 
void LoadSyst (string syst)
 
void LoadUncerts (string syst)
 
void LoadStatUncert (string nom)
 
void AddQuadrature (TH1 *dest, TH1 *hist)
 
void AddQuadrature (TH1 *dest, vector< TH1 * > hists)
 
void AddQuadrature2D (TH2 *dest, TH2 *hist)
 
void AddQuadrature2D (TH2 *dest, vector< TH2 * > hists)
 
TH1 * GetMaxUpDown (TH1 *upHist, TH1 *dwHist)
 
TH2 * GetMaxUpDown2D (TH2 *upHist, TH2 *dwHist)
 
TH1 * GetAbsolute (TH1 *hist)
 
TH2 * GetAbsolute2D (TH2 *hist)
 
TH1 * GetSliceHist (TH2 *hist, int ibinx)
 
TH2 * GetTotalBeamUncert2D ()
 
TH1 * GetTotalBeamUncertENu ()
 
TH1 * GetTotalBeamUncertQ2 ()
 
void LimitTMuBins (TH1 *hist, int ibiny)
 
void MakeSystematicPlots (vector< string > systs, string systFilename)
 
void PlotXSec (string filename, string outputFolder)
 

Variables

map< string, stringsystTitles
 
vector< stringstandardSysts
 
map< string, pair< string, string > > standardSystsUpDown
 
const map< string, stringotherNomSysts
 
vector< stringbeamSysts
 
map< string, Color_t > systColors
 
map< int, double > maxTMuBins
 
const string var3DName = "EAvail"
 
string uncertLabel = "_uncert"
 
const vector< stringvars3D
 
TH3D * nom3DXSec
 
TH2D * nom2DXSec
 
string defaultNominal = "cv"
 
map< string, TH3D * > nominalXSecs3D
 
map< string, TH2D * > nominalXSecs2D
 
map< string, TH1D * > nominalXSecsENu
 
map< string, TH1D * > nominalXSecsQ2
 
map< string, TH3D * > all3DXSecs
 
map< string, TH2D * > xSecs2D
 
map< string, TH1D * > xSecsENu
 
map< string, TH1D * > xSecsQ2
 
map< string, TH2D * > uncerts2D
 
map< string, TH1D * > uncertsENu
 
map< string, TH1D * > uncertsQ2
 
TFile * inputFile
 
TFile * outFile
 
TCanvas * c1
 
const unsigned int N_GENIE = 1000
 
const unsigned int N_PPFX = 100
 

Function Documentation

void AddQuadrature ( TH1 *  dest,
TH1 *  hist 
)

Definition at line 437 of file PlotXSec.C.

References ana::assert(), make_syst_table_plots::ibin, cet::pow(), and std::sqrt().

Referenced by AddQuadrature(), and MakeSystematicPlots().

438 {
439  assert(dest && hist && dest->GetNbinsX() == hist->GetNbinsX());
440  for(int ibin = 0; ibin <= hist->GetNbinsX() + 1; ++ibin)
441  dest->SetBinContent(ibin, sqrt(pow(dest->GetBinContent(ibin), 2) + pow(hist->GetBinContent(ibin), 2)));
442 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
assert(nhit_max >=nhit_nbins)
void AddQuadrature ( TH1 *  dest,
vector< TH1 * >  hists 
)

Definition at line 444 of file PlotXSec.C.

References AddQuadrature(), ana::assert(), and analysePickle::hist.

445 {
446  assert(!hists.empty());
447  for (TH1 * hist : hists)
449 }
TString hists[nhists]
Definition: bdt_com.C:3
void AddQuadrature(TH1 *dest, TH1 *hist)
Definition: PlotXSec.C:437
assert(nhit_max >=nhit_nbins)
void AddQuadrature2D ( TH2 *  dest,
TH2 *  hist 
)

Definition at line 451 of file PlotXSec.C.

References ana::assert(), cet::pow(), and std::sqrt().

Referenced by AddQuadrature2D(), and MakeSystematicPlots().

452 {
453  assert(dest && hist && dest->GetNbinsX() == hist->GetNbinsX() && dest->GetNbinsY() == hist->GetNbinsY());
454  for(int ibinx = 0; ibinx <= dest->GetNbinsX() + 1; ++ibinx)
455  for(int ibiny = 0; ibiny <= dest->GetNbinsY() + 1; ++ibiny)
456  dest->SetBinContent(ibinx, ibiny, sqrt(pow(dest->GetBinContent(ibinx, ibiny), 2) + pow(hist->GetBinContent(ibinx, ibiny), 2)));
457  // for(int ibinx = 0; ibinx <= dest->GetNbinsX() + 1; ++ibinx)
458  // for (int ibiny = 0; ibiny <= dest->GetNbinsY() + 1; ++ibiny)
459  // {
460  // double addErr = hist->GetBinContent(ibinx, ibiny);
461  // double currentErr = dest->GetBinContent(ibinx, ibiny);
462  // dest->SetBinContent(ibinx, ibiny, sqrt( pow(currentErr, 2) + pow(addErr, 2)));
463  // }
464 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
assert(nhit_max >=nhit_nbins)
void AddQuadrature2D ( TH2 *  dest,
vector< TH2 * >  hists 
)

Definition at line 466 of file PlotXSec.C.

References AddQuadrature2D(), ana::assert(), and analysePickle::hist.

467 {
468  assert(!hists.empty());
469  for (TH2 * hist : hists)
471 }
TString hists[nhists]
Definition: bdt_com.C:3
assert(nhit_max >=nhit_nbins)
void AddQuadrature2D(TH2 *dest, TH2 *hist)
Definition: PlotXSec.C:451
void AddUpDwUncerts ( TH1 *  destHist,
TH1 *  upHist,
TH1 *  dwHist 
)

Definition at line 177 of file PlotXSec.C.

References abs(), ana::assert(), cet::sqlite::max(), cet::pow(), and std::sqrt().

Referenced by GetTotalBeamUncertENu(), and GetTotalBeamUncertQ2().

178 {
179  assert(destHist->GetNbinsX() == upHist->GetNbinsX() && destHist->GetNbinsX() == dwHist->GetNbinsX());
180  for(int ibinx = 0; ibinx <= destHist->GetNbinsX() + 1; ++ibinx)
181  {
182  double upErr = abs(upHist->GetBinContent(ibinx));
183  double dwErr = abs(dwHist->GetBinContent(ibinx));
184  double addErr = max(upErr, dwErr);
185  double currentErr = destHist->GetBinContent(ibinx);
186  destHist->SetBinContent(ibinx, sqrt( pow(currentErr, 2) + pow(addErr, 2)));
187  }
188 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void abs(TH1 *hist)
assert(nhit_max >=nhit_nbins)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void AddUpDwUncerts2D ( TH2 *  destHist,
TH2 *  upHist,
TH2 *  dwHist 
)

Definition at line 190 of file PlotXSec.C.

References abs(), ana::assert(), cet::sqlite::max(), cet::pow(), and std::sqrt().

Referenced by GetTotalBeamUncert2D().

191 {
192  assert(destHist->GetNbinsX() == upHist->GetNbinsX() && destHist->GetNbinsX() == dwHist->GetNbinsX());
193  assert(destHist->GetNbinsY() == upHist->GetNbinsY() && destHist->GetNbinsY() == dwHist->GetNbinsY());
194  for(int ibinx = 0; ibinx <= destHist->GetNbinsX() + 1; ++ibinx)
195  for (int ibiny = 0; ibiny <= destHist->GetNbinsY() + 1; ++ibiny)
196  {
197  double upErr = abs(upHist->GetBinContent(ibinx, ibiny));
198  double dwErr = abs(dwHist->GetBinContent(ibinx, ibiny));
199  double addErr = max(upErr, dwErr);
200  double currentErr = destHist->GetBinContent(ibinx, ibiny);
201  destHist->SetBinContent(ibinx, ibiny, sqrt( pow(currentErr, 2) + pow(addErr, 2)));
202  }
203 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void abs(TH1 *hist)
assert(nhit_max >=nhit_nbins)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const TAxis angleAxis ( ,
0.  5,
1.  0 
)

Referenced by FixAxisRange2D().

const TAxis enuAxis ( ,
0.  0,
5.  0 
)

Referenced by FixAxisRangeENu().

void FixAxisRange2D ( TH2 *  hist)

Definition at line 168 of file PlotXSec.C.

References angleAxis(), and muoneAxis().

Referenced by XSec2D(), and XSec2DUncert().

169 {
170  hist->GetXaxis()->SetRangeUser(angleAxis.GetXmin(), angleAxis.GetXmax());
171  hist->GetYaxis()->SetRangeUser(muoneAxis.GetXmin(), muoneAxis.GetXmax());
172  // string histName = hist->GetName();
173  // if (histName.find(uncertLabel) != string::npos)
174  // hist->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
175 }
const TAxis angleAxis(1, 0.5, 1.0)
const TAxis muoneAxis(1, 0.5, 2.5)
void FixAxisRangeENu ( TH1 *  hist)

Definition at line 159 of file PlotXSec.C.

References enuAxis().

Referenced by XSecENu(), and XSecENuUncert().

160 {
161  hist->GetXaxis()->SetRangeUser(enuAxis.GetXmin(), enuAxis.GetXmax());
162  // string histName = hist->GetName();
163  // if (histName.find(uncertLabel) != string::npos)
164  // hist->GetYaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
165 }
const TAxis enuAxis(1, 0.0, 5.0)
void FixAxisRangeQ2 ( TH1 *  hist)

Definition at line 151 of file PlotXSec.C.

References q2Axis().

Referenced by XSecQ2(), and XSecQ2Uncert().

152 {
153  hist->GetXaxis()->SetRangeUser(q2Axis.GetXmin(), q2Axis.GetXmax());
154  // string histName = hist->GetName();
155  // if (histName.find(uncertLabel) != string::npos)
156  // hist->GetYaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
157 }
const TAxis q2Axis(1, 0.0, 2.8)
TH1* GetAbsolute ( TH1 *  hist)

Definition at line 493 of file PlotXSec.C.

References abs(), make_syst_table_plots::ibin, and ana::UniqueName().

Referenced by GetSliceHist(), and MakeSystematicPlots().

494 {
495  TH1 * histAbs = (TH1*) hist->Clone(ana::UniqueName().c_str());
496  for(int ibin = 0; ibin <= hist->GetNbinsX(); ++ibin)
497  histAbs->SetBinContent(ibin, abs(hist->GetBinContent(ibin)));
498  return histAbs;
499 }
void abs(TH1 *hist)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH2* GetAbsolute2D ( TH2 *  hist)

Definition at line 502 of file PlotXSec.C.

References abs(), and ana::UniqueName().

Referenced by MakeSystematicPlots().

503 {
504  TH2 * histAbs = (TH2*) hist->Clone(ana::UniqueName().c_str());
505  for(int ibinx = 0; ibinx <= hist->GetNbinsX(); ++ibinx)
506  for(int ibiny = 0; ibiny <= hist->GetNbinsY(); ++ibiny)
507  histAbs->SetBinContent(ibinx, ibiny, abs(hist->GetBinContent(ibinx, ibiny)));
508  return histAbs;
509 }
void abs(TH1 *hist)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH1* GetMaxUpDown ( TH1 *  upHist,
TH1 *  dwHist 
)

Definition at line 473 of file PlotXSec.C.

References abs(), make_syst_table_plots::ibin, cet::sqlite::max(), and ana::UniqueName().

Referenced by MakeSystematicPlots().

474 {
475  TH1 * upDwError = (TH1*) upHist->Clone(ana::UniqueName().c_str());
476  upDwError->Reset();
477  for (int ibin = 0; ibin <= upHist->GetNbinsX() + 1; ++ibin)
478  upDwError->SetBinContent(ibin, max(abs(upHist->GetBinContent(ibin)), abs(dwHist->GetBinContent(ibin))));
479  return upDwError;
480 }
void abs(TH1 *hist)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH2* GetMaxUpDown2D ( TH2 *  upHist,
TH2 *  dwHist 
)

Definition at line 482 of file PlotXSec.C.

References abs(), cet::sqlite::max(), and ana::UniqueName().

Referenced by MakeSystematicPlots().

483 {
484  TH2 * upDwError = (TH2*) upHist->Clone(ana::UniqueName().c_str());
485  upDwError->Reset();
486  for (int ibinx = 0; ibinx <= upHist->GetNbinsX() + 1; ++ibinx)
487  for (int ibiny = 0; ibiny <= upHist->GetNbinsY() + 1; ++ibiny)
488  upDwError->SetBinContent(ibinx, ibiny, max(abs(upHist->GetBinContent(ibinx, ibiny)), abs(dwHist->GetBinContent(ibinx, ibiny))));
489  return upDwError;
490 }
void abs(TH1 *hist)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH1* GetSliceHist ( TH2 *  hist,
int  ibinx 
)

Definition at line 512 of file PlotXSec.C.

References GetAbsolute(), string, art::to_string(), and ana::UniqueName().

Referenced by MakeSystematicPlots().

513 {
514  // TH1 * histSlice = hist->ProjectionY((hist->GetName() + string("_") + to_string(ibinx)).c_str(), ibinx, ibinx);
515  TH1 * histSlice = hist->ProjectionY((ana::UniqueName()).c_str(), ibinx, ibinx);
516  TH1 * absSlice = GetAbsolute(histSlice);
517  absSlice->SetName((histSlice->GetName() + string("_abs")).c_str());
518  absSlice->SetTitle((to_string(histSlice->GetXaxis()->GetXmin()) + " < cos(#theta_{#mu}) < " + to_string(histSlice->GetXaxis()->GetXmax())).c_str());
519  return absSlice;
520 }
TH1 * GetAbsolute(TH1 *hist)
Definition: PlotXSec.C:493
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string
TH2* GetTotalBeamUncert2D ( )

Definition at line 523 of file PlotXSec.C.

References AddUpDwUncerts2D(), beamSysts, plot_validation_datamc::Clone(), nominalXSecs2D, and uncerts2D.

Referenced by MakeSystematicPlots().

524 {
525  TH2D * totalBeamUncert = (TH2D*) nominalXSecs2D["cv"]->Clone("totalBeamUncert2D");
526  totalBeamUncert->Reset();
527  for (string beamSyst : beamSysts)
528  {
529  TH2 * upSyst = uncerts2D[beamSyst + "_Up"];
530  TH2 * dwSyst = uncerts2D[beamSyst + "_Dw"];
531  AddUpDwUncerts2D(totalBeamUncert, upSyst, dwSyst);
532  }
533  uncerts2D["totalBeam"] = totalBeamUncert;
534  return totalBeamUncert;
535 }
void AddUpDwUncerts2D(TH2 *destHist, TH2 *upHist, TH2 *dwHist)
Definition: PlotXSec.C:190
map< string, TH2D * > uncerts2D
Definition: PlotXSec.C:129
vector< string > beamSysts
Definition: PlotXSec.C:65
map< string, TH2D * > nominalXSecs2D
Definition: PlotXSec.C:120
TH1* GetTotalBeamUncertENu ( )

Definition at line 537 of file PlotXSec.C.

References AddUpDwUncerts(), beamSysts, plot_validation_datamc::Clone(), nominalXSecsENu, and uncertsENu.

Referenced by MakeSystematicPlots().

538 {
539  TH1D * totalBeamUncert = (TH1D*) nominalXSecsENu["cv"]->Clone("totalBeamUncertENu");
540  totalBeamUncert->Reset();
541  for (string beamSyst : beamSysts)
542  {
543  TH1 * upSyst = uncertsENu[beamSyst + "_Up"];
544  TH1 * dwSyst = uncertsENu[beamSyst + "_Dw"];
545  AddUpDwUncerts(totalBeamUncert, upSyst, dwSyst);
546  }
547  uncertsENu["totalBeam"] = totalBeamUncert;
548  return totalBeamUncert;
549 }
map< string, TH1D * > uncertsENu
Definition: PlotXSec.C:130
vector< string > beamSysts
Definition: PlotXSec.C:65
void AddUpDwUncerts(TH1 *destHist, TH1 *upHist, TH1 *dwHist)
Definition: PlotXSec.C:177
map< string, TH1D * > nominalXSecsENu
Definition: PlotXSec.C:121
TH1* GetTotalBeamUncertQ2 ( )

Definition at line 551 of file PlotXSec.C.

References AddUpDwUncerts(), beamSysts, plot_validation_datamc::Clone(), nominalXSecsQ2, and uncertsQ2.

Referenced by MakeSystematicPlots().

552 {
553  TH1D * totalBeamUncert = (TH1D*) nominalXSecsQ2["cv"]->Clone("totalBeamUncertQ2");
554  totalBeamUncert->Reset();
555  for (string beamSyst : beamSysts)
556  {
557  TH1 * upSyst = uncertsQ2[beamSyst + "_Up"];
558  TH1 * dwSyst = uncertsQ2[beamSyst + "_Dw"];
559  AddUpDwUncerts(totalBeamUncert, upSyst, dwSyst);
560  }
561  uncertsQ2["totalBeam"] = totalBeamUncert;
562  return totalBeamUncert;
563 }
map< string, TH1D * > nominalXSecsQ2
Definition: PlotXSec.C:122
vector< string > beamSysts
Definition: PlotXSec.C:65
void AddUpDwUncerts(TH1 *destHist, TH1 *upHist, TH1 *dwHist)
Definition: PlotXSec.C:177
map< string, TH1D * > uncertsQ2
Definition: PlotXSec.C:131
void LimitTMuBins ( TH1 *  hist,
int  ibiny 
)

Definition at line 566 of file PlotXSec.C.

References make_syst_table_plots::ibin, and maxTMuBins.

Referenced by MakeSystematicPlots().

567 {
568  if (maxTMuBins.find(ibiny) != maxTMuBins.end())
569  for (int ibin = 1; ibin <= hist->GetNbinsX(); ++ibin)
570  if (hist->GetXaxis()->GetBinCenter(ibin) >= maxTMuBins.at(ibiny)){
571  hist->SetBinContent(ibin, 0);
572  hist->SetBinError(ibin, 0);
573  }
574 }
map< int, double > maxTMuBins
Definition: PlotXSec.C:94
TH3D* Load3DXSec ( string  syst)

Definition at line 206 of file PlotXSec.C.

References all3DXSecs, inputFile, systTitles, and var3DName.

Referenced by LoadNominal(), LoadSyst(), and XSec2D().

207 {
208  if (all3DXSecs.find(syst) != all3DXSecs.end())
209  return all3DXSecs.at(syst);
210 
211  TH3D* xsec3D = (TH3D*) inputFile->Get((syst + "/" + var3DName + "/xsec_" + syst + "_" + var3DName).c_str());
212  all3DXSecs[syst] = xsec3D;
213  if (!xsec3D)
214  return NULL;
215  xsec3D->SetTitle((systTitles[syst] + " XSec").c_str());
216  xsec3D->Write();
217  return xsec3D;
218 }
TFile * inputFile
Definition: PlotXSec.C:134
map< string, TH3D * > all3DXSecs
Definition: PlotXSec.C:125
map< string, string > systTitles
Definition: PlotXSec.C:23
const string var3DName
Definition: PlotXSec.C:106
void LoadNominal ( string  syst)

Definition at line 348 of file PlotXSec.C.

References dir, Load3DXSec(), nominalXSecs2D, nominalXSecs3D, nominalXSecsENu, nominalXSecsQ2, outFile, XSec2D(), XSecENu(), and XSecQ2().

Referenced by PlotXSec().

349 {
350  if (nominalXSecs3D.find(syst) != nominalXSecs3D.end())
351  return;
352 
353  TDirectory * dir = outFile->GetDirectory(syst.c_str());
354  if (!dir)
355  dir = outFile->mkdir(syst.c_str());
356  dir->cd();
357 
358  TH3D * nomXSec = Load3DXSec(syst);
359  TH2D * nomXSec2D = XSec2D(syst);
360  TH1D * nomXSecENu = XSecENu(syst);
361  TH1D * nomXSecQ2 = XSecQ2(syst);
362 
363  nominalXSecs3D[syst] = nomXSec;
364  nominalXSecs2D[syst] = nomXSec2D;
365  nominalXSecsENu[syst] = nomXSecENu;
366  nominalXSecsQ2 [syst] = nomXSecQ2;
367 
368  outFile->cd();
369 }
TH3D * Load3DXSec(string syst)
Definition: PlotXSec.C:206
TH2D * XSec2D(string syst)
Definition: PlotXSec.C:221
map< string, TH1D * > nominalXSecsQ2
Definition: PlotXSec.C:122
TFile * outFile
Definition: PlotXSec.C:135
TH1D * XSecENu(string syst)
Definition: PlotXSec.C:244
map< string, TH2D * > nominalXSecs2D
Definition: PlotXSec.C:120
TDirectory * dir
Definition: macro.C:5
map< string, TH3D * > nominalXSecs3D
Definition: PlotXSec.C:119
TH1D * XSecQ2(string syst)
Definition: PlotXSec.C:258
map< string, TH1D * > nominalXSecsENu
Definition: PlotXSec.C:121
void LoadStatUncert ( string  nom)

Definition at line 398 of file PlotXSec.C.

References plot_validation_datamc::Clone(), nominalXSecs2D, nominalXSecsENu, nominalXSecsQ2, uncerts2D, uncertsENu, and uncertsQ2.

Referenced by PlotXSec().

399 {
400  if (uncerts2D. find("stat") != uncerts2D.end())
401  return; // Already loaded
402 
403  // Create Stat Uncerts
404  TH2D * statUncert2D = NULL;
405  if (nominalXSecs2D.find(nom) != nominalXSecs2D.end() && nominalXSecs2D[nom])
406  statUncert2D = (TH2D*) nominalXSecs2D[nom]->Clone("stat_uncert_2D");
407  for (int ibinx = 0; statUncert2D && ibinx <= statUncert2D->GetNbinsX() + 1; ++ibinx)
408  for (int ibiny = 0; statUncert2D && ibiny <= statUncert2D->GetNbinsY() + 1; ++ibiny){
409  double content = statUncert2D->GetBinContent(ibinx, ibiny);
410  double uncert = statUncert2D->GetBinError (ibinx, ibiny);
411  statUncert2D->SetTitle("Stat Uncert;#cos{#theta_{#mu}};T_{#mu};;");
412  statUncert2D->SetBinContent(ibinx, ibiny, content == 0.0 ? 0.0 : uncert / content);
413  statUncert2D->SetBinError(ibinx, ibiny, 0);
414  }
415  TH1D * statUncertENu = (TH1D*) nominalXSecsENu[nom]->Clone("stat_uncert_ENu");
416  for (int ibinx = 0; ibinx <= statUncertENu->GetNbinsX() + 1; ++ibinx){
417  double content = statUncertENu->GetBinContent(ibinx);
418  double uncert = statUncertENu->GetBinError(ibinx);
419  statUncertENu->SetTitle("Statistical Uncertainty;E_{nu} (GeV);;;");
420  statUncertENu->SetBinContent(ibinx, content == 0.0 ? 0.0 : uncert / content);
421  statUncertENu->SetBinError(ibinx, 0);
422  }
423  TH1D * statUncertQ2 = (TH1D*) nominalXSecsQ2[nom]->Clone("stat_uncert_Q2");
424  for (int ibinx = 0; ibinx <= statUncertQ2->GetNbinsX() + 1; ++ibinx){
425  double content = statUncertQ2->GetBinContent(ibinx);
426  double uncert = statUncertQ2->GetBinError(ibinx);
427  statUncertQ2->SetTitle("Statistical Uncertainty;Q^{2} (GeV^{2});;;");
428  statUncertQ2->SetBinContent(ibinx, content == 0.0 ? 0.0 : uncert / content);
429  statUncertQ2->SetBinError(ibinx, 0);
430  }
431  if (statUncert2D)
432  uncerts2D ["stat"] = statUncert2D;
433  uncertsENu["stat"] = statUncertENu;
434  uncertsQ2 ["stat"] = statUncertQ2;
435 }
map< string, TH1D * > uncertsENu
Definition: PlotXSec.C:130
map< string, TH1D * > nominalXSecsQ2
Definition: PlotXSec.C:122
map< string, TH2D * > uncerts2D
Definition: PlotXSec.C:129
map< string, TH1D * > uncertsQ2
Definition: PlotXSec.C:131
map< string, TH2D * > nominalXSecs2D
Definition: PlotXSec.C:120
map< string, TH1D * > nominalXSecsENu
Definition: PlotXSec.C:121
void LoadSyst ( string  syst)

Definition at line 371 of file PlotXSec.C.

References dir, Load3DXSec(), outFile, XSec2D(), XSecENu(), and XSecQ2().

Referenced by PlotXSec().

372 {
373  TDirectory * dir = outFile->GetDirectory(syst.c_str());
374  if (!dir)
375  dir = outFile->mkdir(syst.c_str());
376  dir->cd();
377  TH3D * xSec3D = Load3DXSec(syst);
378  TH2D * xSec2D = XSec2D(syst);
379  TH1D * xSecENu = XSecENu(syst);
380  TH1D * xSecQ2 = XSecQ2(syst);
381 
382  outFile->cd();
383 }
TH3D * Load3DXSec(string syst)
Definition: PlotXSec.C:206
TH2D * XSec2D(string syst)
Definition: PlotXSec.C:221
TFile * outFile
Definition: PlotXSec.C:135
TH1D * XSecENu(string syst)
Definition: PlotXSec.C:244
TDirectory * dir
Definition: macro.C:5
TH1D * XSecQ2(string syst)
Definition: PlotXSec.C:258
void LoadUncerts ( string  syst)

Definition at line 385 of file PlotXSec.C.

References dir, outFile, XSec2DUncert(), XSecENuUncert(), and XSecQ2Uncert().

Referenced by PlotXSec().

386 {
387  TDirectory * dir = outFile->GetDirectory(syst.c_str());
388  if (!dir)
389  dir = outFile->mkdir(syst.c_str());
390  dir->cd();
391  TH2 * uncert2D = XSec2DUncert(syst);
392  TH1 * uncertENu = XSecENuUncert(syst);
393  TH1 * uncertQ2 = XSecQ2Uncert(syst);
394 
395  outFile->cd();
396 }
TH2D * XSec2DUncert(string syst)
Definition: PlotXSec.C:272
TH1D * XSecQ2Uncert(string syst)
Definition: PlotXSec.C:324
TFile * outFile
Definition: PlotXSec.C:135
TDirectory * dir
Definition: macro.C:5
TH1D * XSecENuUncert(string syst)
Definition: PlotXSec.C:300
void MakeSystematicPlots ( vector< string systs,
string  systFilename 
)

Definition at line 579 of file PlotXSec.C.

References AddQuadrature(), AddQuadrature2D(), ana::AutoPlaceLegend(), plot_validation_datamc::Clone(), om::cout, allTimeWatchdog::endl, GetAbsolute(), GetAbsolute2D(), GetMaxUpDown(), GetMaxUpDown2D(), ana::Multiverse::GetMinusOneSigmaShift(), ana::Multiverse::GetPlusOneSigmaShift(), GetSliceHist(), GetTotalBeamUncert2D(), GetTotalBeamUncertENu(), GetTotalBeamUncertQ2(), analysePickle::hist, hists, MECModelEnuComparisons::leg, LimitTMuBins(), N_GENIE, N_PPFX, nominalXSecs2D, nominalXSecsENu, nominalXSecsQ2, push_back(), standardSysts, standardSystsUpDown, systColors, systTitles, art::to_string(), uncerts2D, uncertsENu, uncertsQ2, xSecs2D, xSecsENu, and xSecsQ2.

Referenced by PlotXSec().

580 {
581  TH2 * totalUncert2D = NULL;
582  if (nominalXSecs2D.find("cv") != nominalXSecs2D.end() && nominalXSecs2D["cv"])
583  totalUncert2D = (TH2*) nominalXSecs2D["cv"] ->Clone("totalUncert2D");
584  TH1 * totalUncertENu = (TH1*) nominalXSecsENu["cv"]->Clone("totalUncertENu");
585  TH1 * totalUncertQ2 = (TH1*) nominalXSecsQ2["cv"] ->Clone("totalUncertQ2");
586 
587  // Switch the
588  if (totalUncert2D)
589  totalUncert2D->Reset();
590  totalUncertENu ->Reset();
591  totalUncertQ2 ->Reset();
592 
593  map<string, TH2*> uncert2DPlots;
594  map<string, TH1*> uncertENuPlots;
595  map<string, TH1*> uncertQ2Plots;
596 
597  for (string syst : systs)
598  {
599  TH2 * uncert2D = NULL;
600  TH1 * uncertENu = NULL;
601  TH1 * uncertQ2 = NULL;
602  if(std::find(standardSysts.begin(), standardSysts.end(), syst) != standardSysts.end() || syst == "stat")
603  {
604  if (uncerts2D.find(syst) != uncerts2D.end() && uncerts2D[syst]){
605  cout << "Found 2D uncert for: " << syst << endl;
606  uncert2D = GetAbsolute2D(uncerts2D[syst]);
607  }
608  uncertENu = GetAbsolute(uncertsENu[syst]);
609  uncertQ2 = GetAbsolute(uncertsQ2[syst]);
610  }
611  else if (standardSystsUpDown.find(syst) != standardSystsUpDown.end())
612  {
613  string upSyst = standardSystsUpDown[syst].first;
614  string dwSyst = standardSystsUpDown[syst].second;
615  if (uncerts2D.find(upSyst) != uncerts2D.end() && uncerts2D.find(dwSyst) != uncerts2D.end() && uncerts2D[upSyst] && uncerts2D[dwSyst])
616  uncert2D = GetMaxUpDown2D(uncerts2D[upSyst], uncerts2D[dwSyst]);
617  uncertENu = GetMaxUpDown(uncertsENu[upSyst], uncertsENu[dwSyst]);
618  uncertQ2 = GetMaxUpDown(uncertsQ2[upSyst], uncertsQ2[dwSyst]);
619  }
620  else if (syst == "beam")
621  {
622  if (uncerts2D.find("2kA_Up") != uncerts2D.end() && uncerts2D["2kA_Up"])
623  uncert2D = GetTotalBeamUncert2D();
624  uncertENu = GetTotalBeamUncertENu();
625  uncertQ2 = GetTotalBeamUncertQ2();
626  }
627  else if (syst == "genie")
628  {
629  vector<TH1*> genieUnivs2D;
630  vector<TH1*> genieUnivsENu;
631  vector<TH1*> genieUnivsQ2;
632  for (unsigned int iuniv = 0; iuniv < N_GENIE; ++iuniv)
633  {
634  string syst = "genie" + to_string(iuniv);
635  if(xSecs2D.find(syst) != xSecs2D.end() && uncerts2D[syst])
636  genieUnivs2D. push_back(xSecs2D[syst]);
637  genieUnivsENu.push_back(xSecsENu[syst]);
638  genieUnivsQ2. push_back(xSecsQ2[syst]);
639  }
640  ana::Multiverse * genieVerse2D = NULL;
641  if (!genieUnivs2D.empty())
642  genieVerse2D = new ana::Multiverse(genieUnivs2D, 2);
643  ana::Multiverse * genieVerseENu = new ana::Multiverse(genieUnivsENu, 1);
644  ana::Multiverse * genieVerseQ2 = new ana::Multiverse(genieUnivsQ2, 1);
645  if (genieVerse2D)
646  uncert2D = GetMaxUpDown2D(
647  new TH2D(*(TH2D*)genieVerse2D->GetPlusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]),
648  new TH2D(*(TH2D*)genieVerse2D->GetMinusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]));
649  uncertENu = GetMaxUpDown(
650  new TH1D(*(TH1D*)genieVerseENu->GetPlusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]),
651  new TH1D(*(TH1D*)genieVerseENu->GetMinusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]));
652  uncertQ2 = GetMaxUpDown(
653  new TH1D(*(TH1D*)genieVerseQ2->GetPlusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]),
654  new TH1D(*(TH1D*)genieVerseQ2->GetMinusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]));
655 
656  if (uncert2D)
657  uncert2D->Divide(nominalXSecs2D["cv"]);
658  uncertENu->Divide(nominalXSecsENu["cv"]);
659  uncertQ2->Divide(nominalXSecsQ2["cv"]);
660  }
661  else if (syst == "ppfx")
662  {
663  vector<TH1*> ppfxUnivs2D;
664  vector<TH1*> ppfxUnivsENu;
665  vector<TH1*> ppfxUnivsQ2;
666  for (unsigned int iuniv = 0; iuniv < N_PPFX; ++iuniv)
667  {
668  string syst = "ppfx" + to_string(iuniv);
669  if(xSecs2D.find(syst) != xSecs2D.end() && uncerts2D[syst])
670  ppfxUnivs2D. push_back(xSecs2D[syst]);
671  ppfxUnivsENu.push_back(xSecsENu[syst]);
672  ppfxUnivsQ2. push_back(xSecsQ2[syst]);
673  }
674  ana::Multiverse * ppfxVerse2D = NULL;
675  if (!ppfxUnivs2D.empty())
676  ppfxVerse2D = new ana::Multiverse(ppfxUnivs2D, 2);
677  ana::Multiverse * ppfxVerseENu = new ana::Multiverse(ppfxUnivsENu, 1);
678  ana::Multiverse * ppfxVerseQ2 = new ana::Multiverse(ppfxUnivsQ2, 1);
679  if (ppfxVerse2D)
680  uncert2D = GetMaxUpDown2D(
681  new TH2D(*(TH2D*)ppfxVerse2D->GetPlusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]),
682  new TH2D(*(TH2D*)ppfxVerse2D->GetMinusOneSigmaShift(nominalXSecs2D["cv"]) - *nominalXSecs2D["cv"]));
683  uncertENu = GetMaxUpDown(
684  new TH1D(*(TH1D*)ppfxVerseENu->GetPlusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]),
685  new TH1D(*(TH1D*)ppfxVerseENu->GetMinusOneSigmaShift(nominalXSecsENu["cv"]) - *nominalXSecsENu["cv"]));
686  uncertQ2 = GetMaxUpDown(
687  new TH1D(*(TH1D*)ppfxVerseQ2->GetPlusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]),
688  new TH1D(*(TH1D*)ppfxVerseQ2->GetMinusOneSigmaShift(nominalXSecsQ2["cv"]) - *nominalXSecsQ2["cv"]));
689 
690  if (uncert2D)
691  uncert2D->Divide(nominalXSecs2D["cv"]);
692  uncertENu->Divide(nominalXSecsENu["cv"]);
693  uncertQ2->Divide(nominalXSecsQ2["cv"]);
694  }
695  if (totalUncert2D && uncert2D)
696  AddQuadrature2D(totalUncert2D, uncert2D);
697  AddQuadrature(totalUncertENu, uncertENu);
698  AddQuadrature(totalUncertQ2, uncertQ2);
699  if (uncert2D) uncert2D ->SetLineColor(systColors[syst]);
700  if (uncertENu) uncertENu->SetLineColor(systColors[syst]);
701  if (uncertQ2) uncertQ2 ->SetLineColor(systColors[syst]);
702  uncert2DPlots[syst] = uncert2D;
703  uncertENuPlots[syst] = uncertENu;
704  uncertQ2Plots[syst] = uncertQ2;
705  }
706 
707  // Plotting
708  TCanvas * systCanv = new TCanvas();
709  systCanv->SetLeftMargin(0.12);
710  systCanv->SetBottomMargin(0.12);
711  systCanv->SaveAs((systFilename + "[").c_str());
712  TLegend * leg = new TLegend();
713 
714  // 2D Slices
715  if (totalUncert2D){
716  for (int ibinx = 2; ibinx <= totalUncert2D->GetNbinsX(); ++ibinx)
717  {
718  TH1 * totalUncertSlice = GetSliceHist(totalUncert2D, ibinx);
719  LimitTMuBins(totalUncertSlice, ibinx);
720  totalUncertSlice->SetTitle(("Uncertainties: " +
721  to_string(totalUncert2D->GetXaxis()->GetBinLowEdge(ibinx)) + " < cos(#theta_{#mu}) < " +
722  to_string(totalUncert2D->GetXaxis()->GetBinUpEdge(ibinx))).c_str());
723  totalUncertSlice->GetYaxis()->SetRangeUser(0.0, 0.38);
724  totalUncertSlice->Draw("hist");
725  vector<pair<TH1*, string> > histsForLeg;
726  for (string syst : systs){
727  TH1 * histSlice = GetSliceHist(uncert2DPlots[syst], ibinx);
728  LimitTMuBins(histSlice, ibinx);
729  histSlice->SetLineColor(systColors[syst]);
730  histSlice->Draw("hist same");
731  histsForLeg.push_back({histSlice, systTitles[syst]});
732  }
733  leg = ana::AutoPlaceLegend(0.2, 0.35);
734  leg->AddEntry(totalUncertENu, "Total Uncert");
735  for (pair<TH1*, string> hists : histsForLeg)
736  leg->AddEntry(hists.first, hists.second.c_str());
737  leg->Draw();
738  systCanv->SaveAs(systFilename.c_str());
739  systCanv->SaveAs(("plots/SliceSyst_" + to_string(ibinx) + ".png").c_str());
740  leg->Clear();
741  delete leg;
742  }
743  }
744  else{
745  // ENu Plots
746  totalUncertENu->SetTitle("Uncertainties - ENu");
747  totalUncertENu->GetYaxis()->SetTitle("#frac{#delta #sigma}{#sigma}");
748  totalUncertENu->GetYaxis()->CenterTitle();
749  totalUncertENu->Draw("hist");
750  for (string syst : systs){
751  TH1 * hist = uncertENuPlots[syst];
752  hist->SetLineColor(systColors[syst]);
753  hist->Draw("hist same");
754  }
755  leg = ana::AutoPlaceLegend(0.2, 0.35);
756  leg->AddEntry(totalUncertENu, "Total Uncert");
757  for (string syst : systs)
758  leg->AddEntry(uncertENuPlots[syst], systTitles[syst].c_str());
759  leg->Draw();
760  systCanv->SaveAs(systFilename.c_str());
761  leg->Clear();
762  delete leg;
763 
764  // Q2 Plots
765  totalUncertQ2->SetTitle("Uncertainties - Q2");
766  totalUncertQ2->GetYaxis()->SetTitle("#frac{#delta #sigma}{#sigma}");
767  totalUncertQ2->GetYaxis()->CenterTitle();
768  totalUncertQ2->Draw("hist");
769  for (string syst : systs){
770  TH1 * hist = uncertQ2Plots[syst];
771  hist->SetLineColor(systColors[syst]);
772  hist->Draw("hist same");
773  }
774  leg = ana::AutoPlaceLegend(0.2, 0.35);
775  leg->AddEntry(totalUncertQ2, "Total Uncert");
776  for (string syst : systs)
777  leg->AddEntry(uncertQ2Plots[syst], systTitles[syst].c_str());
778  leg->Draw();
779  systCanv->SaveAs(systFilename.c_str());
780  leg->Clear();
781  delete leg;
782  }
783 
784  systCanv->SaveAs((systFilename + "]").c_str());
785 
786  systCanv->Clear();
787  delete systCanv;
788 }
TH2 * GetMaxUpDown2D(TH2 *upHist, TH2 *dwHist)
Definition: PlotXSec.C:482
ratio_hxv Divide(hxv, goal_hxv)
map< string, TH1D * > uncertsENu
Definition: PlotXSec.C:130
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:443
map< string, TH1D * > nominalXSecsQ2
Definition: PlotXSec.C:122
TH1 * GetMaxUpDown(TH1 *upHist, TH1 *dwHist)
Definition: PlotXSec.C:473
vector< string > standardSysts
Definition: PlotXSec.C:41
TH2 * GetAbsolute2D(TH2 *hist)
Definition: PlotXSec.C:502
map< string, TH1D * > xSecsENu
Definition: PlotXSec.C:127
TH2 * GetTotalBeamUncert2D()
Definition: PlotXSec.C:523
TH1 * GetAbsolute(TH1 *hist)
Definition: PlotXSec.C:493
void LimitTMuBins(TH1 *hist, int ibiny)
Definition: PlotXSec.C:566
TString hists[nhists]
Definition: bdt_com.C:3
TH1 * GetTotalBeamUncertQ2()
Definition: PlotXSec.C:551
TH1 * GetTotalBeamUncertENu()
Definition: PlotXSec.C:537
map< string, string > systTitles
Definition: PlotXSec.C:23
if(dump)
map< string, TH2D * > uncerts2D
Definition: PlotXSec.C:129
const std::string cv[Ncv]
const unsigned int N_PPFX
Definition: PlotXSec.C:149
hmean SetLineColor(4)
base_types push_back(int_type())
const unsigned int N_GENIE
Definition: PlotXSec.C:148
map< string, Color_t > systColors
Definition: PlotXSec.C:80
map< string, TH1D * > uncertsQ2
Definition: PlotXSec.C:131
OStream cout
Definition: OStream.cxx:6
map< string, TH2D * > nominalXSecs2D
Definition: PlotXSec.C:120
void AddQuadrature(TH1 *dest, TH1 *hist)
Definition: PlotXSec.C:437
TH1 * GetSliceHist(TH2 *hist, int ibinx)
Definition: PlotXSec.C:512
map< string, TH2D * > xSecs2D
Definition: PlotXSec.C:126
map< string, TH1D * > xSecsQ2
Definition: PlotXSec.C:128
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:717
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
cosmicTree SaveAs("cosmicTree.root")
map< string, pair< string, string > > standardSystsUpDown
Definition: PlotXSec.C:48
void AddQuadrature2D(TH2 *dest, TH2 *hist)
Definition: PlotXSec.C:451
map< string, TH1D * > nominalXSecsENu
Definition: PlotXSec.C:121
const TAxis muoneAxis ( ,
0.  5,
2.  5 
)

Referenced by FixAxisRange2D().

void PlotXSec ( string  filename,
string  outputFolder 
)

Definition at line 796 of file PlotXSec.C.

References beamSysts, c1, make_mec_shifts_plots::GetMaximum(), analysePickle::hist, inputFile, kBlue, kRed, LoadNominal(), LoadStatUncert(), LoadSyst(), LoadUncerts(), MakeSystematicPlots(), N_GENIE, N_PPFX, nominalXSecs2D, nominalXSecsENu, nominalXSecsQ2, otherNomSysts, outFile, standardSysts, standardSystsUpDown, art::to_string(), xSecsENu, and xSecsQ2.

797 {
798  inputFile = new TFile(filename.c_str(), "READ");
799  outFile = new TFile((outputFolder + "/xsecUncertResults.root").c_str(), "RECREATE");
800 
801  c1 = new TCanvas();
802  string outputPDF = outputFolder + "/allUncertainties.pdf";
803  c1->SetTopMargin(0.15);
804  c1->SetBottomMargin(0.15);
805  c1->SetRightMargin(0.15);
806  c1->SaveAs((outputPDF + "[").c_str());
807 
808  string output1D = outputFolder + "/1DXSecs.pdf";
809 
810  string doubleDiffTitle = "#frac{d^{2}#sigma}{dcos(#theta)dTmu} (cm^{2} / GeV / nucleon)";
811  string uncertTitle = "#frac{#delta#sigma}{#sigma}";
812 
813  LoadNominal("cv");
814  for(pair<string, string> nonNominal : otherNomSysts)
815  {
816  string otherNom = nonNominal.second;
817  if (otherNomSysts.find(otherNom) == otherNomSysts.end())
818  LoadNominal(otherNom); // Don't repeat for the same systematic
819  }
820  if (nominalXSecs2D["cv"]){
821  nominalXSecs2D["cv"]->Draw("COLZ");
822  c1->SaveAs(outputPDF.c_str());
823  }
824  LoadStatUncert("cv");
825 
826  // Standard XSecs
827  for (string syst : standardSysts)
828  {
829  LoadSyst(syst);
830  LoadUncerts(syst);
831  }
832 
833  // Standard XSecs
834  for (pair<string, pair<string, string> > systDef : standardSystsUpDown)
835  {
836  string systBase = systDef.first;
837  string systUp = systDef.second.first;
838  string systDw = systDef.second.second;
839 
840  LoadSyst(systUp);
841  LoadUncerts(systUp);
842 
843  LoadSyst(systDw);
844  LoadUncerts(systDw);
845  }
846 
847  // Beam XSecs
848  for (string beamSyst : beamSysts)
849  {
850  string systUp = beamSyst + "_Up";
851  string systDw = beamSyst + "_Dw";
852 
853  LoadSyst(systUp);
854  LoadUncerts(systUp);
855 
856  LoadSyst(systDw);
857  LoadUncerts(systDw);
858  }
859 
860  // Handle Genie and PPFX
861  for (unsigned int igenie = 0; igenie < N_GENIE; ++igenie)
862  {
863  string syst = "genie" + to_string(igenie);
864  LoadSyst(syst);
865  LoadUncerts(syst);
866  }
867  for (unsigned int ippfx = 0; ippfx < N_PPFX; ++ippfx)
868  {
869  string syst = "ppfx" + to_string(ippfx);
870  LoadSyst(syst);
871  LoadUncerts(syst);
872  }
873 
874  c1->SaveAs((outputPDF + "]").c_str());
875 
877  {"muone", "angle", "neutron", "genie", "ppfx", "beam", "stat"},
878  outputFolder + "/systematicPlots.pdf");
879  c1->Clear();
880 
881  nominalXSecsENu["cv"]->GetYaxis()->SetRangeUser(0.0, nominalXSecsENu["cv"]->GetMaximum() * 1.5);
882  nominalXSecsENu["cv"]->SetLineColor(kRed);
883  nominalXSecsENu["cv"]->Draw("hist");
884  for (unsigned int igenie = 0; igenie < N_GENIE; ++igenie){
885  TH1 * hist = xSecsENu["genie" + to_string(igenie)];
886  hist->SetLineColor(kBlue);
887  hist->Draw("hist same");
888  }
889  nominalXSecsENu["cv"]->Draw("hist same");
890  c1->SaveAs((outputFolder + "/enuGenie.png").c_str());
891 
892  nominalXSecsQ2["cv"]->GetYaxis()->SetRangeUser(0.0, nominalXSecsQ2["cv"]->GetMaximum() * 1.5);
893  nominalXSecsQ2["cv"]->SetLineColor(kRed);
894  nominalXSecsQ2["cv"]->Draw("hist");
895  for (unsigned int ippfx = 0; ippfx < N_PPFX; ++ippfx){
896  TH1 * hist = xSecsQ2["ppfx" + to_string(ippfx)];
897  hist->SetLineColor(kBlue);
898  hist->Draw("hist same");
899  }
900  nominalXSecsQ2["cv"]->Draw("hist same");
901  c1->SaveAs((outputFolder + "/q2Genie.png").c_str());
902 
903  outFile->Close();
904 }
TFile * inputFile
Definition: PlotXSec.C:134
enum BeamMode kRed
map< string, TH1D * > nominalXSecsQ2
Definition: PlotXSec.C:122
void MakeSystematicPlots(vector< string > systs, string systFilename)
Definition: PlotXSec.C:579
vector< string > standardSysts
Definition: PlotXSec.C:41
map< string, TH1D * > xSecsENu
Definition: PlotXSec.C:127
string filename
Definition: shutoffs.py:106
void LoadSyst(string syst)
Definition: PlotXSec.C:371
void LoadNominal(string syst)
Definition: PlotXSec.C:348
TCanvas * c1
Definition: PlotXSec.C:138
TFile * outFile
Definition: PlotXSec.C:135
const unsigned int N_PPFX
Definition: PlotXSec.C:149
vector< string > beamSysts
Definition: PlotXSec.C:65
const unsigned int N_GENIE
Definition: PlotXSec.C:148
map< string, TH2D * > nominalXSecs2D
Definition: PlotXSec.C:120
const map< string, string > otherNomSysts
Definition: PlotXSec.C:57
void LoadUncerts(string syst)
Definition: PlotXSec.C:385
map< string, TH1D * > xSecsQ2
Definition: PlotXSec.C:128
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void LoadStatUncert(string nom)
Definition: PlotXSec.C:398
enum BeamMode kBlue
map< string, pair< string, string > > standardSystsUpDown
Definition: PlotXSec.C:48
map< string, TH1D * > nominalXSecsENu
Definition: PlotXSec.C:121
const TAxis q2Axis ( ,
0.  0,
2.  8 
)

Referenced by FixAxisRangeQ2().

const TAxis uncertAxis ( ,
-0.  2,
0.  2 
)
TH2D* XSec2D ( string  syst)

Definition at line 221 of file PlotXSec.C.

References FixAxisRange2D(), Load3DXSec(), systTitles, and xSecs2D.

Referenced by LoadNominal(), and LoadSyst().

222 {
223  if (xSecs2D.find(syst) != xSecs2D.end())
224  return xSecs2D.at(syst);
225 
226  TH3D * syst3DXSec = Load3DXSec(syst);
227  TH2D * syst2DXSec;
228  if (!syst3DXSec)
229  {
230  xSecs2D[syst] = NULL;
231  return NULL;
232  }
233  syst2DXSec = (TH2D*) syst3DXSec->Project3D("yx");
234  syst2DXSec->Scale(1, "width"); // Area-normalize the xsec
235  FixAxisRange2D(syst2DXSec);
236  xSecs2D[syst] = syst2DXSec;
237  syst2DXSec->SetName(("xsec_" + syst + "_2DProj").c_str());
238  syst2DXSec->SetTitle((systTitles[syst] + " XSec").c_str());
239  syst2DXSec->Write();
240  return syst2DXSec;
241 }
TH3D * Load3DXSec(string syst)
Definition: PlotXSec.C:206
map< string, string > systTitles
Definition: PlotXSec.C:23
void FixAxisRange2D(TH2 *hist)
Definition: PlotXSec.C:168
map< string, TH2D * > xSecs2D
Definition: PlotXSec.C:126
TH2D* XSec2DUncert ( string  syst)

Definition at line 272 of file PlotXSec.C.

References ana::assert(), om::cout, allTimeWatchdog::endl, FixAxisRange2D(), nominalXSecs2D, otherNomSysts, uncerts2D, and xSecs2D.

Referenced by LoadUncerts().

273 {
274  if (uncerts2D.find(syst) != uncerts2D.end())
275  return uncerts2D.at(syst);
276 
277  assert(xSecs2D.find(syst) != xSecs2D.end());
278  TH2D * systXSec = xSecs2D.at(syst);
279  if (!systXSec){
280  uncerts2D[syst] = NULL;
281  return NULL;
282  }
283 
284  TH2D * thisNom = nominalXSecs2D["cv"];
285  if (otherNomSysts.find(syst) != otherNomSysts.end()){
286  cout << "Using non-standard nominal : " << otherNomSysts.at(syst) << " for syst " << syst << endl;
287  thisNom = nominalXSecs2D[otherNomSysts.at(syst)];
288  }
289  TH2D * syst2DUncert = new TH2D(*systXSec - *thisNom);
290  syst2DUncert->SetName((syst + "_2d_uncert").c_str());
291  syst2DUncert->SetTitle((syst + " Uncert").c_str());
292  syst2DUncert->Divide(thisNom);
293  FixAxisRange2D(syst2DUncert);
294  // syst2DUncert->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
295  uncerts2D[syst] = syst2DUncert;
296  syst2DUncert->Write();
297  return syst2DUncert;
298 }
map< string, TH2D * > uncerts2D
Definition: PlotXSec.C:129
OStream cout
Definition: OStream.cxx:6
map< string, TH2D * > nominalXSecs2D
Definition: PlotXSec.C:120
const map< string, string > otherNomSysts
Definition: PlotXSec.C:57
void FixAxisRange2D(TH2 *hist)
Definition: PlotXSec.C:168
map< string, TH2D * > xSecs2D
Definition: PlotXSec.C:126
assert(nhit_max >=nhit_nbins)
TH1D* XSecENu ( string  syst)

Definition at line 244 of file PlotXSec.C.

References FixAxisRangeENu(), inputFile, systTitles, and xSecsENu.

Referenced by LoadNominal(), and LoadSyst().

245 {
246  if (xSecsENu.find(syst) != xSecsENu.end())
247  return xSecsENu.at(syst);
248 
249  TH1D * xSecENu = (TH1D*) inputFile->Get((syst + "/ENu/xsec_" + syst + "_ENu").c_str());
250  FixAxisRangeENu(xSecENu);
251  xSecsENu[syst] = xSecENu;
252  xSecENu->SetTitle((systTitles[syst] + " XSec").c_str());
253  xSecENu->Write();
254  return xSecENu;
255 }
TFile * inputFile
Definition: PlotXSec.C:134
map< string, TH1D * > xSecsENu
Definition: PlotXSec.C:127
map< string, string > systTitles
Definition: PlotXSec.C:23
void FixAxisRangeENu(TH1 *hist)
Definition: PlotXSec.C:159
TH1D* XSecENuUncert ( string  syst)

Definition at line 300 of file PlotXSec.C.

References ana::assert(), om::cout, allTimeWatchdog::endl, FixAxisRangeENu(), nominalXSecsENu, otherNomSysts, uncertsENu, and xSecsENu.

Referenced by LoadUncerts().

301 {
302  if (uncertsENu.find(syst) != uncertsENu.end())
303  return uncertsENu.at(syst);
304 
305  assert(xSecsENu.find(syst) != xSecsENu.end());
306  TH1D * systXSec = xSecsENu.at(syst);
307 
308  TH1D * thisNom = nominalXSecsENu["cv"];
309  if (otherNomSysts.find(syst) != otherNomSysts.end()){
310  cout << "Using non-standard nominal : " << otherNomSysts.at(syst) << " for syst " << syst << endl;
311  thisNom = nominalXSecsENu[otherNomSysts.at(syst)];
312  }
313  TH1D * systUncert = new TH1D(*systXSec - *thisNom);
314  systUncert->SetName((syst + "_enu_uncert").c_str());
315  systUncert->SetTitle((syst + " Uncert").c_str());
316  systUncert->Divide(thisNom);
317  FixAxisRangeENu(systUncert);
318  // systUncert->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
319  uncertsENu[syst] = systUncert;
320  systUncert->Write();
321  return systUncert;
322 }
map< string, TH1D * > uncertsENu
Definition: PlotXSec.C:130
map< string, TH1D * > xSecsENu
Definition: PlotXSec.C:127
OStream cout
Definition: OStream.cxx:6
const map< string, string > otherNomSysts
Definition: PlotXSec.C:57
assert(nhit_max >=nhit_nbins)
void FixAxisRangeENu(TH1 *hist)
Definition: PlotXSec.C:159
map< string, TH1D * > nominalXSecsENu
Definition: PlotXSec.C:121
TH1D* XSecQ2 ( string  syst)

Definition at line 258 of file PlotXSec.C.

References FixAxisRangeQ2(), inputFile, systTitles, and xSecsQ2.

Referenced by LoadNominal(), and LoadSyst().

259 {
260  if (xSecsQ2.find(syst) != xSecsQ2.end())
261  return xSecsQ2.at(syst);
262 
263  TH1D * xSecQ2 = (TH1D*) inputFile->Get((syst + "/Q2/xsec_" + syst + "_Q2").c_str());
264  xSecQ2->Scale(1, "width");
265  FixAxisRangeQ2(xSecQ2);
266  xSecsQ2[syst] = xSecQ2;
267  xSecQ2->SetTitle((systTitles[syst] + " XSec").c_str());
268  xSecQ2->Write();
269  return xSecQ2;
270 }
TFile * inputFile
Definition: PlotXSec.C:134
map< string, string > systTitles
Definition: PlotXSec.C:23
map< string, TH1D * > xSecsQ2
Definition: PlotXSec.C:128
void FixAxisRangeQ2(TH1 *hist)
Definition: PlotXSec.C:151
TH1D* XSecQ2Uncert ( string  syst)

Definition at line 324 of file PlotXSec.C.

References ana::assert(), om::cout, allTimeWatchdog::endl, FixAxisRangeQ2(), nominalXSecsQ2, otherNomSysts, uncertsQ2, and xSecsQ2.

Referenced by LoadUncerts().

325 {
326  if (uncertsQ2.find(syst) != uncertsQ2.end())
327  return uncertsQ2.at(syst);
328 
329  assert(xSecsQ2.find(syst) != xSecsQ2.end());
330  TH1D * systXSec = xSecsQ2.at(syst);
331 
332  TH1D * thisNom = nominalXSecsQ2["cv"];
333  if (otherNomSysts.find(syst) != otherNomSysts.end()){
334  cout << "Using non-standard nominal : " << otherNomSysts.at(syst) << " for syst " << syst << endl;
335  thisNom = nominalXSecsQ2[otherNomSysts.at(syst)];
336  }
337  TH1D * systUncert = new TH1D(*systXSec - *thisNom);
338  systUncert->SetName((syst + "_q2_uncert").c_str());
339  systUncert->SetTitle((syst + " Uncert").c_str());
340  systUncert->Divide(thisNom);
341  FixAxisRangeQ2(systUncert);
342  // systUncert->GetZaxis()->SetRangeUser(uncertAxis.GetXmin(), uncertAxis.GetXmax());
343  uncertsQ2[syst] = systUncert;
344  systUncert->Write();
345  return systUncert;
346 }
map< string, TH1D * > nominalXSecsQ2
Definition: PlotXSec.C:122
map< string, TH1D * > uncertsQ2
Definition: PlotXSec.C:131
OStream cout
Definition: OStream.cxx:6
const map< string, string > otherNomSysts
Definition: PlotXSec.C:57
assert(nhit_max >=nhit_nbins)
map< string, TH1D * > xSecsQ2
Definition: PlotXSec.C:128
void FixAxisRangeQ2(TH1 *hist)
Definition: PlotXSec.C:151

Variable Documentation

map<string, TH3D*> all3DXSecs

Definition at line 125 of file PlotXSec.C.

Referenced by Load3DXSec().

vector<string> beamSysts
Initial value:
{
"2kA",
"02mmBeamSpotSize",
"1mmBeamShiftX",
"1mmBeamShiftY",
"3mmHorn1X",
"3mmHorn1Y",
"3mmHorn2X",
"3mmHorn2Y",
"7mmTargetZ",
"MagneticFieldinDecayPipe",
"1mmHornWater"
}

Definition at line 65 of file PlotXSec.C.

Referenced by GetTotalBeamUncert2D(), GetTotalBeamUncertENu(), GetTotalBeamUncertQ2(), and PlotXSec().

TCanvas* c1

Definition at line 138 of file PlotXSec.C.

Referenced by PlotXSec().

string defaultNominal = "cv"

Definition at line 118 of file PlotXSec.C.

map<int, double> maxTMuBins
Initial value:
{
{2, 1.1},
{3, 1.1},
{4, 1.1},
{5, 1.4},
{6, 1.4},
{7, 1.4},
{8, 1.8},
{9, 1.9},
{10, 2.2},
}

Definition at line 94 of file PlotXSec.C.

Referenced by LimitTMuBins().

const unsigned int N_GENIE = 1000

Definition at line 148 of file PlotXSec.C.

Referenced by MakeSystematicPlots(), and PlotXSec().

const unsigned int N_PPFX = 100

Definition at line 149 of file PlotXSec.C.

Referenced by MakeSystematicPlots(), and PlotXSec().

TH2D* nom2DXSec

Definition at line 115 of file PlotXSec.C.

TH3D* nom3DXSec

Definition at line 114 of file PlotXSec.C.

map<string, TH2D*> nominalXSecs2D
map<string, TH3D*> nominalXSecs3D

Definition at line 119 of file PlotXSec.C.

Referenced by LoadNominal().

map<string, TH1D*> nominalXSecsENu
map<string, TH1D*> nominalXSecsQ2
const map<string, string> otherNomSysts
Initial value:
{
{"neutron", "neutron_nom"},
{"light", "nom_good_seventh"},
{"calib", "nom_good_seventh"},
{"calibshift", "nom_good_seventh"}
}

Definition at line 57 of file PlotXSec.C.

Referenced by PlotXSec(), XSec2DUncert(), XSecENuUncert(), and XSecQ2Uncert().

TFile* outFile

Definition at line 135 of file PlotXSec.C.

Referenced by accum_pot_equivalent_fhc(), NDPredictionHandler::AddVar(), FDPredictionHandler::AddVar(), Analyse_Data2DataComp(), Analyse_Data2DataComp_HigherEnergyCuts(), Analyse_Data2DataComp_kNumu2020ND(), AnalyzeNus18Pred(), AnalyzeNus18Systs(), AnalyzeSysts(), ApplyOscillations(), CalcCutVals(), calculateWrongSignNue(), calculateWrongSignNumuQ1(), calculateWrongSignNumuQ2(), calculateWrongSignNumuQ3(), calculateWrongSignNumuQ4(), CalculateXSec(), combineFiles(), Compare_NoExtrap(), Compare_Spectra(), ComparisonPlots_Data(), ComparisonPlots_MC(), CosmicPred(), demo2p5a(), DoTraining(), DrawExtrapSurface(), Evaluate_BDTMLP_Algs_PredNoExtrap(), Evaluate_BDTMLP_Algs_Spectra(), eventListToTextFile(), example_macro(), FillSpectra(), FillTrainingTrees(), FOMPlot(), GenerateCovMx(), get_numi_data_histogram(), LoadNominal(), LoadSyst(), LoadUncerts(), make_dst_cosrejbdttrain(), make_nue_ana2018_pot(), make_nue_ana2019_epoch7d_pot(), make_nue_ana2019_epoch8b_pot(), make_nue_ana2019_fhc_ub_pot(), make_nue_ana2019_rhc_ub_pot(), make_nue_thirdana_pot(), make_numu_thirdana_pot(), make_nus17_pot(), make_rhcpred_2017(), make_rockpred(), make_rockpred_2017(), MakeCosmicSpectra(), MakeCovMx(), MakeExtrapSurface(), MakeISysts(), MakeNumuCovMx(), MakeNus18ExtrapPred(), MakeNus18SidebandPred(), MakeNus18Systs(), makePIDCutTuning(), makePrediction(), MakePrediction(), MakePredictionNoOsc_FHC_FD(), MakePredictionNoOsc_FHC_ND(), MakePredictionNoOsc_RHC_FD(), MakePredictionNoOsc_RHC_ND(), MakeSurface(), MakeSurfaceBinningStudy(), MakeSurfaceLLTest(), mcTruthPredictions(), meanWeight_macro(), MergePredictions(), MergeSurface(), MergeSurfaces(), MichelDecompTest(), MRDiFStudy_FHC_Step2(), MRDiFStudy_RHC_Step2(), multiverse_efficiency_macro(), NoExtrap(), nuCrystalBall(), nue_michelDataMC(), NueExtrap(), NuMu2020_TrimCAFs(), NumuCosmic(), NumuExtrap(), NuSCalculateCorr(), NuSTreeMaker(), PlotCovariance(), PlotNus18Sideband(), PlotXSec(), RunCalibration(), runTwoSampleDecomp(), SaveSpectra(), SystematicComp(), test_micheldecomp(), ThrowFakeData(), TrimCAFs(), UnfoldXSec(), FDPredictionHandler::~FDPredictionHandler(), and NDPredictionHandler::~NDPredictionHandler().

vector<string> standardSysts
Initial value:
{
"calibshift",
"ckv"
}

Definition at line 41 of file PlotXSec.C.

Referenced by MakeSystematicPlots(), and PlotXSec().

map<string, pair<string, string> > standardSystsUpDown
Initial value:
{
{"light", {"lightup", "lightdown"}},
{"calib", {"calibpos", "calibneg"}},
{"muone", {"muoneup", "muonedw"}},
{"angle", {"angleup", "angledw"}},
{"neutron", {"neutron_Up", "neutron_Dw"}}
}

Definition at line 48 of file PlotXSec.C.

Referenced by MakeSystematicPlots(), and PlotXSec().

map<string, Color_t> systColors
Initial value:
{
{"total", kBlack},
{"beam", kBlue},
{"ckv", kOrange},
{"neutron", kGreen},
{"angle", kRed},
{"muone", kViolet},
{"stat", kGray},
{"genie", kCyan},
{"ppfx", kOrange}
}
enum BeamMode kOrange
enum BeamMode kRed
enum BeamMode kViolet
enum BeamMode kGreen
enum BeamMode kBlue

Definition at line 80 of file PlotXSec.C.

Referenced by MakeSystematicPlots().

map<string, string> systTitles
Initial value:
{
{"cv", "Central Value"},
{"light", "Light Level"},
{"calib", "Calibration"},
{"calibshift", "Calib Shape"},
{"ckv", "Cherenkov"},
{"muone", "Muon Energy"},
{"neutron", "Neutron Syst"},
{"angle", "Muon Angle"},
{"neutron_nom","Nominal (Neutron syst)"},
{"nom_good_seventh", "Nominal (calib / light levels)"},
{"genie", "XSec and FSI"},
{"ppfx", "Beam - HP"},
{"stat", "Statistical"},
{"beam", "Beam - Focusing"}
}

Definition at line 23 of file PlotXSec.C.

Referenced by Load3DXSec(), MakeSystematicPlots(), XSec2D(), XSecENu(), and XSecQ2().

string uncertLabel = "_uncert"

Definition at line 107 of file PlotXSec.C.

map<string, TH2D*> uncerts2D
map<string, TH1D*> uncertsENu
map<string, TH1D*> uncertsQ2
const string var3DName = "EAvail"

Definition at line 106 of file PlotXSec.C.

Referenced by Load3DXSec().

const vector<string> vars3D
Initial value:
{
"EAvail"
}

Definition at line 109 of file PlotXSec.C.

map<string, TH2D*> xSecs2D

Definition at line 126 of file PlotXSec.C.

Referenced by MakeSystematicPlots(), XSec2D(), and XSec2DUncert().

map<string, TH1D*> xSecsENu

Definition at line 127 of file PlotXSec.C.

Referenced by MakeSystematicPlots(), PlotXSec(), XSecENu(), and XSecENuUncert().

map<string, TH1D*> xSecsQ2

Definition at line 128 of file PlotXSec.C.

Referenced by MakeSystematicPlots(), PlotXSec(), XSecQ2(), and XSecQ2Uncert().