EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
charm_jet_cross_section.py
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file charm_jet_cross_section.py
1 import matplotlib as mpl
2 import uproot
3 import matplotlib.pyplot as plt
4 import scipy
5 import numpy as np
6 import math
7 import pandas as pd
8 import seaborn as sns
9 import mplhep as hep
10 #import zfit
11 import inspect
12 import sys
13 
14 from concurrent.futures import ThreadPoolExecutor
15 
16 plt.style.use(hep.style.ATLAS)
17 
18 plt.rcParams.update({'font.sans-serif': "Arial",
19  'font.family': "sans-serif",
20  'font.size': 30,
21  'mathtext.fontset': 'custom',
22  'mathtext.rm': 'Arial',
23  })
24 
25 import EICAnalysisTools as eat
26 
27 
28 # Units
29 global u_fb, u_mb, xsection
30 u_fb = 1
31 u_mb = 1e12*u_fb
32 #xsection = 2.408e-08*u_mb
33 
34 ## dataset dictionary (cross-sections, etc.)
35 #global xs_map
36 #xs_map = {}
37 #xs_map['CC_DIS_e18_p275_B15'] = 2.408e-08*u_mb # from PYTHIA8
38 #xs_map['NC_DIS_e18_p275_B15'] = 3.671e-06*u_mb # from PYTHIA8
39 
40 
41 print("Loading data...")
42 
43 branchlist = ["Jet.PT", "Jet.Eta", "Jet.Flavor", "Jet.BTag", "GenJet.PT", "GenJet.Eta", "GenJet.Flavor"]
44 
45 
46 df_ct18nnlo = eat.UprootLoad([f"../CC_DIS_e10_p275_B15_dR5_maxIP3mm_trkpt10_22sigmin_CT18NNLO/*/out.root"], "Delphes",
47  branches=branchlist)
48 #df_ct18nnlo_20rs2 = eat.UprootLoad([f"../CC_DIS_e10_p275_B15_dR5_maxIP3mm_trkpt10_22sigmin_lha_20Rs2/*/out.root"], "Delphes",
49 # branches=branchlist)
50 #df_ct18nnlo_21rs2 = eat.UprootLoad([f"../CC_DIS_e10_p275_B15_dR5_maxIP3mm_trkpt10_22sigmin_lha_21Rs2/*/out.root"], "Delphes",
51 # branches=branchlist)
52 
53 
54 #df_ct18nnlo = df_ct18nnlo[:1000]
55 #df_ct18nnlo_20rs2 = df_ct18nnlo_20rs2[:1000]
56 #df_ct18nnlo_21rs2 = df_ct18nnlo_21rs2[:1000]
57 
58 
59 #df_ct18nlo = UprootLoad(f"../CC_DIS_e18_p275_B15_dR5_maxIP3mm_CT18NLO/*/out.root")
60 #df_ct18Annlo = UprootLoad(f"../CC_DIS_e18_p275_B15_dR5_maxIP3mm_CT18ANNLO/*/out.root")
61 #df_ct18Anlo = UprootLoad(f"../CC_DIS_e18_p275_B15_dR5_maxIP3mm_CT18ANLO/*/out.root")
62 
63 #df_ct18nnlo_nc = UprootLoad(f"../NC_DIS_e18_p275_B15_dR5_maxIP3mm_CT18NNLO/*/out.root")
64 
65 
66 print("Done loading data for the study!")
67 
68 
69 
70 def DifferentialXS(df, x='GenJet.PT', xrange=[10,50], xbins=[10,12.5,15,20,25,30,40,50], which='all', process='CC_DIS_e18_p275'):
71  global u_fb, u_mb, n_gen
72  n_gen = len(df)
73  print(f"n_gen = {n_gen}")
74  jet_pt = np.concatenate(df['GenJet.PT'].to_numpy()).ravel()
75  jet_eta = np.concatenate(df['GenJet.Eta'].to_numpy()).ravel()
76  jet_flavor = np.concatenate(df['GenJet.Flavor'].to_numpy()).ravel()
77 
78  jet_x = np.concatenate(df[x].to_numpy()).ravel()
79 
80 
81  #jet_basics = (0.01 < jet_eta) & (jet_eta < 0.9)
82  jet_basics = (-5.0 < jet_eta) & (jet_eta < 5.0)
83 
84  all_flavor = ( jet_flavor > -999.0 ) & (jet_basics)
85  light_flavor = ( jet_flavor < 4 ) & (jet_basics)
86  charm_flavor = ( jet_flavor == 4 ) & (jet_basics)
87 
88  selection = all_flavor
89  if which == 'charm':
90  selection = charm_flavor
91 
92 
93  (counts, bins) = np.histogram(jet_x[ selection ], range=xrange,
94  bins=xbins)
95  #bins=40)
96 
97  bin_widths = np.diff(bins)
98  bin_centers = bins[:-1] + bin_widths/2
99 
100  errors = np.sqrt(counts)
101  rel_errors = errors/counts
102 
103 
104  # convert counts to dsigma/dpT * 100/fb
105  dsigma_dx = counts * eat.TotalXSection(process) * 100*u_fb**(-1) / n_gen / bin_widths
106  dsigma_dx_errors = rel_errors * dsigma_dx
107 
108  return (bin_centers, bin_widths, dsigma_dx, dsigma_dx_errors)
109 
110 
111 
112 
113 
114 draw_config_eta={}
115 
116 draw_config_eta['xvar'] = 'GenJet.Eta'
117 draw_config_eta['xrange'] = [-4.0, 4.0]
118 draw_config_eta['xbins'] = [-4.0, -3.5, -2.5, -1.0, 1.0, 2.5, 3.5, 4.0]
119 draw_config_eta['ylimits'] = [10, 1e7]
120 draw_config_eta['xlimits'] = [-5,5]
121 draw_config_eta['yunits'] = ''
122 draw_config_eta['xunits'] = ''
123 
124 draw_config_pt={}
125 draw_config_pt['xvar'] = 'GenJet.PT'
126 draw_config_pt['xrange'] = [10,50]
127 draw_config_pt['xbins'] = [10,12.5,15,20,25,30,40,50]
128 draw_config_pt['ylimits'] = [10, 2e5]
129 draw_config_pt['xlimits'] = [0,60]
130 draw_config_pt['yunits'] = '[$\\mathrm{GeV^{-1}}$]'
131 draw_config_pt['xunits'] = '[GeV]'
132 
133 
134 def DrawDiffXSPlot(draw_config={}, process='CC_DIS_e18_p275'):
135  xvar=draw_config['xvar']
136  xrange=draw_config['xrange']
137  xbins=draw_config['xbins']
138  ylimits=draw_config['ylimits']
139  xlimits=draw_config['xlimits']
140  yunits=draw_config['yunits']
141  xunits=draw_config['xunits']
142 
143  fig, ax = plt.subplots(figsize=(12,12))
144  plt.axis('off')
145 
146 
147 
148  if process.find('CC') != -1:
149  gridspec = fig.add_gridspec(ncols=1, nrows=1, width_ratios=[1], height_ratios=[1])
150 
151 
152  # Log plot of differential cross-sections
153  ax1 = fig.add_subplot(gridspec[0, 0])
154 
155  ax1.grid(which='both', axis='both')
156  (bins, bin_widths, xs, xs_error) = DifferentialXS(df_ct18nnlo, x=xvar, xrange=xrange, xbins=xbins, which='all', process='CC_DIS_e10_p275_CT18NNLO')
157  ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='all jets (CT18NNLO)', marker='o', ms=10, ls='none', linewidth=2, color='red')
158 
159  #(bins, bin_widths, xs, xs_error) = DifferentialXS(df_ct18nnlo_20rs2, x=xvar, xrange=xrange, xbins=xbins, which='all', process='CC_DIS_e10_p275_CT1820Rs2')
160  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='all jets (CT18NNLO 20Rs2)', marker='o', ms=10, ls='none', linewidth=2, color='saddlebrown')
161 
162  #(bins, bin_widths, xs, xs_error) = DifferentialXS(df_ct18nnlo_21rs2, x=xvar, xrange=xrange, xbins=xbins, which='all', process='CC_DIS_e10_p275_CT1821Rs2')
163  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='all jets (CT18NNLO 21Rs2)', marker='o', ms=10, ls='none', linewidth=2, color='magenta')
164 
165  #(bins, bin_widths, xs, xs_error) = DifferentialXS(df_ct18Annlo, x=xvar, xrange=xrange, xbins=xbins, which='all')
166  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='all jets (CT18ANNLO)', marker='s', ms=10, ls='none', linewidth=2, color='saddlebrown')
167 
168  #(bins, bin_widths, xs, xs_error) = DifferentialXS(df_ct18nlo, x=xvar, xrange=xrange, xbins=xbins, which='all')
169  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='all jets (CT18NLO)', marker='v', ms=10, ls='none', linewidth=2, color='magenta')
170 
171  #(bins, bin_widths, xs, xs_error) = DifferentialXS(df_ct18Anlo, x=xvar, xrange=xrange, xbins=xbins, which='all')
172  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='all jets (CT18ANLO)', marker='D', ms=10, ls='none', linewidth=2, color='coral')
173 
174  charm_ct18nnlo = DifferentialXS(df_ct18nnlo, x=xvar, xrange=xrange, xbins=xbins, which='charm')
175  (bins, bin_widths, xs, xs_error) = charm_ct18nnlo
176  ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18NNLO)', marker='D', ms=10, ls='none', linewidth=2, fillstyle='none', color='blue')
177 
178  #charm_ct18nnlo_20rs2 = DifferentialXS(df_ct18nnlo_20rs2, x=xvar, xrange=xrange, xbins=xbins, which='charm')
179  #(bins, bin_widths, xs, xs_error) = charm_ct18nnlo_20rs2
180  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18NNLO 20Rs2)', marker='D', ms=10, ls='none', linewidth=2, fillstyle='none', color='saddlebrown')
181 
182  #charm_ct18nnlo_21rs2 = DifferentialXS(df_ct18nnlo_21rs2, x=xvar, xrange=xrange, xbins=xbins, which='charm')
183  #(bins, bin_widths, xs, xs_error) = charm_ct18nnlo_21rs2
184  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18NNLO 21rs2)', marker='D', ms=10, ls='none', linewidth=2, fillstyle='none', color='magenta')
185 
186  #charm_ct18Annlo = DifferentialXS(df_ct18Annlo, x=xvar, xrange=xrange, xbins=xbins, which='charm')
187  #(bins, bin_widths, xs, xs_error) = charm_ct18Annlo
188  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18ANNLO)', marker='s', ms=10, ls='none', linewidth=2, fillstyle='none', color='saddlebrown')
189 
190  #charm_ct18nlo = DifferentialXS(df_ct18nlo, x=xvar, xrange=xrange, xbins=xbins, which='charm')
191  #(bins, bin_widths, xs, xs_error) = charm_ct18nlo
192  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18NLO)', marker='v', ms=10, ls='none', linewidth=2, fillstyle='none', color='magenta')
193 
194  #charm_ct18Anlo = DifferentialXS(df_ct18Anlo, x=xvar, xrange=xrange, xbins=xbins, which='charm')
195  #(bins, bin_widths, xs, xs_error) = charm_ct18Anlo
196  #ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18ANLO)', marker='D', ms=10, ls='none', linewidth=2, fillstyle='none', color='coral')
197 
198 
199  xvar_symbol="p_T"
200  if xvar.find('Eta') != -1:
201  xvar_symbol='\\eta'
202 
203  plt.ylabel(f'$d\\sigma/d{xvar_symbol} \\times 100\\mathrm{{fb^{{-1}}}}$ {yunits}')
204  plt.xlabel(f'Generated Jet ${xvar_symbol}$ [GeV]')
205  plt.yscale('log')
206 
207  plt.title("CC-DIS, 10x275GeV, $Q^2>100\\mathrm{GeV^2}$", fontsize=20)
208 
209  ax1.set_ylim(ylimits)
210  ax1.set_xlim(xlimits)
211 
212  ax1.legend(fontsize=18)
213 
214 
215  # ratio plot for charm production
216  #ax2 = fig.add_subplot(gridspec[1, 0])
217 
218  #ax2.grid(which='both', axis='both')
219 
220  #ratio_charm_ct18nnlo_ct18nnlo = charm_ct18nnlo[2]/charm_ct18nnlo[2]
221  #ratio_charm_ct18nnlo20rs2_ct18nnlo = charm_ct18nnlo_20rs2[2]/charm_ct18nnlo[2]
222  #ratio_charm_ct18nnlo21rs2_ct18nnlo = charm_ct18nnlo_21rs2[2]/charm_ct18nnlo[2]
223  #ratio_charm_ct18Annlo_ct18nnlo = charm_ct18Annlo[2]/charm_ct18nnlo[2]
224  #ratio_charm_ct18nlo_ct18nnlo = charm_ct18nlo[2]/charm_ct18nnlo[2]
225  #ratio_charm_ct18Anlo_ct18nnlo = charm_ct18Anlo[2]/charm_ct18nnlo[2]
226 
227  #ratioerr_charm_ct18nnlo_ct18nnlo = ratio_charm_ct18nnlo_ct18nnlo*np.sqrt((charm_ct18nnlo[3]/charm_ct18nnlo[2])**2)
228  #ratioerr_charm_ct18nnlo20rs2_ct18nnlo = ratio_charm_ct18nnlo20rs2_ct18nnlo*np.sqrt((charm_ct18nnlo_20rs2[3]/charm_ct18nnlo_20rs2[2])**2)
229  #ratioerr_charm_ct18nnlo21rs2_ct18nnlo = ratio_charm_ct18nnlo21rs2_ct18nnlo*np.sqrt((charm_ct18nnlo_21rs2[3]/charm_ct18nnlo_21rs2[2])**2)
230  #ratioerr_charm_ct18Annlo_ct18nnlo = ratio_charm_ct18Annlo_ct18nnlo*np.sqrt((charm_ct18Annlo[3]/charm_ct18Annlo[2])**2 + (charm_ct18nnlo[3]/charm_ct18nnlo[2])**2)
231  #ratioerr_charm_ct18nlo_ct18nnlo = ratio_charm_ct18nlo_ct18nnlo*np.sqrt((charm_ct18nlo[3]/charm_ct18nlo[2])**2 + (charm_ct18nnlo[3]/charm_ct18nnlo[2])**2)
232  #ratioerr_charm_ct18Anlo_ct18nnlo = ratio_charm_ct18Anlo_ct18nnlo*np.sqrt((charm_ct18Anlo[3]/charm_ct18Anlo[2])**2 + (charm_ct18nnlo[3]/charm_ct18nnlo[2])**2)
233 
234  #print(ratio_charm_ct18nnlo20rs2_ct18nnlo)
235  #print(ratio_charm_ct18nnlo21rs2_ct18nnlo)
236  #print(ratio_charm_ct18Annlo_ct18nnlo)
237  #print(ratio_charm_ct18nlo_ct18nnlo)
238  #print(ratio_charm_ct18Anlo_ct18nnlo)
239 
240  #ax2.axhline(1.0, color='k', linestyle='-', linewidth=2)
241  #ax2.errorbar(bins, ratio_charm_ct18nnlo_ct18nnlo, xerr = bin_widths/2, yerr=ratioerr_charm_ct18nnlo_ct18nnlo, marker='s', ms=10, ls='none', linewidth=2, fillstyle='none', color='saddlebrown')
242  #errorboxes = [mpl.patches.Rectangle((x - xe, y - ye), 2*xe, 2*ye)
243  # for x, y, xe, ye in zip(bins, ratio_charm_ct18nnlo_ct18nnlo, (bin_widths/2), ratioerr_charm_ct18nnlo_ct18nnlo)]
244 
245  # Create patch collection with specified colour/alpha
246  #pc = mpl.collections.PatchCollection(errorboxes, facecolor='k', alpha=0.25)
247 
248  # Add collection to axes
249  #ax2.add_collection(pc)
250 
251  #ax2.errorbar(bins, ratio_charm_ct18nnlo20rs2_ct18nnlo, xerr = bin_widths/2, marker='s', ms=10, ls='none', linewidth=2, fillstyle='none', color='saddlebrown')
252  #ax2.errorbar(bins, ratio_charm_ct18nnlo21rs2_ct18nnlo, xerr = bin_widths/2, marker='s', ms=10, ls='none', linewidth=2, fillstyle='none', color='magenta')
253  #ax2.errorbar(bins, ratio_charm_ct18Annlo_ct18nnlo, xerr = bin_widths/2, marker='s', ms=10, ls='none', linewidth=2, fillstyle='none', color='saddlebrown')
254  #ax2.errorbar(bins, ratio_charm_ct18nlo_ct18nnlo, xerr = bin_widths/2, marker='v', ms=10, ls='none', linewidth=2, fillstyle='none', color='magenta')
255  #ax2.errorbar(bins, ratio_charm_ct18Anlo_ct18nnlo, xerr = bin_widths/2, marker='D', ms=10, ls='none', linewidth=2, fillstyle='none', color='coral')
256 
257  #ax2.set_ylim([0.25, 1.75])
258  #ax2.set_xlim(xlimits)
259  #plt.ylabel('Ratio to\nCT18NNLO', fontsize=18)
260  elif process.find('NC') != -1:
261  gridspec = fig.add_gridspec(ncols=1, nrows=1, width_ratios=[1], height_ratios=[1])
262 
263 
264  # Log plot of differential cross-sections
265  ax1 = fig.add_subplot(gridspec[0, 0])
266 
267  ax1.grid(which='both', axis='both')
268  (bins, bin_widths, xs, xs_error) = DifferentialXS(df_ct18nnlo_nc, x=xvar, xrange=xrange, xbins=xbins, which='all', process=process)
269  ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='all jets (CT18NNLO)', marker='o', ms=10, ls='none', linewidth=2, color='red')
270 
271  charm_ct18nnlo = DifferentialXS(df_ct18nnlo_nc, x=xvar, xrange=xrange, xbins=xbins, which='charm', process=process)
272  (bins, bin_widths, xs, xs_error) = charm_ct18nnlo
273  ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18NNLO)', marker='D', ms=10, ls='none', linewidth=2, fillstyle='none', color='blue')
274 
275  xvar_symbol="p_T"
276  if xvar.find('Eta') != -1:
277  xvar_symbol='\\eta'
278 
279  plt.ylabel(f'$d\\sigma/d{xvar_symbol} \\times 100\\mathrm{{fb^{{-1}}}}$ {yunits}')
280  plt.xlabel(f'Generated Jet ${xvar_symbol}$ [GeV]')
281  plt.yscale('log')
282 
283  plt.title("NC-DIS, 18x275GeV, $Q^2>100\\mathrm{GeV^2}$", fontsize=20)
284 
285  ax1.set_ylim(ylimits)
286  ax1.set_xlim(xlimits)
287 
288  ax1.legend(fontsize=18)
289 
290 
291  plt.tight_layout()
292 
293  plt.savefig(f"jet_differential_xs_{xvar}_{process}.png")
294  plt.savefig(f"jet_differential_xs_{xvar}_{process}.pdf")
295 
296 ### Charged Current Scattering
297 DrawDiffXSPlot(draw_config_pt, process='CC_DIS_e10_p275')
298 DrawDiffXSPlot(draw_config_eta, process='CC_DIS_e10_p275')
299 
300 sys.exit()
301 
302 ### Neutral Current Scattering
303 draw_config_pt={}
304 draw_config_pt['xvar'] = 'GenJet.PT'
305 draw_config_pt['xrange'] = [10,50]
306 draw_config_pt['xbins'] = [10,12.5,15,20,25,30,40,50]
307 draw_config_pt['ylimits'] = [1e2, 5e8]
308 draw_config_pt['xlimits'] = [0,60]
309 draw_config_pt['yunits'] = '[$\\mathrm{GeV^{-1}}$]'
310 draw_config_pt['xunits'] = '[GeV]'
311 #DrawDiffXSPlot(draw_config_pt, process='NC_DIS_e18_p275')
312 
313 
314 # Flavour tagging study code
315 def DifferentialTaggingYield(df, x='Jet.PT', xrange=[10,50], xbins=[10,12.5,15,20,25,30,40,50], which='all', process='CC_DIS_e18_p275_B15'):
316  global u_fb, u_mb, xs_map, n_gen
317  n_gen = len(df)
318  print(f"n_gen = {n_gen}")
319  jet_pt = np.concatenate(df['Jet.PT'].to_numpy()).ravel()
320  jet_eta = np.concatenate(df['Jet.Eta'].to_numpy()).ravel()
321  jet_flavor = np.concatenate(df['Jet.Flavor'].to_numpy()).ravel()
322  jet_tag = np.concatenate(df['Jet.BTag'].to_numpy()).ravel()
323 
324  jet_x = np.concatenate(df[x].to_numpy()).ravel()
325 
326 
327  #jet_basics = (0.01 < jet_eta) & (jet_eta < 0.9)
328  jet_basics = (jet_tag == 1)
329 
330  all_flavor = ( jet_flavor > -999.0 ) & (jet_basics)
331  light_flavor = ( jet_flavor < 4 ) & (jet_basics)
332  charm_flavor = ( jet_flavor == 4 ) & (jet_basics)
333 
334  selection = all_flavor
335  if which == 'charm':
336  selection = charm_flavor
337  elif which == 'light':
338  selection = light_flavor
339 
340  (counts, bins) = np.histogram(jet_x[ selection ], range=xrange,
341  bins=xbins)
342  #bins=40)
343 
344  bin_widths = np.diff(bins)
345  bin_centers = bins[:-1] + bin_widths/2
346 
347  errors = np.sqrt(counts)
348  rel_errors = errors/counts
349 
350 
351  # convert counts to dsigma/dpT * 100/fb
352  dsigma_dx = counts * eat.TotalXSection(process) * 100*u_fb**(-1) / n_gen / bin_widths
353  dsigma_dx_errors = rel_errors * dsigma_dx
354 
355  return (bin_centers, bin_widths, dsigma_dx, dsigma_dx_errors)
356 
357 
358 
359 # Inputs
360 
361 draw_config_eta={}
362 
363 draw_config_eta['xvar'] = 'Jet.Eta'
364 draw_config_eta['xrange'] = [-4.0, 4.0]
365 draw_config_eta['xbins'] = [-4.0, -3.5, -2.5, -1.0, 1.0, 2.5, 3.5, 4.0]
366 draw_config_eta['ylimits'] = [1, 1e3]
367 draw_config_eta['xlimits'] = [-5,5]
368 draw_config_eta['yunits'] = ''
369 draw_config_eta['xunits'] = ''
370 draw_config_eta['data'] = df_ct18nnlo
371 
372 draw_config_pt={}
373 draw_config_pt['xvar'] = 'Jet.PT'
374 draw_config_pt['xrange'] = [10,50]
375 draw_config_pt['xbins'] = [10,12.5,15,20,25,30,40,50]
376 draw_config_pt['ylimits'] = [1, 1e3]
377 draw_config_pt['xlimits'] = [0,60]
378 draw_config_pt['yunits'] = '[$\\mathrm{GeV^{-1}}$]'
379 draw_config_pt['xunits'] = '[GeV]'
380 draw_config_pt['data'] = df_ct18nnlo
381 
382 
383 
384 
385 def DrawDiffTagYieldPlot(draw_config={}, process='CC_DIS_e18_p275'):
386  xvar=draw_config['xvar']
387  xrange=draw_config['xrange']
388  xbins=draw_config['xbins']
389  ylimits=draw_config['ylimits']
390  xlimits=draw_config['xlimits']
391  yunits=draw_config['yunits']
392  xunits=draw_config['xunits']
393 
394  fig, ax = plt.subplots(figsize=(12,8))
395  plt.axis('off')
396 
397 
398  gridspec = fig.add_gridspec(ncols=1, nrows=1, width_ratios=[1], height_ratios=[1])
399 
400 
401  # Log plot of differential cross-sections
402  ax1 = fig.add_subplot(gridspec[0, 0])
403 
404  ax1.grid(which='both', axis='both')
405 
406  (bins, bin_widths, xs, xs_error) = DifferentialTaggingYield(draw_config['data'], x=xvar, xrange=xrange, xbins=xbins, which='light', process=process)
407  ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='light jets (CT18NNLO)', marker='o', ms=10, ls='none', linewidth=2, color='red')
408 
409 
410  charm_ct18nnlo = DifferentialTaggingYield(draw_config['data'], x=xvar, xrange=xrange, xbins=xbins, which='charm', process=process)
411  (bins, bin_widths, xs, xs_error) = charm_ct18nnlo
412  ax1.errorbar(bins, xs, xerr = bin_widths/2, yerr=xs_error, label='charm jets (CT18NNLO)', marker='D', ms=10, ls='none', linewidth=2, fillstyle='none', color='blue')
413 
414 
415 
416  xvar_symbol="p_T"
417  if xvar.find('Eta') != -1:
418  xvar_symbol='\\eta'
419 
420  plt.ylabel(f'$\\varepsilon_{{tag}} \\times d\\sigma/d{xvar_symbol} \\times 100\\mathrm{{fb^{{-1}}}}$ {yunits}')
421  plt.xlabel(f'Reconstructed Jet ${xvar_symbol}$ {xunits}')
422  plt.yscale('log')
423 
424  process_name = 'CC-DIS'
425  if process.find('NC') != -1:
426  process_name = 'NC-DIS'
427  plt.title(f"{process_name}, 18x275GeV, $Q^2>100\\mathrm{{GeV^2}}$", fontsize=20)
428 
429  ax1.set_ylim(ylimits)
430  ax1.set_xlim(xlimits)
431 
432  ax1.legend(fontsize=18)
433 
434 
435  plt.tight_layout()
436 
437  plt.savefig(f"jet_differential_tagyield_{xvar}_{process}.png")
438  plt.savefig(f"jet_differential_tagyield_{xvar}_{process}.pdf")
439 
440 
441 #DrawDiffTagYieldPlot(draw_config_pt)
442 #DrawDiffTagYieldPlot(draw_config_eta)
443 
444 # Neutral Current Process
445 draw_config_pt={}
446 draw_config_pt['xvar'] = 'Jet.PT'
447 draw_config_pt['xrange'] = [10,50]
448 draw_config_pt['xbins'] = [10,12.5,15,20,25,30,40,50]
449 draw_config_pt['ylimits'] = [1e-1, 1e6]
450 draw_config_pt['xlimits'] = [0,60]
451 draw_config_pt['yunits'] = '[$\\mathrm{GeV^{-1}}$]'
452 draw_config_pt['xunits'] = '[GeV]'
453 draw_config_pt['data'] = df_ct18nnlo_nc
454 #DrawDiffTagYieldPlot(draw_config_pt, process='NC_DIS_e18_p275_B15')
455 
456 
457 ### Tagging Efficiency
458 # Flavour tagging study code
459 def DifferentialTaggingEfficiency(df, x='Jet.PT', xrange=[10,50], xbins=[10,12.5,15,20,25,30,40,50], which='all'):
460  n_gen = len(df)
461  print(f"n_gen = {n_gen}")
462  jet_pt = np.concatenate(df['Jet.PT'].to_numpy()).ravel()
463  jet_eta = np.concatenate(df['Jet.Eta'].to_numpy()).ravel()
464  jet_flavor = np.concatenate(df['Jet.Flavor'].to_numpy()).ravel()
465  jet_tag = np.concatenate(df['Jet.BTag'].to_numpy()).ravel()
466 
467  jet_x = np.concatenate(df[x].to_numpy()).ravel()
468 
469 
470  jet_tagged = (jet_tag == 1)
471 
472  all_flavor = ( jet_flavor > -999.0 )
473  light_flavor = ( jet_flavor < 4 )
474  charm_flavor = ( jet_flavor == 4 )
475 
476  selection = all_flavor
477  if which == 'charm':
478  selection = charm_flavor
479  elif which == 'light':
480  selection = light_flavor
481 
482  (tru_counts, bins) = np.histogram(jet_x[ selection ], range=xrange,
483  bins=xbins)
484 
485  (tag_counts, bins) = np.histogram(jet_x[ (selection) & (jet_tagged) ], range=xrange,
486  bins=xbins)
487 
488 
489  eff = tag_counts/tru_counts
490  n_err = scipy.stats.binom.interval(1.0-math.exp(-1), tru_counts, eff)
491 
492  eff_err=[np.fabs(n_err[0]/tru_counts-eff), np.fabs(n_err[1]/tru_counts-eff)]
493 
494  bin_widths = np.diff(bins)
495  bin_centers = bins[:-1] + bin_widths/2
496 
497 
498  return (bin_centers, bin_widths, eff, eff_err)
499 
500 
501 # Inputs
502 draw_config_eta={}
503 
504 draw_config_eta['xvar'] = 'Jet.Eta'
505 draw_config_eta['xrange'] = [-4.0, 4.0]
506 draw_config_eta['xbins'] = [-4.0, -3.5, -2.5, -1.0, 1.0, 2.5, 3.5, 4.0]
507 draw_config_eta['ylimits'] = [1e-5, 1]
508 draw_config_eta['xlimits'] = [-5,5]
509 draw_config_eta['yunits'] = ''
510 draw_config_eta['xunits'] = ''
511 
512 draw_config_pt={}
513 draw_config_pt['xvar'] = 'Jet.PT'
514 draw_config_pt['xrange'] = [10,50]
515 draw_config_pt['xbins'] = [10,12.5,15,20,25,30,40,50]
516 draw_config_pt['ylimits'] = [1e-5, 1]
517 draw_config_pt['xlimits'] = [0,60]
518 draw_config_pt['yunits'] = '[$\\mathrm{GeV^{-1}}$]'
519 draw_config_pt['xunits'] = '[GeV]'
520 
521 
522 
523 def DrawDiffTagEfficiencyPlot(draw_config={}):
524  xvar=draw_config['xvar']
525  xrange=draw_config['xrange']
526  xbins=draw_config['xbins']
527  ylimits=draw_config['ylimits']
528  xlimits=draw_config['xlimits']
529  yunits=draw_config['yunits']
530  xunits=draw_config['xunits']
531 
532 
533  fig, ax = plt.subplots(figsize=(12,8))
534  plt.axis('off')
535 
536 
537  gridspec = fig.add_gridspec(ncols=1, nrows=1, width_ratios=[1], height_ratios=[1])
538 
539 
540  # Log plot of differential cross-sections
541  ax1 = fig.add_subplot(gridspec[0, 0])
542 
543  ax1.grid(which='both', axis='both')
544 
545  (bins, bin_widths, eff, eff_error) = DifferentialTaggingEfficiency(df_ct18nnlo, x=xvar, xrange=xrange, xbins=xbins, which='light')
546  ax1.errorbar(bins, eff, xerr = bin_widths/2, yerr=eff_error, label='light jets (CT18NNLO)', marker='o', ms=10, ls='none', linewidth=2, color='red')
547 
548 
549  (bins, bin_widths, eff, eff_error) = DifferentialTaggingEfficiency(df_ct18nnlo, x=xvar, xrange=xrange, xbins=xbins, which='charm')
550  ax1.errorbar(bins, eff, xerr = bin_widths/2, yerr=eff_error, label='charm jets (CT18NNLO)', marker='D', ms=10, ls='none', linewidth=2, color='blue')
551 
552 
553 
554  xvar_symbol="p_T"
555  if xvar.find('Eta') != -1:
556  xvar_symbol='\\eta'
557 
558  plt.ylabel(f'$\\varepsilon_{{tag}}$ {yunits}')
559  plt.xlabel(f'Reconstructed Jet ${xvar_symbol}$ {xunits}')
560 
561  plt.title("CC-DIS, 18x275GeV, $Q^2>100\\mathrm{GeV^2}$", fontsize=20)
562 
563  ax1.set_ylim(ylimits)
564  ax1.set_xlim(xlimits)
565 
566  ax1.legend(fontsize=18)
567 
568  plt.yscale('log')
569 
570  plt.tight_layout()
571 
572  plt.savefig(f"jet_differential_tagefficiency_{xvar}.png")
573  plt.savefig(f"jet_differential_tagefficiency_{xvar}.pdf")
574 
575 #DrawDiffTagEfficiencyPlot(draw_config_pt)
576 #DrawDiffTagEfficiencyPlot(draw_config_eta)