BcVegPy

GitLab Repository ./enable BcVegPy/HEAD

Authors: Philip Ilten ; Jake Pfaller

Keywords: B_c mesons, heavy quarkonium, Les Houches Accord, Fortran

Plugins: LHAupBcVegPy (LHAup)

Dependencies: PYTHIA8>=8.318, GFORTRAN

Pythia Versions: 8.318

PythiaParallel: No

XML Documentation: No

Test Image: pythia8/dev:test

License: GNU GPL v2 or later

---

This package interfaces the BCVEGPY Fortran generator for the hadronic production of \(B_c\) mesons with Pythia. BCVEGPY events are passed to Pythia through the Les Houches Accord user-process class LHAupBcVegPy, after which Pythia performs the showering, hadronization, and decays. The S-wave, P-wave, color-octet, and 2S \(B_c\) states of BCVEGPY are available.

Requirements

• Pythia 8.318.
• gfortran, used to build the BCVEGPY Fortran source and to link against libgfortran.

Building

The package uses the standard PYTHIA-CONTRIB build system. The pythia8-config script must either be along PATH or passed to configure.

./configure --with-pythia8-config=/path/to/pythia8-config
make -j

This first builds the Fortran source into lib/libBcVegPyFortran.a, and then the plugin library lib/libBcVegPy.so.

Usage

The BCVEGPY grid files are located through the BCVEGPYDATAROOT environment variable, which must point to share/BcVegPy of this package.

export BCVEGPYDATAROOT=$PWD/share/BcVegPy

The generator is then loaded in Pythia with the following setting. Note that loading the plugin sets Beams:frameType = 5, adds the 2S states 100541 (\(B_c(2S)^+\)) and 100543 (\(B_c^*(2S)^+\)), with \(B_c^{(*)+}\pi\pi\) decays, to the particle data if they are not already defined, and updates the P-wave states 10543, 10541, 20543, and 545 (states 3 to 6 below) if they have not been changed. These are given the absolute values of the BcVegPy:masses of states 3 to 6, which are below the \(B+D\) threshold, as observed by LHCb, zero widths, and radiative decays: \(^1P_1 \to B_c^+\gamma\), and \(^3P_J \to B_c^{*+}\gamma\). If any of these particles is changed, e.g., its mass, it is not updated, and so its decays should also be set.

Init:plugins = {libBcVegPy.so::LHAupBcVegPy}

The center-of-mass energy, in GeV, must be set with the BCVEGPY command upcom ecm, e.g., BcVegPy:cmnds += {upcom ecm 14000}. Pre-generated VEGAS grids, share/BcVegPy/data/<energy>/grade*.dat, are provided for 7, 8, 13, 13.6, and 14 TeV, and are read by default (loggrade ivegasopen 0); other energies are not supported by the provided grids. Setting loggrade ivegasopen 1 instead recalculates the grid for the requested state with VEGAS and rewrites the corresponding existing grid file, while with loggrade iusecurdir 1 the grid files are taken from the current directory rather than BCVEGPYDATAROOT. A minimal test that loads and initializes the plugin is provided in share/BcVegPy/examples.

cd share/BcVegPy/examples
make testPlugins
./testPlugins LHAupBcVegPy.cmnd

Settings

The following Pythia settings are defined by the package.

• BcVegPy:mesonIDs (default {541, 543, 10543, 10541, 20543, 545, 541, 541, 100541, 100543}): PDG IDs of the ten \(B_c\) states listed below. Exactly ten values must be given.
• BcVegPy:massRatios (default 3.733 for each state): the ratio \(m_b/m_c\) for each of the ten \(B_c\) states. Exactly ten values must be given.
• BcVegPy:masses (default {-6.27447, -6.3333, -6.741, -6.706, -6.750, -6.768, -6.27447, -6.27447, -6.8721, -6.900}): the mass, in GeV, of each of the ten \(B_c\) states. Exactly ten values must be given, and the resulting masses must be positive. A positive mass is used as is. For a negative mass, the mass of the corresponding BcVegPy:mesonIDs entry in the Pythia particle data is used first, followed by the absolute value of the given mass if that meson is not defined.
• BcVegPy:weightMax (default 0): the maximum event weight, in pb, for unweighted events, i.e., XMAXUP, above which BCVEGPY limits the event weights, see below. For 0, the BCVEGPY maximum is used, i.e., 200000 pb, or the maximum found by VEGAS with loggrade ivegasopen 1.
• BcVegPy:defaultCmnds: list of default commands used to configure BCVEGPY, see below. The defaults are given in src/Package.cc, and normally should not need to be changed.
• BcVegPy:cmndFiles (default empty): list of files with BCVEGPY commands, one command per line. Empty lines, and lines not beginning with a letter or digit, e.g., ! a comment, are ignored.
• BcVegPy:cmnds (default empty): list of additional BCVEGPY commands, which can be appended with +=.

