make_syst_table_plots.py

Go to the documentation of this file.

import os, sys, shlex, subprocess
from optparse import OptionParser
from ROOT import *

gROOT.SetBatch( kTRUE )

import root_style

gROOT.SetBatch( kTRUE )
root_style.SetRootEnv()


systs_xsec = (
"CCQE z-exp EV shift #1",
"CCQE z-exp EV shift #2",
"CCQE z-exp EV shift #3",
"CCQE z-exp EV shift #4",
"MaCCRES",
"MvCCRES",
"MaNCRES",
"MvNCRES",
"ZNormCCQE",
"RPA shape: higher-Q^{2} enhancement (2020)",
"RPA shape: low-Q^{2} suppression (2020)",
"RES low-Q^2 suppression",
"DIS vnCC1pi",
"hN FSI mean free path",
"hN FSI fate fraction eigenvector #1",
"MEC E_{#nu} shape, neutrinos",
"MEC E_{#nu} shape, antineutrinos",
"MEC 2020 (q_{0}, |#vec{q}|) response, neutrinos",
"MEC 2020 (q_{0}, |#vec{q}|) response, antineutrinos",
"MEC initial state np fraction, neutrinos",
"MEC initial state np fraction, antineutrinos",
"Radiative corrections for #nu_{e}",
"Radiative corrections for #bar{#nu}_{e}",
"Second class currents",
"Genie PC 0",
"Genie PC 1",
"Genie PC 2",
"Genie PC 3",
"Genie PC 4",
"#nu_{#tau} Scale",
)

systs_other = (
"Flux Component 00",
"Flux Component 01",
"Flux Component 02",
"Flux Component 03",
"Flux Component 04",
"AbsCalib",
"RelCalib",
"CalibShape",
"CalibDrift",
"Lightlevel",
"Cherenkov",
"Uncorr ND Mu Energy Scale 2020",
"Uncorr MuCat Mu Energy 2020",
"Neutron Pile-up 2020",
"Corr Mu Energy Scale 2020",
"Uncorr FD Mu Energy Scale 2020",
"Lepton Angle Syst ND XZ",
"Lepton Angle Syst ND YZ",
"Lepton Angle Syst FD XZ",
"Lepton Angle Syst FD YZ",
"Neutron visible energy systematic 2018",
"Acceptance ND to FD Kinematics Signal FHC 2020",
"Acceptance ND to FD Kinematics Signal RHC 2020",
"Michel Electrons Tagging Uncertainty",
)


