prettyPlotComps.py

Go to the documentation of this file.

import ROOT,sys
from array import array

ROOT.TGaxis.SetMaxDigits(3)

plotvars = ["flux","wrongsign"]

#Draw
canvas = ROOT.TCanvas("canvas","canvas");

for i in range(0,1):

  plotvar = plotvars[i]

  # Get the histo files
  f_1 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/OptimizedEngineeredNov2017/neutrino/flux//histos_g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_neutrino_LBNEFD_fastmc.root")  
#  f_2 = ROOT.TFile("/dune/data2/users/ljf26/fluxfiles/g4lbne/v3r3p7/QGSP_BERT/NuMI_Opt_2015CDR_Be_80GeV_WideDecayPipe/neutrino/flux//histos_g4lbne_v3r3p7_QGSP_BERT_NuMI_Opt_2015CDR_Be_80GeV_WideDecayPipe_neutrino_LBNEFD_fastmc.root")
#  f_3 = ROOT.TFile("/dune/data2/users/ljf26/fluxfiles/g4lbne/v3r3p7/QGSP_BERT/NuMI_Opt_2015CDR_Be_80GeV_LongDecayPipe/neutrino/flux//histos_g4lbne_v3r3p7_QGSP_BERT_NuMI_Opt_2015CDR_Be_80GeV_LongDecayPipe_neutrino_LBNEFD_fastmc.root")
#  f_4 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r3p6/QGSP_BERT/CP_run5_9116_80GeV/neutrino/flux//histos_g4lbne_v3r3p6_QGSP_BERT_CP_run5_9116_80GeV_neutrino_LBNEFD_fastmc.root")
#  f_5 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r3p6/QGSP_BERT/CP_run5_9116_80GeV_LongDecayPipe/neutrino/flux//histos_g4lbne_v3r3p6_QGSP_BERT_CP_run5_9116_80GeV_LongDecayPipe_neutrino_LBNEFD_fastmc.root")
#  f_6 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/Perfect_80GeV_LongDecayPipe/neutrino/flux//histos_g4lbne_v3r5p4_QGSP_BERT_Perfect_80GeV_LongDecayPipe_neutrino_LBNEFD_fastmc.root")

  """
  if plotvar=="wrongsign":
    f_1 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/NuMI_Opt_2015CDR_80GeV/antineutrino/flux//histos_g4lbne_v3r5p4_QGSP_BERT_NuMI_Opt_2015CDR_80GeV_antineutrino_LBNEFD_fastmc.root")  
    f_2 = ROOT.TFile("/dune/data2/users/ljf26/fluxfiles/g4lbne/v3r3p7/QGSP_BERT/NuMI_Opt_2015CDR_Be_80GeV_WideDecayPipe/antineutrino/flux//histos_g4lbne_v3r3p7_QGSP_BERT_NuMI_Opt_2015CDR_Be_80GeV_WideDecayPipe_antineutrino_LBNEFD_fastmc.root")
    f_3 = ROOT.TFile("/dune/data2/users/ljf26/fluxfiles/g4lbne/v3r3p7/QGSP_BERT/NuMI_Opt_2015CDR_Be_80GeV_LongDecayPipe/antineutrino/flux//histos_g4lbne_v3r3p7_QGSP_BERT_NuMI_Opt_2015CDR_Be_80GeV_LongDecayPipe_antineutrino_LBNEFD_fastmc.root")
    f_4 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r3p6/QGSP_BERT/CP_run5_9116_80GeV/antineutrino/flux//histos_g4lbne_v3r3p6_QGSP_BERT_CP_run5_9116_80GeV_antineutrino_LBNEFD_fastmc.root")
    f_5 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r3p6/QGSP_BERT/CP_run5_9116_80GeV_LongDecayPipe/antineutrino/flux//histos_g4lbne_v3r3p6_QGSP_BERT_CP_run5_9116_80GeV_LongDecayPipe_antineutrino_LBNEFD_fastmc.root")
    f_6 = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r5p4/QGSP_BERT/Perfect_80GeV_LongDecayPipe/antineutrino/flux//histos_g4lbne_v3r3p6_QGSP_BERT_Perfect_80GeV_LongDecayPipe_antineutrino_LBNEFD_fastmc.root")
  """

  histo_title = "numu_flux"
  
  h_1 = f_1.Get(histo_title)
#  h_2 = f_2.Get(histo_title)
#  h_3 = f_3.Get(histo_title)
#  h_4 = f_4.Get(histo_title)
#  h_5 = f_5.Get(histo_title)
#  h_6 = f_6.Get(histo_title)

  if(plotvar=="flux" or plotvar=="wrongsign"):
    # Rebin and convert to /GeV units
    newbins = [0.0]
    for i in range(61):
      newbins.append((i+1)*0.125)
