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 30 of file scatteractionsfinder.h.

Inheritance diagram for smash::ScatterActionsFinder:

Collaboration 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
std::unique_ptr< StringProcess >	string_process_interface_
	Class that deals with strings, interfacing Pythia. More...
const CollisionCriterion	coll_crit_
	Specifies which collision criterion is used. More...
const double	elastic_parameter_
	Elastic cross section parameter (in mb). More...
const int	testparticles_
	Number of test particles. More...
const bool	isotropic_
	Do all collisions isotropically. More...
const bool	two_to_one_
	Enable 2->1 processes. More...
const ReactionsBitSet	incl_set_
	List of included 2<->2 reactions. More...
const MultiParticleReactionsBitSet	incl_multi_set_
	List of included multi-particle reactions. More...
const double	scale_xs_
	Factor by which all (partial) cross sections are scaled. More...
const double	additional_el_xs_
	Additional constant elastic cross section. More...
const double	low_snn_cut_
	Elastic collsions between two nucleons with sqrt_s below low_snn_cut_ are excluded. More...
const bool	strings_switch_
	Switch to turn off string excitation. More...
const bool	use_AQM_
	Switch to control whether to use AQM or not. More...
const bool	strings_with_probability_
	Decide whether to implement string fragmentation based on a probability. More...
const NNbarTreatment	nnbar_treatment_
	Switch for NNbar reactions. 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...
const double	maximum_cross_section_
const bool	allow_first_collisions_within_nucleus_
	If particles within nucleus are allowed to collide for their first time. More...
const bool	only_warn_for_high_prob_
	Switch to turn off throwing an exception for collision probabilities larger than 1. More...

Constructor & Destructor Documentation

ScatterActionsFinder()

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

Constructor of the finder with the given parameters.

Parameters

