PlotUnfolding.C

Go to the documentation of this file.

// PlotUnfolding.C

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TStyle.h"

#include <string>
#include <vector>

unsigned int nuniv = 0;
unsigned int nppfx = 0;

const std::string xs_mode0 = "xs_nominal";
const std::string xs_mode1 = "xs_bkgcons";

const std::string xs_mode = xs_mode0;

const std::string yxaxistitle = ";Cos(#theta_{#mu});Muon Kinetic Energy (GeV)";
const std::string zaxistitle = "Available Energy";

const std::vector<std::string> vars3D{"EAvail"};
const std::vector<std::string> vars1D
{
  "ENu",
  "Q2"
};

std::map<std::string, std::string> systNames = {
  // {"lightup",         "syst_specs_lightup.root"},
  // {"lightdown",       "syst_specs_lightdw.root"},
  // {"calibpos",        "syst_specs_calibpos.root"},
  // {"calibneg",        "syst_specs_calibneg.root"},
  // {"MuESDw", "Muon Energy Down"},
  // {"MuESUp", "Muon Energy Up"},
  {"angleup",         "Angle Up"},
  {"angledw",         "Angle Down"},
  {"ckv",             "Cherenkov"},
  {"neutron_Up",      "Neutron Up"},
  {"neutron_Dw",      "Neutron Down"}
};

std::map<std::string, std::string> systUpDownPairs = {
  {"neutron", "Neutron"}
};

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

TAxis* xaxis = new TAxis(1, 0.5, 1.0);
TAxis* yaxis = new TAxis(1, 0.5, 2.5);
TAxis* zaxis = new TAxis(1, -0.25, 0.25);
TAxis* eaxis = new TAxis(1, 0.0, 2.5);

std::vector<std::string> updownlabel = {"_Up", "_Dw"};

TH2F * makeSystematic2D(TH2F* syst, TH2F* nominal){
  TH2F * rat = new TH2F(*syst - *nominal);
  rat->Divide(nominal);
  return rat;
}

TH1F * makeSystematic1D(TH1F* syst, TH1F* nominal){
  TH1F * rat = new TH1F(*syst - *nominal);
  rat->Divide(nominal);
  return rat;
}

