make_true_q0q3_plots.py

Go to the documentation of this file.

import os, sys
from array import array
from ROOT import *
from math import *
from optparse import OptionParser

import root_style
#import canMan

gROOT.SetBatch( kTRUE )
root_style.SetRootEnv()

def RebinHistQ0Q3( hist ) :

  hist.RebinX( 2 )
  hist.RebinY( 2 )

def SetRangeQ0Q3( hist ) :

  hist.GetXaxis().SetRangeUser( 0.0, 1.0 )
  hist.GetYaxis().SetRangeUser( 0.0, 0.8 )
  #hist.GetYaxis().SetLimits( 1.0e-3, 0.8 / hist.GetYaxis().GetXmax() )
  #print hist.GetYaxis().GetXmax()

def GetMaximum( hist ) :

  return hist.GetBinContent( hist.GetMaximumBin() )

def SetMECWeightsPalette() :

  # A colour palette that goes blue->white->red, useful for correlation matrices

  n_color_contours = 999
  #n_color_contours = 25
  NRGBs = 3;
  initialized = False
  colors = array( "i" )

  stops = array( "d", [ 0.00, 0.20, 1.00 ] )
  red   = array( "d", [ 0.00, 1.00, 1.00 ] )
  green = array( "d", [ 0.00, 1.00, 0.00 ] )
  blue  = array( "d", [ 1.00, 1.00, 0.00 ] )
  colmin = TColor.CreateGradientColorTable( NRGBs, stops, red, green, blue, n_color_contours )
  for i in range( n_color_contours ) : colors.insert( i, colmin + i )

  gStyle.SetNumberContours( n_color_contours )
  gStyle.SetPalette( n_color_contours, colors )

def SetTitleQ0( hist ) :

  hist.GetXaxis().SetTitle( "True #left|#vec{q}#right| (GeV/c)" )

def SetDecimalsAllAxes( hist ) :

  hist.GetXaxis().SetDecimals( kTRUE )
  hist.GetYaxis().SetDecimals( kTRUE )
  hist.GetZaxis().SetDecimals( kTRUE )

def SetHistLineAtt( hist, width, color ) :

  hist.SetLineWidth( width )
  hist.SetLineColor( color )