[in]	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 308 of file scatteractionsfinder.cc.

    : coll_crit_(parameters.coll_crit),
      elastic_parameter_(
          config.take({"Collision_Term", "Elastic_Cross_Section"}, -1.)),
      testparticles_(parameters.testparticles),
      isotropic_(config.take({"Collision_Term", "Isotropic"}, false)),
      two_to_one_(parameters.two_to_one),
      incl_set_(parameters.included_2to2),
      incl_multi_set_(parameters.included_multi),
      scale_xs_(parameters.scale_xs),
      additional_el_xs_(parameters.additional_el_xs),
      low_snn_cut_(parameters.low_snn_cut),
      strings_switch_(parameters.strings_switch),
      use_AQM_(parameters.use_AQM),
      strings_with_probability_(parameters.strings_with_probability),
      nnbar_treatment_(parameters.nnbar_treatment),
      box_length_(parameters.box_length),
      string_formation_time_(config.take(
          {"Collision_Term", "String_Parameters", "Formation_Time"}, 1.)),
      maximum_cross_section_(parameters.maximum_cross_section),
      allow_first_collisions_within_nucleus_(
          parameters.allow_collisions_within_nucleus),
      only_warn_for_high_prob_(config.take(
          {"Collision_Term", "Only_Warn_For_High_Probability"}, false)) {
  if (is_constant_elastic_isotropic()) {
    logg[LFindScatter].info(
        "Constant elastic isotropic cross-section mode:", " using ",
        elastic_parameter_, " mb as maximal cross-section.");
  }
  if (incl_multi_set_.any() && 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 (incl_multi_set_[IncludedMultiParticleReactions::Deuteron_3to2] == 1 &&
      (incl_set_[IncludedReactions::PiDeuteron_to_pidprime] == 1 ||
       incl_set_[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 (incl_multi_set_[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 ((nnbar_treatment_ == NNbarTreatment::TwoToFive &&
       incl_multi_set_[IncludedMultiParticleReactions::NNbar_5to2] != 1) ||
      (incl_multi_set_[IncludedMultiParticleReactions::NNbar_5to2] == 1 &&
       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 (nnbar_treatment_ == NNbarTreatment::Resonances &&
      incl_set_[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 (strings_switch_) {
    auto subconfig = config["Collision_Term"]["String_Parameters"];
    string_process_interface_ = 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));
  }
}

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/c]. 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 [23] (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 [4] (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/c]

Definition at line 65 of file scatteractionsfinder.h.

                                                        {
    if (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 (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);
      }
    }
  }

Here is the call graph for this function:

Here is the caller graph for this function:

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 616 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 (incl_multi_set_.any()) {
        // Also, check for 3 particle scatterings with stochastic criterion
        for (const ParticleData& p3 : search_list) {
          if (incl_multi_set_[IncludedMultiParticleReactions::Deuteron_3to2] ==
                  1 ||
              incl_multi_set_[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 (incl_multi_set_
                    [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 (incl_multi_set_[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;
}

Here is the call graph for this function:

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/c]
[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 681 of file scatteractionsfinder.cc.

                                                                 {
  std::vector<ActionPtr> actions;
  if (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;
}

Here is the call graph for this function:

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/c]
[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 702 of file scatteractionsfinder.cc.

                                                                 {
  std::vector<ActionPtr> actions;
  if (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;
}

Here is the call graph for this function:

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 197 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 212 of file scatteractionsfinder.h.

                                                    {
    return ParticleType::list_all().size() == 1 && !two_to_one_ && isotropic_ &&
           elastic_parameter_ > 0.;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

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 229 of file scatteractionsfinder.h.

                                                              {
    return (is_constant_elastic_isotropic() ? elastic_parameter_
                                            : maximum_cross_section_) /
           testparticles * fm2_mb * M_1_PI;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

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 730 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 = make_unique<ScatterAction>(
                A, B, time, isotropic_, string_formation_time_);
            if (strings_switch_) {
              act->set_string_interface(string_process_interface_.get());
            }
            act->add_all_scatterings(
                elastic_parameter_, two_to_one_, incl_set_, incl_multi_set_,
                low_snn_cut_, strings_switch_, use_AQM_,
                strings_with_probability_, nnbar_treatment_, scale_xs_,
                additional_el_xs_);
            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;
      }
    }
  }
}

Here is the call graph for this function:

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 1118 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;
  /*
  // Round to output precision.
  for (auto& p : plab) {
    p = std::floor((p * 100000) + 0.5) / 100000;
  }
  */
  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 = make_unique<ScatterAction>(
        a_data, b_data, 0., isotropic_, string_formation_time_);
    if (strings_switch_) {
      act->set_string_interface(string_process_interface_.get());
    }
    act->add_all_scatterings(elastic_parameter_, two_to_one_, incl_set_,
                             incl_multi_set_, low_snn_cut_, strings_switch_,
                             use_AQM_, strings_with_probability_,
                             nnbar_treatment_, scale_xs_, additional_el_xs_);
    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");
  }
}

Here is the call graph for this function:

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 262 of file scatteractionsfinder.h.

                                          {
    if (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 [4] (3.27). 2. A stochastic collision criterion as introduced in Staudenmaier:2021lrg [50].

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 408 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 (!allow_first_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 (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;
  }
 
  // Create ScatterAction object.
  ScatterActionPtr act = make_unique<ScatterAction>(
      data_a, data_b, time_until_collision, isotropic_, string_formation_time_,
      box_length_);
 
  if (coll_crit_ == CollisionCriterion::Stochastic) {
    act->set_stochastic_pos_idx();
  }
 
  if (strings_switch_) {
    act->set_string_interface(string_process_interface_.get());
  }
 
  // Distance squared calculation not needed for stochastic criterion
  const double distance_squared =
      (coll_crit_ == CollisionCriterion::Geometric)
          ? act->transverse_distance_sqr()
          : (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 (coll_crit_ != CollisionCriterion::Stochastic &&
      distance_squared >= max_transverse_distance_sqr(testparticles_)) {
    return nullptr;
  }
 
  // Add various subprocesses.
  act->add_all_scatterings(elastic_parameter_, two_to_one_, incl_set_,
                           incl_multi_set_, low_snn_cut_, strings_switch_,
                           use_AQM_, strings_with_probability_,
                           nnbar_treatment_, scale_xs_, additional_el_xs_);
 
  double xs =
      act->cross_section() * fm2_mb / static_cast<double>(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 (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 = ", testparticles_);
 
    if (prob > 1.) {
      std::stringstream err;
      err << "Probability larger than 1 for stochastic rates. ( P_22 = " << prob
          << " )\nConsider using smaller timesteps.";
      if (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 (coll_crit_ == CollisionCriterion::Geometric ||
             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);
  }
 
  // Using std::move here is redundant with newer compilers, but required for
  // supporting GCC 4.8. Once we drop this support, std::move should be removed.
  return std::move(act);
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 [50]. 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 544 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 (!allow_first_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 =
      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, incl_multi_set_);
 
  /* 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(testparticles_, plist.size() - 1);
 
  // 5. Check that probability is smaller than one
  if (prob > 1.) {
    std::stringstream err;
    err << "Probability larger than 1 for stochastic rates. ( P_nm = " << prob
        << " )\nConsider using smaller timesteps.";
    if (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 std::move(act);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

string_process_interface_

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

private

Class that deals with strings, interfacing Pythia.

Definition at line 317 of file scatteractionsfinder.h.

coll_crit_

const CollisionCriterion smash::ScatterActionsFinder::coll_crit_

private

Specifies which collision criterion is used.

Definition at line 319 of file scatteractionsfinder.h.

elastic_parameter_

const double smash::ScatterActionsFinder::elastic_parameter_

private

Elastic cross section parameter (in mb).

Definition at line 321 of file scatteractionsfinder.h.

testparticles_

const int smash::ScatterActionsFinder::testparticles_

private

Number of test particles.

Definition at line 323 of file scatteractionsfinder.h.

isotropic_

const bool smash::ScatterActionsFinder::isotropic_

private

Do all collisions isotropically.

Definition at line 325 of file scatteractionsfinder.h.

two_to_one_

const bool smash::ScatterActionsFinder::two_to_one_

private

Enable 2->1 processes.

Definition at line 327 of file scatteractionsfinder.h.

incl_set_

const ReactionsBitSet smash::ScatterActionsFinder::incl_set_

private

List of included 2<->2 reactions.

Definition at line 329 of file scatteractionsfinder.h.

incl_multi_set_

const MultiParticleReactionsBitSet smash::ScatterActionsFinder::incl_multi_set_

private

List of included multi-particle reactions.

Definition at line 331 of file scatteractionsfinder.h.

scale_xs_

const double smash::ScatterActionsFinder::scale_xs_

private

Factor by which all (partial) cross sections are scaled.

Definition at line 333 of file scatteractionsfinder.h.

additional_el_xs_

const double smash::ScatterActionsFinder::additional_el_xs_

private

Additional constant elastic cross section.

Definition at line 335 of file scatteractionsfinder.h.

low_snn_cut_

const double smash::ScatterActionsFinder::low_snn_cut_

private

Elastic collsions between two nucleons with sqrt_s below low_snn_cut_ are excluded.

Definition at line 340 of file scatteractionsfinder.h.

strings_switch_

const bool smash::ScatterActionsFinder::strings_switch_

private

Switch to turn off string excitation.

Definition at line 342 of file scatteractionsfinder.h.

use_AQM_

const bool smash::ScatterActionsFinder::use_AQM_

private

Switch to control whether to use AQM or not.

Definition at line 344 of file scatteractionsfinder.h.

strings_with_probability_

const bool smash::ScatterActionsFinder::strings_with_probability_

private

Decide whether to implement string fragmentation based on a probability.

Definition at line 346 of file scatteractionsfinder.h.

nnbar_treatment_

const NNbarTreatment smash::ScatterActionsFinder::nnbar_treatment_

private

Switch for NNbar reactions.

Definition at line 348 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 354 of file scatteractionsfinder.h.

string_formation_time_

const double smash::ScatterActionsFinder::string_formation_time_

private

Parameter for formation time.

Definition at line 356 of file scatteractionsfinder.h.

maximum_cross_section_

const double smash::ScatterActionsFinder::maximum_cross_section_

private

See also

input_collision_term_

Definition at line 358 of file scatteractionsfinder.h.

allow_first_collisions_within_nucleus_

const bool smash::ScatterActionsFinder::allow_first_collisions_within_nucleus_

private

If particles within nucleus are allowed to collide for their first time.

Definition at line 360 of file scatteractionsfinder.h.

only_warn_for_high_prob_

const bool smash::ScatterActionsFinder::only_warn_for_high_prob_

private

Switch to turn off throwing an exception for collision probabilities larger than 1.

In larger production runs it is ok, if the probability rarely slips over 1.

Definition at line 366 of file scatteractionsfinder.h.

---

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

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