nucleon.h

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file nucleon.h

// TRENTO: Reduced Thickness Event-by-event Nuclear Topology
// Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
// TRENTO3D: Three-dimensional extension of TRENTO by Weiyao Ke
// MIT License

#ifndef NUCLEON_H
#define NUCLEON_H

#include <boost/math/constants/constants.hpp>

#include "fast_exp.h"
#include "fwd_decl.h"
#include "random.h"

namespace trento {

class Nucleon;

class NucleonProfile {
 public:
  explicit NucleonProfile(const VarMap& var_map);

  double radius() const;

  double max_impact() const;

  void fluctuate();

  double thickness(double distance_sqr) const;

  double deterministic_thickness(double distance_sqr) const;

  double norm_Tpp(double bpp_sqr) const;

  bool participate(Nucleon& A, Nucleon& B) const;

 private:
  const double width_sqr_;

  const double trunc_radius_sqr_;

  const double max_impact_sqr_;

  const double neg_one_div_two_width_sqr_;

  const double neg_one_div_four_width_sqr_;

  const double one_div_four_pi_;

  const double cross_sec_param_;

  const FastExp<double> fast_exp_;

  std::gamma_distribution<double> fluct_dist_;

  double prefactor_;
 
  bool with_ncoll_;
};

class Nucleon {
 public:
  Nucleon() = default;

  double x() const;

  double y() const;

  double z() const;

  bool is_participant() const;

 private:
  friend class Nucleus;

  friend bool NucleonProfile::participate(Nucleon&, Nucleon&) const;

  void set_position(double x, double y, double z);

  void set_participant();

  double x_, y_, z_;

  bool participant_;
};

// These functions are short, called very often, and account for a large
// fraction of the total computation time, so request inlining.

// Nucleon inline member functions

inline double Nucleon::x() const {
  return x_;
}

inline double Nucleon::y() const {
  return y_;
}

inline double Nucleon::z() const {
  return z_;
}

inline bool Nucleon::is_participant() const {
  return participant_;
}

inline void Nucleon::set_position(double x, double y, double z) {
  x_ = x;
  y_ = y;
  z_ = z;
  participant_ = false;
}

inline void Nucleon::set_participant() {
  participant_ = true;
}

// NucleonProfile inline member functions

inline double NucleonProfile::radius() const {
  return std::sqrt(trunc_radius_sqr_);
}

inline double NucleonProfile::max_impact() const {
  return std::sqrt(max_impact_sqr_);
}

inline void NucleonProfile::fluctuate() {
  prefactor_ = fluct_dist_(random::engine) *
     math::double_constants::one_div_two_pi / width_sqr_;
}

inline double NucleonProfile::thickness(double distance_sqr) const {
  if (distance_sqr > trunc_radius_sqr_)
    return 0.;
  return prefactor_ * fast_exp_(neg_one_div_two_width_sqr_*distance_sqr);
}

// WK
inline double NucleonProfile::deterministic_thickness(double distance_sqr) const {
  if (distance_sqr > trunc_radius_sqr_)
    return 0.;
  return math::double_constants::one_div_two_pi / width_sqr_ 
                * fast_exp_(neg_one_div_two_width_sqr_*distance_sqr);
}

// WK
inline double NucleonProfile::norm_Tpp(double bpp_sqr) const  {
  return one_div_four_pi_ / width_sqr_ 
                * fast_exp_(neg_one_div_four_width_sqr_*bpp_sqr);
}

inline bool NucleonProfile::participate(Nucleon& A, Nucleon& B) const {
  // If both nucleons are already participants, there's nothing to do, unless
  // in Ncoll mode
  if (A.is_participant() && B.is_participant() && (! with_ncoll_))
    return true;

  double dx = A.x() - B.x();
  double dy = A.y() - B.y();
  double distance_sqr = dx*dx + dy*dy;

  // Check if nucleons are out of range.
  if (distance_sqr > max_impact_sqr_)
    return false;

  // The probability is
  //   P = 1 - exp(...)
  // which we could sample as
  //   P > U
  // where U is a standard uniform (0, 1) random number.  We can also compute
  //   1 - P = exp(...)
  // and then sample
  //   (1 - P) > (1 - U)
  // or equivalently
  //   (1 - P) < U
  auto one_minus_prob = std::exp(
      -std::exp(cross_sec_param_ - .25*distance_sqr/width_sqr_));

  // Sample one random number and decide if this pair participates.
  if (one_minus_prob < random::canonical<double>()) {
    A.set_participant();
    B.set_participant();
    return true;
  }

  return false;
}

}  // namespace trento

#endif  // NUCLEON_H