diff --git a/README.md b/README.md
index f003872..d07a5ab 100644
--- a/README.md
+++ b/README.md
@@ -101,7 +101,7 @@ Particle settings should be defined after the corresponding `@#` tag using the s
 * `energy`:
   * Sets the pp collision energy (in TeV) used to get the correct parent 
     kinematics
-  * Supported: `7`, `8`, `13` and `14`
+  * Supported: `7`, `8`, `13`, `13.6` and `14`
   * more types may be added in $RAPIDSIM_ROOT/rootfiles/fonll or $RAPIDSIM_CONFIG/rootfiles/fonll
 
 * `parent`:
diff --git a/rootfiles/fonll/LHCb13.6.root b/rootfiles/fonll/LHCb13.6.root
new file mode 100644
index 0000000..b88c0bc
Binary files /dev/null and b/rootfiles/fonll/LHCb13.6.root differ
diff --git a/rootfiles/fonll/LHCc13.6.root b/rootfiles/fonll/LHCc13.6.root
new file mode 100644
index 0000000..a15fe1a
Binary files /dev/null and b/rootfiles/fonll/LHCc13.6.root differ
diff --git a/rootfiles/fonll/README.md b/rootfiles/fonll/README.md
new file mode 100644
index 0000000..c8abbbb
--- /dev/null
+++ b/rootfiles/fonll/README.md
@@ -0,0 +1,3 @@
+# FONLL Heavy Quark Production
+
+Files in this directory contain histograms for dsigma/dpT and dsigma/deta with data obtained from https://www.lpthe.jussieu.fr/~cacciari/fonll/fonllform.html.  It uses CTEQ6.6 PDFs except 13.6 TeV which uses NNPDF30nlo_as0118 set.  For b quark, it is cross-section to B hadron rather than bare quark while for c quark it uses mixture of 70% of D0 and 30% of D+.  The cross-sections are done in the region 0.125 GeV <= pT <= 299.875 and -7.98 <= eta <= 7.98 with 1200 points for each variable. Other parameters are kept at default values.
diff --git a/utils/makeFONLLHistograms.py b/utils/makeFONLLHistograms.py
new file mode 100644
index 0000000..299e165
--- /dev/null
+++ b/utils/makeFONLLHistograms.py
@@ -0,0 +1,39 @@
+'''Script to create histograms for FONLL cross-sections used to generate
+kinematics. Cross-sections are obtained from
+https://www.lpthe.jussieu.fr/~cacciari/fonll/fonllform.html and this script
+needs two of them, as a pT and eta which are inputs and script outputs root
+file with necessary histograms. Details about individual settings are in
+separate README.md in corresponding directory.'''
+
+import argparse
+import numpy
+import yaml
+import ROOT
+
+parser = argparse.ArgumentParser( prog='FONLL histogram maker',
+                    description='' )
+
+parser.add_argument('-pt')
+parser.add_argument('-eta')
+parser.add_argument('-out')
+
+args = parser.parse_args()
+
+print(args)
+
+datapt = numpy.loadtxt( args.pt )
+dataeta = numpy.loadtxt( args.eta )
+
+outFile = ROOT.TFile.Open( args.out, 'RECREATE' )
+
+hheta = ROOT.TH1D( 'eta', '', 1200, -8, 8 )
+hhpt = ROOT.TH1D( 'pT', '', 1200, 0, 300 )
+
+for ii in range(1, 1200):
+  xbin = hheta.FindBin(dataeta[ii-1][0])
+  hheta.SetBinContent(xbin, dataeta[ii-1][1])
+  xbin = hhpt.FindBin(datapt[ii-1][0])
+  hhpt.SetBinContent(xbin, datapt[ii-1][1])
+
+outFile.Write()
+outFile.Close()