14 #include "Pythia8/Pythia.h" 34 bool isotropic,
double string_formation_time)
35 :
Action({in_part_a, in_part_b}, time),
36 total_cross_section_(0.),
37 isotropic_(isotropic),
38 string_formation_time_(string_formation_time) {}
50 const auto &log = logger<LogArea::ScatterAction>();
90 "ScatterAction::generate_final_state: Invalid process type " +
91 std::to_string(static_cast<int>(
process_type_)) +
" was requested. " +
98 new_particle.boost_momentum(
102 new_particle.set_4position(middle_point);
108 double elastic_parameter,
bool two_to_one,
ReactionsBitSet included_2to2,
109 double low_snn_cut,
bool strings_switch,
bool use_AQM,
114 elastic_parameter, two_to_one, included_2to2, low_snn_cut, strings_switch,
126 if (!strings_with_probability &&
151 return (1. / std::sqrt(1.0 -
beta_cm().sqr()));
181 const auto &log = logger<LogArea::ScatterAction>();
183 " position difference [fm]: ", pos_diff,
184 ", momentum difference [GeV]: ", mom_diff);
186 const double dp2 = mom_diff.
sqr();
187 const double dr2 = pos_diff.
sqr();
192 const double dpdr = pos_diff * mom_diff;
202 return dr2 - dpdr * dpdr / dp2;
221 return 5.5 * p8 / (7.7 + p8);
223 return std::min(9.0, 5.334 + 0.67 * (plab - 2.));
240 }
else if (plab < 0.6) {
241 return 16.53 * (plab - 0.225);
242 }
else if (plab < 1.6) {
243 return -1.63 * plab + 7.16;
250 double kinetic_energy_cm) {
259 const auto &log = logger<LogArea::ScatterAction>();
274 const double mass_a = masses.first;
275 const double mass_b = masses.second;
277 const std::array<double, 2> t_range = get_t_range<double>(
278 kinetic_energy_cm, mass_in_a, mass_in_b, mass_a, mass_b);
284 double mandelstam_s_new = 0.;
300 double bb, a, plab =
plab_from_s(mandelstam_s_new);
307 a = (plab < 0.8) ? 1. : 0.64 / (plab * plab);
317 t = t_range[0] + t_range[1] - t;
321 1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
334 t = t_range[0] + t_range[1] - t;
337 1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
341 const std::array<double, 4>
p{1.46434, 5.80311, -6.89358, 1.94302};
342 const double a =
p[0] + mass_a * (
p[1] + mass_a * (
p[2] + mass_a *
p[3]));
346 double t = t_range[0];
351 t = t_range[0] + t_range[1] - t;
354 1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
357 phitheta.distribute_isotropically();
367 const double p_f =
pCM(kinetic_energy_cm, mass_a, mass_b);
369 log.warn(
"Particle: ", p_a->
pdgcode(),
" radial momentum: ", p_f);
370 log.warn(
"Etot: ", kinetic_energy_cm,
" m_a: ", mass_a,
" m_b: ", mass_b);
378 log.debug(
"p_a: ", *p_a,
"\np_b: ", *p_b);
400 const size_t index_tmax = (t0 > t1) ? 0 : 1;
401 const double form_time_begin =
425 const auto &log = logger<LogArea::ScatterAction>();
429 "resonance_formation: " 430 "Incorrect number of particles in final state: ";
445 const size_t index_tmax = (t0 > t1) ? 0 : 1;
446 const double begin_form_time =
458 log.debug(
"Momentum of the new particle: ",
467 const auto &log = logger<LogArea::Pythia>();
474 bool success =
false;
476 const int ntry_max = 10000;
477 while (!success && ntry < ntry_max) {
504 log.error(
"Unknown string process required.");
508 if (ntry == ntry_max) {
545 if (tform_in > tform_out) {
554 out_mom += data.momentum();
557 log.debug(
"Outgoing momenta string:", out_mom);
565 out <<
" (not performed)";
void add_collisions(CollisionBranchList pv)
Add several new collision channels at once.
void add_collision(CollisionBranchPtr p)
Add a new collision channel.
double string_probability(bool strings_switch, bool use_transition_probability, bool use_AQM, bool treat_nnbar_with_strings) const
double high_energy() const
Determine the parametrized total cross section at high energies for the given collision, which is non-zero for Baryon-Baryon and Nucleon-Pion scatterings currently.
bool is_Deltastar() const
double partial_cross_section_
Partial cross-section to the chosen outgoing channel.
bool next_NDiffSoft()
Soft Non-diffractive process is modelled in accordance with dual-topological approach Capella:1978ig...
int antiparticle_sign() const
The ThreeVector class represents a physical three-vector with the components .
void generate_final_state() override
Generate the final-state of the scattering process.
T expo(T A, T x1, T x2)
Draws a random number x from an exponential distribution exp(A*x), where A is assumed to be positive...
double diffractive. Two strings are formed, one from A and one from B.
FourVector get_interaction_point() const
Get the interaction point.
constexpr double really_small
Numerical error tolerance.
ThreeVector velocity() const
Get the velocity (3-vector divided by zero component).
bool is_string_soft_process(ProcessType p)
Check if a given process type is a soft string excitation.
double cm_momentum() const
Get the momentum of the center of mass of the incoming particles in the calculation frame...
ProcessType process_type_
type of process
double gamma_cm() const
Get the gamma factor corresponding to a boost to the center of mass frame of the colliding particles...
Thrown when ScatterAction is called to perform with unknown ProcessType.
ThreeVector beta_cm() const
Get the velocity of the center of mass of the scattering particles in the calculation frame...
Guard type that safely disables floating point traps for the scope in which it is placed...
Collection of useful constants that are known at compile time.
a special case of baryon-antibaryon annihilation.
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
StringProcess * string_process_
Pointer to interface class for strings.
FourVector total_momentum() const
Sum of 4-momenta of incoming particles.
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept
constexpr T pow_int(const T base, unsigned const exponent)
Efficient template for calculating integer powers using squaring.
2->2 inelastic scattering
const FourVector & momentum() const
Get the particle's 4-momentum.
constexpr double nucleon_mass
Nucleon mass in GeV.
double get_partial_weight() const override
Get the partial cross section of the chosen channel.
CollisionBranchList string_excitation(double total_string_xs, StringProcess *string_process, bool use_AQM) const
Determine the cross section for string excitations, which is given by the difference between the para...
ProcessType get_type() const override
void resonance_formation()
Perform a 2->1 resonance-formation process.
NNbarTreatment
Treatment of N Nbar Annihilation.
void string_excitation()
Todo(ryu): document better - it is not really UrQMD-based, isn't it? Perform the UrQMD-based string e...
void boost(const ThreeVector &v)
Apply a full Lorentz boost of momentum and position.
static double Cugnon_bnp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic np...
double plab_from_s(double mandelstam_s, double mass)
Convert Mandelstam-s to p_lab in a fixed-target collision.
double sqr() const
calculate the square of the vector (which is a scalar)
ScatterAction(const ParticleData &in_part1, const ParticleData &in_part2, double time, bool isotropic=false, double string_formation_time=1.0)
Construct a ScatterAction object.
double total_cross_section_
Total hadronic cross section.
void rotate_z_axis_to(ThreeVector &r)
Rotate the z-axis onto the vector r.
double mandelstam_s() const
Determine the Mandelstam s variable,.
ThreeVector threevec() const
ParticleList particle_list() const
elastic scattering: particles remain the same, only momenta change
void init(const ParticleList &incoming, double tcoll)
initialization feed intial particles, time of collision and gamma factor of the center of mass...
int charge() const
The charge of the particle.
CollisionBranch is a derivative of ProcessBranch, which is used to represent particular final-state c...
bool next_SDiff(bool is_AB_to_AX)
Single-diffractive process is based on single pomeron exchange described in Ingelman:1984ns.
const FourVector & position() const
Get the particle's position in Minkowski space.
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
bool next_DDiff()
Double-diffractive process ( A + B -> X + X ) is similar to the single-diffractive process...
void set_4momentum(const FourVector &momentum_vector)
Set the particle's 4-momentum directly.
ParticleList incoming_particles_
List with data of incoming particles.
hard string process involving 2->2 QCD process by PYTHIA.
std::pair< FourVector, FourVector > get_potential_at_interaction_point() const
Get the skyrme and asymmetry potential at the interaction point.
CollisionBranchList generate_collision_list(double elastic_parameter, bool two_to_one_switch, ReactionsBitSet included_2to2, double low_snn_cut, bool strings_switch, bool use_AQM, bool strings_with_probability, NNbarTreatment nnbar_treatment, StringProcess *string_process) const
Generate a list of all possible collisions between the incoming particles with the given c...
FourVector total_momentum_of_outgoing_particles() const
Calculate the total kinetic momentum of the outgoing particles.
double cm_momentum_squared() const
Get the squared momentum of the center of mass of the incoming particles in the calculation frame...
void add_all_scatterings(double elastic_parameter, bool two_to_one, ReactionsBitSet included_2to2, double low_snn_cut, bool strings_switch, bool use_AQM, bool strings_with_probability, NNbarTreatment nnbar_treatment)
Add all possible scattering subprocesses for this action object.
void elastic_scattering()
Perform an elastic two-body scattering, i.e. just exchange momentum.
std::bitset< 6 > ReactionsBitSet
Container for the 2 to 2 reactions in the code.
Action is the base class for a generic process that takes a number of incoming particles and transfor...
bool isotropic_
Do this collision isotropically?
ParticleList get_final_state()
a function to get the final state particle list which is called after the collision ...
Use string fragmentation.
static double Cugnon_bpp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic pp...
T power(T n, T xMin, T xMax)
Draws a random number according to a power-law distribution ~ x^n.
The cross section class assembels everything that is needed to calculate the cross section and return...
const ParticleType & type() const
Get the type of the particle.
void sample_2body_phasespace()
Sample the full 2-body phase-space (masses, momenta, angles) in the center-of-mass frame for the fina...
PdgCode pdgcode() const
Get the pdgcode of the particle.
resonance formation (2->1)
CollisionBranchList collision_channels_
List of possible collisions.
double get_total_weight() const override
Get the total cross section of scattering particles.
single diffractive AB->XB.
Angles provides a common interface for generating directions: i.e., two angles that should be interpr...
virtual double cross_section() const
Get the total cross section of the scattering particles.
bool next_BBbarAnn()
Baryon-antibaryon annihilation process Based on what UrQMD Bass:1998ca, Bleicher:1999xi does...
void sample_angles(std::pair< double, double > masses, double kinetic_energy_cm) override
Sample final-state angles in a 2->2 collision (possibly anisotropic).
non-diffractive. Two strings are formed both have ends in A and B.
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
The FourVector class holds relevant values in Minkowski spacetime with (+, −, −, −) metric signature.
bool next_NDiffHard()
Hard Non-diffractive process is based on PYTHIA 8 with partonic showers and interactions.
double transverse_distance_sqr() const
Calculate the transverse distance of the two incoming particles in their local rest frame...
ParticleData contains the dynamic information of a certain particle.
(41-45) soft string excitations.
const double time_of_execution_
Time at which the action is supposed to be performed (absolute time in the lab frame in fm/c)...
void format_debug_output(std::ostream &out) const override
Writes information about this scatter action to the out stream.
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
void inelastic_scattering()
Perform an inelastic two-body scattering, i.e. new particles are formed.