The commands are applied in the order BcVegPy:defaultCmnds, then the files of BcVegPy:cmndFiles, and finally BcVegPy:cmnds, where a later command overrides an earlier one.

BCVEGPY Commands

BCVEGPY is configured by setting variables in its Fortran common blocks. Each command has the form

<common block> <variable> <value> [<index>]

and is case insensitive.

• <value>: must be a complete number of the variable's type, e.g., 4.5 is invalid for an integer variable. Complex values are given as <real>:<imaginary>.
• <index>: is an optional index of an array, and is given as a colon-separated list in Fortran order, i.e., the first index varies fastest. Each index is a zero-based offset from the lower bound of its dimension. Without indices, the first element is set.

If an invalid command is given, e.g., an unknown variable, an index outside its dimension, or a malformed value, it is reported as an error and Pythia initialization fails.

The \(B_c\) state to generate is selected with counter ibcstate <state>, while mixevnt imix 1 mixes the states.

State	Description
1	color-singlet \(^1S_0\)
2	color-singlet \(^3S_1\)
3	color-singlet \(^1P_1\)
4	color-singlet \(^3P_0\)
5	color-singlet \(^3P_1\)
6	color-singlet \(^3P_2\)
7	color-octet \(^1S_0\)
8	color-octet \(^3S_1\)
9	color-singlet \(2S\) \(^1S_0\)
10	color-singlet \(2S\) \(^3S_1\)

The \(B_c\) mass of each state is set by BcVegPy:masses. By default, all of these masses are negative, so the masses are taken from the Pythia particle data, which for the P-wave and 2S states are set by the plugin as described above, and can be changed with, e.g., 541:m0 = 6.2745, while the default BCVEGPY masses are used for any meson that is not defined. The 2S states added to the particle data by the plugin use the absolute values of the masses of states 9 and 10. The quark masses and wave-function parameters (pmbc, pmb, pmc, fbc, and fbcc of the upcom and wavezero common blocks) are recomputed from these masses, the selected state, and BcVegPy:massRatios after all commands have been read, so these should not be set directly. Note that the provided grids were calculated with the default BCVEGPY masses. The grids are still valid for other masses, but are less efficient, and the cross sections stored with them, which determine the fraction of each state when the states are mixed, are not updated.

By default, unweighted events are generated (usertran idpp 1), for which BCVEGPY limits the event weights to its maximum of 200000 pb (XMAXUP). With the Pythia PDFs and \(\alpha_s\), this underestimates the cross section, e.g., by a factor of about 2 for the S-wave states at 13 TeV, and biases the kinematic distributions. Each event with a limited weight is reported as a warning, so their number is given by Pythia::stat(). The maximum can be raised with BcVegPy:weightMax, which reduces this bias, but also the unweighting efficiency, e.g., at 13 TeV with the provided grids, a maximum of 1e7 pb limits the S-wave cross sections by less than 1%, with an efficiency of about 1%. Weighted events, using the Les Houches strategy 4, are generated with usertran idpp 4, where the event weight, in pb, is given by Info::weight(). This avoids the limit, but with the provided grids the weights vary strongly, so that many more events are needed for the same statistical precision. Only the values 1 and 4 of idpp are supported, and other values fail the initialization.

Parameters