def ReplaceSystName( syst_name ) :

  syst_name = syst_name.replace( "antineutrinos", "#bar{#nu}" )
  syst_name = syst_name.replace( "neutrinos", "#nu" )
  syst_name = syst_name.replace( "CCQE z-exp EV shift", "CCQE Z-Exp EV" )
  syst_name = syst_name.replace( "ZNormCCQE", "CCQE Z-Exp Norm." )
  syst_name = syst_name.replace( "RPA shape: higher-Q^{2} enhancement (2020)", "CCQE RPA High-Q^{2} Enh." )
  syst_name = syst_name.replace( "RPA shape: low-Q^{2} suppression (2020)", "CCQE RPA Low-Q^{2} Supp." )
  syst_name = syst_name.replace( "RES low-Q^2 suppression", "RES Low-Q^{2} Supp." )
  syst_name = syst_name.replace( "hN FSI mean free path", "hN FSI MFP" )
  syst_name = syst_name.replace( "hN FSI fate fraction eigenvector #1", "hN FSI Fate Frac. EV #1" )
  syst_name = syst_name.replace( "MEC E_{#nu} shape", "MEC E_{#nu} Shape" )
  syst_name = syst_name.replace( "MEC 2020 (q_{0}, |#vec{q}|) response", "MEC (q_{0},#left|#vec{q}#right|)" )
  syst_name = syst_name.replace( "Radiative corrections for", "Rad. Corr." )
  syst_name = syst_name.replace( "MEC initial state np fraction", "MEC np Frac." )
  syst_name = syst_name.replace( "AbsCalib", "Abs. Calib." )
  syst_name = syst_name.replace( "RelCalib", "Rel. Calib." )
  syst_name = syst_name.replace( "CalibShape", "Calib. Shape" )
  syst_name = syst_name.replace( "CalibDrift", "Calib. Drift" )
  syst_name = syst_name.replace( "Lightlevel ", "Light Level" )
  syst_name = syst_name.replace( "Neutron visible energy systematic 2018", "Neutron Vis. E" )
  syst_name = syst_name.replace( "Uncorr ND Mu Energy Scale 2020", "E_{#mu} Scale ND" )
  syst_name = syst_name.replace( "Uncorr FD Mu Energy Scale 2020", "E_{#mu} Scale FD" )
  syst_name = syst_name.replace( "Corr Mu Energy Scale 2020", "E_{#mu} Scale Corr" )
  syst_name = syst_name.replace( "Uncorr MuCat Mu Energy 2020", "E_{#mu} #mu Catcher" )
  syst_name = syst_name.replace( "Neutron Pile-up 2020", "E_{#mu} Neutron Pile-up ND" )
  syst_name = syst_name.replace( "Acceptance ND to FD Kinematics Signal FHC 2020", "Nue Signal Accept. FHC" )
  syst_name = syst_name.replace( "Acceptance ND to FD Kinematics Signal RHC 2020", "Nue Signal Accept. RHC" )
  syst_name = syst_name.replace( "Michel Electrons Tagging Uncertainty", "Michel Tagging" )
  syst_name = syst_name.replace( "Flux Component 0", "Flux PC " )

  return syst_name

