diff --git a/meson.build b/meson.build index d2e0bfb6..dbf6f02c 100644 --- a/meson.build +++ b/meson.build @@ -80,7 +80,7 @@ common_syms = ['chromo_openlogfile','chromo_closelogfile','npyrng'] sibyll_syms_base = [ 'sibyll','sibyll_ini','sib_sigma_hp','sib_sigma_hair','sib_sigma_hnuc', 'int_nuc','decsib','decpar','sibini','sibhep','sib_list', - 'isib_pid2pdg','isib_pdg2pid','pdg_ini', 'sibnuc', 'sigma_nuc_nuc'] + 'isib_pid2pdg','isib_pdg2pid','pdg_ini', 'sibnuc', 'sigma_nuc_nuc', 'int_h_nuc'] qgs_syms_base = ['cqgsini','qgsect','qgini','qgconf','qgreg','chepevt','qgcrossc','cqgshh_ha_cs'] urqmd_syms = [ 'urqmd','init','uinit','set0','params','uounit','strini','loginit','loadwtab', diff --git a/src/chromo/models/sibyll.py b/src/chromo/models/sibyll.py index 540dcbbb..683bb71e 100644 --- a/src/chromo/models/sibyll.py +++ b/src/chromo/models/sibyll.py @@ -111,10 +111,15 @@ def _get_charge(self, npart): return self._lib.schg.ichg[:npart] def _get_impact_parameter(self): - return self._lib.cnucms.b + return self._lib.cnucms.b if self._lib.cnucms.na > 0 else self._lib.s_cncm0.b def _get_n_wounded(self): - return self._lib.cnucms.na, self._lib.cnucms.nb + na = self._lib.cnucms.na if self._lib.cnucms.na > 0 else 1 + nb = ( + self._lib.cnucms.nb if self._lib.cnucms.nb > 0 else self._lib.s_cncm0.na + ) # Handle hadron-Nucleus case + + return na, nb def _history_zero_indexing(self): # Sibyll has only mothers diff --git a/tests/test_sibyll.py b/tests/test_sibyll.py index dc03de8b..c450168e 100644 --- a/tests/test_sibyll.py +++ b/tests/test_sibyll.py @@ -120,3 +120,96 @@ def test_cross_section(model): c.non_diffractive, c.inelastic - c.diffractive_xb - c.diffractive_ax - c.diffractive_xx, ) + + +def check_wounded(model, kin): + gen = model(kin) + gen.set_stable(111) + + wounded_na_list = [] + wounded_nb_list = [] + wounded_b_list = [] + cnucms_na_list = [] + cnucms_nb_list = [] + cnucms_b_list = [] + s_cncm0_na_list = [] + s_cncm0_b_list = [] + + for event in gen(10): + + # Check wounded candidates [Main check] + wounded_na, wounded_nb = event.n_wounded[0], event.n_wounded[1] + wounded_b = event.impact_parameter + wounded_na_list.append(wounded_na != 0) + wounded_nb_list.append(wounded_nb != 0) + wounded_b_list.append(wounded_b != 0) + + # Check cnucms status [Detail check] + cnucms_na, cnucms_nb = event._lib.cnucms.na, event._lib.cnucms.nb + cnucms_b = event.impact_parameter + cnucms_na_list.append(cnucms_na != 0) + cnucms_nb_list.append(cnucms_nb != 0) + cnucms_b_list.append(cnucms_b != 0) + + # Check s_cncm0 status [Detail check] + s_cncm0_na_list.append(event._lib.s_cncm0.na != 0) + s_cncm0_b_list.append(event._lib.s_cncm0.b != 0) + + main_errors = [] + detail_errors = [] + + if not any(wounded_na_list): + main_errors.append(f"{gen.label=} | wounded_na_list is empty") + if not any(wounded_nb_list): + main_errors.append(f"{gen.label=} | wounded_nb_list is empty") + if not any(wounded_b_list): + main_errors.append(f"{gen.label=} | wounded_b_list is empty") + + if not any(cnucms_na_list): + detail_errors.append(f"{gen.label=} | cnucms_na_list is empty") + if not any(cnucms_nb_list): + detail_errors.append(f"{gen.label=} | cnucms_nb_list is empty") + if not any(cnucms_b_list): + detail_errors.append(f"{gen.label=} | cnucms_b_list is empty") + if not any(s_cncm0_na_list): + detail_errors.append(f"{gen.label=} | s_cncm0_na_list is empty") + if not any(s_cncm0_b_list): + detail_errors.append(f"{gen.label=} | s_cncm0_b_list is empty") + + return main_errors, detail_errors + + +@pytest.mark.parametrize("model", get_sibylls()) +def test_wounded_proton_proton(model): + pytest.xfail( + reason="No impact parameter in proton-proton collision, neither in " + "event._lib.s_cncm0.b or event._lib.cnucms.b" + ) + main_errors, detail_errors = run_in_separate_process( + check_wounded, model, CenterOfMass(5 * TeV, "p", "p") + ) + assert not main_errors + + +@pytest.mark.parametrize("model", get_sibylls()) +def test_wounded_proton_nucleus(model): + main_errors, detail_errors = run_in_separate_process( + check_wounded, model, CenterOfMass(5 * TeV, "p", (16, 8)) + ) + assert not main_errors + + +@pytest.mark.parametrize("model", get_sibylls()) +def test_wounded_nucleus_proton(model): + main_errors, detail_errors = run_in_separate_process( + check_wounded, model, CenterOfMass(5 * TeV, (16, 8), "p") + ) + assert not main_errors + + +@pytest.mark.parametrize("model", get_sibylls()) +def test_wounded_nucleus_nucleus(model): + main_errors, detail_errors = run_in_separate_process( + check_wounded, model, CenterOfMass(5 * TeV, (16, 8), (16, 8)) + ) + assert not main_errors