MECModelEnuComparisons.py

Go to the documentation of this file.

"""
  MECModelEnuComparisons.py:
    Comparisons of the Enu distributions for various MEC models.
    Origin of plots used in Fig. 14 of Doc 23264 and in the xsec tuning paper.

  Original author: J. Wolcott <jwolcott@fnal.gov>
             Date: April 2017 (modified for paper, Oct 2019)

"""

import os.path
import sys

import ROOT

COLORS = {
    "DefaultPlusMECWithNC_MEC": ROOT.kOrange-1,
    "LocalFGNievesQEAndMEC_MEC": ROOT.kRed,
    "Martini":                   ROOT.kBlue,
    "SuSA":                      ROOT.kViolet,
}

NAMES = {
    "DefaultPlusMECWithNC_MEC": "Empirical MEC",
    "LocalFGNievesQEAndMEC_MEC": "Val#grave{e}ncia MEC", ##"Nieves et al. MEC (GENIE)",
    "Martini":                   "Lyon MEC", #"Martini and Ericson MEC (PRC 80, 065501)",
    "SuSA":                      "SuSA MEC", #"Megias et al. MEC (PRD 94, 093004)"
}

DRAW_ORDER = ["DefaultPlusMECWithNC_MEC", "LocalFGNievesQEAndMEC_MEC", "Martini", "SuSA"]

OUT_DIR = "/nova/ana/users/jwolcott/2p2h_tuning/paper"
OUT_TYPES = ("png", "pdf")

###############

def SaveCanvas(c, name_stub, out_types=OUT_TYPES):
        for out_type in out_types:
                c.SaveAs("%s.%s" % (name_stub, out_type))

###############

SPLINE_FILES = {
        "DefaultPlusMECWithNC_MEC": "/cvmfs/nova.opensciencegrid.org/externals/genie_xsec/v2_12_0/NULL/DefaultPlusMECWithNC/data/xsec_graphs.root",
        "LocalFGNievesQEAndMEC_MEC": "/cvmfs/nova.opensciencegrid.org/externals/genie_xsec/v2_12_0/NULL/LocalFGNievesQEAndMEC/data/xsec_graphs.root",

        # these two were digitized from the plots in the papers
        "Martini": "/nova/ana/xsec_tuning/2019_paper/Martini_MEC_sigma.root",
        "SuSA": "/nova/ana/xsec_tuning/2019_paper/SuSA_MEC_sigma.root"
}

graphs = {}
for model, fname in SPLINE_FILES.iteritems():
        f = ROOT.TFile(fname)
        g = graphs[model] = f.Get("nu_mu_C12/mec_cc")
        if not g:
                graphs[model] = f.Get("xsec_graph")

c = ROOT.TCanvas()
leg = ROOT.TLegend(0.35, 0.65, 0.85, 0.9)
leg.SetFillStyle(0)
leg.SetTextFont(42)
i = 0
for gname in DRAW_ORDER:
        if gname not in graphs:
                continue
        g = graphs[gname]
        g.SetLineColor(COLORS[gname])
        g.SetLineWidth(2)
        leg.AddEntry(g, NAMES[gname], "l")
        opt = "l"
        if i == 0:
                opt = "a" + opt
        g.Draw(opt)

        i += 1
leg.Draw()

g = graphs["DefaultPlusMECWithNC_MEC"]
g.GetXaxis().SetRangeUser(0, 10)
g.GetYaxis().SetRangeUser(0, 3)
g.SetTitle(";E_{#nu} (GeV);Total cross section (10^{-39} cm^{2})")
SaveCanvas(c, os.path.join(OUT_DIR, "sigma_comparison"))

# make them max out at the same place
x = ROOT.Double()
y = ROOT.Double()
EmpMECmax = 0
for pt_num in range(graphs["DefaultPlusMECWithNC_MEC"].GetN()):
        graphs["DefaultPlusMECWithNC_MEC"].GetPoint(pt_num, x, y)
        if y > EmpMECmax:
                EmpMECmax = float(y)

