nucleus.h

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file nucleus.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 NUCLEUS_H
#define NUCLEUS_H

#include <memory>
#include <random>
#include <string>
#include <vector>

#include "fwd_decl.h"
#include "nucleon.h"

#ifdef TRENTO_HDF5
// forward declaration for std::unique_ptr<H5::DataSet> member of ManualNucleus
namespace H5 {
class DataSet;
}
#endif

namespace trento {

// Alias for a smart pointer to a Nucleus.
using NucleusPtr = std::unique_ptr<Nucleus>;

class Nucleus {
 public:
  static NucleusPtr create(const std::string& species, double nucleon_width, double nucleon_dmin = 0);

  virtual ~Nucleus() = default;

  virtual double radius() const = 0;

  void sample_nucleons(double offset);

  using size_type = std::vector<Nucleon>::size_type;
  using iterator = std::vector<Nucleon>::iterator;
  using const_iterator = std::vector<Nucleon>::const_iterator;

  // size = number of nucleons
  size_type size() const noexcept
  { return nucleons_.size(); }

  // non-const overload
  iterator begin() noexcept
  { return nucleons_.begin(); }
  iterator end() noexcept
  { return nucleons_.end(); }

  // const overload
  const_iterator begin() const noexcept
  { return nucleons_.begin(); }
  const_iterator end() const noexcept
  { return nucleons_.end(); }

  // forced const
  const_iterator cbegin() const noexcept
  { return nucleons_.cbegin(); }
  const_iterator cend() const noexcept
  { return nucleons_.cend(); }

 protected:
  explicit Nucleus(std::size_t A);

  void set_nucleon_position(iterator nucleon, double x, double y, double z);

 private:
  virtual void sample_nucleons_impl() = 0;

  std::vector<Nucleon> nucleons_;

  double offset_;
};

// Now declare Nucleus subclasses.
// Each subclass must do the following:
//   - Define a public constructor which (at minimum) initializes a Nucleus
//     with the required number of nucleons.
//   - Override and implement the pure virtual functions
//     radius() and sample_nucleons_impl().

class Proton : public Nucleus {
 public:
  Proton();

  virtual double radius() const override;

 private:
  virtual void sample_nucleons_impl() override;
};

class Deuteron : public Nucleus {
 public:
  Deuteron(double a = 0.457, double b = 2.35);

  virtual double radius() const override;

 private:
  virtual void sample_nucleons_impl() override;

  const double a_, b_;
};

class MinDistNucleus : public Nucleus {
 protected:
  MinDistNucleus(std::size_t A, double dmin = 0);

  bool is_too_close(const_iterator nucleon) const;

 private:
  const double dminsq_;
};

class WoodsSaxonNucleus : public MinDistNucleus {
 public:
  WoodsSaxonNucleus(std::size_t A, double R, double a, double dmin = 0);

  virtual double radius() const override;

 private:
  virtual void sample_nucleons_impl() override;

  const double R_, a_;

  mutable std::piecewise_linear_distribution<double> woods_saxon_dist_;
};

class DeformedWoodsSaxonNucleus : public MinDistNucleus {
 public:
  DeformedWoodsSaxonNucleus(std::size_t A, double R, double a,
                            double beta2, double beta4, double dmin = 0);

  virtual double radius() const override;

 private:
  virtual void sample_nucleons_impl() override;

  double deformed_woods_saxon_dist(double r, double cos_theta) const;

  const double R_, a_, beta2_, beta4_;

  const double rmax_;
};

#ifdef TRENTO_HDF5

class ManualNucleus : public Nucleus {
 public:
  static std::unique_ptr<ManualNucleus> create(const std::string& path);

  virtual ~ManualNucleus() override;

  virtual double radius() const override;

 private:
  ManualNucleus(std::unique_ptr<H5::DataSet> dataset,
                std::size_t nconfigs, std::size_t A, double rmax);

  virtual void sample_nucleons_impl() override;

  const std::unique_ptr<H5::DataSet> dataset_;

  const double rmax_;

  std::uniform_int_distribution<std::size_t> index_dist_;
};

#endif  // TRENTO_HDF5

}  // namespace trento

#endif  // NUCLEUS_H