plot_nd_spectra_2018.C File Reference

#include "CAFAna/Analysis/Exposures.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "Utilities/rootlogon.C"
#include "TFile.h"
#include "TH1.h"
#include "THStack.h"
#include "TCanvas.h"
#include "TGaxis.h"
#include "TGraph.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TLine.h"
#include "TSystem.h"
#include "TPaveLabel.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <cstddef>

Go to the source code of this file.

Functions
void	SetAxis (TH1 *h, bool visible, int color, int style=1, bool scale=true)
TString	POTstr (double pot)
double	DerivePowScale (TH1 *h)
TH1 *	MakeRatio (TH1 *num, TH1 *denom, int color=1)
void	One ()
void	Zero ()
void	Legend (TH1 *h_to_clone)
void	RatioLegend ()
void	NeutrinoMode (std::string mode)
void	plot_nd_spectra_2018 (std::string mode="fhc")

Function Documentation

double DerivePowScale

(

TH1 *

h

)

Definition at line 209 of file plot_nd_spectra_2018.C.

References makeTrainCVSamples::int, and std::log10().

Referenced by SetAxis().

                              {
  // "Next smallest power of 1000 than maximum bin" - from plot_validation_datamc.py
  return 3*(int(log10(h->GetMaximum())))/3;
}

void Legend

(

TH1 *

h_to_clone

)

Definition at line 247 of file plot_nd_spectra_2018.C.

References make_syst_table_plots::h, MECModelEnuComparisons::i, makeTrainCVSamples::int, kGreen, kRed, MECModelEnuComparisons::leg, and moon_position_table_new3::second.

Referenced by plot_nd_spectra_2018().

                            {
  std::vector <std::pair<int, TString>> entries = {{kBlack, "ND data"},
                                                   {kRed, "Total MC"},
                                                   {kAzure, "NC"},
                                                   {kGreen+2, "#nu_{#mu} CC"},
                                                   {kMagenta, "#nu_{e} CC"},
                                                   {kGreen-3, "#bar{#nu}_{#mu} CC"},
                                                   {kMagenta-2, "#bar{#nu}_{e} CC"},
                                                   {kBlack, "Uncorr. MC"},
                                                   {17, "Uncorr. MC"},
                                                   {1, "Decomp MC"}};
                                                       

  TLegend* leg = new TLegend(.78, .4, .95, .89);

  leg->SetFillColor(kWhite);
  leg->SetBorderSize(1);
  leg->SetMargin(.22);

  TH1* h = (TH1F*)h_to_clone->Clone();
  leg->SetTextSize(0.55*h->GetXaxis()->GetTitleSize());
  h->SetLineStyle(1);
  h->SetLineColor(14);
  h->SetMarkerStyle(kFullCircle);
  
  leg->AddEntry(h, entries[0].second, "lep");
  for(int i = 1; i < (int) entries.size()-2; i++){
    h->SetLineColor(entries[i].first);
    if(i == 7) h->SetLineStyle(7); // Uncorr MC
    else h->SetLineStyle(1);
    if(i == 8){
      h->SetMarkerColor(entries[i].first); // Ratio to Uncorr MC
      h->SetMarkerStyle(kFullCircle);
      leg->AddEntry(h->Clone(), entries[i].second, "lep");
      continue;
    }
    if(i == 9) {
      h->SetMarkerColor(entries[i].first);
      h->SetMarkerStyle(kFullCircle);
      leg->AddEntry(h->Clone(), entries[i].second, "lep");
      continue;
    }
    leg->AddEntry(h->Clone(), entries[i].second, "l");                   
  }
  leg->SetTextSize(1.3*leg->GetTextSize());
  leg->Draw();

}

TH1 * MakeRatio	(	TH1 *	num,
		TH1 *	denom,
		int	color = 1
	)

Definition at line 214 of file plot_nd_spectra_2018.C.

References runNovaSAM::ret.

Referenced by plot_nd_spectra_2018().

                                                {
  TH1* ret = (TH1*)num->Clone();
  ret->Divide(denom);

  ret->SetMarkerStyle(kFullCircle);
  ret->SetMarkerColor(color);
  ret->SetLineColor(color);
  ret->GetYaxis()->SetRangeUser(.5, 1.5);
  ret->GetYaxis()->SetTitle("Data / MC");
  ret->GetXaxis()->SetNdivisions(1900, kFALSE);
  //  ret->GetXaxis()->SetTickLength(1.);
  return ret;
}

