#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:

Collaboration diagram for smash::ScatterActionsFinder:

Public Member Functions
	ScatterActionsFinder (Configuration config, const ExperimentParameters &parameters, const std::vector< bool > &nucleon_has_interacted, int N_tot, int N_proj)
	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 cell_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, 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 (const ParticleData &data_a, const ParticleData &data_b, double dt, const std::vector< FourVector > &beam_momentum={}, const double cell_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...

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 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 std::vector< bool > &	nucleon_has_interacted_
	Parameter to record whether the nucleon has experienced a collision or not. More...

const int	N_tot_
	Record the total number of the nucleons in the two colliding nuclei. More...

const int	N_proj_
	Record the number of the nucleons in the projectile. More...

const double	string_formation_time_
	Parameter for formation time. More...

Constructor & Destructor Documentation

smash::ScatterActionsFinder::ScatterActionsFinder	(	Configuration	config,
		const ExperimentParameters &	parameters,
		const std::vector< bool > &	nucleon_has_interacted,
		int	N_tot,
		int	N_proj
	)

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.
[in]	nucleon_has_interacted	Flags to record whether an initial nucleon has interacted with another particle not from the same nucleus. The flags are used if we want to exclude the first collisions among the nucleons within the same nucleus.
[in]	N_tot	Total number of the initial nucleons. This number, as well as the next parameter, will be used to determine whether two intial nucleons are within the same nucleus if we'd like to exclude the first collisions among them.
[in]	N_proj	Total projectile number

Definition at line 218 of file scatteractionsfinder.cc.

     : coll_crit_(config.take({"Collision_Term", "Collision_Criterion"},
                              CollisionCriterion::Geometric)),
       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),
       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),
       nucleon_has_interacted_(nucleon_has_interacted),
       N_tot_(N_tot),
       N_proj_(N_proj),
       string_formation_time_(config.take(
           {"Collision_Term", "String_Parameters", "Formation_Time"}, 1.)) {
   if (coll_crit_ == CollisionCriterion::Stochastic &&
       !(is_constant_elastic_isotropic())) {
     throw std::invalid_argument(
         "The stochastic collision criterion is only supported for elastic (and "
         "isotropic)\n2-to-2 reactions of one particle species. Change you "
         "config accordingly.");
   }
   if (is_constant_elastic_isotropic()) {
     const auto& log = logger<LogArea::FindScatter>();
     log.info("Constant elastic isotropic cross-section mode:", " using ",
              elastic_parameter_, " mb as maximal cross-section.");
   }
   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));
   }
 }

Here is the call graph for this function:

Member Function Documentation

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 as in Xu:2004mz.

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.

UrQMD collision time in computational frame, see Bass:1998ca (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 78 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();
 
       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:

ActionList smash::ScatterActionsFinder::find_actions_in_cell	(	const ParticleList &	search_list,
		double	dt,
		const double	cell_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]	cell_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 399 of file scatteractionsfinder.cc.

                                                       {
   std::vector<ActionPtr> actions;
   for (const ParticleData& p1 : search_list) {
     for (const ParticleData& p2 : search_list) {
       if (p1.id() < p2.id()) {
         // Check if a collision is possible.
         ActionPtr act = check_collision(p1, p2, dt, beam_momentum, cell_vol);
         if (act) {
           actions.push_back(std::move(act));
         }
       }
     }
   }
   return actions;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

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 417 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(p1, p2, dt, beam_momentum);
       if (act) {
         actions.push_back(std::move(act));
       }
     }
   }
   return actions;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

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 438 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(p1, p2, dt, beam_momentum);
       if (act) {
         actions.push_back(std::move(act));
       }
     }
   }
   return actions;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

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

                                                                {
     return ActionList();
   }

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 201 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:

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

inline

The maximal distance over which particles can interact, 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.

Definition at line 216 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:

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 466 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_,
                                      low_snn_cut_, strings_switch_, use_AQM_,
                                      strings_with_probability_,
                                      nnbar_treatment_);
             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:

