error in Sybill2.3d strange D meson decays?

[kcarloni]

Hi, I think the D_s decays in Sybill2.3d are missing the tau+nu_tau decay mode, which has a branching fraction of roughly 6%:

In [1]: import chromo
   ...: 

In [2]: import numpy as np
   ...: 

In [3]: kinematics = chromo.kinematics.CenterOfMass(
   ...:     13 * chromo.constants.TeV,
   ...:     "proton", "proton"
   ...: )

In [4]: generator = chromo.models.Sibyll23d(kinematics)
   ...: 
 ====================================================
 |                                                  |
 |                 S I B Y L L  2.3d                |
 |                                                  |
 |         HADRONIC INTERACTION MONTE CARLO         |
 |                        BY                        |
 |            Eun-Joo AHN, Felix RIEHN              |
 |      R. ENGEL, A. FEDYNITCH, R.S. FLETCHER,      |
 |       T.K. GAISSER, P. LIPARI, T. STANEV         |
 |                                                  |
 | Publication to be cited when using this program: |
 | Eun-Joo AHN et al., Phys.Rev. D80 (2009) 094003  |
 | F. RIEHN et al., hep-ph: 1912.03300              |
 | last modifications: F. Riehn (05/20/2020)        |
 ====================================================

 SIG_AIR_INI: initializing target: (i,A)           1           0 air..
 SIG_AIR_INI: initializing target: (i,A)           2          14 nit..
 SIG_AIR_INI: initializing target: (i,A)           3          16 oxy..

In [5]: n_try = 10_000
   ...: num_taus = 0
   ...: num_D_strange = 0

In [6]: for evt in generator(n_try):
   ...:     num_taus += np.sum( evt.pid == 15 )
   ...:     num_taus += np.sum( evt.pid == -15 )
   ...:     num_D_strange += np.sum( evt.pid == 431 )
   ...:     num_D_strange += np.sum( evt.pid == -431 )
   ...: 

In [7]: num_D_strange
   ...: 
Out[7]: 469

In [8]: num_taus
   ...: 
Out[8]: 0

Thanks very much for this package -- it's been amazing to use!

[afedynitch]

Hi @kcarloni.

It might be a feature of SIBYLL for this particular decay channel. AFAIK the decay routines are from the predecessor of PYTHIA6 and may not be perfectly up to date. To workaround this issue, one possible solution is to "route" the desired particles to PYTHIA8 for decays. We have a class called Pythia8DecayHandler, but it's not active by default for SIBYLL. Unfortunately, the code must be somewhat modified for taking custom settings about particles should be forwarded to PYTHIA. I'll try working on it whenever time permits.

[kcarloni]

Hi, thanks for your reply! I think I successfully rerouted it by turning off charmed meson decays + then manually passing events to the Pythia8DecayHandler, e.g.:

k_default = chromo.kinematics.CenterOfMass( 13 * chromo.constants.TeV, "proton", "proton" )
g = chromo.models.Sibyll23d( k_default )

use_pythia_decay = [ 411, 421, 431 ]
for p in use_pythia_decay:
        g.set_stable(p)
        
stable_pids = [ 11, 12, 14, 16, 22, 2212 ]
decay_h = chromo.decay_handler.Pythia8DecayHandler(stable_pids)

for evt in g(10):
        decay_h(evt)
        # ....

Just wanted to let the team know / make an issue in case others notice this.