MECModelEnuComparisons.py
Go to the documentation of this file.
1 """
2  MECModelEnuComparisons.py:
3  Comparisons of the Enu distributions for various MEC models.
4  Origin of plots used in Fig. 14 of Doc 23264 and in the xsec tuning paper.
5 
6  Original author: J. Wolcott <jwolcott@fnal.gov>
7  Date: April 2017 (modified for paper, Oct 2019)
8 
9 """
10 
11 import os.path
12 import sys
13 
14 import ROOT
15 
16 COLORS = {
17  "DefaultPlusMECWithNC_MEC": ROOT.kOrange-1,
18  "LocalFGNievesQEAndMEC_MEC": ROOT.kRed,
19  "Martini": ROOT.kBlue,
20  "SuSA": ROOT.kViolet,
21 }
22 
23 NAMES = {
24  "DefaultPlusMECWithNC_MEC": "Empirical MEC",
25  "LocalFGNievesQEAndMEC_MEC": "Val#grave{e}ncia MEC", ##"Nieves et al. MEC (GENIE)",
26  "Martini": "Martini MEC", #"Martini et al. MEC (PRC 80, 065501)",
27  "SuSA": "SuSA MEC", #"Megias et al. MEC (PRD 94, 093004)"
28 }
29 
30 DRAW_ORDER = ["DefaultPlusMECWithNC_MEC", "LocalFGNievesQEAndMEC_MEC", "Martini", "SuSA"]
31 
32 OUT_DIR = "/nova/ana/users/jwolcott/2p2h_tuning/paper"
33 OUT_TYPES = ("png", "pdf")
34 
35 ###############
36 
37 def SaveCanvas(c, name_stub, out_types=OUT_TYPES):
38  for out_type in out_types:
39  c.SaveAs("%s.%s" % (name_stub, out_type))
40 
41 ###############
42 
43 SPLINE_FILES = {
44  "DefaultPlusMECWithNC_MEC": "/cvmfs/nova.opensciencegrid.org/externals/genie_xsec/v2_12_0/NULL/DefaultPlusMECWithNC/data/xsec_graphs.root",
45  "LocalFGNievesQEAndMEC_MEC": "/cvmfs/nova.opensciencegrid.org/externals/genie_xsec/v2_12_0/NULL/LocalFGNievesQEAndMEC/data/xsec_graphs.root",
46 
47  # these two were digitized from the plots in the papers
48  "Martini": "/nova/ana/xsec_tuning/2019_paper/Martini_MEC_sigma.root",
49  "SuSA": "/nova/ana/xsec_tuning/2019_paper/SuSA_MEC_sigma.root"
50 }
51 
52 graphs = {}
53 for model, fname in SPLINE_FILES.iteritems():
54  f = ROOT.TFile(fname)
55  g = graphs[model] = f.Get("nu_mu_C12/mec_cc")
56  if not g:
57  graphs[model] = f.Get("xsec_graph")
58 
59 c = ROOT.TCanvas()
60 leg = ROOT.TLegend(0.35, 0.65, 0.85, 0.9)
61 leg.SetFillStyle(0)
62 leg.SetTextFont(42)
63 i = 0
64 for gname in DRAW_ORDER:
65  if gname not in graphs:
66  continue
67  g = graphs[gname]
68  g.SetLineColor(COLORS[gname])
69  g.SetLineWidth(2)
70  leg.AddEntry(g, NAMES[gname], "l")
71  opt = "l"
72  if i == 0:
73  opt = "a" + opt
74  g.Draw(opt)
75 
76  i += 1
77 leg.Draw()
78 
79 g = graphs["DefaultPlusMECWithNC_MEC"]
80 g.GetXaxis().SetRangeUser(0, 10)
81 g.GetYaxis().SetRangeUser(0, 3)
82 g.SetTitle(";E_{#nu} (GeV);Total cross section (10^{-39} cm^{2})")
83 SaveCanvas(c, os.path.join(OUT_DIR, "sigma_comparison"))
84 
85 # make them max out at the same place
86 x = ROOT.Double()
87 y = ROOT.Double()
88 EmpMECmax = 0
89 for pt_num in range(graphs["DefaultPlusMECWithNC_MEC"].GetN()):
90  graphs["DefaultPlusMECWithNC_MEC"].GetPoint(pt_num, x, y)
91  if y > EmpMECmax:
92  EmpMECmax = float(y)
93 
94 EmpMEC3GeV = graphs["DefaultPlusMECWithNC_MEC"].Eval(3)
95 EmpMEC10GeV = graphs["DefaultPlusMECWithNC_MEC"].Eval(10)
96 
97 altgraphs = {}
98 for gname in DRAW_ORDER:
99  if gname not in graphs:
100  continue
101  g = graphs[gname]
102 
103  g.GetPoint(g.GetN() - 1, x, y)
104  # print gname, x, y
105  if x > 10:
106  norm = EmpMEC10GeV / g.Eval(10)
107  elif x > 2:
108  norm = EmpMEC3GeV / g.Eval(3)
109  else:
110  maximum = 0
111  for pt_num in range(g.GetN()):
112  g.GetPoint(pt_num, x, y)
113  if y > maximum:
114  maximum = y
115  norm = EmpMECmax / maximum
116  gp = altgraphs[gname] = ROOT.TGraph()
117  for pt_num in range(g.GetN()):
118  g.GetPoint(pt_num, x, y)
119  gp.SetPoint(pt_num, x, y * norm)
120 
121 c.Clear()
122 altmp = ROOT.TMultiGraph()
123 for gname, g in altgraphs.iteritems():
124  g.SetLineColor(COLORS[gname])
125  g.SetLineWidth(2)
126  altmp.Add(g)
127 altmp.Draw("al")
128 altmp.GetXaxis().SetRangeUser(0, 10)
129 altmp.GetXaxis().CenterTitle()
130 altmp.GetYaxis().SetRangeUser(0, 5)
131 altmp.GetYaxis().CenterTitle()
132 altmp.GetXaxis().SetTitle("E_{#nu} (GeV)")
133 altmp.GetYaxis().SetTitle("Renormalized #sigma_{MEC} (10^{-39} cm^{2})")
134 leg.Draw()
135 SaveCanvas(c, os.path.join(OUT_DIR, "sigma_comparison_renormalized"))
136 
137 rat_graphs = {}
138 for gname, g in altgraphs.iteritems():
139  if gname == "DefaultPlusMECWithNC_MEC":
140  continue
141  gp = rat_graphs[gname] = ROOT.TGraph()
142  for pt_num in range(g.GetN()):
143  g.GetPoint(pt_num, x, y)
144  EmpMECVal = graphs["DefaultPlusMECWithNC_MEC"].Eval(x)
145  gp.SetPoint(pt_num, x, y / EmpMECVal if EmpMECVal > 0 else 0)
146 
147 c.Clear()
148 ratmp = ROOT.TMultiGraph()
149 for gname, g in rat_graphs.iteritems():
150  g.SetLineColor(COLORS[gname])
151  g.SetLineWidth(2)
152  ratmp.Add(g)
153 ratmp.Draw("al")
154 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")}
155 for fn in fns.itervalues():
156  fn.SetLineColor(ROOT.kGreen + 2)
157  fn.SetLineStyle(2)
158  fn.Draw("same")
159 ratmp.GetXaxis().SetRangeUser(0, 10)
160 ratmp.GetXaxis().CenterTitle()
161 ratmp.GetYaxis().SetRangeUser(0, 2.5)
162 ratmp.GetXaxis().SetTitle("E_{#nu} (GeV)")
163 ratmp.GetYaxis().SetTitle("Shape ratio to Empirical MEC")
164 ratmp.GetYaxis().CenterTitle()
165 leg.GetListOfPrimitives().RemoveFirst()
166 leg.AddEntry(fns["pos"], "Uncertainty envelope", "l")
167 leg.Draw()
168 SaveCanvas(c, os.path.join(OUT_DIR, "sigma_comparison_renormalized_ratio"))
169 
170 
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
def SaveCanvas(c, name_stub, out_types=OUT_TYPES)