EmpMEC3GeV = graphs["DefaultPlusMECWithNC_MEC"].Eval(3)
EmpMEC10GeV = graphs["DefaultPlusMECWithNC_MEC"].Eval(10)

altgraphs = {}
for gname in DRAW_ORDER:
        if gname not in graphs:
                continue
        g = graphs[gname]

        g.GetPoint(g.GetN() - 1, x, y)
        #    print gname, x, y
        if x > 10:
                norm = EmpMEC10GeV / g.Eval(10)
        elif x > 2:
                norm = EmpMEC3GeV / g.Eval(3)
        else:
                maximum = 0
                for pt_num in range(g.GetN()):
                        g.GetPoint(pt_num, x, y)
                        if y > maximum:
                                maximum = y
                norm = EmpMECmax / maximum
        gp = altgraphs[gname] = ROOT.TGraph()
        for pt_num in range(g.GetN()):
                g.GetPoint(pt_num, x, y)
                gp.SetPoint(pt_num, x, y * norm)

c.Clear()
altmp = ROOT.TMultiGraph()
for gname, g in altgraphs.iteritems():
        g.SetLineColor(COLORS[gname])
        g.SetLineWidth(2)
        altmp.Add(g)
altmp.Draw("al")
altmp.GetXaxis().SetRangeUser(0, 10)
altmp.GetXaxis().CenterTitle()
altmp.GetYaxis().SetRangeUser(0, 5)
altmp.GetYaxis().CenterTitle()
altmp.GetXaxis().SetTitle("E_{#nu} (GeV)")
altmp.GetYaxis().SetTitle("Renormalized #sigma_{MEC} (10^{-39} cm^{2})")
altmp.GetYaxis().SetDecimals()
leg.Draw()
SaveCanvas(c, os.path.join(OUT_DIR, "sigma_comparison_renormalized"))

rat_graphs = {}
for gname, g in altgraphs.iteritems():
        if gname == "DefaultPlusMECWithNC_MEC":
                continue
        gp = rat_graphs[gname] = ROOT.TGraph()
        for pt_num in range(g.GetN()):
                g.GetPoint(pt_num, x, y)
                EmpMECVal = graphs["DefaultPlusMECWithNC_MEC"].Eval(x)
                gp.SetPoint(pt_num, x, y / EmpMECVal if EmpMECVal > 0 else 0)

c.Clear()
ratmp = ROOT.TMultiGraph()
for gname, g in rat_graphs.iteritems():
        g.SetLineColor(COLORS[gname])
        g.SetLineWidth(2)
        ratmp.Add(g)
ratmp.Draw("al")
fns = {direction: ROOT.TF1("fn_%s" % direction, "1 %s 1./(1.+2.5*x)" % ("+" if direction == "pos" else "-"), ratmp.GetXaxis().GetXmin(), ratmp.GetXaxis().GetXmax()) for direction in ("pos", "neg")}
for fn in fns.itervalues():
        fn.SetLineColor(ROOT.kGreen + 2)
        fn.SetLineStyle(2)
        fn.Draw("same")
ratmp.GetXaxis().SetRangeUser(0, 10)
ratmp.GetXaxis().CenterTitle()
ratmp.GetYaxis().SetRangeUser(0, 2.5)
ratmp.GetXaxis().SetTitle("E_{#nu} (GeV)")
ratmp.GetYaxis().SetTitle("Shape ratio to Empirical MEC")
ratmp.GetYaxis().CenterTitle()
ratmp.GetYaxis().SetDecimals()
leg.GetListOfPrimitives().RemoveFirst()
leg.AddEntry(fns["pos"], "Uncertainty envelope", "l")
leg.Draw()
SaveCanvas(c, os.path.join(OUT_DIR, "sigma_comparison_renormalized_ratio"))