scatteractionsfinder.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2022
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_SCATTERACTIONSFINDER_H_
#define SRC_INCLUDE_SMASH_SCATTERACTIONSFINDER_H_
 
#include <memory>
#include <set>
#include <vector>
 
#include "action.h"
#include "actionfinderfactory.h"
#include "configuration.h"
#include "scatteraction.h"
#include "scatteractionsfinderparameters.h"
 
namespace smash {
 
class ScatterActionsFinder : public ActionFinderInterface {
 public:
  ScatterActionsFinder(Configuration &config,
                       const ExperimentParameters &parameters);
 
  inline double collision_time(
      const ParticleData &p1, const ParticleData &p2, double dt,
      const std::vector<FourVector> &beam_momentum) const {
    if (finder_parameters_.coll_crit == CollisionCriterion::Stochastic) {
      return dt * random::uniform(0., 1.);
    } else {
      /*
       * For frozen Fermi motion:
       * If particles have not yet interacted and are the initial nucleons,
       * perform action finding with beam momentum instead of Fermi motion
       * corrected momentum. That is because the particles are propagated with
       * the beam momentum until they interact.
       */
      if (p1.id() < 0 || p2.id() < 0) {
        throw std::runtime_error("Invalid particle ID for Fermi motion");
      }
      const bool p1_has_no_prior_interactions =
          (static_cast<uint64_t>(p1.id()) <                 // particle from
           static_cast<uint64_t>(beam_momentum.size())) &&  // initial nucleus
          (p1.get_history().collisions_per_particle == 0);
 
      const bool p2_has_no_prior_interactions =
          (static_cast<uint64_t>(p2.id()) <                 // particle from
           static_cast<uint64_t>(beam_momentum.size())) &&  // initial nucleus
          (p2.get_history().collisions_per_particle == 0);
 
      const FourVector p1_mom = (p1_has_no_prior_interactions)
                                    ? beam_momentum[p1.id()]
                                    : p1.momentum();
      const FourVector p2_mom = (p2_has_no_prior_interactions)
                                    ? beam_momentum[p2.id()]
                                    : p2.momentum();
      if (finder_parameters_.coll_crit == CollisionCriterion::Covariant) {
        const FourVector delta_x = p1.position() - p2.position();
        const double p1_sqr = p1_mom.sqr();
        const double p2_sqr = p2_mom.sqr();
        const double p1_dot_x = p1_mom.Dot(delta_x);
        const double p2_dot_x = p2_mom.Dot(delta_x);
        const double p1_dot_p2 = p1_mom.Dot(p2_mom);
        const double denominator = std::pow(p1_dot_p2, 2) - p1_sqr * p2_sqr;
        if (unlikely(std::abs(denominator) < really_small * really_small)) {
          return -1.0;
        }
 
        const double time_1 = (p2_sqr * p1_dot_x - p1_dot_p2 * p2_dot_x) *
                              p1_mom.x0() / denominator;
        const double time_2 = -(p1_sqr * p2_dot_x - p1_dot_p2 * p1_dot_x) *
                              p2_mom.x0() / denominator;
        return (time_1 + time_2) / 2;
      } else {
        const ThreeVector dv_times_e1e2 =
            p1_mom.threevec() * p2_mom.x0() - p2_mom.threevec() * p1_mom.x0();
        const double dv_times_e1e2_sqr = dv_times_e1e2.sqr();
        if (dv_times_e1e2_sqr < really_small) {
          return -1.0;
        }
        const ThreeVector dr =
            p1.position().threevec() - p2.position().threevec();
        return -(dr * dv_times_e1e2) *
               (p1_mom.x0() * p2_mom.x0() / dv_times_e1e2_sqr);
      }
    }
  }
 
  ActionList find_actions_in_cell(
      const ParticleList &search_list, double dt, const double gcell_vol,
      const std::vector<FourVector> &beam_momentum) const override;
 
  ActionList find_actions_with_neighbors(
      const ParticleList &search_list, const ParticleList &neighbors_list,
      double dt, const std::vector<FourVector> &beam_momentum) const override;
 
  ActionList find_actions_with_surrounding_particles(
      const ParticleList &search_list, const Particles &surrounding_list,
      double dt, const std::vector<FourVector> &beam_momentum) const override;
 
  ActionList find_final_actions(const Particles & /*search_list*/,
                                bool /*only_res*/ = false) const override {
    return ActionList();
  }
 
  inline bool is_constant_elastic_isotropic() const {
    return ParticleType::list_all().size() == 1 &&
           !finder_parameters_.two_to_one && isotropic_ &&
           finder_parameters_.elastic_parameter > 0.;
  }
 
  double max_transverse_distance_sqr(int testparticles) const {
    return (is_constant_elastic_isotropic()
                ? finder_parameters_.elastic_parameter
                : finder_parameters_.maximum_cross_section) /
           testparticles * fm2_mb * M_1_PI;
  }
 
  void dump_reactions() const;
 
  void dump_cross_sections(const ParticleType &a, const ParticleType &b,
                           double m_a, double m_b, bool final_state,
                           std::vector<double> &plab) const;
 
  StringProcess *get_process_string_ptr() {
    if (finder_parameters_.strings_switch) {
      return string_process_interface_.get();
    } else {
      return NULL;
    }
  }
 
 private:
  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;
 
  ActionPtr check_collision_multi_part(const ParticleList &plist, double dt,
                                       const double gcell_vol) const;
 
  ScatterActionsFinderParameters finder_parameters_;
  std::unique_ptr<StringProcess> string_process_interface_;
  const bool isotropic_;
  const double box_length_;
  const double string_formation_time_;
};
 
ScatterActionsFinderParameters create_finder_parameters(
    Configuration &config, const ExperimentParameters &parameters);
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_SCATTERACTIONSFINDER_H_