smash::ScatterActionsFinder Class Reference

#include <scatteractionsfinder.h>

A simple scatter finder: Just loops through all particles and checks each pair for a collision.

It supports two collision criteria: a geometric and stochastic criterion.

Definition at line 31 of file scatteractionsfinder.h.

Inheritance diagram for smash::ScatterActionsFinder:

Public Member Functions
	ScatterActionsFinder (Configuration &config, const ExperimentParameters &parameters)
	Constructor of the finder with the given parameters. More...
double	collision_time (const ParticleData &p1, const ParticleData &p2, double dt, const std::vector< FourVector > &beam_momentum) const
	Determine the collision time of the two particles. More...
ActionList	find_actions_in_cell (const ParticleList &search_list, double dt, const double gcell_vol, const std::vector< FourVector > &beam_momentum) const override
	Search for all the possible collisions within one cell. More...
ActionList	find_actions_with_neighbors (const ParticleList &search_list, const ParticleList &neighbors_list, double dt, const std::vector< FourVector > &beam_momentum) const override
	Search for all the possible collisions among the neighboring cells. More...
ActionList	find_actions_with_surrounding_particles (const ParticleList &search_list, const Particles &surrounding_list, double dt, const std::vector< FourVector > &beam_momentum) const override
	Search for all the possible secondary collisions between the outgoing particles and the rest. More...
ActionList	find_final_actions (const Particles &, bool=false) const override
	Find some final collisions at the end of the simulation. More...
bool	is_constant_elastic_isotropic () const
	If there is only one particle sort, no decays (only elastic scatterings are possible), scatterings are isotropic and cross-section fixed to elastic_parameter_ independently on momenta, then maximal cross-section is elastic_parameter_. More...
double	max_transverse_distance_sqr (int testparticles) const
	The maximal distance over which particles can interact in case of the geometric criterion, related to the number of test particles and the maximal cross section. More...
void	dump_reactions () const
	Prints out all the 2-> n (n > 1) reactions with non-zero cross-sections between all possible pairs of particle types. More...
void	dump_cross_sections (const ParticleType &a, const ParticleType &b, double m_a, double m_b, bool final_state, std::vector< double > &plab) const
	Print out partial cross-sections of all processes that can occur in the collision of a(mass = m_a) and b(mass = m_b). More...
StringProcess *	get_process_string_ptr ()
Public Member Functions inherited from smash::ActionFinderInterface
virtual	~ActionFinderInterface ()=default

Private Member Functions
ActionPtr	check_collision_two_part (const ParticleData &data_a, const ParticleData &data_b, double dt, const std::vector< FourVector > &beam_momentum={}, const double gcell_vol=0.0) const
	Check for a single pair of particles (id_a, id_b) if a collision will happen in the next timestep and create a corresponding Action object in that case. More...
ActionPtr	check_collision_multi_part (const ParticleList &plist, double dt, const double gcell_vol) const
	Check for multiple i.e. More...

Private Attributes
ScatterActionsFinderParameters	finder_parameters_
	Struct collecting several parameters. More...
std::unique_ptr< StringProcess >	string_process_interface_
	Class that deals with strings, interfacing Pythia. More...
const bool	isotropic_
	Do all collisions isotropically. More...
const double	box_length_
	Box length: needed to determine coordinates of collision correctly in case of collision through the wall. More...
const double	string_formation_time_
	Parameter for formation time. More...

Constructor & Destructor Documentation

ScatterActionsFinder()

smash::ScatterActionsFinder::ScatterActionsFinder	(	Configuration &	config,
		const ExperimentParameters &	parameters
	)

Constructor of the finder with the given parameters.

Parameters

[in,out]	config	Configuration of smash from which we take: 1) A global elastic cross section [mb]. It will be used regardless of the species of the colliding particles. It won't be used if the value is negative. 2) An option determining whether all the scatterings are isotropic 3) Parameters of the string process
[in]	parameters	Struct of parameters determining whether to exclude some certain types of scatterings and switching among the methods to treat with the NNbar collisions.