def GetSystTable( fname, systs, addLepAng = False ) :

  table = []
  for line in open( fname ) :

    if "Source of Uncertainty" in line : continue
    if "hline" in line : continue

    line = line.strip()
    if not "Statistical" in line : line = line[ 0 : line.find( "\\" ) ]
    line = line.split( " & " )
    syst_name = line[0]

    if   systs == "xsec"  and syst_name not in systs_xsec  : continue
    elif systs == "other" and syst_name not in systs_other : continue

    table.append( line )

    # Ugly hack to add zero-filled lepton angle uncertainty for no-pT extrap
    if addLepAng and "Detector Response" in syst_name :
      table.append( ( "Lepton Angle", "+0 / -0", "+0 / -0", "+0 / -0" ) )
    if addLepAng and "Uncorr FD Mu Energy Scale 2020" in syst_name :
      table.append( ( "Lepton Angle Syst ND XZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )
      table.append( ( "Lepton Angle Syst ND YZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )
      table.append( ( "Lepton Angle Syst FD XZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )
      table.append( ( "Lepton Angle Syst FD YZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )

  return table


if __name__ == "__main__" :

  ###############################
  # Parse Options
  ###############################

  systs_options = ( "summary", "xsec", "other" )

  parser = OptionParser()
  parser.add_option( "--systs"  , dest = "systs"  , type = "str"         , default = None , help = "Systs to plot: " + ", ".join( systs_options ) )
  parser.add_option( "--compare", dest = "compare", action = "store_true", default = False, help = "Plot tables for pT vs. no pT extrap comparison" )
  parser.add_option( "--pt"     , dest = "pt"     , action = "store_true", default = False, help = "Plot tables for pT extrap" )

  if len(sys.argv) < 2 : parser.parse_args( "--help".split() )

  ( opts, args ) = parser.parse_args()

  if not opts.systs : sys.exit( "Specify systs" )
  if opts.systs not in systs_options : sys.exit( "Wrong systs" )

  systs = opts.systs

  isSummary = systs == "summary"
  isXSec    = systs == "xsec"
  isOther   = systs == "other"

  table_type = "full"
  if isSummary : table_type = "short"

  box_colors = ( kRed+1, kOrange+1, kGreen+1, kBlue )

  systs_tables = []
  out_tag = ""
  if opts.pt :
    out_tag = "pt"
    systs_tables.append( GetSystTable( "tables/syst_table_ptExtrap_%s.txt" % table_type, systs ) )
  elif opts.compare :
    out_tag = "pt_vs_no_pt"
    systs_tables.append( GetSystTable( "tables/syst_table_%s.txt" % table_type, systs, True ) )
    systs_tables.append( GetSystTable( "tables/syst_table_ptExtrap_%s.txt" % table_type, systs ) )
  else :
    sys.exit( "Please specify \"pt\" or \"compare\"" )

  plots_dir = "./syst_plots"
  if not os.path.isdir( plots_dir ) : os.makedirs( plots_dir )

  ntables = len( systs_tables )
  nsysts = len( systs_tables[0] )

  if ntables > len( box_colors ) : sys.exit( "Fewer colors than tables" )

  c = TCanvas()
  canvas_width = 1000
  canvas_height = 800
  if isXSec :
    canvas_height = 1200
  c.SetCanvasSize( canvas_width, canvas_height )
  c.SetLeftMargin( 0.3 )
  c.SetRightMargin( 0.1 )
  c.SetTopMargin( 0.1 )
  c.SetGridy( 1 )

  osc_vars = ( "s23", "dCP", "m32" )
  for var in osc_vars :

    var_index = osc_vars.index( var ) + 1

    xtitle = ""
    xmax = 1.0
    ndiv = 502
    if   var == "s23" :
      xtitle = "sin^{2}#theta_{23}"
      if   isXSec :
        xmax = 0.004
        ndiv = 402
      elif isOther :
        xmax = 0.02
        ndiv = 402
      else :
        xmax = 0.04
        ndiv = 402
    elif var == "dCP" :
      xtitle = "#delta_{CP} / #pi"
      if   isXSec :
        xmax = 0.05
        ndiv = 502
      elif isOther :
        xmax = 0.04
        ndiv = 402
      else :
        xmax = 0.3
        ndiv = 302
    elif var == "m32" :
      xtitle = "#Delta m^{2}_{32} (#times 10^{-3} eV^{2})"
      if   isXSec :
        xmax = 0.01
      elif isOther :
        xmax = 0.02
        ndiv = 502
      else :
        xmax = 0.05
        ndiv = 502

    h = TH2F( var, "", 10, -xmax, +xmax, nsysts, 0, nsysts )

    boxes = []

    for syst_idx in range( nsysts ) :

      ibin = nsysts - syst_idx

      syst_name = systs_tables[0][syst_idx][0]
      syst_name = ReplaceSystName( syst_name )

      h.GetYaxis().SetBinLabel( ibin, syst_name )
      
      for table_idx in range( len( systs_tables ) ) :

        errors = systs_tables[table_idx][syst_idx][var_index].split( "/" )

        xlo = float( errors[1] )
        xhi = float( errors[0] )
        space = 0.25
        dy = ( 1.0 - space ) / ntables
        yhi = ibin - ( space / 2 ) - table_idx * dy
        ylo = yhi - dy
        boxes.append( TBox( xlo, ylo, xhi, yhi ) )
        boxes[-1].SetFillColor( box_colors[ table_idx ] )
        boxes[-1].SetLineWidth( 0 )


    line = TLine()
    line.SetLineWidth( 1 )
    line.SetLineStyle( 3 )
    line.SetLineColor( kBlack )

    root_style.ApplyAxisStyleAll( h ) 
    h.GetXaxis().SetTitle( "Uncertainty in " + xtitle )
    if isXSec : h.GetYaxis().SetLabelSize( 0.035 )
    h.GetXaxis().SetNdivisions( ndiv )
    h.GetXaxis().SetNoExponent( kTRUE )
    h.Draw()
    for b in boxes : b.Draw()
    line.DrawLine( 0, 0, 0, nsysts ) 
    for outType in ( "pdf", "C" ) :
      c.SaveAs( os.path.join( plots_dir, "syst_table_%s_%s_%s.%s" % ( var, out_tag, systs, outType ) ) )