The BCVEGPY parameters that can be set with commands are listed below, with the defaults of BcVegPy:defaultCmnds. The descriptions are quoted from the BCVEGPY papers, BCVEGPY1.0 hep-ph/0309120, BCVEGPY2.0 hep-ph/0504017, and BCVEGPY2.1 hep-ph/0604238, which describe the stand-alone program with PYTHIA 6. Notes on the use of each parameter with this package, where it differs, follow the quotes.

• upcom ecm (no default, must be set): "total energy for the hadron collision (in units GeV)" [1.0]. Grids are provided for 7, 8, 13, 13.6, and 14 TeV, see above.
• counter ibcstate (default 1): "state of the \(B_c\) (\(B_c^*\)) meson" [1.0], where "IBCSTATE=1, for \(B_c[^1S_0]\); IBCSTATE=2, for \(B_c^*[^3S_1]\); IBCSTATE=3, 4, ..., 6, for the four P-wave states \(B_c^*[^1P_1]\), \(B_c^*[^3P_0]\), \(B_c^*[^3P_1]\) and \(B_c^*[^3P_2]\) respectively" [2.0]. This version also provides the color-octet S-wave states 7 and 8, and the 2S states 9 and 10, see the table above.
• mixevnt imix (default 1) and mixevnt imixtype (default 1): "A useful trick for generating the mixed type of \(B_c\) events is suggested (three types of mixed events are provided in the generator, e.g. by setting the parameter IMIXTYPE=2, one can generate the mixed events for the two color-singlet \((c\bar{b})\)-quarkonium states \((c\bar{b})_1(^1S_0)\) and \((c\bar{b})_1(^3S_1)\))" [2.1]. With imix 1, each event is of one state, selected with the probabilities given by the cross sections stored in the grid files: all ten states for imixtype 1, states 1 and 2 for imixtype 2, and states 3 to 8 for imixtype 3. With imix 0, only the state counter ibcstate is generated. Note that the event weight of the selected state is also proportional to its cross section, so that each state is effectively weighted by its cross section squared, and that the stored cross sections are those of the setup used for the grids. For predictions, the states should therefore be generated separately with imix 0.
• confine ptcut (default 0.0): "\(p_T\) cut for the \(B_c\) (\(B_c^*\)) meson (in units GeV, value can be freely selected by users)" [1.0]. The matrix element is zero below the cut, so grids generated with the cut are more efficient.
• confine etacut (default 1000000000.0): "Y cut for the \(B_c\) (\(B_c^*\)) meson (value can be freely selected by the users)" [1.0], i.e., a cut on the absolute rapidity.
• funtrans nq2 (default 3): "choice of \(Q^2\), the type of the characteristic energy scale squared in the production (in units GeV\(^2\)). Here seven choices are available: i). NQ2=1: \(Q^2 = \hat{s}/4\) (\(\hat{s}\) is the squared center-of-mass energy of the subprocess); ii). NQ2=2: \(Q^2 = \hat{s}\); iii). NQ2=3: \(Q^2 = p^2_{TB_c} + m^2_{B_c}\); iv). NQ2=4: \(Q^2 = (\sqrt{p^2_{TB_c} + m^2_{B_c}} + \sqrt{p^2_{Tb} + m^2_b} + \sqrt{p^2_{Tc} + m^2_c})^2\); v). NQ2=5: \(Q^2 = (\sqrt{p^2_{TB_c} + m^2_{B_c}} + \sqrt{p^2_{Tb} + m^2_b} + \sqrt{p^2_{Tc} + m^2_c})^2/9\); vi). NQ2=6: \(Q^2 = p^2_{Tb} + m^2_b\) for the \(\alpha_s\) in parton distribution functions and in the coupling to the b-quark line, and \(Q^2 = 4m^2_c\) for the \(\alpha_s\) in the coupling to the \(\bar{c}\)-quark line; vii). NQ2=7: \(Q^2 = p^2_{Tb} + m^2_b\)" [1.0]. This version also provides nq2 8, with the scale of nq2 3 and the two \(\alpha_s\) scales of nq2 6, and nq2 9, with \(Q^2 = 4m_b^2\).
• funtrans npdfu (default 2): "choice of the collision type of hadrons ... NPDFU=1 for \(p-\bar{p}\) and NPDFU=2 for \(p-p\)" [1.0].
• loggrade igenerate (default 1): "whether to generate complete events. IGENERATE=0, when users wish the simulation to stop after the generation of the 'final state' containing the \(B_c\) meson, b-jet and \(\bar{c}\)-jet of the subprocess \(gg \to B_c + b + \bar{c}\); IGENERATE=1, when users wish that complete events including the \(B_c\) production are to be generated" [1.0]. This should be left at 1, since for unweighted events igenerate 0 would set the maximum event weight to zero, so it fails the initialization.
• loggrade ivegasopen (default 0): "whether switch on/off the VEGAS subroutine: IVEGASOPEN=1 for using VEGAS; IVEGASOPEN=0 for not using VEGAS" [1.0]. "With the option IVEGASOPEN=1, one can use the VEGAS program to achieve an important sampling-function, which may be used to increase the efficiency of the Monte Carlo simulation" [2.0]. The grid file of the state is then overwritten, see Generating Grids.
• loggrade igrade (default 1): "whether to use the grade generated by previous running VEGAS when IVEGASOPEN=0; IGRADE=1 means to use; IGRADE=0 means not to use" [1.0]. "By setting IVEGASOPEN=0 and IGRADE=1, one may generate the mixed \(B_c\) events just by reading the existent importance sampling function and total cross-section from 'previous' VEGAS runs ... with setting IVEGASOPEN=0 and IGRADE=0, VEGAS or the importance sampling function (grade function) generated in previous runs will not be used at all" [2.1]. Without the grids, the generation is very inefficient, so igrade 0 is only useful for tests.
• loggrade iusecurdir (default 0): this parameter is not described in the papers. The grid files are read from the current directory for iusecurdir 1, rather than from $BCVEGPYDATAROOT/data/<energy>.
• loggrade ievntdis (default 0): "INUMEVNT=: indicates whether to switch on or off to keep the event number distributions other than the differential distributions of \(p_T\), \(y\) and etc." [2.0]. In this version, it only changes a printout.
• subopen isubonly (default 0): "whether to keep the information only of the sub-process \(gg \to B_c(B_c^*) + b + \bar{c}\); ISUBONLY=0 for the full hadronic production, i.e. the structure functions are connected; ISUBONLY=1 for the subprocess only" [1.0]. This is only intended for checks of the subprocess, and has not been tested with this package.
• subopen subenergy (default 100.0): "the energy (in units GeV) of the sub-process \(gg \to B_c(B_c^*) + b + \bar{c}\). It is needed only when ISUBONLY=1" [1.0].
• usertran ishower (default 1): "whether to switch on/off the showers, including initial and final states, multiple interactions, hadronization" [1.0]. In this version, it only sets PYTHIA 6 switches, and so has no effect. The showers, multiparton interactions, and hadronization are instead set by Pythia, e.g., with PartonLevel:ISR, PartonLevel:FSR, PartonLevel:MPI, and HadronLevel:all.
• usertran idpp (default 1): the Les Houches strategy IDWTUP, the "master switch dictating how the event weights and the cross-sections should be interpreted ... when IDWTUP=1, parton-level events have a weight at the input to PYTHIA. Events are then either accepted or rejected, so that fully generated events at the output have a common weight" [1.0]. Only the values 1 and 4 are supported, see above.
• vegasinf number (default 1000000): "total number of times for calling the integrand (VEGAS parameter, see VEGAS manual). The parameter is needed only when IVEGASOPEN=1 in each iteration" [1.0].
• vegasinf nitmx (default 20): "upper limit for the number of iterations (VEGAS parameter, see VEGAS manual). The parameter is needed only when IVEGASOPEN=1" [1.0].
• vegcross iveggrade (default 0): "indicates whether to switch on or off to use the existed grade (importance sampling function) before running VEGAS. IVEGGRADE=1 means to switch on and IVEGGRADE=0 means to switch off. With this parameter, one can improve the efficiency of the importance sampling function" [2.0]. In this version, the existing grid is not read for iveggrade 1, so this should be left at 0.
• qqbar iqqbar (default 0): "indicates whether to switch on or off to use the less important light quark-antiquark annihilation subprocess \(q + \bar{q} \to (c\bar{b}) + \bar{c} + b\) to generate the \((c\bar{b})\)-quarkonium events, where \((c\bar{b})\)-quarkonium is in \(^1S_0\) or \(^3S_1\) state. IQQBAR=1 means to switch on and IQQBAR=0 means to switch off. When IQQBAR=0, it will use the default gluon-gluon fusion subprocess to generate the \((c\bar{b})\)-quarkonium events" [2.0]. "Numerical results of the calculation show that its contributions to the production are smaller than those from the gluon-gluon fusion (~ 1%)" [2.0]. For the states other than 1 and 2, where BCVEGPY would stop the program, iqqbar 1 fails the initialization.
• qqbar iqcode (default 2): "indicates the species of the light quark for the light quark-antiquark annihilation mechanism. Under the option IQQBAR=1, IQCODE=1 means \(u\) and \(\bar{u}\), IQCODE=2 means \(d\) and \(\bar{d}\), while IQCODE=3 means \(s\) and \(\bar{s}\) used in the quark-antiquark annihilation subprocess" [2.0].
• outpdf ioutpdf (default 0): "indicates whether to switch on or off to use the three latest version of PDFs: CTEQ6L, GRV98L and MRST2001L offered in the program for the hadronic production. IOUTPDF=1 means to switch on, while IOUTPDF=0 means to switch off and then the inner PDFs of PYTHIA is used to generate the events" [2.0]. With ioutpdf 0, the PDFs and \(\alpha_s\) are those of Pythia, e.g., set with PDF:pSet, SigmaProcess:alphaSvalue, and SigmaProcess:alphaSorder. The PDF tables needed for ioutpdf 1 are not provided with this package: the CTEQ tables are read from $BCVEGPYDATAROOT/data, and the GRV table from the current directory.
• outpdf ipdfnum (default 300): "indicates which one of the three latest version of PDFs: CTEQ6L, GRV98L and MRST2001L is used for the hadronic production. It comes into operation only under the option IOUTPDF=1. IPDFNUM=100 means to use GRV98L; IPDFNUM=200 means to use MRST2001L; IPDFNUM=300 means to use CTEQ6L" [2.0]. In this version, ipdfnum 300 uses CTEQ6L1, from cteq6l1.tbl, and no \(\alpha_s\) is provided for ipdfnum 200.
• coloct ioctet (default 0): "indicates whether to switch on or off to use the subprocess \(gg \to (c\bar{b})_8 + b + \bar{c}\) to generate the color-octet components \(|(c\bar{b})_8(^1S_0)g\rangle\) and \(|(c\bar{b})_8(^3S_1)g\rangle\). IOCTET=1 means to switch on and IOCTET=0 means to switch off" [2.0]. In this version, the color-octet states are instead selected with counter ibcstate 7 and 8, which set ioctet 1, while the other states set ioctet 0. For the P-wave and 2S states, or with iqqbar 1, where BCVEGPY would stop the program, ioctet 1 fails the initialization.
• octmatrix coeoct (default 0.1): "the value for \(\Delta_S(v)\)" [2.0], where "\(\Delta_S(v)\) is of the order \(v^2\) or so, and one can take it to be within the region of 0.10 to 0.30" [2.0]. In this version, the states 7 and 8 set \(\Delta_S(v) = 0.2\), so this setting has no effect. Since the color-octet cross sections are proportional to \(\Delta_S(v)^2\), they can be scaled to other values.