Definition at line 34 of file scatteractionsfinder.cc.

    : finder_parameters_(create_finder_parameters(config, parameters)),
      isotropic_(config.take({"Collision_Term", "Isotropic"}, false)),
      box_length_(parameters.box_length),
      string_formation_time_(config.take(
          {"Collision_Term", "String_Parameters", "Formation_Time"}, 1.)) {
  if (is_constant_elastic_isotropic()) {
    logg[LFindScatter].info(
        "Constant elastic isotropic cross-section mode:", " using ",
        finder_parameters_.elastic_parameter, " mb as maximal cross-section.");
  }
  if (finder_parameters_.included_multi.any() &&
      finder_parameters_.coll_crit != CollisionCriterion::Stochastic) {
    throw std::invalid_argument(
        "Multi-body reactions (like e.g. 3->1 or 3->2) are only possible with "
        "the stochastic "
        "collision "
        "criterion. Change your config accordingly.");
  }
 
  if (finder_parameters_
              .included_multi[IncludedMultiParticleReactions::Deuteron_3to2] ==
          1 &&
      (finder_parameters_
               .included_2to2[IncludedReactions::PiDeuteron_to_pidprime] == 1 ||
       finder_parameters_
               .included_2to2[IncludedReactions::NDeuteron_to_Ndprime] == 1)) {
    throw std::invalid_argument(
        "To prevent double counting it is not possible to enable deuteron 3->2 "
        "reactions\nand reactions involving the d' at the same time\ni.e. to "
        "include \"Deuteron_3to2\" in `Multi_Particle_Reactions` and\n "
        "\"PiDeuteron_to_pidprime\" "
        "or \"NDeuteron_to_Ndprime\" in `Included_2to2` at the same time.\n"
        "Change your config accordingly.");
  }
 
  if (finder_parameters_
              .included_multi[IncludedMultiParticleReactions::Deuteron_3to2] ==
          1 &&
      ParticleType::try_find(pdg::dprime)) {
    throw std::invalid_argument(
        "Do not use the d' resonance and enable \"Deuteron_3to2\" "
        "`Multi_Particle_Reactions` at the same time. Either use the direct "
        "3-to-2 reactions or the d' together with \"PiDeuteron_to_pidprime\" "
        "and \"NDeuteron_to_Ndprime\" in `Included_2to2`. Otherwise the "
        "deuteron 3-to-2 reactions would be double counted.");
  }
 
  if ((finder_parameters_.nnbar_treatment == NNbarTreatment::TwoToFive &&
       finder_parameters_
               .included_multi[IncludedMultiParticleReactions::NNbar_5to2] !=
           1) ||
      (finder_parameters_
               .included_multi[IncludedMultiParticleReactions::NNbar_5to2] ==
           1 &&
       finder_parameters_.nnbar_treatment != NNbarTreatment::TwoToFive)) {
    throw std::invalid_argument(
        "In order to conserve detailed balance, when \"NNbar_5to2\" is "
        "included in\n`Multi_Particle_Reactions`, the `NNbarTreatment` has to "
        "be set to \"two to five\" and vice versa.");
  }
 
  if (finder_parameters_.nnbar_treatment == NNbarTreatment::Resonances &&
      finder_parameters_.included_2to2[IncludedReactions::NNbar] != 1) {
    throw std::invalid_argument(
        "'NNbar' has to be in the list of allowed 2 to 2 processes "
        "to enable annihilation to go through resonances");
  }
 
  if (finder_parameters_.strings_switch) {
    auto subconfig = config.extract_sub_configuration(
        {"Collision_Term", "String_Parameters"}, Configuration::GetEmpty::Yes);
    string_process_interface_ = std::make_unique<StringProcess>(
        subconfig.take({"String_Tension"}, 1.0), string_formation_time_,
        subconfig.take({"Gluon_Beta"}, 0.5),
        subconfig.take({"Gluon_Pmin"}, 0.001),
        subconfig.take({"Quark_Alpha"}, 2.0),
        subconfig.take({"Quark_Beta"}, 7.0),
        subconfig.take({"Strange_Supp"}, 0.16),
        subconfig.take({"Diquark_Supp"}, 0.036),
        subconfig.take({"Sigma_Perp"}, 0.42),
        subconfig.take({"StringZ_A_Leading"}, 0.2),
        subconfig.take({"StringZ_B_Leading"}, 2.0),
        subconfig.take({"StringZ_A"}, 2.0), subconfig.take({"StringZ_B"}, 0.55),
        subconfig.take({"String_Sigma_T"}, 0.5),
        subconfig.take({"Form_Time_Factor"}, 1.0),
        subconfig.take({"Mass_Dependent_Formation_Times"}, false),
        subconfig.take({"Prob_proton_to_d_uu"}, 1. / 3.),
        subconfig.take({"Separate_Fragment_Baryon"}, true),
        subconfig.take({"Popcorn_Rate"}, 0.15),
        subconfig.take({"Use_Monash_Tune"},
                       parameters.use_monash_tune_default.value()));
  }
}

Member Function Documentation

collision_time()

double smash::ScatterActionsFinder::collision_time ( const ParticleData &  p1, const ParticleData &  p2, double  dt, const std::vector< FourVector > &  beam_momentum  ) const

inline

Determine the collision time of the two particles.

Time of the closest approach is taken as collision time, if the geometric collision criterion is used. For stochastic criterion the time is distributed uniformly within the timestep.

Parameters

[in]	p1	First incoming particle
[in]	p2	Second incoming particle
[in]	dt	The maximum time interval at the current time step [fm]
[in]	beam_momentum	[GeV] List of beam momenta for each particle; only necessary for frozen Fermi motion

Returns

How long does it take for the two incoming particles to propagate before scattering [fm]. It's set equal to -1 if the two particles are not moving relative to each other.

JAM collision times from the closest approach in the two-particle center-of-mass-framem, see Hirano:2012yy [26] (5.13) and (5.14). The scatteraction is performed at the mean of these two times.

UrQMD collision time in computational frame, see Bass:1998ca [5] (3.28): position of particle 1: \(r_1\) [fm] position of particle 2: \(r_2\) [fm] velocity of particle 1: \(v_1\) velocity of particle 1: \(v_2\)

\[t_{coll} = - (r_1 - r_2) . (v_1 - v_2) / (v_1 - v_2)^2\]

[fm]

