From 6938f677dbd036041079af73e66128a137946f26 Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Thu, 3 Oct 2024 14:08:58 -0400 Subject: [PATCH] fix runtime error in get_decay_momenta... function, and add a check for on-shell mediators to ensure non-zero coupling --- src/DarkNews/integrands.py | 14 ++++++++------ src/DarkNews/processes.py | 4 ++-- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/DarkNews/integrands.py b/src/DarkNews/integrands.py index dba33db..fd7af91 100755 --- a/src/DarkNews/integrands.py +++ b/src/DarkNews/integrands.py @@ -675,7 +675,7 @@ def get_momenta_from_vegas_samples(vsamples, MC_case): return four_momenta -def get_decay_momenta_from_vegas_samples(vsamples, decay_case, PN_LAB): +def get_decay_momenta_from_vegas_samples(vsamples, decay_case, PN_LAB, rng=None): """ Construct the four momenta of all final state particles in the decay process from the vegas weights. @@ -684,10 +684,12 @@ def get_decay_momenta_from_vegas_samples(vsamples, decay_case, PN_LAB): vsamples (np.ndarray): integration samples obtained from vegas as hypercube coordinates. Always in the interval [0,1]. - MC_case (DarkNews.process.dec_case): the decay class of DarkNews + decay_case (DarkNews.process.dec_case): the decay class of DarkNews PN_LAB (np.ndarray): four-momentum of the upscattered N in the lab frame: [E, pX, pY, pZ] + rng: random number generating function that can be passed to phase_space function calls + Returns: dict: each key corresponds to a set of four momenta for a given particle involved, so the values are 2D np.ndarrays with each row a different event and each column a different @@ -732,7 +734,7 @@ def get_decay_momenta_from_vegas_samples(vsamples, decay_case, PN_LAB): "m3": m_mediator, # Z' } # Phnl, Phnl_daughter, Pz' - P1LAB_decay, P2LAB_decay, P3LAB_decay = phase_space.two_body_decay(N_decay_samples, boost=boost_scattered_N, **masses_decay, rng=MC_case.rng) + P1LAB_decay, P2LAB_decay, P3LAB_decay = phase_space.two_body_decay(N_decay_samples, boost=boost_scattered_N, **masses_decay) # Z' boost parameters boost_Z = { @@ -751,7 +753,7 @@ def get_decay_momenta_from_vegas_samples(vsamples, decay_case, PN_LAB): "m3": mm, # \ell- } # PZ', pe-, pe+ - P1LAB_decayZ, P2LAB_decayZ, P3LAB_decayZ = phase_space.two_body_decay(Z_decay_samples, boost=boost_Z, **masses_decay, rng=MC_case.rng) + P1LAB_decayZ, P2LAB_decayZ, P3LAB_decayZ = phase_space.two_body_decay(Z_decay_samples, boost=boost_Z, **masses_decay) four_momenta["P_decay_N_parent"] = P1LAB_decay four_momenta["P_decay_N_daughter"] = P2LAB_decay @@ -782,7 +784,7 @@ def get_decay_momenta_from_vegas_samples(vsamples, decay_case, PN_LAB): P2LAB_decay, P3LAB_decay, P4LAB_decay, - ) = phase_space.three_body_decay(N_decay_samples, boost=boost_scattered_N, **masses_decay, rng=MC_case.rng) + ) = phase_space.three_body_decay(N_decay_samples, boost=boost_scattered_N, **masses_decay) four_momenta["P_decay_N_parent"] = P1LAB_decay four_momenta["P_decay_ell_minus"] = P2LAB_decay @@ -804,7 +806,7 @@ def get_decay_momenta_from_vegas_samples(vsamples, decay_case, PN_LAB): "m3": 0.0, # gamma } # Phnl, Phnl', Pgamma - P1LAB_decay, P2LAB_decay, P3LAB_decay = phase_space.two_body_decay(N_decay_samples, boost=boost_scattered_N, **masses_decay, rng=MC_case.rng) + P1LAB_decay, P2LAB_decay, P3LAB_decay = phase_space.two_body_decay(N_decay_samples, boost=boost_scattered_N, **masses_decay) four_momenta["P_decay_N_parent"] = P1LAB_decay four_momenta["P_decay_N_daughter"] = P2LAB_decay diff --git a/src/DarkNews/processes.py b/src/DarkNews/processes.py index 6edee13..90478cf 100644 --- a/src/DarkNews/processes.py +++ b/src/DarkNews/processes.py @@ -279,12 +279,12 @@ def __init__(self, nu_parent, nu_daughter, final_lepton1, final_lepton2, TheoryM ## Is the mediator on shell? self.vector_on_shell = ( - (TheoryModel.mzprime is not None) and (self.m_parent - self.m_daughter > TheoryModel.mzprime) and (TheoryModel.mzprime > self.mm + self.mp) + (TheoryModel.mzprime is not None) and (self.m_parent - self.m_daughter > TheoryModel.mzprime) and (TheoryModel.mzprime > self.mm + self.mp) and (self.Dih>0) ) self.vector_off_shell = not self.vector_on_shell self.scalar_on_shell = ( - (TheoryModel.mhprime is not None) and (self.m_parent - self.m_daughter > TheoryModel.mhprime) and (TheoryModel.mhprime > self.mm + self.mp) + (TheoryModel.mhprime is not None) and (self.m_parent - self.m_daughter > TheoryModel.mhprime) and (TheoryModel.mhprime > self.mm + self.mp) and (self.Sih>0) ) self.scalar_off_shell = not self.scalar_on_shell