PlotNus18Sideband.C

Go to the documentation of this file.

#ifdef __CINT__
void PlotNus18Sideband(std::string inFile)
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "OscLib/OscCalcSterile.h"

#include "NuXAna/macros/Nus18/Nus18Utilities.h"

#include "TH1.h"
#include "THStack.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TLegend.h"

#include <utility>

using namespace ana;

void PlotNus18Sideband(std::string inFile) {

  TH1::AddDirectory(0);

  // Read in the data objects
  // IPrediction* pred_genie_e = LoadFromFile<PredictionNoExtrap>(inFile.c_str(), "pred_genie").release();
  // Spectrum* spec_data_e     = LoadFromFile<Spectrum>(inFile.c_str(), "spec_data").release();
  // Spectrum* spec_cos_e      = LoadFromFile<Spectrum>(inFile.c_str(), "spec_cos").release();

  IPrediction* pred_genie_e = LoadFromFile<PredictionNoExtrap>(inFile.c_str(), "pred_genie_e").release();
  Spectrum* spec_data_e     = LoadFromFile<Spectrum>(inFile.c_str(), "spec_data_e").release();
  Spectrum* spec_cos_e      = LoadFromFile<Spectrum>(inFile.c_str(), "spec_cos_e").release();

  IPrediction* pred_genie_cvn = LoadFromFile<PredictionNoExtrap>(inFile.c_str(), "pred_genie_cvn").release();
  Spectrum* spec_data_cvn     = LoadFromFile<Spectrum>(inFile.c_str(), "spec_data_cvn").release();
  Spectrum* spec_cos_cvn      = LoadFromFile<Spectrum>(inFile.c_str(), "spec_cos_cvn").release();

  IPrediction* pred_genie_nhits = LoadFromFile<PredictionNoExtrap>(inFile.c_str(), "pred_genie_nhits").release();
  Spectrum* spec_data_nhits     = LoadFromFile<Spectrum>(inFile.c_str(), "spec_data_nhits").release();
  Spectrum* spec_cos_nhits      = LoadFromFile<Spectrum>(inFile.c_str(), "spec_cos_nhits").release();

  // Set up output
  std::string fileName = inFile.substr(0, inFile.find(".root"));
  std::string saveName = fileName + "_ana.root";
  TFile* outFile = new TFile(saveName.c_str(), "RECREATE");
  TCanvas* canv = new TCanvas("canv","",800,600);
  
  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc = DefaultSterileCalc(3);
  SetAngles(calc);

  TH1* hDataE   = spec_data_e->ToTH1(spec_data_e->POT());
  TH1* hSigE    = pred_genie_e->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(spec_data_e->POT());
  TH1* hBeamBgE = pred_genie_e->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(spec_data_e->POT());
  TH1* hCosE    = spec_cos_e->ToTH1(spec_data_e->Livetime(), kLivetime);

  TH1* hDataCVN   = spec_data_cvn->ToTH1(spec_data_cvn->POT());
  TH1* hSigCVN    = pred_genie_cvn->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(spec_data_cvn->POT());
  TH1* hBeamBgCVN = pred_genie_cvn->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(spec_data_cvn->POT());
  TH1* hCosCVN    = spec_cos_cvn->ToTH1(spec_data_cvn->Livetime(), kLivetime);

  TH1* hDataNHits   = spec_data_nhits->ToTH1(spec_data_nhits->POT());
  TH1* hSigNHits    = pred_genie_nhits->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(spec_data_nhits->POT());
  TH1* hBeamBgNHits = pred_genie_nhits->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(spec_data_nhits->POT());
  TH1* hCosNHits    = spec_cos_nhits->ToTH1(spec_data_nhits->Livetime(), kLivetime);

  // Plot everything!
  TLegend* leg = new TLegend(0.5,0.6,0.88,0.88);

  // energy
  std::cout << "Energy:" << std::endl;
  std::cout << "  Data: " << hDataE->Integral() << std::endl;
  std::cout << "  NC: " << hSigE->Integral() << std::endl;
  std::cout << "  Beam background: " << hBeamBgE->Integral() << std::endl;
  std::cout << "  Cosmics: " << hCosE->Integral() << std::endl;

  hDataE->SetMarkerStyle(8);
  hDataE->SetMarkerSize(1);
  hDataE->SetMarkerColor(1);
  hDataE->SetLineColor(1);
  hSigE->SetLineColor(kAzure+2);
  hBeamBgE->SetLineColor(kGreen+2);
  hCosE->SetLineColor(kOrange+2);

  leg->Clear();
  leg->AddEntry(hDataE, "Data", "p");
  leg->AddEntry(hSigE, "Signal", "l");
  leg->AddEntry(hBeamBgE, "Beam background", "l");
  leg->AddEntry(hCosE, "Cosmic background", "l");

  THStack* stackE = new THStack("stackE", "");
  stackE->Add(hCosE);
  stackE->Add(hBeamBgE);
  stackE->Add(hSigE);

  canv->cd();
  canv->Clear();
  hDataE->Draw();
  stackE->Draw("hist nostack same");
  leg->Draw();
  outFile->cd();
  canv->Write("SidebandDataMCEnergy");

  canv->cd();
  canv->Clear();
  hSigE->SetFillColor(kAzure+2);
  hBeamBgE->SetFillColor(kGreen+2);
  hCosE->SetFillColor(kOrange+2);
  stackE->SetMaximum(200);
  stackE->Draw("hist");
  stackE->GetXaxis()->SetTitle("Energy Deposited (GeV)");
  hDataE->Draw("same");
  leg->Draw();
  outFile->cd();
  canv->Write("SidebandDataMCEnergyStack");

  // cvn
  hDataCVN->SetMarkerStyle(8);
  hDataCVN->SetMarkerSize(1);
  hDataCVN->SetMarkerColor(1);
  hDataCVN->SetLineColor(1);
  hSigCVN->SetLineColor(kAzure+2);
  hBeamBgCVN->SetLineColor(kGreen+2);
  hCosCVN->SetLineColor(kOrange+2);

  leg->Clear();
  leg->AddEntry(hDataCVN, "Data", "p");
  leg->AddEntry(hSigCVN, "Signal", "l");
  leg->AddEntry(hBeamBgCVN, "Beam background", "l");
  leg->AddEntry(hCosCVN, "Cosmic background", "l");

  THStack* stackCVN = new THStack("stackCVN", "");
  stackCVN->Add(hCosCVN);
  stackCVN->Add(hBeamBgCVN);
  stackCVN->Add(hSigCVN);

  canv->cd();
  canv->Clear();
  hDataCVN->Draw();
  stackCVN->Draw("hist nostack same");
  leg->Draw();
  outFile->cd();
  canv->Write("SidebandDataMCCVN");

  canv->cd();
  canv->Clear();
  hSigCVN->SetFillColor(kAzure+2);
  hBeamBgCVN->SetFillColor(kGreen+2);
  hCosCVN->SetFillColor(kOrange+2);
  stackCVN->SetMaximum(60);
  stackCVN->Draw("hist");
  stackCVN->GetXaxis()->SetTitle("CVNnc Identifier");
  hDataCVN->Draw("same");
  leg->Draw();
  outFile->cd();
  canv->Write("SidebandDataMCCVNncStack");

  // nhits
  hDataNHits->SetMarkerStyle(8);
  hDataNHits->SetMarkerSize(1);
  hDataNHits->SetMarkerColor(1);
  hDataNHits->SetLineColor(1);
  hSigNHits->SetLineColor(kAzure+2);
  hBeamBgNHits->SetLineColor(kGreen+2);
  hCosNHits->SetLineColor(kOrange+2);

  leg->Clear();
  leg->AddEntry(hDataNHits, "Data", "p");
  leg->AddEntry(hSigNHits, "Signal", "l");
  leg->AddEntry(hBeamBgNHits, "Beam background", "l");
  leg->AddEntry(hCosNHits, "Cosmic background", "l");

  THStack* stackNHits = new THStack("stackNHits", "");
  stackNHits->Add(hCosNHits);
  stackNHits->Add(hBeamBgNHits);
  stackNHits->Add(hSigNHits);

  canv->cd();
  canv->Clear();
  hDataNHits->Draw();
  stackNHits->Draw("hist nostack same");
  leg->Draw();
  outFile->cd();
  canv->Write("SidebandDataMCNHits");

  canv->cd();
  canv->Clear();
  hSigNHits->SetFillColor(kAzure+2);
  hBeamBgNHits->SetFillColor(kGreen+2);
  hCosNHits->SetFillColor(kOrange+2);
  stackNHits->SetMaximum(60);
  stackNHits->Draw("hist");
  stackNHits->GetXaxis()->SetTitle("Number of Hits");
  hDataNHits->Draw("same");
  leg->Draw();
  outFile->cd();
  canv->Write("SidebandDataMCNHitsStack");

  outFile->Close();
  delete canv;
  delete outFile;

  return;

}

#endif