Definition at line 66 of file scatteractionsfinder.h.

                                                        {
    if (finder_parameters_.coll_crit == CollisionCriterion::Stochastic) {
      return dt * random::uniform(0., 1.);
    } else {
      /*
       * For frozen Fermi motion:
       * If particles have not yet interacted and are the initial nucleons,
       * perform action finding with beam momentum instead of Fermi motion
       * corrected momentum. That is because the particles are propagated with
       * the beam momentum until they interact.
       */
      if (p1.id() < 0 || p2.id() < 0) {
        throw std::runtime_error("Invalid particle ID for Fermi motion");
      }
      const bool p1_has_no_prior_interactions =
          (static_cast<uint64_t>(p1.id()) <                 // particle from
           static_cast<uint64_t>(beam_momentum.size())) &&  // initial nucleus
          (p1.get_history().collisions_per_particle == 0);
 
      const bool p2_has_no_prior_interactions =
          (static_cast<uint64_t>(p2.id()) <                 // particle from
           static_cast<uint64_t>(beam_momentum.size())) &&  // initial nucleus
          (p2.get_history().collisions_per_particle == 0);
 
      const FourVector p1_mom = (p1_has_no_prior_interactions)
                                    ? beam_momentum[p1.id()]
                                    : p1.momentum();
      const FourVector p2_mom = (p2_has_no_prior_interactions)
                                    ? beam_momentum[p2.id()]
                                    : p2.momentum();
      if (finder_parameters_.coll_crit == CollisionCriterion::Covariant) {
        const FourVector delta_x = p1.position() - p2.position();
        const double p1_sqr = p1_mom.sqr();
        const double p2_sqr = p2_mom.sqr();
        const double p1_dot_x = p1_mom.Dot(delta_x);
        const double p2_dot_x = p2_mom.Dot(delta_x);
        const double p1_dot_p2 = p1_mom.Dot(p2_mom);
        const double denominator = std::pow(p1_dot_p2, 2) - p1_sqr * p2_sqr;
        if (unlikely(std::abs(denominator) < really_small * really_small)) {
          return -1.0;
        }
 
        const double time_1 = (p2_sqr * p1_dot_x - p1_dot_p2 * p2_dot_x) *
                              p1_mom.x0() / denominator;
        const double time_2 = -(p1_sqr * p2_dot_x - p1_dot_p2 * p1_dot_x) *
                              p2_mom.x0() / denominator;
        return (time_1 + time_2) / 2;
      } else {
        const ThreeVector dv_times_e1e2 =
            p1_mom.threevec() * p2_mom.x0() - p2_mom.threevec() * p1_mom.x0();
        const double dv_times_e1e2_sqr = dv_times_e1e2.sqr();
        if (dv_times_e1e2_sqr < really_small) {
          return -1.0;
        }
        const ThreeVector dr =
            p1.position().threevec() - p2.position().threevec();
        return -(dr * dv_times_e1e2) *
               (p1_mom.x0() * p2_mom.x0() / dv_times_e1e2_sqr);
      }
    }
  }

find_actions_in_cell()

ActionList smash::ScatterActionsFinder::find_actions_in_cell ( const ParticleList &  search_list, double  dt, const double  gcell_vol, const std::vector< FourVector > &  beam_momentum  ) const

overridevirtual

Search for all the possible collisions within one cell.