if __name__ == "__main__":

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

  #default_indir  = "./"
  default_indir = "/pnfs/nova/persistent/users/mislivec/xsec_tuning_paper/"
  default_outdir = './plots/true_q0q3'

  #default_indir  = "/pnfs/nova/persistent/users/mislivec/xs_tuning_hists/rpa_fix"
  #default_outdir = "/nova/app/users/mislivec/scripts_xs_tuning_blessed_plots/plots_true_q0q3_rpa_fix"

  parser = OptionParser()
  parser.add_option( '--beam'  , dest = 'beam'  , default = None          , help = 'Beam (FHC or RHC)' )
  parser.add_option( '--sim', dest = 'sim', action = "store_true", default = False, help = "Add \"NOvA Simulation\" label" )
  #parser.add_option( '--outdir', dest = 'outdir', default = default_outdir, help = 'Output plots directory:\n%s' % default_outdir )
  #parser.add_option( '--indir' , dest = 'indir' , default = default_indir , help = 'Input plots directory:\n%s' % default_indir )

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

  ( opts, args ) = parser.parse_args()

  if not opts.beam : sys.exit( 'Beam (FHC or RHC) requred' )

  #hists_file_name_data = os.path.join( opts.indir, "hists_2018_%s_data_full.root" % opts.beam )
  #hists_file_name_mc   = os.path.join( opts.indir, "hists_2018_%s_mc_full.root"   % opts.beam ) 

  hists_file_name_data = os.path.join( default_indir, "hists_2018_%s_data_full.root" % opts.beam )
  hists_file_name_mc   = os.path.join( default_indir, "hists_2018_%s_mc_full.root"   % opts.beam ) 

  if not os.path.exists( hists_file_name_data ) : sys.exit( 'File does not exist:\n%s' % hists_file_name_data )
  if not os.path.exists( hists_file_name_mc   ) : sys.exit( 'File does not exist:\n%s' % hists_file_name_mc )
  print 'Input MC histogram file:\n', hists_file_name_data
  print 'Input MC histogram file:\n', hists_file_name_mc

  hists_file_data = TFile( hists_file_name_data )
  hists_file_mc   = TFile( hists_file_name_mc   )

  tree_pot_data = hists_file_data.Get( 'pot' )
  pot_data = 0.0
  for i in range( tree_pot_data.GetEntries() ) :
    tree_pot_data.GetEntry( i )
    pot_data += tree_pot_data.pot_data

  tree_pot_mc = hists_file_mc.Get( 'pot' )
  pot_mc = 0.0
  for i in range( tree_pot_mc.GetEntries() ) :
    tree_pot_mc.GetEntry( i )
    pot_mc += tree_pot_mc.pot_mc

  mc_pot_scale = pot_data / pot_mc

  print 'pot_data = %.2e' % pot_data
  print 'pot_mc   = %.2e' % pot_mc
  print 'mc_pot_scale = %.3f' % mc_pot_scale

  event_scale_power_1d = 4
  event_scale_1d = 1.0 / pow( 10.0, event_scale_power_1d )
  text_events_1d = '10^{%d} Events' % event_scale_power_1d
  scale_hist_1d = mc_pot_scale * event_scale_1d 

  for k in hists_file_mc.GetListOfKeys() :
    name = k.GetName()
    #print name
    #if 'TrueQ0_vs_TrueQ3' in name : print name
    #if 'TrueQ0_vs_TrueQ3_MECShapeUp_MEC_' in name : print name
    #if "TrueQ0_vs_TrueQ3_MECShape" in name : print name

  isRHC = opts.beam.upper() == 'RHC'

  tag_neutrino_mec_weights = 'Nu'
  if isRHC : tag_neutrino_mec_weights = 'Nubar'

  weights_file_name = os.path.join( os.environ[ 'SRT_PUBLIC_CONTEXT' ], "MCReweight/data/xs/rw_empiricalMEC2018_RPAfix_%s.root" % tag_neutrino_mec_weights.lower() )
  if not os.path.exists( weights_file_name ) : sys.exit( 'File does not exist:\n%s' % weights_file_name )

  print 'Input MEC weights file:'
  print weights_file_name

  #sys.exit()

  weights_file = TFile( weights_file_name )
  for k in weights_file.GetListOfKeys() : print k.GetName()

  #plots_dir = opts.outdir
  plots_dir = default_outdir
  if not os.path.isdir( plots_dir ) : os.makedirs( plots_dir )

  tag_neutrino_mec_hists = 'Numu'
  if isRHC : tag_neutrino_mec_hists = 'Numubar'

  hist_nominal    = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_Nominal_%s'      % tag_neutrino_mec_hists )
  hist_tuned      = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_%s'              % tag_neutrino_mec_hists )
  hist_valencia   = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_Valencia_%s'     % tag_neutrino_mec_hists )
  hist_minerva    = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_MINERvA_%s'      % tag_neutrino_mec_hists )
  hist_shift_up   = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeUp_MEC_%s'   % tag_neutrino_mec_hists )
  hist_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeDown_MEC_%s' % tag_neutrino_mec_hists )
  hist_weights    = weights_file .Get( '%s_mec_weights_smoothed'              % tag_neutrino_mec_hists.lower() )

  '''
  hist_all_mec            = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC' )
  hist_all_mec_shift_up   = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeUp_MEC' )
  hist_all_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeDown_MEC' )
  hist_non_mec            = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_TotMC' )
  hist_non_mec_shift_up   = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeUp_TotMC' )
  hist_non_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeDown_TotMC' )
  hist_non_mec           .Add( hist_all_mec, -1.0 )
  hist_non_mec_shift_up  .Add( hist_all_mec_shift_up, -1.0 )
  hist_non_mec_shift_down.Add( hist_all_mec_shift_down, -1.0 )
  '''

  #hist_all_mec            = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC' )
  #hist_all_mec_shift_up   = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsQELike_MEC' )
  #hist_all_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsRESLike_MEC' )
  hist_non_mec            = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_TotMC' )
  hist_non_mec_shift_up   = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsQELike_TotMC' )
  hist_non_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsRESLike_TotMC' )
  hist_non_mec           .Add( hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC' ), -1.0 )
  hist_non_mec_shift_up  .Add( hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsQELike_MEC' ), -1.0 )
  hist_non_mec_shift_down.Add( hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsRESLike_MEC' ), -1.0 )


  root_style.ApplyAxisStyleAll( hist_nominal )
  root_style.ApplyAxisStyleAll( hist_tuned )
  root_style.ApplyAxisStyleAll( hist_valencia )
  root_style.ApplyAxisStyleAll( hist_minerva )
  root_style.ApplyAxisStyleAll( hist_shift_up )
  root_style.ApplyAxisStyleAll( hist_shift_down )
  root_style.ApplyAxisStyleAll( hist_weights )
  root_style.ApplyAxisStyleAll( hist_non_mec )
  root_style.ApplyAxisStyleAll( hist_non_mec_shift_up )
  root_style.ApplyAxisStyleAll( hist_non_mec_shift_down )

  '''
  hist_nominal           .Scale( mc_pot_scale * event_scale_1d )
  hist_tuned             .Scale( mc_pot_scale * event_scale_1d )
  hist_valencia          .Scale( mc_pot_scale * event_scale_1d )
  hist_minerva           .Scale( mc_pot_scale * event_scale_1d )
  hist_shift_up          .Scale( mc_pot_scale * event_scale_1d )
  hist_shift_down        .Scale( mc_pot_scale * event_scale_1d )
  hist_non_mec           .Scale( mc_pot_scale * event_scale_1d )
  hist_non_mec_shift_up  .Scale( mc_pot_scale * event_scale_1d )
  hist_non_mec_shift_down.Scale( mc_pot_scale * event_scale_1d )
  '''

  RebinHistQ0Q3( hist_nominal )
  RebinHistQ0Q3( hist_tuned )
  RebinHistQ0Q3( hist_valencia )
  RebinHistQ0Q3( hist_minerva )
  RebinHistQ0Q3( hist_shift_up )
  RebinHistQ0Q3( hist_shift_down )
  RebinHistQ0Q3( hist_non_mec )
  RebinHistQ0Q3( hist_non_mec_shift_up )
  RebinHistQ0Q3( hist_non_mec_shift_down )

  SetTitleQ0( hist_nominal )
  SetTitleQ0( hist_tuned )
  SetTitleQ0( hist_valencia )
  SetTitleQ0( hist_minerva )
  SetTitleQ0( hist_shift_up )
  SetTitleQ0( hist_shift_down )
  SetTitleQ0( hist_non_mec )
  SetTitleQ0( hist_non_mec_shift_up )
  SetTitleQ0( hist_non_mec_shift_down )
  SetTitleQ0( hist_weights ) 

  SetDecimalsAllAxes( hist_nominal )
  SetDecimalsAllAxes( hist_tuned )
  SetDecimalsAllAxes( hist_weights )
  #SetDecimalsAllAxes(  )

  margin_left  = gStyle.GetPadLeftMargin()
  margin_right = gStyle.GetPadRightMargin()
  margin_top   = gStyle.GetPadTopMargin()
  text_size = 0.05
  text_font = 42
  text_angle = 0.0
  setNDC = True

  xpos_label_mec     = 0.37
  xpos_label_no_mec  = 0.37
  xpos_legend_mec    = xpos_label_mec - 0.02
  xpos_legend_no_mec = xpos_label_no_mec - 0.02
  ypos_label_top     = 1.0 - gStyle.GetPadTopMargin() - 0.02
  ypos_label_mec     = ypos_label_top - text_size
  ypos_label_no_mec  = ypos_label_top - text_size
  if not isRHC :
    ypos_label_mec    -= 0.005
    #ypos_label_no_mec -= 0.005
  ypos_legend_mec    = ypos_label_mec    - text_size - 0.01
  ypos_legend_no_mec = ypos_label_no_mec - text_size - 0.02

  text_mec_neutrino = '#nu_{#mu}'
  if isRHC : text_mec_neutrino = '#bar{#nu}_{#mu}'

  text_beam = "Neutrino Beam"
  if isRHC : text_beam = "Antineutrino Beam"

  text_nova          = "NOvA" #"NO#nuA"
  text_minerva       = "MINERvA" #"MINER#nuA"
  #text_mec           = "%s CC MEC" % text_mec_neutrino
  text_numu_cc       = "%s CC" % text_mec_neutrino
  text_no_2p2h       = "#nu_{#mu} + #bar{#nu}_{#mu} CC No 2p2h"
  text_empirical_mec = "Empirical MEC"
  text_valencia_mec  = "Valencia MEC"
  text_nova_tune     = "2018 %s #nu + #bar{#nu} Tune" % text_nova
  text_nova_2p2h     = "%s 2p2h" % text_nova
  text_minerva_2p2h  = "%s 2p2h" % text_minerva
  text_simulation    = "%s Simulation" % text_nova

  label_beam       = root_style.GetLabel( text_beam   , xpos_label_mec   , ypos_label_top   , text_size, kBlack, text_font, 13, text_angle, setNDC )
  #label_mec        = root_style.GetLabel( text_mec    , xpos_label_mec   , ypos_label_mec   , text_size, kBlack, text_font, 13, text_angle, setNDC )
  label_numu_cc    = root_style.GetLabel( text_numu_cc, xpos_label_mec   , ypos_label_mec   , text_size, kBlack, text_font, 13, text_angle, setNDC )
  label_no_mec     = root_style.GetLabel( text_no_2p2h, xpos_label_no_mec, ypos_label_no_mec, text_size, kBlack, text_font, 13, text_angle, setNDC )
  label_simulation = root_style.GetLabel( text_simulation, 1.0 - margin_right, 1.0 - margin_top * 0.7, 2/30., kBlue , text_font, 31, text_angle, setNDC )

  canvas = TCanvas()

  for var in ( 'TrueQ0', 'TrueQ3' ) :

    xmax = 1.0
    if var == 'TrueQ0' :
      hist_1D_nominal    = hist_nominal   .ProjectionY( '%s_Empirical' % var )
      hist_1D_tuned      = hist_tuned     .ProjectionY( '%s_2018CV' % var )
      hist_1D_valencia   = hist_valencia  .ProjectionY( '%s_Valencia' % var )
      hist_1D_minerva    = hist_minerva   .ProjectionY( '%s_MINERvA' % var )
      hist_1D_shift_up   = hist_shift_up  .ProjectionY( '%s_ShiftUp' % var )
      hist_1D_shift_down = hist_shift_down.ProjectionY( '%s_ShiftDown' % var )

      hist_1D_non_mec            = hist_non_mec           .ProjectionY( '%s_NonMEC' % var )
      hist_1D_non_mec_shift_up   = hist_non_mec_shift_up  .ProjectionY( '%s_NonMEC_ShiftUp' % var )
      hist_1D_non_mec_shift_down = hist_non_mec_shift_down.ProjectionY( '%s_NonMEC_ShiftDown' % var )

      xmax = 1.0
    else :
      hist_1D_nominal    = hist_nominal   .ProjectionX( '%s_Empirical' % var )
      hist_1D_tuned      = hist_tuned     .ProjectionX( '%s_2018CV' % var )
      hist_1D_valencia   = hist_valencia  .ProjectionX( '%s_Valencia' % var )
      hist_1D_minerva    = hist_minerva   .ProjectionX( '%s_MINERvA' % var )
      hist_1D_shift_up   = hist_shift_up  .ProjectionX( '%s_ShiftUp' % var )
      hist_1D_shift_down = hist_shift_down.ProjectionX( '%s_ShiftDown' % var )

      hist_1D_non_mec            = hist_non_mec           .ProjectionX( '%s_NonMEC' % var )
      hist_1D_non_mec_shift_up   = hist_non_mec_shift_up  .ProjectionX( '%s_NonMEC_ShiftUp' % var )
      hist_1D_non_mec_shift_down = hist_non_mec_shift_down.ProjectionX( '%s_NonMEC_ShiftDown' % var )

      xmax = 1.2

    root_style.ApplyAxisStyleAll( hist_1D_nominal )
    root_style.ApplyAxisStyleAll( hist_1D_tuned )
    root_style.ApplyAxisStyleAll( hist_1D_valencia )
    root_style.ApplyAxisStyleAll( hist_1D_minerva )
    root_style.ApplyAxisStyleAll( hist_1D_shift_up )
    root_style.ApplyAxisStyleAll( hist_1D_shift_down )
    root_style.ApplyAxisStyleAll( hist_1D_non_mec )
    root_style.ApplyAxisStyleAll( hist_1D_non_mec_shift_up )
    root_style.ApplyAxisStyleAll( hist_1D_non_mec_shift_down )

    hist_1D_nominal   .GetXaxis().SetRangeUser( 0, xmax )
    hist_1D_tuned     .GetXaxis().SetRangeUser( 0, xmax )
    hist_1D_valencia  .GetXaxis().SetRangeUser( 0, xmax )
    hist_1D_minerva   .GetXaxis().SetRangeUser( 0, xmax )
    hist_1D_shift_up  .GetXaxis().SetRangeUser( 0, xmax )
    hist_1D_shift_down.GetXaxis().SetRangeUser( 0, xmax )
    #hist_1D_non_mec           .GetXaxis().SetRangeUser( 0, xmax )
    #hist_1D_non_mec_shift_up  .GetXaxis().SetRangeUser( 0, xmax )
    #hist_1D_non_mec_shift_down.GetXaxis().SetRangeUser( 0, xmax )

    hist_1D_nominal           .Scale( scale_hist_1d )
    hist_1D_tuned             .Scale( scale_hist_1d )
    hist_1D_valencia          .Scale( scale_hist_1d )
    hist_1D_minerva           .Scale( scale_hist_1d )
    hist_1D_shift_up          .Scale( scale_hist_1d )
    hist_1D_shift_down        .Scale( scale_hist_1d )
    hist_1D_non_mec           .Scale( scale_hist_1d )
    hist_1D_non_mec_shift_up  .Scale( scale_hist_1d )
    hist_1D_non_mec_shift_down.Scale( scale_hist_1d )

    hist_1D_nominal           .GetYaxis().SetTitle( text_events_1d )
    hist_1D_tuned             .GetYaxis().SetTitle( text_events_1d )
    hist_1D_valencia          .GetYaxis().SetTitle( text_events_1d )
    hist_1D_minerva           .GetYaxis().SetTitle( text_events_1d )
    hist_1D_shift_up          .GetYaxis().SetTitle( text_events_1d )
    hist_1D_shift_down        .GetYaxis().SetTitle( text_events_1d )
    hist_1D_non_mec           .GetYaxis().SetTitle( text_events_1d )
    hist_1D_non_mec_shift_up  .GetYaxis().SetTitle( text_events_1d )
    hist_1D_non_mec_shift_down.GetYaxis().SetTitle( text_events_1d )

    hist_line_width = 3
    SetHistLineAtt( hist_1D_nominal   , hist_line_width, kGreen+2 )
    SetHistLineAtt( hist_1D_tuned     , hist_line_width, kBlack )
    SetHistLineAtt( hist_1D_valencia  , hist_line_width, kRed )
    SetHistLineAtt( hist_1D_minerva   , hist_line_width, kBlue )
    SetHistLineAtt( hist_1D_shift_up  , hist_line_width, kBlue )
    SetHistLineAtt( hist_1D_shift_down, hist_line_width, kRed )
    SetHistLineAtt( hist_1D_non_mec           , hist_line_width, kBlack )
    SetHistLineAtt( hist_1D_non_mec_shift_up  , hist_line_width, kBlue )
    SetHistLineAtt( hist_1D_non_mec_shift_down, hist_line_width, kRed )

    max_nominal    = GetMaximum( hist_1D_nominal )
    max_tuned      = GetMaximum( hist_1D_tuned )
    max_valencia   = GetMaximum( hist_1D_valencia )
    max_minerva    = GetMaximum( hist_1D_minerva )
    max_shift_up   = GetMaximum( hist_1D_shift_up )
    max_shift_down = GetMaximum( hist_1D_shift_down )
    max_non_mec            = GetMaximum( hist_1D_non_mec )
    max_non_mec_shift_up   = GetMaximum( hist_1D_non_mec_shift_up )
    max_non_mec_shift_down = GetMaximum( hist_1D_non_mec_shift_down )

    legend1 = root_style.GetLegend()
    legend1.SetTextSize( text_size )
    legend1.SetTextFont( text_font )
    #legend1.SetHeader( '%s CC MEC' % text_mec_neutrino )
    legend1.AddEntry( hist_1D_tuned  , text_nova_2p2h  , 'L' )
    legend1.AddEntry( hist_1D_nominal, text_empirical_mec, 'L' )
    legend1.SetX1( xpos_legend_mec )
    legend1.SetY2( ypos_legend_mec )
    root_style.SetLegendPosY1( legend1 )

    canvas.Clear()

    hist_1D_tuned.SetMaximum( 1.8 * max( max_nominal, max_tuned ) )
    hist_1D_tuned.Draw( 'hist' )
    hist_1D_nominal.Draw( 'hist same' )
    hist_1D_tuned.Draw( 'hist same' )
    legend1.Draw()
    if opts.sim : label_simulation.Draw()
    label_beam.Draw()
    label_numu_cc.Draw()
    root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Empirical' % ( tag_neutrino_mec_hists, var ) ) )

    text_variable = 'blah'
    if var == 'TrueQ0' : 
      text_variable = 'true energy transfer, q_0,'
    else :
      text_variable = 'true three-momentum transfer magnitude, |q|,'

    text_caption = 'The %s near detector %s distribution for true %s CC MEC events that pass the 2018 numu selection.  ' % ( opts.beam.upper(), text_variable, tag_neutrino_mec_hists.lower() )

    caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Empirical.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
    caption.write( text_caption + 'Shown are the 2018 CV MEC and GENIE Empirical MEC.' )
    caption.close()

    legend2 = root_style.GetLegend()
    legend2.SetTextSize( text_size )
    legend2.SetTextFont( text_font )
    #legend2.SetHeader( '%s CC MEC' % text_mec_neutrino )
    legend2.AddEntry( hist_1D_tuned   , text_nova_2p2h   , 'L' )
    legend2.AddEntry( hist_1D_valencia, text_valencia_mec, 'L' )
    legend2.SetX1( xpos_legend_mec )
    legend2.SetY2( ypos_legend_mec )
    root_style.SetLegendPosY1( legend2 )

    canvas.Clear()

    hist_1D_tuned.SetMaximum( 1.8 * max( max_valencia, max_tuned ) )
    hist_1D_tuned.Draw( 'hist' )
    hist_1D_valencia.Draw( 'hist same' )
    hist_1D_tuned.Draw( 'hist same' )
    legend2.Draw()
    if opts.sim : label_simulation.Draw()
    label_beam.Draw()
    label_numu_cc.Draw()
    root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Valencia' % ( tag_neutrino_mec_hists, var ) ) )

    caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Valencia.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
    caption.write( text_caption + 'Shown are the 2018 CV MEC and Valencia MEC.' )
    caption.close()

    legend3 = root_style.GetLegend()
    legend3.SetTextSize( text_size )
    legend3.SetTextFont( text_font )
    #legend3.SetHeader( '%s CC MEC' % text_mec_neutrino )
    legend3.AddEntry( hist_1D_tuned  , text_nova_2p2h  , 'L' )
    legend3.AddEntry( hist_1D_minerva, text_minerva_2p2h, 'L' )
    legend3.SetX1( xpos_legend_mec )
    legend3.SetY2( ypos_legend_mec )
    root_style.SetLegendPosY1( legend3 )

    canvas.Clear()

    hist_1D_tuned.SetMaximum( 1.8 * max( max_minerva, max_tuned ) )
    hist_1D_tuned.Draw( 'hist' )
    hist_1D_minerva.Draw( 'hist same' )
    hist_1D_tuned.Draw( 'hist same' )
    legend3.Draw()
    if opts.sim : label_simulation.Draw()
    label_beam.Draw()
    label_numu_cc.Draw()
    root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_MINERvA' % ( tag_neutrino_mec_hists, var ) ) )

    caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_MINERvA.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
    caption.write( text_caption + 'Shown are the 2018 CV MEC and the MINERvA MEC tune.' )
    caption.close()

    legend4 = root_style.GetLegend()
    legend4.SetTextSize( text_size )
    legend4.SetTextFont( text_font )
    #legend4.SetHeader( '%s CC MEC' % text_mec_neutrino )
    legend4.AddEntry( hist_1D_tuned   , text_nova_2p2h   , 'L' )
    legend4.AddEntry( hist_1D_minerva , text_minerva_2p2h, 'L' )
    legend4.AddEntry( hist_1D_valencia, text_valencia_mec, 'L' )
    legend4.SetX1( xpos_legend_mec )
    legend4.SetY2( ypos_legend_mec )
    root_style.SetLegendPosY1( legend4 )

    canvas.Clear()

    hist_1D_tuned.SetMaximum( 2.0 * max( max_valencia, max_minerva, max_tuned ) )
    hist_1D_tuned.Draw( 'hist' )
    hist_1D_minerva.Draw( 'hist same' )
    hist_1D_valencia.Draw( 'hist same' )
    hist_1D_tuned.Draw( 'hist same' )
    legend4.Draw()
    if opts.sim : label_simulation.Draw()
    label_beam.Draw()
    label_numu_cc.Draw()
    root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_All' % ( tag_neutrino_mec_hists, var ) ) )

    caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_All.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
    caption.write( text_caption + 'Shown are the 2018 CV MEC, Valencia MEC, and the MINERvA MEC tune.' )
    caption.close()

    legend5 = root_style.GetLegend()
    legend5.SetTextSize( text_size )
    legend5.SetTextFont( text_font )
    #legend5.SetHeader( '%s CC MEC' % text_mec_neutrino )
    legend5.AddEntry( hist_1D_tuned   ,   text_nova_2p2h    , 'L' )
    legend5.AddEntry( hist_1D_shift_down, text_nova_2p2h + ' -1 #sigma', 'L' )
    legend5.AddEntry( hist_1D_shift_up,   text_nova_2p2h + ' +1 #sigma', 'L' )
    #legend5.AddEntry( hist_1D_shift_down, '-1 #sigma', 'L' )
    #legend5.AddEntry( hist_1D_shift_up,   '+1 #sigma', 'L' )
    legend5.SetX1( xpos_legend_mec )
    legend5.SetY2( ypos_legend_mec )
    root_style.SetLegendPosY1( legend5 )

    canvas.Clear()

    hist_1D_tuned.SetMaximum( 2.0 * max( max_shift_up, max_shift_down, max_tuned ) )
    if var == 'TrueQ0' : hist_1D_tuned.GetXaxis().SetRangeUser( 0.0, 0.8 )
    hist_1D_tuned.Draw( 'hist' )
    hist_1D_shift_down.Draw( 'hist same' )
    hist_1D_shift_up.Draw( 'hist same' )
    hist_1D_tuned.Draw( 'hist same' )
    legend5.Draw()
    if opts.sim : label_simulation.Draw()
    label_beam.Draw()
    label_numu_cc.Draw()
    root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_Shifts' % ( tag_neutrino_mec_hists, var ) ) )

    caption = open( os.path.join( plots_dir, '%s_%s_MEC_Shifts.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
    caption.write( text_caption + 'Shown are the 2018 CV MEC and 2018 MEC (q_0,|q|) +/-1 sigma shifts.  The shifts affect both shape and normalization.' )
    caption.close()

    legend6 = root_style.GetLegend()
    legend6.SetTextSize( text_size )
    legend6.SetTextFont( text_font )
    legend6.AddEntry( hist_1D_non_mec           , text_nova_tune, 'L' )
    legend6.AddEntry( hist_1D_non_mec_shift_up  , 'QE-Like Shifts', 'L' )
    legend6.AddEntry( hist_1D_non_mec_shift_down, 'RES-Like Shifts', 'L' )
    legend6.SetX1( xpos_legend_no_mec )
    legend6.SetY2( ypos_legend_no_mec )
    root_style.SetLegendPosY1( legend6 )

    canvas.Clear()

    ymax_non_mec = 2.0
    if var == 'TrueQ0' : ymax_non_mec = 1.3

    hist_1D_non_mec.SetMaximum( ymax_non_mec * max( max_non_mec_shift_up, max_non_mec_shift_down, max_non_mec ) )
    hist_1D_non_mec.GetXaxis().SetDecimals( kTRUE )
    hist_1D_non_mec.Draw( 'hist' )
    hist_1D_non_mec_shift_down.Draw( 'hist same' )
    hist_1D_non_mec_shift_up.Draw( 'hist same' )
    hist_1D_non_mec.Draw( 'hist same' )
    legend6.Draw()
    if opts.sim : label_simulation.Draw()
    label_beam.Draw()
    label_no_mec.Draw()
    root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_NonMEC_Shifts' % ( opts.beam.upper(), var ) ) )

    caption = open( os.path.join( plots_dir, '%s_%s_NonMEC_Shifts.txt' % ( opts.beam.upper(), var ) ), 'w' )
    caption.write( 'The %s near detector %s distribution for events that pass the 2018 numu selection and are not true numu/numubar CC MEC.  ' % ( opts.beam.upper(), text_variable ) )
    caption.write( 'Shown are the 2018 CV and the "QE-Like" and "RES-Like" shifts used to fit the 2018 MEC (q_0,|q|) +/-1 sigma shifts.  ' )
    caption.write( 'See the 2018 cross-section tuning tech note (DocDB 27755) for more details.' )
    caption.close()


  SetRangeQ0Q3( hist_nominal )
  SetRangeQ0Q3( hist_tuned )

  event_scale_power_2d = 3
  event_scale_2d = 1.0 / pow( 10.0, event_scale_power_2d )
  text_events_2d = '10^{%d} Events' % event_scale_power_2d
  scale_hist_2d  = mc_pot_scale * event_scale_2d

  hist_nominal.Scale( scale_hist_2d )
  hist_tuned  .Scale( scale_hist_2d )

  hist_nominal.GetZaxis().SetTitle( text_events_2d )
  hist_tuned  .GetZaxis().SetTitle( text_events_2d )
  hist_weights.GetZaxis().SetTitle( 'Weight' )

  canvas.SetLeftMargin(0.13)
  canvas.SetRightMargin(0.17)
  #gROOT.ForceStyle()

  text_weights_neutrino = '#nu'
  if isRHC : text_weights_neutrino = '#bar{#nu}'

  text_weights_1 = "Weights Applied to"
  text_weights_2 = "%s Empirical MEC" % text_weights_neutrino
  text_empirical = "%s CC Empirical MEC" % text_mec_neutrino
  text_nova_2p2h = "%s CC %s 2p2h" % ( text_mec_neutrino, text_nova )

  xpos_label_2d        = canvas.GetLeftMargin() + 0.03 # text_size
  ypos_label_weights_1 = ypos_label_top
  ypos_label_weights_2 = ypos_label_top - text_size
  ypos_label_nova_2p2h = ypos_label_top - text_size
  ypos_label_empirical = ypos_label_top - text_size - 0.005
  #ypos_label_mec_2d    = ypos_label_top - text_size * 2 - 0.01
  #if not isRHC : ypos_label_mec_2d -= 0.005

  label_weights_1 = root_style.GetLabel( text_weights_1, xpos_label_2d, ypos_label_weights_1, text_size, kWhite, text_font, 13, text_angle, setNDC )
  label_weights_2 = root_style.GetLabel( text_weights_2, xpos_label_2d, ypos_label_weights_2, text_size, kWhite, text_font, 13, text_angle, setNDC )
  label_empirical = root_style.GetLabel( text_empirical, xpos_label_2d, ypos_label_empirical, text_size, kBlack, text_font, 13, text_angle, setNDC )
  label_nova_2p2h = root_style.GetLabel( text_nova_2p2h, xpos_label_2d, ypos_label_nova_2p2h, text_size, kBlack, text_font, 13, text_angle, setNDC )
  label_beam.SetX( xpos_label_2d )
  #label_mec .SetX( xpos_label_2d )
  #label_mec .SetY( ypos_label_mec_2d )

  for i in range( 1, hist_weights.GetNbinsX() + 1 ) :
    for j in range( 1, hist_weights.GetNbinsY() + 1 ) :
      if hist_weights.GetBinContent( i, j ) < 0 :
        hist_weights.SetBinContent( i, j, 0.0 )

  #hist_weights.GetZaxis().SetRangeUser( 0, 5.0 )
  hist_weights.GetZaxis().SetRangeUser( -1.0e-6, 5.0 )

  SetMECWeightsPalette()

  hist_weights.Draw( 'colz' )
  label_weights_1.Draw()
  label_weights_2.Draw()
  if opts.sim : label_simulation.Draw()
  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_Weights' % tag_neutrino_mec_weights ) )

  gStyle.SetPalette( kBird )

  text_neutrino = 'neutrino'
  if isRHC : text_neutrino = 'antineutrino'

  caption = open( os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_Weights.txt' % tag_neutrino_mec_weights ), 'w' )
  caption.write( 'The 2018 %s MEC CV weights in true energy transfer, q_0, vs. true three-momentum transfer magnitude, |q|, bins.  ' % text_neutrino )
  caption.write( 'Non-zero weights in the kinematically forbidden region, q_0 > |q|, are an artifact of smoothing the weights.' )
  caption.close()

  zaxis_tag = "_FixedZ"
  zmax = 10
  if isRHC : zmax = 3

  hist_nominal.Draw( 'colz' )
  label_beam.Draw()
  label_empirical.Draw()
  if opts.sim : label_simulation.Draw()
  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_Empirical" ) )

  hist_nominal.GetZaxis().SetRangeUser( 0, zmax )
  hist_nominal.Draw( 'colz' )
  label_beam.Draw()
  label_empirical.Draw()
  if opts.sim : label_simulation.Draw()
  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_Empirical" + zaxis_tag ) )

  text_caption  = 'The %s near detector true energy transfer, q_0, vs. true three-momentum transfer magnitude, |q|, ' % opts.beam.upper()
  text_caption += 'distribution for true %s CC MEC events that pass the 2018 numu selection.  ' % tag_neutrino_mec_hists.lower()

  caption = open( os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_Empirical.txt' % tag_neutrino_mec_hists ), 'w' )
  caption.write( text_caption + 'Shown is the GENIE Empirical MEC.' )
  caption.close()
 
  hist_tuned.Draw( 'colz' )
  label_beam.Draw()
  label_nova_2p2h.Draw()
  if opts.sim : label_simulation.Draw()
  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_2018CV" ) )

  hist_tuned.GetZaxis().SetRangeUser( 0, zmax )
  hist_tuned.Draw( 'colz' )
  label_beam.Draw()
  label_nova_2p2h.Draw()
  if opts.sim : label_simulation.Draw()
  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_2018CV" + zaxis_tag ) )

  caption = open( os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_2018CV.txt' % tag_neutrino_mec_hists ), 'w' )
  caption.write( text_caption + 'Shown is the 2018 CV MEC.' )
  caption.close()