The other parameters of the papers are set as follows with this package. The masses PMBC, PMB, and PMC are set by BcVegPy:masses and BcVegPy:massRatios. The wave function at the origin FBC, which is "the radial wave function \(R(0)\) for S-wave production or the first derivative of the radial wave function \(R'(0)\) for P-wave production" [2.0], is fixed to \(R(0) = 1.241\) GeV\(^{3/2}\) for the 1S states, \(R(0) = 0.991\) GeV\(^{3/2}\) for the 2S states, and \(R'(0) = 0.44833\) GeV\(^{5/2}\) for the P-wave states. The number of events NEV is set by the Pythia event loop. The PYTHIA 6 parameters MSTP(51), MSTU(111), and PARU(111), "choice of the proton parton-distribution set", "order of \(\alpha_s\) in the evaluation in PYALPS", and "constant value of \(\alpha_s\)" [1.0], are replaced by PDF:pSet, SigmaProcess:alphaSorder, and SigmaProcess:alphaSvalue.

Generating Grids

New grid files, e.g., after changing the PDFs or BcVegPy:masses, or for a new energy, can be generated with share/BcVegPy/examples/grids.cc, which follows scripts/grids.py of the LHCb BcVegPyData package.

cd share/BcVegPy/examples
make grids
./grids 14000 grids14TeV -t          # quick test of the setup only
./grids 14000 grids14TeV [FILE ...]  # full grids, which can take days

The program runs the VEGAS integration of BCVEGPY (loggrade ivegasopen 1) for each of the ten states at the given energy, in GeV, and writes the grid files to the given directory, without generating any events. The grid files in the directory are overwritten. The -t option uses only a few VEGAS points, to test the setup, and the resulting grids should not be used. Optional Pythia command files can set, e.g., the PDFs, BcVegPy:masses, or the number of VEGAS points per iteration and iterations, with the defaults BcVegPy:cmnds += {vegasinf number 1000000, vegasinf nitmx 20}. For the energies known to BCVEGPY, i.e., 7, 8, 10, 13, 13.6, and 14 TeV, the grid files can then be copied to share/BcVegPy/data/<energy>, e.g., data/14TeV. Otherwise, the grid files can be used by running from their directory with BcVegPy:cmnds += {loggrade iusecurdir 1}.

LHCb:2019tea Comparison

The directory share/BcVegPy/validate/LHCb:2019tea compares Pythia + BcVegPy with the LHCb measurement of the \(B_c^-\) production fraction, $f_c/(f_u + f_d), and production asymmetry at 7 and 13 TeV, arXiv:1910.13404.

cd share/BcVegPy/validate/LHCb:2019tea
make
./generate 13000 [NBC] [NB] [OUTPUT] [FILE ...]  # writes pythia13000.dat
./generate 7000                                  # writes pythia7000.dat
python3 compare.py                               # writes compare.pdf

The generate program generates NBC weighted events (by default 10000) for each of the ten BCVEGPY states and NB inclusive non-diffractive events (by default 1000000) for the \(B^0\) and \(B^+\) mesons, and writes the binned predictions to OUTPUT. Optional Pythia command files are read for both samples. The compare.py script requires numpy and matplotlib, and plots the predictions against the LHCb fit versus \(p_T\), the average, and the production asymmetries. Note that BCVEGPY generates both charges equally, so the predicted asymmetry is zero.

Chang:2005bf Comparison

The directory share/BcVegPy/validate/Chang:2005bf compares BcVegPy with the BCVEGPY calculation of the color-singlet S-wave and P-wave, and color-octet S-wave, \(B_c\) states at the LHC (14 TeV) and the Tevatron (1.96 TeV), hep-ph/0502155.

cd share/BcVegPy/validate/Chang:2005bf
make
./generate 14000 [N] [OUTPUT] [FILE ...]  # writes bcvegpy14000.dat
./generate 1960                           # writes bcvegpy1960.dat
python3 compare.py                        # writes compare.pdf

The generate program uses the parameters of the paper, generates N weighted events (by default 200000) at the parton level for each of the eight states and each of the pT cuts of the paper, 0, 5, 20, 35, and 50 GeV, and writes the cross sections and distributions to OUTPUT. The VEGAS grids are first generated in grids<ECM>, and are reused by later runs. With the defaults, each energy takes roughly 1.5 hours. The compare.py script requires numpy and matplotlib, compares the total cross sections (Table II) and the color-octet to color-singlet ratios versus the pT cut (Table III), and plots the distributions as in Figs. 1 to 4 of the paper.

Development

The C++ bindings to the Fortran procedures and common blocks, include/BcVegPy/Wrapper.h and src/Wrapper.cc, are generated automatically and should not be edited by hand. They can be regenerated from the Fortran source, which requires gfortran, as follows.

cd fortran
python3 wrapper.py

References

BCVEGPY is described in hep-ph/0309120, hep-ph/0502155, hep-ph/0504017, and hep-ph/0604238.