This function is only used for counting the primary collisions at the beginning of each time step. (Although it's also called afterwards for searching the secondary collisions among the outgoing particles, no new actions will be found since the scattered pairs cannot scatter again.)

Parameters

[in]	search_list	A list of particles within one cell
[in]	dt	The maximum time interval at the current time step [fm]
[in]	gcell_vol	Volume of searched grid cell [fm^3]
[in]	beam_momentum	[GeV] List of beam momenta for each particle; only necessary for frozen Fermi motion

Returns

A list of possible scatter actions

Implements smash::ActionFinderInterface.

Definition at line 445 of file scatteractionsfinder.cc.

                                                      {
  std::vector<ActionPtr> actions;
  for (const ParticleData& p1 : search_list) {
    for (const ParticleData& p2 : search_list) {
      // Check for 2 particle scattering
      if (p1.id() < p2.id()) {
        ActionPtr act =
            check_collision_two_part(p1, p2, dt, beam_momentum, gcell_vol);
        if (act) {
          actions.push_back(std::move(act));
        }
      }
      if (finder_parameters_.included_multi.any()) {
        // Also, check for 3 particle scatterings with stochastic criterion
        for (const ParticleData& p3 : search_list) {
          if (finder_parameters_.included_multi
                      [IncludedMultiParticleReactions::Deuteron_3to2] == 1 ||
              finder_parameters_.included_multi
                      [IncludedMultiParticleReactions::Meson_3to1] == 1) {
            if (p1.id() < p2.id() && p2.id() < p3.id()) {
              ActionPtr act =
                  check_collision_multi_part({p1, p2, p3}, dt, gcell_vol);
              if (act) {
                actions.push_back(std::move(act));
              }
            }
          }
          for (const ParticleData& p4 : search_list) {
            if (finder_parameters_.included_multi
                    [IncludedMultiParticleReactions::A3_Nuclei_4to2]) {
              if (p1.id() < p2.id() && p2.id() < p3.id() && p3.id() < p4.id()) {
                ActionPtr act =
                    check_collision_multi_part({p1, p2, p3, p4}, dt, gcell_vol);
                if (act) {
                  actions.push_back(std::move(act));
                }
              }
            }
            if (finder_parameters_.included_multi
                        [IncludedMultiParticleReactions::NNbar_5to2] == 1 &&
                search_list.size() >= 5) {
              for (const ParticleData& p5 : search_list) {
                if ((p1.id() < p2.id() && p2.id() < p3.id() &&
                     p3.id() < p4.id() && p4.id() < p5.id()) &&
                    (p1.is_pion() && p2.is_pion() && p3.is_pion() &&
                     p4.is_pion() && p5.is_pion())) {
                  // at the moment only pure pion 5-body reactions
                  ActionPtr act = check_collision_multi_part(
                      {p1, p2, p3, p4, p5}, dt, gcell_vol);
                  if (act) {
                    actions.push_back(std::move(act));
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  return actions;
}

find_actions_with_neighbors()

ActionList smash::ScatterActionsFinder::find_actions_with_neighbors ( const ParticleList &  search_list, const ParticleList &  neighbors_list, double  dt, const std::vector< FourVector > &  beam_momentum  ) const

overridevirtual

Search for all the possible collisions among the neighboring cells.

This function is only used for counting the primary collisions at the beginning of each time step.

Parameters

[in]	search_list	A list of particles within the current cell
[in]	neighbors_list	A list of particles within the neighboring cell
[in]	dt	The maximum time interval at the current time step [fm]
[in]	beam_momentum	[GeV] List of beam momenta for each particle; only necessary for frozen Fermi motion

Returns

A list of possible scatter actions

Implements smash::ActionFinderInterface.

Definition at line 510 of file scatteractionsfinder.cc.

                                                                 {
  std::vector<ActionPtr> actions;
  if (finder_parameters_.coll_crit == CollisionCriterion::Stochastic) {
    // Only search in cells
    return actions;
  }
  for (const ParticleData& p1 : search_list) {
    for (const ParticleData& p2 : neighbors_list) {
      assert(p1.id() != p2.id());
      // Check if a collision is possible.
      ActionPtr act = check_collision_two_part(p1, p2, dt, beam_momentum);
      if (act) {
        actions.push_back(std::move(act));
      }
    }
  }
  return actions;
}

find_actions_with_surrounding_particles()

ActionList smash::ScatterActionsFinder::find_actions_with_surrounding_particles ( const ParticleList &  search_list, const Particles &  surrounding_list, double  dt, const std::vector< FourVector > &  beam_momentum  ) const

overridevirtual

Search for all the possible secondary collisions between the outgoing particles and the rest.

Parameters

[in]	search_list	A list of particles within the current cell
[in]	surrounding_list	The whole particle list
[in]	dt	The maximum time interval at the current time step [fm]
[in]	beam_momentum	[GeV] List of beam momenta for each particle; only necessary for frozen Fermi motion

Returns

A list of possible scatter actions

Implements smash::ActionFinderInterface.

Definition at line 531 of file scatteractionsfinder.cc.

                                                                 {
  std::vector<ActionPtr> actions;
  if (finder_parameters_.coll_crit == CollisionCriterion::Stochastic) {
    // Only search in cells
    return actions;
  }
  for (const ParticleData& p2 : surrounding_list) {
    /* don't look for collisions if the particle from the surrounding list is
     * also in the search list */
    auto result = std::find_if(
        search_list.begin(), search_list.end(),
        [&p2](const ParticleData& p) { return p.id() == p2.id(); });
    if (result != search_list.end()) {
      continue;
    }
    for (const ParticleData& p1 : search_list) {
      // Check if a collision is possible.
      ActionPtr act = check_collision_two_part(p1, p2, dt, beam_momentum);
      if (act) {
        actions.push_back(std::move(act));
      }
    }
  }
  return actions;
}

find_final_actions()

ActionList smash::ScatterActionsFinder::find_final_actions ( const Particles &  , bool  = false  ) const

inlineoverridevirtual

Find some final collisions at the end of the simulation.

Todo:

Seems to do nothing.

Implements smash::ActionFinderInterface.

Definition at line 198 of file scatteractionsfinder.h.

                                                               {
    return ActionList();
  }

is_constant_elastic_isotropic()

bool smash::ScatterActionsFinder::is_constant_elastic_isotropic ( ) const

inline

If there is only one particle sort, no decays (only elastic scatterings are possible), scatterings are isotropic and cross-section fixed to elastic_parameter_ independently on momenta, then maximal cross-section is elastic_parameter_.

This knowledge can be used for improving performance.

Returns

A boolean indicating whether all the scatterings are elastic and isotropic

Definition at line 213 of file scatteractionsfinder.h.

                                                    {
    return ParticleType::list_all().size() == 1 &&
           !finder_parameters_.two_to_one && isotropic_ &&
           finder_parameters_.elastic_parameter > 0.;
  }

max_transverse_distance_sqr()

double smash::ScatterActionsFinder::max_transverse_distance_sqr ( int  testparticles ) const

inline

The maximal distance over which particles can interact in case of the geometric criterion, related to the number of test particles and the maximal cross section.

Parameters

[in]

testparticles

Number of test particles.

Returns

Maximal transverse distance squared. [fm \(^{2}\)] Particle pairs whose transverse distance is larger than this are not checked for collisions in case of the default geometric collision criterion.

Definition at line 231 of file scatteractionsfinder.h.

                                                              {
    return (is_constant_elastic_isotropic()
                ? finder_parameters_.elastic_parameter
                : finder_parameters_.maximum_cross_section) /
           testparticles * fm2_mb * M_1_PI;
  }

dump_reactions()

void smash::ScatterActionsFinder::dump_reactions

(

)

const

Prints out all the 2-> n (n > 1) reactions with non-zero cross-sections between all possible pairs of particle types.

Definition at line 559 of file scatteractionsfinder.cc.

                                                {
  constexpr double time = 0.0;
 
  const size_t N_isotypes = IsoParticleType::list_all().size();
  const size_t N_pairs = N_isotypes * (N_isotypes - 1) / 2;
 
  std::cout << N_isotypes << " iso-particle types." << std::endl;
  std::cout << "They can make " << N_pairs << " pairs." << std::endl;
  std::vector<double> momentum_scan_list = {0.1, 0.3, 0.5, 1.0,
                                            2.0, 3.0, 5.0, 10.0};
  for (const IsoParticleType& A_isotype : IsoParticleType::list_all()) {
    for (const IsoParticleType& B_isotype : IsoParticleType::list_all()) {
      if (&A_isotype > &B_isotype) {
        continue;
      }
      bool any_nonzero_cs = false;
      std::vector<std::string> r_list;
      for (const ParticleTypePtr A_type : A_isotype.get_states()) {
        for (const ParticleTypePtr B_type : B_isotype.get_states()) {
          if (A_type > B_type) {
            continue;
          }
          ParticleData A(*A_type), B(*B_type);
          for (auto mom : momentum_scan_list) {
            A.set_4momentum(A.pole_mass(), mom, 0.0, 0.0);
            B.set_4momentum(B.pole_mass(), -mom, 0.0, 0.0);
            ScatterActionPtr act = std::make_unique<ScatterAction>(
                A, B, time, isotropic_, string_formation_time_, -1, false);
            if (finder_parameters_.strings_switch) {
              act->set_string_interface(string_process_interface_.get());
            }
            act->add_all_scatterings(finder_parameters_);
            const double total_cs = act->cross_section();
            if (total_cs <= 0.0) {
              continue;
            }
            any_nonzero_cs = true;
            for (const auto& channel : act->collision_channels()) {
              const auto type = channel->get_type();
              std::string r;
              if (is_string_soft_process(type) ||
                  type == ProcessType::StringHard) {
                r = A_type->name() + B_type->name() + std::string(" → strings");
              } else {
                std::string r_type =
                    (type == ProcessType::Elastic) ? std::string(" (el)")
                    : (channel->get_type() == ProcessType::TwoToTwo)
                        ? std::string(" (inel)")
                        : std::string(" (?)");
                r = A_type->name() + B_type->name() + std::string(" → ") +
                    channel->particle_types()[0]->name() +
                    channel->particle_types()[1]->name() + r_type;
              }
              isoclean(r);
              r_list.push_back(r);
            }
          }
        }
      }
      std::sort(r_list.begin(), r_list.end());
      r_list.erase(std::unique(r_list.begin(), r_list.end()), r_list.end());
      if (any_nonzero_cs) {
        for (auto r : r_list) {
          std::cout << r;
          if (r_list.back() != r) {
            std::cout << ", ";
          }
        }
        std::cout << std::endl;
      }
    }
  }
}

dump_cross_sections()

void smash::ScatterActionsFinder::dump_cross_sections	(	const ParticleType &	a,
		const ParticleType &	b,
		double	m_a,
		double	m_b,
		bool	final_state,
		std::vector< double > &	plab
	)		const

Print out partial cross-sections of all processes that can occur in the collision of a(mass = m_a) and b(mass = m_b).

Parameters

[in]	a	The specie of the first incoming particle.
[in]	b	The specie of the second incoming particle.
[in]	m_a	Mass of species a [GeV].
[in]	m_b	Mass of species b [GeV].
[in]	final_state	Whether the final state cross sections should be printed.
[in]	plab	Optional momenta in lab frame to be evaluated [GeV]. Ignored if empty.

Definition at line 942 of file scatteractionsfinder.cc.

                                                     {
  typedef std::vector<std::pair<double, double>> xs_saver;
  std::map<std::string, xs_saver> xs_dump;
  std::map<std::string, double> outgoing_total_mass;
 
  ParticleData a_data(a), b_data(b);
  int n_momentum_points = 200;
  constexpr double momentum_step = 0.02;
  if (plab.size() > 0) {
    n_momentum_points = plab.size();
    // Remove duplicates.
    std::sort(plab.begin(), plab.end());
    plab.erase(std::unique(plab.begin(), plab.end()), plab.end());
  }
  for (int i = 0; i < n_momentum_points; i++) {
    double momentum;
    if (plab.size() > 0) {
      momentum = pCM_from_s(s_from_plab(plab.at(i), m_a, m_b), m_a, m_b);
    } else {
      momentum = momentum_step * (i + 1);
    }
    a_data.set_4momentum(m_a, momentum, 0.0, 0.0);
    b_data.set_4momentum(m_b, -momentum, 0.0, 0.0);
    const double sqrts = (a_data.momentum() + b_data.momentum()).abs();
    const ParticleList incoming = {a_data, b_data};
    ScatterActionPtr act = std::make_unique<ScatterAction>(
        a_data, b_data, 0., isotropic_, string_formation_time_, -1, false);
    if (finder_parameters_.strings_switch) {
      act->set_string_interface(string_process_interface_.get());
    }
    act->add_all_scatterings(finder_parameters_);
    decaytree::Node tree(a.name() + b.name(), act->cross_section(), {&a, &b},
                         {&a, &b}, {&a, &b}, {});
    const CollisionBranchList& processes = act->collision_channels();
    for (const auto& process : processes) {
      const double xs = process->weight();
      if (xs <= 0.0) {
        continue;
      }
      if (!final_state) {
        std::stringstream process_description_stream;
        process_description_stream << *process;
        const std::string& description = process_description_stream.str();
        double m_tot = 0.0;
        for (const auto& ptype : process->particle_types()) {
          m_tot += ptype->mass();
        }
        outgoing_total_mass[description] = m_tot;
        if (!xs_dump[description].empty() &&
            std::abs(xs_dump[description].back().first - sqrts) <
                really_small) {
          xs_dump[description].back().second += xs;
        } else {
          xs_dump[description].push_back(std::make_pair(sqrts, xs));
        }
      } else {
        std::stringstream process_description_stream;
        process_description_stream << *process;
        const std::string& description = process_description_stream.str();
        ParticleTypePtrList initial_particles = {&a, &b};
        ParticleTypePtrList final_particles = process->particle_types();
        auto& process_node =
            tree.add_action(description, xs, std::move(initial_particles),
                            std::move(final_particles));
        decaytree::add_decays(process_node, sqrts);
      }
    }
    xs_dump["total"].push_back(std::make_pair(sqrts, act->cross_section()));
    // Total cross-section should be the first in the list -> negative mass
    outgoing_total_mass["total"] = -1.0;
    if (final_state) {
      // tree.print();
      auto final_state_xs = tree.final_state_cross_sections();
      deduplicate(final_state_xs);
      for (const auto& p : final_state_xs) {
        // Don't print empty columns.
        //
        // FIXME(steinberg): The better fix would be to not have them in the
        // first place.
        if (p.name_ == "") {
          continue;
        }
        outgoing_total_mass[p.name_] = p.mass_;
        xs_dump[p.name_].push_back(std::make_pair(sqrts, p.cross_section_));
      }
    }
  }
  // Get rid of cross sections that are zero.
  // (This only happens if their is a resonance in the final state that cannot
  // decay with our simplified assumptions.)
  for (auto it = begin(xs_dump); it != end(xs_dump);) {
    // Sum cross section over all energies.
    const xs_saver& xs = (*it).second;
    double sum = 0;
    for (const auto& p : xs) {
      sum += p.second;
    }
    if (sum == 0.) {
      it = xs_dump.erase(it);
    } else {
      ++it;
    }
  }
 
  // Nice ordering of channels by summed pole mass of products
  std::vector<std::string> all_channels;
  for (const auto& channel : xs_dump) {
    all_channels.push_back(channel.first);
  }
  std::sort(all_channels.begin(), all_channels.end(),
            [&](const std::string& str_a, const std::string& str_b) {
              return outgoing_total_mass[str_a] < outgoing_total_mass[str_b];
            });
 
  // Print header
  std::cout << "# Dumping partial " << a.name() << b.name()
            << " cross-sections in mb, energies in GeV" << std::endl;
  std::cout << "   sqrt_s";
  // Align everything to 16 unicode characters.
  // This should be enough for the longest channel name (7 final-state
  // particles).
  for (const auto& channel : all_channels) {
    std::cout << utf8::fill_left(channel, 16, ' ');
  }
  std::cout << std::endl;
 
  // Print out all partial cross-sections in mb
  for (int i = 0; i < n_momentum_points; i++) {
    double momentum;
    if (plab.size() > 0) {
      momentum = pCM_from_s(s_from_plab(plab.at(i), m_a, m_b), m_a, m_b);
    } else {
      momentum = momentum_step * (i + 1);
    }
    a_data.set_4momentum(m_a, momentum, 0.0, 0.0);
    b_data.set_4momentum(m_b, -momentum, 0.0, 0.0);
    const double sqrts = (a_data.momentum() + b_data.momentum()).abs();
    std::printf("%9.6f", sqrts);
    for (const auto& channel : all_channels) {
      const xs_saver energy_and_xs = xs_dump[channel];
      size_t j = 0;
      for (; j < energy_and_xs.size() && energy_and_xs[j].first < sqrts; j++) {
      }
      double xs = 0.0;
      if (j < energy_and_xs.size() &&
          std::abs(energy_and_xs[j].first - sqrts) < really_small) {
        xs = energy_and_xs[j].second;
      }
      std::printf("%16.6f", xs);  // Same alignment as in the header.
    }
    std::printf("\n");
  }
}

get_process_string_ptr()

StringProcess* smash::ScatterActionsFinder::get_process_string_ptr ( )

inline

Returns

Pointer to the string process class object. If string is turned off, the null pointer is returned.

Definition at line 265 of file scatteractionsfinder.h.

                                          {
    if (finder_parameters_.strings_switch) {
      return string_process_interface_.get();
    } else {
      return NULL;
    }
  }

check_collision_two_part()

ActionPtr smash::ScatterActionsFinder::check_collision_two_part ( const ParticleData &  data_a, const ParticleData &  data_b, double  dt, const std::vector< FourVector > &  beam_momentum = {}, const double  gcell_vol = 0.0  ) const

private

Check for a single pair of particles (id_a, id_b) if a collision will happen in the next timestep and create a corresponding Action object in that case.

Two criteria for the collision decision are supported: 1. The default geometric criterion from UrQMD Bass:1998ca [5] (3.27). 2. A stochastic collision criterion as introduced in Staudenmaier:2021lrg [56].

Parameters

[in]	data_a	First incoming particle
[in]	data_b	Second incoming particle
[in]	dt	Maximum time interval within which a collision can happen
[in]	beam_momentum	[GeV] List of beam momenta for each particle; only necessary for frozen Fermi motion
[in]	gcell_vol	(optional) volume of grid cell in which the collision is checked

Returns

A null pointer if no collision happens or an action which contains the information of the outgoing particles.

Note: gcell_vol is optional, since only find_actions_in_cell has (and needs) this information for the stochastic collision criterion.

Definition at line 211 of file scatteractionsfinder.cc.

                                  {
  /* If the two particles
   * 1) belong to one of the two colliding nuclei, and
   * 2) both of them have never experienced any collisions,
   * then the collisions between them are banned. */
  if (!finder_parameters_.allow_collisions_within_nucleus) {
    assert(data_a.id() >= 0);
    assert(data_b.id() >= 0);
    bool in_same_nucleus = (data_a.belongs_to() == BelongsTo::Projectile &&
                            data_b.belongs_to() == BelongsTo::Projectile) ||
                           (data_a.belongs_to() == BelongsTo::Target &&
                            data_b.belongs_to() == BelongsTo::Target);
    bool never_interacted_before =
        data_a.get_history().collisions_per_particle == 0 &&
        data_b.get_history().collisions_per_particle == 0;
    if (in_same_nucleus && never_interacted_before) {
      return nullptr;
    }
  }
 
  // No grid or search in cell means no collision for stochastic criterion
  if (finder_parameters_.coll_crit == CollisionCriterion::Stochastic &&
      gcell_vol < really_small) {
    return nullptr;
  }
 
  // Determine time of collision.
  const double time_until_collision =
      collision_time(data_a, data_b, dt, beam_momentum);
 
  // Check that collision happens in this timestep.
  if (time_until_collision < 0. || time_until_collision >= dt) {
    return nullptr;
  }
 
  // Determine which total cross section to use
  bool incoming_parametrized = (finder_parameters_.total_xs_strategy ==
                                TotalCrossSectionStrategy::TopDown);
  if (finder_parameters_.total_xs_strategy ==
      TotalCrossSectionStrategy::TopDownMeasured) {
    const PdgCode& pdg_a = data_a.type().pdgcode();
    const PdgCode& pdg_b = data_b.type().pdgcode();
    incoming_parametrized = parametrization_exists(pdg_a, pdg_b);
  }
 
  // Create ScatterAction object.
  ScatterActionPtr act = std::make_unique<ScatterAction>(
      data_a, data_b, time_until_collision, isotropic_, string_formation_time_,
      box_length_, incoming_parametrized);
 
  if (finder_parameters_.coll_crit == CollisionCriterion::Stochastic) {
    act->set_stochastic_pos_idx();
  }
 
  if (finder_parameters_.strings_switch) {
    act->set_string_interface(string_process_interface_.get());
  }
 
  // Distance squared calculation not needed for stochastic criterion
  const double distance_squared =
      (finder_parameters_.coll_crit == CollisionCriterion::Geometric)
          ? act->transverse_distance_sqr()
      : (finder_parameters_.coll_crit == CollisionCriterion::Covariant)
          ? act->cov_transverse_distance_sqr()
          : 0.0;
 
  // Don't calculate cross section if the particles are very far apart.
  // Not needed for stochastic criterion because of cell structure.
  if (finder_parameters_.coll_crit != CollisionCriterion::Stochastic &&
      distance_squared >=
          max_transverse_distance_sqr(finder_parameters_.testparticles)) {
    return nullptr;
  }
 
  if (incoming_parametrized) {
    act->set_parametrized_total_cross_section(finder_parameters_);
  } else {
    // Add various subprocesses.
    act->add_all_scatterings(finder_parameters_);
  }
 
  double xs = act->cross_section() * fm2_mb /
              static_cast<double>(finder_parameters_.testparticles);
 
  // Take cross section scaling factors into account
  xs *= data_a.xsec_scaling_factor(time_until_collision);
  xs *= data_b.xsec_scaling_factor(time_until_collision);
 
  if (finder_parameters_.coll_crit == CollisionCriterion::Stochastic) {
    const double v_rel = act->relative_velocity();
    /* Collision probability for 2-particle scattering, see
     * \iref{Staudenmaier:2021lrg}. */
    const double prob = xs * v_rel * dt / gcell_vol;
 
    logg[LFindScatter].debug(
        "Stochastic collison criterion parameters (2-particles):\nprob = ",
        prob, ", xs = ", xs, ", v_rel = ", v_rel, ", dt = ", dt,
        ", gcell_vol = ", gcell_vol,
        ", testparticles = ", finder_parameters_.testparticles);
 
    if (prob > 1.) {
      std::stringstream err;
      err << "Probability larger than 1 for stochastic rates. ( P_22 = " << prob
          << " )\n"
          << data_a.type().name() << data_b.type().name() << " with masses "
          << data_a.momentum().abs() << " and " << data_b.momentum().abs()
          << " at sqrts[GeV] = " << act->sqrt_s()
          << " with xs[fm^2]/Ntest = " << xs
          << "\nConsider using smaller timesteps.";
      if (finder_parameters_.only_warn_for_high_prob) {
        logg[LFindScatter].warn(err.str());
      } else {
        throw std::runtime_error(err.str());
      }
    }
 
    // probability criterion
    double random_no = random::uniform(0., 1.);
    if (random_no > prob) {
      return nullptr;
    }
 
  } else if (finder_parameters_.coll_crit == CollisionCriterion::Geometric ||
             finder_parameters_.coll_crit == CollisionCriterion::Covariant) {
    // just collided with this particle
    if (data_a.id_process() > 0 && data_a.id_process() == data_b.id_process()) {
      logg[LFindScatter].debug("Skipping collided particles at time ",
                               data_a.position().x0(), " due to process ",
                               data_a.id_process(), "\n    ", data_a, "\n<-> ",
                               data_b);
 
      return nullptr;
    }
 
    // Cross section for collision criterion
    const double cross_section_criterion = xs * M_1_PI;
 
    // distance criterion according to cross_section
    if (distance_squared >= cross_section_criterion) {
      return nullptr;
    }
 
    logg[LFindScatter].debug("particle distance squared: ", distance_squared,
                             "\n    ", data_a, "\n<-> ", data_b);
  }
 
  // Include possible outgoing branches
  if (incoming_parametrized) {
    act->add_all_scatterings(finder_parameters_);
  }
 
  return act;
}

check_collision_multi_part()

ActionPtr smash::ScatterActionsFinder::check_collision_multi_part ( const ParticleList &  plist, double  dt, const double  gcell_vol  ) const

private

Check for multiple i.e.

more than 2 particles if a collision will happen in the next timestep and create a corresponding Action object in that case.

This is only possible for the stochastic collision criterion, which is introduced in Staudenmaier:2021lrg [56]. Following the same general idea as for the 2-particle scatterings, probabilities for multi-particle scatterings can be derived.

Parameters

[in]	plist	List of incoming particles
[in]	dt	Maximum time interval within which a collision can happen
[in]	gcell_vol	volume of grid cell in which the collision is checked

Returns

A null pointer if no collision happens or an action which contains the information of the outgoing particles.

Definition at line 368 of file scatteractionsfinder.cc.

                                                                        {
  /* If all particles
   * 1) belong to the two colliding nuclei
   * 2) are within the same nucleus
   * 3) have never experienced any collisons,
   * then the collision between them are banned also for multi-particle
   * interactions. */
  if (!finder_parameters_.allow_collisions_within_nucleus) {
    bool all_projectile =
        std::all_of(plist.begin(), plist.end(), [&](const ParticleData& data) {
          return data.belongs_to() == BelongsTo::Projectile;
        });
    bool all_target =
        std::all_of(plist.begin(), plist.end(), [&](const ParticleData& data) {
          return data.belongs_to() == BelongsTo::Target;
        });
    bool none_collided =
        std::all_of(plist.begin(), plist.end(), [&](const ParticleData& data) {
          return data.get_history().collisions_per_particle == 0;
        });
    if ((all_projectile || all_target) && none_collided) {
      return nullptr;
    }
  }
  // No grid or search in cell
  if (gcell_vol < really_small) {
    return nullptr;
  }
 
  /* Optimisation for later: Already check here at the beginning
   * if collision with plist is possible before constructing actions. */
 
  // 1. Determine time of collision.
  const double time_until_collision = dt * random::uniform(0., 1.);
 
  // 2. Create ScatterAction object.
  ScatterActionMultiPtr act =
      std::make_unique<ScatterActionMulti>(plist, time_until_collision);
 
  act->set_stochastic_pos_idx();
 
  // 3. Add possible final states (dt and gcell_vol for probability calculation)
  act->add_possible_reactions(dt, gcell_vol, finder_parameters_.included_multi);
 
  /* 4. Return total collision probability
   *    Scales with 1 over the number of testpartciles to the power of the
   *    number of incoming particles - 1 */
  const double prob =
      act->get_total_weight() /
      std::pow(finder_parameters_.testparticles, plist.size() - 1);
 
  // 5. Check that probability is smaller than one
  if (prob > 1.) {
    std::stringstream err;
    err << "Probability " << prob << " larger than 1 for stochastic rates for ";
    for (const ParticleData& data : plist) {
      err << data.type().name();
    }
    err << " at sqrts[GeV] = " << act->sqrt_s()
        << "\nConsider using smaller timesteps.";
    if (finder_parameters_.only_warn_for_high_prob) {
      logg[LFindScatter].warn(err.str());
    } else {
      throw std::runtime_error(err.str());
    }
  }
 
  // 6. Perform probability decisions
  double random_no = random::uniform(0., 1.);
  if (random_no > prob) {
    return nullptr;
  }
 
  return act;
}

Member Data Documentation

finder_parameters_

ScatterActionsFinderParameters smash::ScatterActionsFinder::finder_parameters_

private

Struct collecting several parameters.

Definition at line 320 of file scatteractionsfinder.h.

string_process_interface_

std::unique_ptr<StringProcess> smash::ScatterActionsFinder::string_process_interface_

private

Class that deals with strings, interfacing Pythia.

Definition at line 322 of file scatteractionsfinder.h.

isotropic_

const bool smash::ScatterActionsFinder::isotropic_

private

Do all collisions isotropically.

Definition at line 324 of file scatteractionsfinder.h.

box_length_

const double smash::ScatterActionsFinder::box_length_

private

Box length: needed to determine coordinates of collision correctly in case of collision through the wall.

Ignored if negative.

Definition at line 330 of file scatteractionsfinder.h.

string_formation_time_

const double smash::ScatterActionsFinder::string_formation_time_

private

Parameter for formation time.

Definition at line 332 of file scatteractionsfinder.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/scatteractionsfinder.h
• src/scatteractionsfinder.cc