void NeutrinoMode

(

std::string

mode

)

Definition at line 312 of file plot_nd_spectra_2018.C.

References string.

Referenced by plot_nd_spectra_2018().

                                 {
  std::string txt = mode == "fhc"? "Neutrino Mode" : "Antineutrino Mode";
  TLatex* lb = new TLatex(.1, .95, txt.c_str());
  lb->SetNDC();
  lb->SetTextAlign(13);
  lb->SetTextColor(kGray+3);
  lb->SetTextSize(2/45.);
  lb->Draw();

}

void One

(

)

Definition at line 238 of file plot_nd_spectra_2018.C.

References MECModelEnuComparisons::g.

Referenced by osc::_PMNS< T >::EvalEqn11(), and plot_nd_spectra_2018().

          {
  TGraph* g = new TGraph();
  g->SetPoint(0, -1e10,1);
  g->SetPoint(1, 1e10, 1);
  g->SetLineWidth(2);
  g->SetLineStyle(7);
  g->Draw("same");
}

void plot_nd_spectra_2018

(

std::string

mode = "fhc"

)

Definition at line 44 of file plot_nd_spectra_2018.C.

References plot_validation_datamc::c, plot_validation_datamc::Clone(), colors, comps, om::cout, plot_validation_datamc_2018::decomp, allTimeWatchdog::endl, plot_validation_datamc::fdata, plot_validation_datamc_2018::fdecomp, fmc, genie::utils::style::Format(), make_mec_shifts_plots::GetMaximum(), GetYaxis(), MECModelEnuComparisons::i, makeTrainCVSamples::int, calib::j, kGreen, kRed, PandAna.Demos.pi0_spectra::labels, Legend(), ana::Spectrum::LoadFrom(), MakeRatio(), maxY, mc, submit_nova_art::mode, NeutrinoMode(), ana::Nue2018ThreeBinAxis(), ana::Nue2018ThreeBinDivisions(), ana::Nue2018ThreeBinLabels(), One(), plot_validation_datamc::p1, plot_validation_datamc::p2, POTstr(), Preliminary(), RatioLegend(), plot_validation_datamc_2018::rdecomp, SetAxis(), string, art::to_string(), PandAna.Demos.pi0_spectra::tot, and Zero().

                                                  {
  
  std::string fdir = "/nova/ana/nu_e_ana/Ana2018/DataMC/";

  TFile* fdata = new TFile((fdir+"20180331_Ana2018_"+mode+"_datamc_decomp_data_validation.root").c_str(), "READ");
  TFile *fmc = new TFile((fdir+ "20180331_Ana2018_"+mode+"_datamc_decomp_mc_validation.root").c_str(), "READ");
  TFile* fdecomp = new TFile((fdir + "20180331_Ana2018_"+mode+"_datamc_decomp_mc_decomp_validation.root").c_str(), "READ");

  if(fmc->IsZombie()) {std::cout << "MC File not found "; return; }
  if(fdata->IsZombie()) {std::cout << "Data File not found "; return; }
  if(fdecomp->IsZombie()) {std::cout << "MC Decomp File not found"; return; }

  // Pulling Data        
  auto spec = Spectrum::LoadFrom(fdata, (TString)"Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}");
  double datapot = spec->POT();
  TH1* hdata = spec->ToTH1(datapot);
 
  std::vector <TH1*> mc;
  std::vector <TH1*> decomp;
  std::vector <std::string> comps = {"0_nue", "1_numu", "2_nc", "3_antinue", "4_antinumu"};
  TH1* tempmc;
  TH1* tempdecomp;
  unique_ptr<Spectrum> tempspec;
  // Load spectra from uncorrected mc and decomp corrected mc files
  for(std::string comp : comps ) {
    tempspec = Spectrum::LoadFrom(fmc, (TString::Format("Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat=%s}",comp.c_str())));
    tempmc = tempspec->ToTH1(datapot);
    
    std::cout << "Loaded MC Spectrum Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat="+comp+"}" << endl;; 
    
    tempspec = Spectrum::LoadFrom(fdecomp, (TString::Format("Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat=%s}",comp.c_str())));
    tempdecomp = tempspec->ToTH1(datapot);
    
    std::cout << "Loaded Decomp Spectrum Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat="+comp+"}" << endl;

    mc.push_back(tempmc);
    decomp.push_back(tempdecomp);
  }
  
  // Add up all of the components to get total MC
  TH1* totalmc = (TH1*)mc[0]->Clone("totalmc");
  TH1* totaldecomp = (TH1*)decomp[0]->Clone("totaldecomp");
  for(int i = 1; i < (int)mc.size(); i++) {
    totalmc->Add( mc[i] );
    totaldecomp->Add( decomp[i] );
    std::cout << "Adding component " + std::to_string(i) << endl;
  }
  mc.push_back(totalmc);
  decomp.push_back(totaldecomp);
  std::cout << "Added total MC and Decomp" << endl;

  int tot = mc.size()-1;
  TH1* rdecomp;
  TH1* rmc;
  // make the histograms for the ratios
  rdecomp = MakeRatio(hdata, decomp[tot]);
  rmc = MakeRatio(hdata, mc[tot], 17);


  // Lets start drawing
  // Split plot
  TCanvas* c = new TCanvas("c");
  double ysplit=0.35;
  TPad* p1 =new TPad("p1", "", 0,0,1,1);
  TPad* p2 =new TPad("p2","",0,0,1,1);
  p1->SetBottomMargin(ysplit);
  p2->SetTopMargin(1-ysplit);
  p1->SetFillStyle(0);
  p2->SetFillStyle(0);

  p1->cd();


  std::vector <int> colors = {kMagenta,
                              kGreen+2,
                              kAzure,
                              kMagenta-2,
                              kGreen-3,
                              kRed};
  std::vector<std::string> labels = {"#nu_{e} CC",
                                     "#nu_{#mu} CC",
                                     "NC",
                                     "#bar{#nu}_{e} CC",
                                     "#bar{#nu}_{#mu} CC"};
                                     
                

  decomp[tot]->GetYaxis()->SetTitle(POTstr(datapot));
  decomp[tot]->GetYaxis()->CenterTitle();

  SetAxis(decomp[tot], true, colors[tot]);
  double maxY = decomp[tot]->GetMaximum();
  decomp[tot]->GetYaxis()->SetRangeUser(0, 1.30*maxY);
  decomp[tot]->Draw("hist");
  Legend(decomp[tot]);

  // Makes it look fancy
  Nue2018ThreeBinDivisions(true);
  Nue2018ThreeBinLabels(1.42*maxY, 0.06, kGray+3, true);

  SetAxis(mc[tot], false, colors[tot], 7);
  mc[tot]->Draw("hist same");
  for(int j = 0; j < (int)mc.size()-1 ; j++) {
    SetAxis(decomp[j], false, colors[j]);
    decomp[j]->Draw("hist same");
    SetAxis(mc[j], false, colors[j], 7);
    mc[j]->Draw("hist same");
  }

  hdata->SetMarkerStyle(kFullCircle);
  SetAxis(hdata, false, 1);
  hdata->Draw("ep same");
  
  // Add fancy axes
  Nue2018ThreeBinAxis(decomp[tot], false, true, true); 
  decomp[tot]->GetXaxis()->SetLabelSize(0.);
  decomp[tot]->GetXaxis()->SetTitleSize(0.);
  decomp[tot]->GetYaxis()->SetTitleSize(0.75*decomp[tot]->GetYaxis()->GetTitleSize());
  decomp[tot]->GetYaxis()->SetRangeUser(0.001, 
                            decomp[tot]->GetMaximum());

  // The empty bins of the histograms still get drawn at y=0. Try to 
  // draw over these colored lines with black
  Zero();

  // Draw the ratio
  p2->cd();

  SetAxis(rdecomp, true, 1, 1, false);
  rdecomp->Draw("ep");

  SetAxis(rmc, false, 17, 1, false);
  rmc->Draw("ep same");
  One();
  Nue2018ThreeBinDivisions(true);
  Nue2018ThreeBinAxis(rdecomp, true, true, true);

  rdecomp->GetYaxis()->SetTitleSize(0.75*rdecomp->GetYaxis()->GetTitleSize());
  rmc->GetYaxis()->SetRangeUser(.5,1.5);
  rdecomp->GetYaxis()->SetRangeUser(.5,1.5);
  RatioLegend();
  Preliminary();
  NeutrinoMode(mode);

  p1->Draw();
  p2->Draw("same");
  
  gPad->Print("Energy_Ana2018_Bin_CVNeLow_split.pdf");
}