#    for i in range(4):
#      newbins.append(8.0 + (i+1)*0.5)
  #for i in range(24):
  #    newbins.append(8.0 + (i+1)*0.5)
  #for i in range(7):
  #    newbins.append(20.0 + (i+1)*2)

    h_1 = h_1.Rebin(len(newbins)-1,"1_flux",array('d',newbins))
    #h_2 = h_2.Rebin(len(newbins)-1,"2_flux",array('d',newbins))
    #h_3 = h_3.Rebin(len(newbins)-1,"3_flux",array('d',newbins))
    #h_4 = h_4.Rebin(len(newbins)-1,"4_flux",array('d',newbins))
    #h_5 = h_5.Rebin(len(newbins)-1,"5_flux",array('d',newbins))
    #h_6 = h_6.Rebin(len(newbins)-1,"6_flux",array('d',newbins))


    h_1.Scale(1,"width") # convert to per GeV from flux per bin
    #h_2.Scale(1,"width") # convert to per GeV from flux per bin
    #h_3.Scale(1,"width") # convert to per GeV from flux per bin
    #h_4.Scale(1,"width") # convert to per GeV from flux per bin
    #h_5.Scale(1,"width") # convert to per GeV from flux per bin
    #h_6.Scale(1,"width") # convert to per GeV from flux per bin

    # Scale to flux per year instead of flux per POT
    scale_factor_120GeV = 1.1e21
    scale_factor_80GeV = 1.47e21
    scale_factor_60GeV = 1.89e21
      
    h_1.Scale(scale_factor_120GeV);
    #h_2.Scale(scale_factor_80GeV);
    #h_3.Scale(scale_factor_80GeV);
    #h_4.Scale(scale_factor_80GeV);
    #h_5.Scale(scale_factor_80GeV);
    #h_6.Scale(scale_factor_80GeV);
  
  # Rename axes appropriately
  y_title = "Unoscillated #nu_{#mu}s /  GeV / m^{2} / Year"
  h_1.GetYaxis().SetTitle(y_title);
  #h_2.GetYaxis().SetTitle(y_title);
  #h_3.GetYaxis().SetTitle(y_title);
  #h_4.GetYaxis().SetTitle(y_title);
  #h_5.GetYaxis().SetTitle(y_title);
  #h_6.GetYaxis().SetTitle(y_title);

  h_1.GetXaxis().SetTitle("Neutrino Energy (GeV)")

  h_1.SetLineColor(2);
  #h_2.SetLineColor(6);
  #h_3.SetLineColor(1);
  #h_4.SetLineColor(2);
  #h_5.SetLineColor(4);
  #h_6.SetLineColor(15);

  #h_2.SetMaximum(100e9)
  
  a_title = "#nu_{#mu} Flux, #nu Mode, Far Detector"
  if(plotvar=="wrongsign"):
    a_title = "#nu_{#mu} Flux, #bar{#nu} Mode"
  ROOT.gStyle.SetTitleY(0.93);
  ROOT.gStyle.SetTitleFont(42,"t");

  h_1.SetTitle(a_title);
  #h_2.SetTitle(a_title);
  #h_3.SetTitle(a_title);
  #h_4.SetTitle(a_title);
  #h_5.SetTitle(a_title);
  #h_6.SetTitle(a_title);

  #maxes = [h_1.GetMaximum(),h_2.GetMaximum(),h_3.GetMaximum(),h_4.GetMaximum(),h_5.GetMaximum(),h_6.GetMaximum()]
  plot_max = h_1.GetMaximum()*1.2
#  if(plotvar=="flux"):
#    plot_max = 70e9


  for histo in [h_1]:
    histo.GetXaxis().CenterTitle()
    histo.GetYaxis().CenterTitle()
    histo.SetMaximum(plot_max)

  h_1.DrawCopy("hist");


  ROOT.gStyle.SetOptTitle(1)
  ROOT.gPad.Update()

  title = ROOT.gPad.GetPrimitive("title");
  title.SetBorderSize(0);
  #title.SetTextSize(0.04);
  #title.SetTextFont(42)
  
  #canvas.SetLogx();
  #h_6.DrawCopy("histsame");
#  h_2.DrawCopy("histsame");
#  h_3.DrawCopy("histsame");
  h_1.DrawCopy("histsame");
#  h_5.DrawCopy("histsame");
  #h_1.DrawCopy("histsame");

  if(plotvar=="flux"):
    leg = ROOT.TLegend(0.5,0.6,0.9,0.9);
    #leg.AddEntry(h_6,"Perfect Focusing, 250x4 m DP","l");
 #   leg.AddEntry(h_5,"Optimized, 241x4 m DP","l");
#    leg.AddEntry(h_4,"Optimized","l");
#    leg.AddEntry(h_3,"Enhanced Reference, 250x4 m DP","l");
#    leg.AddEntry(h_2,"Enhanced Reference","l");
    leg.AddEntry(h_1,"DUNE Flux (Nov 2017)","l");
    leg.SetFillStyle(0);
    leg.SetTextFont(42);
  leg.Draw();
  
  ROOT.gPad.RedrawAxis()
  ROOT.gPad.Update()

"""
  if not plotvar=="flux":
    latex = ROOT.TLatex()
    latex.SetTextSize(0.07);
    latex.SetTextColor(14)
    latex.SetNDC()
    latex.DrawLatex(.2,.3,"PLACEHOLDER PLOT");
"""
ROOT.gEnv.Print()
canvas.Print("DUNE_Numu_Flux_Nov2017.eps")
canvas.Print("DUNE_Numu_Flux_Nov2017.pdg")
canvas.Print("DUNE_Numu_Flux_Nov2017.png")