scatteractionsfinder.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2020
 *      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"
 
namespace smash {
 
class ScatterActionsFinder : public ActionFinderInterface {
 public:
  ScatterActionsFinder(Configuration config,
                       const ExperimentParameters &parameters,
                       const std::vector<bool> &nucleon_has_interacted,
                       int N_tot, int N_proj);
 
  inline double collision_time(
      const ParticleData &p1, const ParticleData &p2, double dt,
      const std::vector<FourVector> &beam_momentum) const {
    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;
 
        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 && !two_to_one_ && isotropic_ &&
           elastic_parameter_ > 0.;
  }
 
  double max_transverse_distance_sqr(int testparticles) const {
    return (is_constant_elastic_isotropic() ? elastic_parameter_
                                            : 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 (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;
 
  std::unique_ptr<StringProcess> string_process_interface_;
  const CollisionCriterion coll_crit_;
  const double elastic_parameter_;
  const int testparticles_;
  const bool isotropic_;
  const bool two_to_one_;
  const ReactionsBitSet incl_set_;
  const MultiParticleReactionsBitSet incl_multi_set_;
  const double scale_xs_;
  const double additional_el_xs_;
  const double low_snn_cut_;
  const bool strings_switch_;
  const bool use_AQM_;
  const bool strings_with_probability_;
  const NNbarTreatment nnbar_treatment_;
  const double box_length_;
  const std::vector<bool> &nucleon_has_interacted_;
  const int N_tot_;
  const int N_proj_;
  const double string_formation_time_;
  const double maximum_cross_section_;
  const bool only_warn_for_high_prob_;
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_SCATTERACTIONSFINDER_H_