TString POTstr

(

double

pot

)

Definition at line 205 of file plot_nd_spectra_2018.C.

References genie::utils::style::Format().

Referenced by plot_nd_spectra_2018().

                           {
  return TString(TString::Format("10^{3} Events / %.3g#times10^{20} POT", pot/1e20)); 
}

void RatioLegend

(

)

Definition at line 296 of file plot_nd_spectra_2018.C.

References make_syst_table_plots::h, and MECModelEnuComparisons::leg.

Referenced by plot_nd_spectra_2018().

                  {
  TLegend* leg = new TLegend(.12, .13, .47, .18);
  //  leg->SetFillStyle(0);
  leg->SetNColumns(2);
  leg->SetBorderSize(1);
  TH1* h = new TH1F("","", 1, 0, 1);
  h->SetMarkerColor(17);
  h->SetLineColor(17);
  h->SetMarkerStyle(kFullCircle);
  leg->AddEntry(h->Clone(), "Uncorr. MC", "lep");
  h->SetMarkerColor(1);
  h->SetLineColor(1);
  leg->AddEntry(h->Clone(), "Decomp MC", "lep");
  leg->Draw();
}

void SetAxis	(	TH1 *	h,
		bool	visible,
		int	color,
		int	style = 1,
		bool	scale = true
	)