Here is the caller graph for this function:

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 852 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_,
                              low_snn_cut_, strings_switch_, use_AQM_,
                              strings_with_probability_, nnbar_treatment_);
     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:

Here is the caller graph for this function:

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

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

Here is the call graph for this function:

ActionPtr smash::ScatterActionsFinder::check_collision	(	const ParticleData &	data_a,
		const ParticleData &	data_b,
		double	dt,
		const std::vector< FourVector > &	beam_momentum = `{}`,
		const double	cell_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 (3.27). 2. A stochastic collision criterion e.g. employed by BAMPS Xu:2004mz (Sec.IIB). Note, the latter is currently only tested for a box with a fixed elastic cross section.

More details on the stochastic collision criterion can be found here:

P. Danielewicz and G. F. Bertsch, Nucl. Phys. A533, 712 (1991).
A. Lang, H. Babovsky, W. Cassing, U. Mosel, H. G. Reusch, and K. Weber, J. Comp. Phys. 106, 391 (1993).
W. Cassing, Nucl. Phys. A700, 618 (2002).

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]	cell_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: cell_vol is optional, since only find_actions_in_cell has (and needs) this information for the stochastic collision criterion.

Definition at line 274 of file scatteractionsfinder.cc.

                                                                              {
   const auto& log = logger<LogArea::FindScatter>();
 
   /* If the two particles
    * 1) belong to the two colliding nuclei
    * 2) are within the same nucleus
    * 3) both of them have never experienced any collisons,
    * then the collision between them are banned. */
   assert(data_a.id() >= 0);
   assert(data_b.id() >= 0);
   if (data_a.id() < N_tot_ && data_b.id() < N_tot_ &&
       ((data_a.id() < N_proj_ && data_b.id() < N_proj_) ||
        (data_a.id() >= N_proj_ && data_b.id() >= N_proj_)) &&
       !(nucleon_has_interacted_[data_a.id()] ||
         nucleon_has_interacted_[data_b.id()])) {
     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_);
 
   if (strings_switch_) {
     act->set_string_interface(string_process_interface_.get());
   }
 
   if (coll_crit_ == CollisionCriterion::Stochastic) {
     // No grid or search in cell
     if (cell_vol < really_small) {
       return nullptr;
     }
 
     // Add various subprocesses.
     act->add_all_scatterings(elastic_parameter_, two_to_one_, incl_set_,
                              low_snn_cut_, strings_switch_, use_AQM_,
                              strings_with_probability_, nnbar_treatment_);
 
     const double xs = act->cross_section() * fm2_mb;
 
     // Relative velocity calculation, see e.g. \iref{Seifert:2017oyb}, eq. (5)
     const double m1 = act->incoming_particles()[0].effective_mass();
     const double m1_sqr = m1 * m1;
     const double m2 = act->incoming_particles()[1].effective_mass();
     const double m2_sqr = m2 * m2;
     const double e1 = act->incoming_particles()[0].momentum().x0();
     const double e2 = act->incoming_particles()[1].momentum().x0();
     const double m_s = act->mandelstam_s();
     const double lambda = (m_s - m1_sqr - m2_sqr) * (m_s - m1_sqr - m2_sqr) -
                           4. * m1_sqr * m2_sqr;
     const double v_rel = std::sqrt(lambda) / (2. * e1 * e2);
 
     // Collision probability, see e.g. \iref{Xu:2004mz}, eq. (11)
     const double p_22 = xs * v_rel * dt / (cell_vol * testparticles_);
 
     log.debug("Stochastic collison criterion parameters:\np_22 = ", p_22,
               ", xs = ", xs, ", v_rel = ", v_rel, ", dt = ", dt,
               ", cell_vol = ", cell_vol, ", testparticles = ", testparticles_);
 
     if (p_22 > 1.) {
       std::stringstream err;
       err << "Probability larger than 1 for stochastic rates. ( P = " << p_22
           << " )\nUse smaller timesteps.";
       throw std::runtime_error(err.str());
     }
 
     // probability criterion
     double random_no = random::uniform(0., 1.);
     if (random_no > p_22) {
       return nullptr;
     }
 
   } else if (coll_crit_ == CollisionCriterion::Geometric) {
     // just collided with this particle
     if (data_a.id_process() > 0 && data_a.id_process() == data_b.id_process()) {
       log.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;
     }
 
     const double distance_squared = act->transverse_distance_sqr();
 
     // Don't calculate cross section if the particles are very far apart.
     if (distance_squared >= max_transverse_distance_sqr(testparticles_)) {
       return nullptr;
     }
 
     // Add various subprocesses.
     act->add_all_scatterings(elastic_parameter_, two_to_one_, incl_set_,
                              low_snn_cut_, strings_switch_, use_AQM_,
                              strings_with_probability_, nnbar_treatment_);
 
     // Cross section for collision criterion
     double cross_section_criterion = act->cross_section() * fm2_mb * M_1_PI /
                                      static_cast<double>(testparticles_);
 
     // Take cross section scaling factors into account
     cross_section_criterion *= data_a.xsec_scaling_factor(time_until_collision);
     cross_section_criterion *= data_b.xsec_scaling_factor(time_until_collision);
 
     // distance criterion according to cross_section
     if (distance_squared >= cross_section_criterion) {
       return nullptr;
     }
 
     log.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:

Member Data Documentation

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

private

Class that deals with strings, interfacing Pythia.

Definition at line 294 of file scatteractionsfinder.h.

const CollisionCriterion smash::ScatterActionsFinder::coll_crit_

private

Specifies which collision criterion is used.

Definition at line 296 of file scatteractionsfinder.h.

const double smash::ScatterActionsFinder::elastic_parameter_

private

Elastic cross section parameter (in mb).

Definition at line 298 of file scatteractionsfinder.h.

const int smash::ScatterActionsFinder::testparticles_

private

Number of test particles.

Definition at line 300 of file scatteractionsfinder.h.

const bool smash::ScatterActionsFinder::isotropic_

private

Do all collisions isotropically.

Definition at line 302 of file scatteractionsfinder.h.

const bool smash::ScatterActionsFinder::two_to_one_

private

Enable 2->1 processes.

Definition at line 304 of file scatteractionsfinder.h.

const ReactionsBitSet smash::ScatterActionsFinder::incl_set_

private

List of included 2<->2 reactions.

Definition at line 306 of file scatteractionsfinder.h.

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

const bool smash::ScatterActionsFinder::strings_switch_

private

Switch to turn off string excitation.

Definition at line 313 of file scatteractionsfinder.h.

const bool smash::ScatterActionsFinder::use_AQM_

private

Switch to control whether to use AQM or not.

Definition at line 315 of file scatteractionsfinder.h.

const bool smash::ScatterActionsFinder::strings_with_probability_

private

Decide whether to implement string fragmentation based on a probability.

Definition at line 317 of file scatteractionsfinder.h.

const NNbarTreatment smash::ScatterActionsFinder::nnbar_treatment_

private

Switch for NNbar reactions.

Definition at line 319 of file scatteractionsfinder.h.

const std::vector<bool>& smash::ScatterActionsFinder::nucleon_has_interacted_

private

Parameter to record whether the nucleon has experienced a collision or not.

Definition at line 323 of file scatteractionsfinder.h.

const int smash::ScatterActionsFinder::N_tot_

private

Record the total number of the nucleons in the two colliding nuclei.

Definition at line 325 of file scatteractionsfinder.h.

const int smash::ScatterActionsFinder::N_proj_

private

Record the number of the nucleons in the projectile.

Definition at line 327 of file scatteractionsfinder.h.

const double smash::ScatterActionsFinder::string_formation_time_

private

Parameter for formation time.

Definition at line 329 of file scatteractionsfinder.h.

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

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

Public Member Functions

Private Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation