diff --git a/src/particles/transformation/ToFixedT.H b/src/particles/transformation/ToFixedT.H index 1515075b5..8c7bcfb9b 100644 --- a/src/particles/transformation/ToFixedT.H +++ b/src/particles/transformation/ToFixedT.H @@ -61,10 +61,13 @@ namespace transformation amrex::ParticleReal const y = p.pos(RealAoS::y); amrex::ParticleReal const t = p.pos(RealAoS::t); + // small tolerance to avoid NaN for pz<0: + amrex::ParticleReal const tol = 1.0e-8_prt; + // compute value of reference pzd = beta*gamma amrex::ParticleReal const argd = -1.0_prt + pow(m_ptd, 2); - AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid pzd arg (<=0)"); - amrex::ParticleReal const pzdf = argd > 0.0_prt ? sqrt(argd) : 0.0_prt; + // AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid pzd arg (<=0)"); + amrex::ParticleReal const pzdf = argd > 0.0_prt ? sqrt(argd) : tol; // transform momenta to dynamic units (eg, so that momenta are // normalized by mc): @@ -74,8 +77,8 @@ namespace transformation // compute value of particle pz = beta*gamma amrex::ParticleReal const arg = -1.0_prt + pow(m_ptd+pt, 2) - pow(px, 2) - pow(py, 2); - AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pz arg (<=0)"); - amrex::ParticleReal const pzf = arg > 0.0_prt ? sqrt(arg) : 0.0_prt; + // AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pz arg (<=0)"); + amrex::ParticleReal const pzf = arg > 0.0_prt ? sqrt(arg) : tol; // transform position and momentum (from fixed s to fixed t) p.pos(RealAoS::x) = x + px*t/(m_ptd+pt);