Definition at line 194 of file plot_nd_spectra_2018.C.

References DerivePowScale(), and norm.

Referenced by plot_nd_spectra_2018().

                                                                     {
  if(scale) h->Scale(.001);
  h->SetLineStyle(style);
  h->SetLineColor(color);
  if(!visible){
    h->SetTickLength(0.);
    h->SetLabelSize(0.);
  }
  double norm = DerivePowScale(h);
}

void Zero

(

)

Definition at line 228 of file plot_nd_spectra_2018.C.

References MECModelEnuComparisons::g.

Referenced by stan::mcmc::base_nuts< Model, diag_e_metric, expl_leapfrog, BaseRNG >::build_tree(), stan::variational::normal_meanfield::calc_grad(), stan::variational::normal_fullrank::calc_grad(), cmf::ChiSqrCalculator::ChiSqrCovMat(), stan::mcmc::mock_hamiltonian< mock_model, rng_t >::dphi_dq(), stan::mcmc::divergent_hamiltonian< mock_model, rng_t >::dphi_dq(), stan::mcmc::divergent_hamiltonian< mock_model, rng_t >::dtau_dp(), stan::mcmc::diag_e_metric< Model, BaseRNG >::dtau_dq(), stan::mcmc::unit_e_metric< stan_model, boost::ecuyer1988 >::dtau_dq(), stan::mcmc::dense_e_metric< Model, BaseRNG >::dtau_dq(), stan::mcmc::mock_hamiltonian< mock_model, rng_t >::dtau_dq(), stan::mcmc::divergent_hamiltonian< mock_model, rng_t >::dtau_dq(), stan::math::exp_action_chain_vd(), osc::OscCalcPMNSOptEigen::FillCache(), ana::LikelihoodCovMxExperiment::GetExpectedSpectrum(), ana::Hist::Initialized(), ana::LikelihoodCovMxExperiment::LikelihoodCovMxExperiment(), stan::math::internal::ops_partials_edge< Dx, std::vector< Eigen::Matrix< fvar< Dx >, R, C > > >::ops_partials_edge(), stan::math::internal::ops_partials_edge< double, std::vector< Eigen::Matrix< var, R, C > > >::ops_partials_edge(), plot_nd_spectra_2018(), stan::variational::normal_meanfield::set_to_zero(), stan::variational::normal_fullrank::set_to_zero(), eta_adapt_big_test::SetUp(), eta_should_fail_test::SetUp(), eta_adapt_small_test::SetUp(), advi_test::SetUp(), eta_adapt_test::SetUp(), stochastic_gradient_ascent_test::SetUp(), TEST(), TEST_F(), stan::mcmc::base_nuts< Model, diag_e_metric, expl_leapfrog, BaseRNG >::transition(), stan::math::wishart_rng(), cmf::CMFDecorrelator::writeResults(), and ana::Hist::Zero().

            {
  TGraph* g = new TGraph();
  g->SetPoint(0, -1e10,0.003);
  g->SetPoint(1, 1e10, 0.003);
  g->SetLineWidth(2);
  g->SetLineStyle(1);
  g->Draw("same");

}