194 const std::vector<bool>& nucleon_has_interacted,
int N_tot,
int N_proj)
195 : elastic_parameter_(
196 config.take({
"Collision_Term",
"Elastic_Cross_Section"}, -1.)),
197 testparticles_(parameters.testparticles),
198 isotropic_(config.take({
"Collision_Term",
"Isotropic"},
false)),
199 two_to_one_(parameters.two_to_one),
200 incl_set_(parameters.included_2to2),
201 low_snn_cut_(parameters.low_snn_cut),
202 strings_switch_(parameters.strings_switch),
203 use_AQM_(parameters.use_AQM),
204 strings_with_probability_(parameters.strings_with_probability),
205 nnbar_treatment_(parameters.nnbar_treatment),
206 nucleon_has_interacted_(nucleon_has_interacted),
209 string_formation_time_(config.take(
210 {
"Collision_Term",
"String_Parameters",
"Formation_Time"}, 1.)) {
211 if (is_constant_elastic_isotropic()) {
212 const auto& log = logger<LogArea::FindScatter>();
213 log.info(
"Constant elastic isotropic cross-section mode:",
" using ",
214 elastic_parameter_,
" mb as maximal cross-section.");
216 if (strings_switch_) {
217 auto subconfig = config[
"Collision_Term"][
"String_Parameters"];
218 string_process_interface_ = make_unique<StringProcess>(
219 subconfig.take({
"String_Tension"}, 1.0), string_formation_time_,
220 subconfig.take({
"Gluon_Beta"}, 0.5),
221 subconfig.take({
"Gluon_Pmin"}, 0.001),
222 subconfig.take({
"Quark_Alpha"}, 2.0),
223 subconfig.take({
"Quark_Beta"}, 5.0),
224 subconfig.take({
"Strange_Supp"}, 0.12),
225 subconfig.take({
"Diquark_Supp"}, 0.03),
226 subconfig.take({
"Sigma_Perp"}, 0.42),
227 subconfig.take({
"Leading_Frag_Mean"}, 1.0),
228 subconfig.take({
"Leading_Frag_Width"}, 0.6),
229 subconfig.take({
"StringZ_A"}, 0.68), subconfig.take({
"StringZ_B"}, 0.3),
230 subconfig.take({
"String_Sigma_T"}, 0.5),
231 subconfig.take({
"Form_Time_Factor"}, 1.0),
232 subconfig.take({
"Use_Yoyo_Model"},
true),
233 subconfig.take({
"Prob_proton_to_d_uu"}, 1. / 3.));
241 const auto& log = logger<LogArea::FindScatter>();
247 log.debug(
"Skipping collided particles at time ", data_a.
position().
x0(),
248 " due to process ", data_a.
id_process(),
"\n ", data_a,
258 assert(data_a.
id() >= 0);
259 assert(data_b.
id() >= 0);
269 const double time_until_collision =
collision_time(data_a, data_b);
272 if (time_until_collision < 0. || time_until_collision >= dt) {
277 ScatterActionPtr act = make_unique<ScatterAction>(
283 const double distance_squared = act->transverse_distance_sqr();
296 double cross_section_criterion = act->cross_section() *
fm2_mb * M_1_PI /
299 cross_section_criterion *= data_a.xsec_scaling_factor(time_until_collision);
300 cross_section_criterion *= data_b.xsec_scaling_factor(time_until_collision);
303 if (distance_squared >= cross_section_criterion) {
308 log.debug(
"particle distance squared: ", distance_squared,
"\n ", data_a,
312 return std::move(act);
316 const ParticleList& search_list,
double dt)
const {
317 std::vector<ActionPtr> actions;
320 if (p1.id() < p2.id()) {
324 actions.push_back(std::move(act));
333 const ParticleList& search_list,
const ParticleList& neighbors_list,
335 std::vector<ActionPtr> actions;
338 assert(p1.id() != p2.id());
342 actions.push_back(std::move(act));
350 const ParticleList& search_list,
const Particles& surrounding_list,
352 std::vector<ActionPtr> actions;
356 auto result = std::find_if(
357 search_list.begin(), search_list.end(),
359 if (result != search_list.end()) {
366 actions.push_back(std::move(act));
374 constexpr
double time = 0.0;
377 const size_t N_pairs = N_isotypes * (N_isotypes - 1) / 2;
379 std::cout << N_isotypes <<
" iso-particle types." << std::endl;
380 std::cout <<
"They can make " << N_pairs <<
" pairs." << std::endl;
381 std::vector<double> momentum_scan_list = {0.1, 0.3, 0.5, 1.0,
382 2.0, 3.0, 5.0, 10.0};
385 if (&A_isotype > &B_isotype) {
388 bool any_nonzero_cs =
false;
389 std::vector<std::string> r_list;
392 if (A_type > B_type) {
396 for (
auto mom : momentum_scan_list) {
397 A.set_4momentum(A.pole_mass(), mom, 0.0, 0.0);
399 ScatterActionPtr act = make_unique<ScatterAction>(
408 const double total_cs = act->cross_section();
409 if (total_cs <= 0.0) {
412 any_nonzero_cs =
true;
413 for (
const auto& channel : act->collision_channels()) {
414 const auto type = channel->get_type();
418 r = A_type->name() + B_type->name() + std::string(
" → strings");
422 ? std::string(
" (el)")
424 ? std::string(
" (inel)")
425 : std::string(
" (?)");
426 r = A_type->name() + B_type->name() + std::string(
" → ") +
427 channel->particle_types()[0]->name() +
428 channel->particle_types()[1]->name() + r_type;
436 std::sort(r_list.begin(), r_list.end());
437 r_list.erase(std::unique(r_list.begin(), r_list.end()), r_list.end());
438 if (any_nonzero_cs) {
439 for (
auto r : r_list) {
441 if (r_list.back() != r) {
445 std::cout << std::endl;
475 namespace decaytree {
523 Node(
const std::string& name,
double weight,
524 ParticleTypePtrList&& initial_particles,
525 ParticleTypePtrList&& final_particles, ParticleTypePtrList&& state,
526 std::vector<Node>&& children)
546 ParticleTypePtrList&& initial_particles,
547 ParticleTypePtrList&& final_particles) {
549 ParticleTypePtrList state(
state_);
550 for (
const auto&
p : initial_particles) {
551 state.erase(std::find(state.begin(), state.end(),
p));
553 for (
const auto&
p : final_particles) {
557 std::sort(state.begin(), state.end(),
559 return a->
name() < b->name();
562 Node new_node(name, weight, std::move(initial_particles),
563 std::move(final_particles), std::move(state), {});
564 children_.emplace_back(std::move(new_node));
575 std::vector<FinalStateCrossSection> result;
588 for (uint64_t i = 0; i < depth; i++) {
593 child.print_helper(depth + 1);
609 uint64_t depth, std::vector<FinalStateCrossSection>& result,
610 const std::string& name,
double weight,
611 bool show_intermediate_states =
false)
const {
619 std::string new_name;
622 if (show_intermediate_states) {
624 if (!new_name.empty()) {
632 for (
const auto& s :
state_) {
633 new_name += s->name();
636 if (show_intermediate_states) {
645 child.final_state_cross_sections_helper(depth + 1, result, new_name,
646 weight, show_intermediate_states);
660 const DecayBranchPtr& decay,
661 ParticleTypePtrList& final_state) {
662 std::stringstream name;
663 name <<
"[" << res_name <<
"->";
664 for (
const auto&
p : decay->particle_types()) {
666 final_state.push_back(
p);
689 uint32_t n_unstable = 0;
696 n_unstable != 0 ? 1. /
static_cast<double>(n_unstable) : 1.;
701 ParticleTypePtrList parts;
703 auto& new_node = node.
add_action(name, norm * decay->weight(), {ptype},
718 static void deduplicate(std::vector<FinalStateCrossSection>& final_state_xs) {
719 std::sort(final_state_xs.begin(), final_state_xs.end(),
722 auto current = final_state_xs.begin();
723 while (current != final_state_xs.end()) {
724 auto adjacent = std::adjacent_find(
725 current, final_state_xs.end(),
727 return a.
name_ == b.name_;
730 if (adjacent != final_state_xs.end()) {
731 adjacent->cross_section_ += (adjacent + 1)->cross_section_;
732 final_state_xs.erase(adjacent + 1);
739 double m_a,
double m_b,
740 bool final_state)
const {
741 typedef std::vector<std::pair<double, double>> xs_saver;
742 std::map<std::string, xs_saver> xs_dump;
743 std::map<std::string, double> outgoing_total_mass;
746 constexpr
int n_momentum_points = 200;
747 constexpr
double momentum_step = 0.02;
748 for (
int i = 1; i < n_momentum_points + 1; i++) {
749 const double momentum = momentum_step * i;
750 a_data.set_4momentum(m_a, momentum, 0.0, 0.0);
752 const double sqrts = (a_data.momentum() + b_data.
momentum()).abs();
753 const ParticleList incoming = {a_data, b_data};
754 ScatterActionPtr act = make_unique<ScatterAction>(
763 {&a, &b}, {&a, &b}, {});
764 const CollisionBranchList& processes = act->collision_channels();
765 for (
const auto& process : processes) {
766 const double xs = process->weight();
771 std::stringstream process_description_stream;
772 process_description_stream << *process;
773 const std::string& description = process_description_stream.str();
775 for (
const auto& ptype : process->particle_types()) {
776 m_tot += ptype->mass();
778 outgoing_total_mass[description] = m_tot;
779 if (!xs_dump[description].empty() &&
780 std::abs(xs_dump[description].back().first - sqrts) <
782 xs_dump[description].back().second += xs;
784 xs_dump[description].push_back(std::make_pair(sqrts, xs));
787 std::stringstream process_description_stream;
788 process_description_stream << *process;
789 const std::string& description = process_description_stream.str();
790 ParticleTypePtrList initial_particles = {&a, &b};
791 ParticleTypePtrList final_particles = process->particle_types();
793 tree.add_action(description, xs, std::move(initial_particles),
794 std::move(final_particles));
798 xs_dump[
"total"].push_back(std::make_pair(sqrts, act->cross_section()));
800 outgoing_total_mass[
"total"] = -1.0;
803 auto final_state_xs = tree.final_state_cross_sections();
805 for (
const auto&
p : final_state_xs) {
806 outgoing_total_mass[
p.name_] =
p.mass_;
807 xs_dump[
p.name_].push_back(std::make_pair(sqrts,
p.cross_section_));
813 std::vector<std::string> all_channels;
814 for (
const auto channel : xs_dump) {
815 all_channels.push_back(channel.first);
817 std::sort(all_channels.begin(), all_channels.end(),
818 [&](
const std::string& str_a,
const std::string& str_b) {
819 return outgoing_total_mass[str_a] < outgoing_total_mass[str_b];
823 std::cout <<
"# Dumping partial " << a.
name() << b.
name()
824 <<
" cross-sections in mb, energies in GeV" << std::endl;
825 std::cout <<
" sqrt_s";
829 for (
const auto channel : all_channels) {
832 std::cout << std::endl;
835 for (
int i = 1; i < n_momentum_points; i++) {
836 const double momentum = momentum_step * i;
837 a_data.set_4momentum(m_a, momentum, 0.0, 0.0);
839 const double sqrts = (a_data.momentum() + b_data.
momentum()).abs();
840 printf(
"%9.5f", sqrts);
841 for (
const auto channel : all_channels) {
842 const xs_saver energy_and_xs = xs_dump[channel];
844 for (; j < energy_and_xs.size() && energy_and_xs[j].first < sqrts; j++) {
847 if (j < energy_and_xs.size() &&
848 std::abs(energy_and_xs[j].first - sqrts) <
really_small) {
849 xs = energy_and_xs[j].second;
851 printf(
"%16.6f", xs);
void print_helper(uint64_t depth) const
Internal helper function for print, to be called recursively to print all nodes.
Node(const Node &)=delete
Cannot be copied.
std::string name_
Name of the final state.
static double collision_time(const ParticleData &p1, const ParticleData &p2)
Determine the collision time of the two particles.
double mass_
Total mass of final state particles.
static void deduplicate(std::vector< FinalStateCrossSection > &final_state_xs)
Deduplicate the final-state cross sections by summing.
const int testparticles_
Number of test particles.
constexpr double really_small
Numerical error tolerance.
const ReactionsBitSet incl_set_
List of included 2<->2 reactions.
bool is_string_soft_process(ProcessType p)
Check if a given process type is a soft string excitation.
const bool two_to_one_
Enable 2->1 processes.
static std::string make_decay_name(const std::string &res_name, const DecayBranchPtr &decay, ParticleTypePtrList &final_state)
Generate name for decay and update final state.
const double low_snn_cut_
Elastic collsions between two nucleons with sqrt_s below low_snn_cut_ are excluded.
const double string_formation_time_
Parameter for formation time.
void isoclean(std::string &s)
Remove ⁺, ⁻, ⁰ from string.
Collection of useful constants that are known at compile time.
double pole_mass() const
Get the particle's pole mass ("on-shell").
const bool strings_with_probability_
Decide whether to implement string fragmentation based on a probability.
double max_transverse_distance_sqr(int testparticles) const
The maximal distance over which particles can interact, related to the number of test particles and t...
ActionPtr check_collision(const ParticleData &data_a, const ParticleData &data_b, double dt) const
Check for a single pair of particles (id_a, id_b) if a collision will happen in the next timestep and...
constexpr double fm2_mb
mb <-> fm^2 conversion factor.
2->2 inelastic scattering
std::vector< Node > children_
Possible actions after this action.
const FourVector & momentum() const
Get the particle's 4-momentum.
const std::vector< bool > & nucleon_has_interacted_
Parameter to record whether the nucleon has experienced a collision or not.
Represent a final-state cross section.
Interface to the SMASH configuration files.
std::string fill_left(const std::string &s, size_t width, char fill=' ')
Fill string with characters to the left until the given width is reached.
std::unique_ptr< StringProcess > string_process_interface_
Class that deals with strings, interfacing Pythia.
const bool strings_switch_
Switch to turn off string excitation.
ParticleTypePtrList initial_particles_
Initial-state particle types in this action.
ActionList find_actions_in_cell(const ParticleList &search_list, double dt) const override
Search for all the possible collisions within one cell.
std::string name_
Name for printing.
elastic scattering: particles remain the same, only momenta change
double weight_
Weight (cross section or branching ratio).
void final_state_cross_sections_helper(uint64_t depth, std::vector< FinalStateCrossSection > &result, const std::string &name, double weight, bool show_intermediate_states=false) const
Internal helper function for final_state_cross_sections, to be called recursively to calculate all fi...
const FourVector & position() const
Get the particle's position in Minkowski space.
Particle type contains the static properties of a particle species.
ScatterActionsFinder(Configuration config, const ExperimentParameters ¶meters, const std::vector< bool > &nucleon_has_interacted, int N_tot, int N_proj)
Constructor of the finder with the given parameters.
void dump_cross_sections(const ParticleType &a, const ParticleType &b, double m_a, double m_b, bool final_state) const
Print out partial cross-sections of all processes that can occur in the collision of a(mass = m_a) an...
void set_4momentum(const FourVector &momentum_vector)
Set the particle's 4-momentum directly.
IsoParticleType is a class to represent isospin multiplets.
const NNbarTreatment nnbar_treatment_
Switch for NNbar reactions.
uint32_t id_process() const
Get the id of the last action.
hard string process involving 2->2 QCD process by PYTHIA.
FinalStateCrossSection(const std::string &name, double cross_section, double mass)
Construct a final-state cross section.
ParticleTypePtrList final_particles_
Final-state particle types in this action.
void dump_reactions() const
Prints out all the 2-> n (n > 1) reactions with non-zero cross-sections between all possible pairs of...
std::vector< FinalStateCrossSection > final_state_cross_sections() const
const DecayModes & decay_modes() const
const bool use_AQM_
Switch to control whether to use AQM or not.
const int N_proj_
Record the number of the nucleons in the projectile.
const bool isotropic_
Do all collisions isotropically.
const int N_tot_
Record the total number of the nucleons in the two colliding nuclei.
void print() const
Print the decay tree starting with this node.
A pointer-like interface to global references to ParticleType objects.
static void add_decays(Node &node)
Add nodes for all decays possible from the given node and all of its children.
Node of a decay tree, representing a possible action (2-to-2 or 1-to-2).
ParticleTypePtrList state_
Particle types corresponding to the global state after this action.
The Particles class abstracts the storage and manipulation of particles.
int id() const
Get the id of the particle.
Node(const std::string &name, double weight, ParticleTypePtrList &&initial_particles, ParticleTypePtrList &&final_particles, ParticleTypePtrList &&state, std::vector< Node > &&children)
double cross_section_
Corresponding cross section in mb.
Node & add_action(const std::string &name, double weight, ParticleTypePtrList &&initial_particles, ParticleTypePtrList &&final_particles)
Add an action to the children of this node.
ActionList find_actions_with_surrounding_particles(const ParticleList &search_list, const Particles &surrounding_list, double dt) const override
Search for all the possible secondary collisions between the outgoing particles and the rest...
Helper structure for Experiment.
ParticleData contains the dynamic information of a certain particle.
const double elastic_parameter_
Elastic cross section parameter (in mb).
static const IsoParticleTypeList & list_all()
Returns a list of all IsoParticleTypes.
ActionList find_actions_with_neighbors(const ParticleList &search_list, const ParticleList &neighbors_list, double dt) const override
Search for all the possible collisions among the neighboring cells.
const DecayBranchList & decay_mode_list() const
const std::string & name() const