void PlotUnfolding(std::string inputFile, std::string outputDir){
  // Load the inputs
  std::string outputFile = outputDir + "/UnfoldingResults.root";
  std::string outputPDF  = outputDir + "/all_uncertainties.pdf";
  std::string outputFormat = ".png";
  TFile * fin            = new TFile(inputFile.c_str());
  TFile * fout           = new TFile(outputFile.c_str(), "RECREATE");
  TDirectory * systDir = fout->mkdir("systematics");
  fin->cd();

  TCanvas * c1 = new TCanvas("c1", "c1", 1600, 1200);
  c1->SaveAs((outputPDF + "[").c_str());
  fout->cd();

  gStyle->SetPaintTextFormat(".3f");

  for (std::string var : vars3D){

    // Load nominal
    TH3F * nom3D            = (TH3F*) fin->Get((xs_mode + "/hXS3Dsec_neutron_nom_" + var).c_str());
    TH2F * nom2D            = (TH2F*) fin->Get((xs_mode + "/hXS2Dsec_neutron_nom_" + var).c_str());
    TH1F * nom1D            = (TH1F*) fin->Get((xs_mode + "/hXS1Dsec_neutron_nom_" + var).c_str());
    TH2F * nomEff           = (TH2F*) fin->Get(("h2DEff_nom_" + var).c_str());

    for (std::pair<std::string, std::string> syst : systNames){
      std::string systName = syst.first;
      std::string systTitle = syst.second;
      TH2F * systHist2D = (TH2F*) fin->Get((xs_mode + "/hXS2Dsec_" + systName + "_" + var).c_str());
      TH2F * systUncert2D = makeSystematic2D(systHist2D, nom2D);
      systUncert2D->SetName((systName + "_fracuncert").c_str());
      systUncert2D->SetTitle(("(" + systTitle + " - nominal) / nominal" + yxaxistitle).c_str());
      systUncert2D->GetXaxis()->SetRangeUser(xaxis->GetXmin(), xaxis->GetXmax());
      systUncert2D->GetYaxis()->SetRangeUser(yaxis->GetXmin(), yaxis->GetXmax());
      systUncert2D->GetZaxis()->SetRangeUser(zaxis->GetXmin(), zaxis->GetXmax());
      systUncert2D->Draw("COLZ");
      systUncert2D->Write();
      c1->SaveAs((outputDir + "/" + systName + "_fracuncert" + outputFormat).c_str());
      c1->SaveAs(outputPDF.c_str());

      TH1F * systHist1D = (TH1F*) fin->Get((xs_mode + "/hXS1Dsec_" + systName + "_" + var).c_str());
      TH1F * systUncert1D = makeSystematic1D(systHist1D, nom1D);
      systUncert1D->SetName((systName + "_fracuncert_1d").c_str());
      systUncert1D->SetTitle(("(" + systTitle + " - nominal) / nominal" + zaxistitle).c_str());
      systUncert1D->GetXaxis()->SetRangeUser(eaxis->GetXmin(), eaxis->GetXmax());
      systUncert1D->Draw("hist");
      systUncert1D->Write();
      c1->SaveAs((outputDir + "/" + systName + "_fracuncert_1d" + outputFormat).c_str());
      c1->SaveAs(outputPDF.c_str());      
    }

    for (std::string syst : beamSysts){
      for (unsigned int iupdown = 0; iupdown < updownlabel.size(); iupdown++){
        std::string systName = syst + updownlabel[iupdown];
        std::string systTitle = syst + updownlabel[iupdown];
        TH2F * systHist2D = (TH2F*) fin->Get((xs_mode + "/hXS2Dsec_" + systName + "_" + var).c_str());
        TH2F * systUncert2D = makeSystematic2D(systHist2D, nom2D);
        systUncert2D->SetName((systName + "_fracuncert").c_str());
        systUncert2D->SetTitle(("(" + systTitle + " - nominal) / nominal" + yxaxistitle).c_str());
        systUncert2D->GetXaxis()->SetRangeUser(xaxis->GetXmin(), xaxis->GetXmax());
        systUncert2D->GetYaxis()->SetRangeUser(yaxis->GetXmin(), yaxis->GetXmax());
        systUncert2D->GetZaxis()->SetRangeUser(zaxis->GetXmin(), zaxis->GetXmax());
        systUncert2D->Draw("COLZ");
        systUncert2D->Write();
        c1->SaveAs((outputDir + "/" + systName + "_fracuncert" + outputFormat).c_str());
        c1->SaveAs(outputPDF.c_str());
      }
    }
    nom2D->SetName("nominal_2Dxsec");
    nom2D->SetTitle("Nominal Cross-section;Cos(#theta_{#mu});Muon Kinetic Energy (GeV)");
    nom2D->GetXaxis()->SetRangeUser(xaxis->GetXmin(), xaxis->GetXmax());
    nom2D->GetYaxis()->SetRangeUser(yaxis->GetXmin(), yaxis->GetXmax());
    gPad->SetRightMargin(0.15);
    nom2D->Draw("COLZ");
    nom2D->Write();
    c1->SaveAs((outputDir + "/nominal_xsec" + outputFormat).c_str());

    nom1D->SetName("nominal_1Dxsec");
    nom1D->SetTitle("Nominal Cross-section;Available Energy (GeV)");
    nom1D->GetXaxis()->SetRangeUser(eaxis->GetXmin(), eaxis->GetXmax());
    nom1D->Draw("hist");
    nom1D->Write();
    c1->SaveAs((outputDir + "/nominal_xsec_1d" + outputFormat).c_str());

    nomEff->SetName("nominal_eff");
    nomEff->SetTitle(("Nominal Efficiency" + yxaxistitle).c_str());
    nomEff->GetXaxis()->SetRangeUser(xaxis->GetXmin(), xaxis->GetXmax());
    nomEff->GetYaxis()->SetRangeUser(yaxis->GetXmin(), yaxis->GetXmax());
    nomEff->GetZaxis()->SetRangeUser(0.0, 1.0);
    nomEff->Draw("COLZ");
    nomEff->Write();
    c1->SaveAs((outputDir + "/nominal_eff" + outputFormat).c_str());

    for (std::pair<std::string, std::string> syst : systNames){
      std::string systName = syst.first;
      std::string systTitle = syst.second;
      TH2F * systEff = (TH2F*) fin->Get(("h2DEff_" + systName + "_" + var).c_str());
      systEff->SetName((systName + "_eff").c_str());
      systEff->SetTitle((systTitle + " Efficiency" + yxaxistitle).c_str());
      systEff->GetXaxis()->SetRangeUser(xaxis->GetXmin(), xaxis->GetXmax());
      systEff->GetYaxis()->SetRangeUser(yaxis->GetXmin(), yaxis->GetXmax());
      systEff->GetZaxis()->SetRangeUser(0.0, 1.0);
      systEff->Draw("COLZ");
      systEff->Write();
      c1->SaveAs((outputDir + "/" + systName + "_eff" + outputFormat).c_str());
    }

    for (std::string syst : beamSysts)
      for (unsigned int iupdown = 0; iupdown < updownlabel.size(); iupdown++){
        std::string systName = syst + updownlabel[iupdown];
        std::string systTitle = syst + updownlabel[iupdown];
        TH2F * systEff = (TH2F*) fin->Get(("h2DEff_" + systName + "_" + var).c_str());
        systEff->SetName((systName + "_eff").c_str());
        systEff->SetTitle((systTitle + " Efficiency" + yxaxistitle).c_str());
        systEff->GetXaxis()->SetRangeUser(xaxis->GetXmin(), xaxis->GetXmax());
        systEff->GetYaxis()->SetRangeUser(yaxis->GetXmin(), yaxis->GetXmax());
        systEff->GetZaxis()->SetRangeUser(0.0, 1.0);
        systEff->Draw("COLZ");
        systEff->Write();
        c1->SaveAs((outputDir + "/" + systName + "_eff" + outputFormat).c_str());
      }
  }

  c1->SaveAs((outputPDF + "]").c_str());
  fout->Close();
}