scatteractionsfinder.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */

#ifndef SRC_INCLUDE_SCATTERACTIONSFINDER_H_
#define SRC_INCLUDE_SCATTERACTIONSFINDER_H_

#include <memory>
#include <set>
#include <vector>

#include "action.h"
#include "actionfinderfactory.h"
#include "configuration.h"
#include "constants.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();

      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 cell_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(const ParticleData &data_a,
                            const ParticleData &data_b, double dt,
                            const std::vector<FourVector> &beam_momentum = {},
                            const double cell_vol = 0.0) 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 double low_snn_cut_;
  const bool strings_switch_;
  const bool use_AQM_;
  const bool strings_with_probability_;
  const NNbarTreatment nnbar_treatment_;
  const std::vector<bool> &nucleon_has_interacted_;
  const int N_tot_;
  const int N_proj_;
  const double string_formation_time_;
};

}  // namespace smash

#endif  // SRC_INCLUDE_SCATTERACTIONSFINDER_H_