customnucleus.cc

Go to the documentation of this file.

/*
 *    Copyright (c) 2014-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
#include <cmath>
#include <fstream>
#include <map>
#include <string>
#include <vector>
 
#include "smash/constants.h"
#include "smash/customnucleus.h"
#include "smash/particledata.h"
#include "smash/particletype.h"
#include "smash/pdgcode.h"
#include "smash/pdgcode_constants.h"
#include "smash/random.h"
 
namespace smash {
static constexpr int LCollider = LogArea::Collider::id;
 
std::unique_ptr<std::ifstream> CustomNucleus::filestream_shared_ = nullptr;
 
CustomNucleus::CustomNucleus(Configuration& config, int testparticles,
                             bool same_file) {
  // Read in file directory from config
  std::string particle_list_file_directory =
      config.take({"Custom", "File_Directory"});
  // Read in file name from config
  std::string particle_list_file_name = config.take({"Custom", "File_Name"});
 
  if (particles_.size() != 0) {
    throw std::runtime_error(
        "Your Particle List is already filled before reading in from the "
        "external file."
        "Something went wrong. Please check your config.");
  }
  /*
   * Counts number of nucleons in one nucleus as it is specialized
   * by the user in the config file.
   * It is needed to read in the proper number of nucleons for one
   * nucleus and to restart at the listreading for the following
   * nucleus as one does not want to read configurations twice.
   */
  std::map<PdgCode, int> particle_list = config.take({"Particles"});
  for (const auto& particle : particle_list)
    number_of_nucleons_ += particle.second * testparticles;
  /*
   * "if" statement makes sure the streams to the file are initialized
   * properly.
   */
  const std::string path =
      file_path(particle_list_file_directory, particle_list_file_name);
  if (same_file && !filestream_shared_) {
    filestream_shared_ = make_unique<std::ifstream>(path);
    used_filestream_ = &filestream_shared_;
  } else if (!same_file) {
    filestream_ = make_unique<std::ifstream>(path);
    used_filestream_ = &filestream_;
  } else {
    used_filestream_ = &filestream_shared_;
  }
 
  custom_nucleus_ = readfile(**used_filestream_, number_of_nucleons_);
  fill_from_list(custom_nucleus_);
  // Inherited from nucleus class (see nucleus.h)
  set_parameters_automatic();
}
 
void CustomNucleus::fill_from_list(const std::vector<Nucleoncustom>& vec) {
  particles_.clear();
  index = 0;
  // checking if particle is proton or neutron
  for (const auto& it : vec) {
    PdgCode pdgcode;
    if (it.isospin == 1) {
      pdgcode = pdg::p;
    } else if (it.isospin == 0) {
      pdgcode = pdg::n;
    } else {
      throw std::runtime_error(
          "Your particles charges are not 1 = proton or 0 = neutron."
          "Check whether your list is correct or there is an error.");
    }
    // setting parameters for the particles in the particlelist in smash
    const ParticleType& current_type = ParticleType::find(pdgcode);
    double current_mass = current_type.mass();
    particles_.emplace_back(current_type);
    particles_.back().set_4momentum(current_mass, 0.0, 0.0, 0.0);
  }
}
 
ThreeVector CustomNucleus::distribute_nucleon() {
  /*
   * As only arrange_nucleons is called at the beginning of every
   * event it is important to have readfile and fill from list
   * called again when a new event starts. The constructor is only
   * called twice to initialize the first target and projectile.
   * Therefore this if statement is implemented.
   */
  if (index >= custom_nucleus_.size()) {
    custom_nucleus_ = readfile(**used_filestream_, number_of_nucleons_);
    fill_from_list(custom_nucleus_);
  }
  const auto& pos = custom_nucleus_.at(index);
  index++;
  ThreeVector nucleon_position(pos.x, pos.y, pos.z);
  // rotate nucleon about euler angle
  nucleon_position.rotate(euler_phi_, euler_theta_, euler_psi_);
 
  return nucleon_position;
}
 
void CustomNucleus::arrange_nucleons() {
  /* Randomly generate Euler angles for rotation everytime a new
   * custom nucleus is initialiezed. Therefore this is done 2 times per
   * event.
   */
  Nucleus::random_euler_angles();
 
  for (auto i = begin(); i != end(); i++) {
    // Initialize momentum
    i->set_4momentum(i->pole_mass(), 0.0, 0.0, 0.0);
    /* Sampling the Woods-Saxon, get the radial
     * position and solid angle for the nucleon. */
    ThreeVector pos = distribute_nucleon();
    // Set the position of the nucleon.
    i->set_4position(FourVector(0.0, pos));
  }
  // Recenter
  align_center();
}
 
void CustomNucleus::generate_fermi_momenta() {
  Nucleus::generate_fermi_momenta();
  logg[LCollider].warn() << "Fermi motion activated with a custom nucleus.\n";
  logg[LCollider].warn() << "Be aware that generating the Fermi momenta\n"
                         << "assumes nucleons distributed according to a\n"
                         << "Woods-Saxon distribution.";
}
 
std::string CustomNucleus::file_path(const std::string& file_directory,
                                     const std::string& file_name) {
  if (file_directory.back() == '/') {
    return file_directory + file_name;
  } else {
    return file_directory + '/' + file_name;
  }
}
 
std::vector<Nucleoncustom> CustomNucleus::readfile(std::ifstream& infile,
                                                   int particle_number) const {
  int A = particle_number;
  std::string line;
  std::vector<Nucleoncustom> custom_nucleus;
  // read in only A particles for one nucleus
  for (int i = 0; i < A; ++i) {
    std::getline(infile, line);
    // make sure the stream goes back to the beginning when it hits end of file
    if (infile.eof()) {
      infile.clear();
      infile.seekg(0, infile.beg);
      std::getline(infile, line);
    }
    Nucleoncustom nucleon;
    std::istringstream iss(line);
    if (!(iss >> nucleon.x >> nucleon.y >> nucleon.z >>
          nucleon.spinprojection >> nucleon.isospin)) {
      throw std::runtime_error(
          "SMASH could not read in a line from your initial nuclei input file."
          "Check if your file has the following format: x y z spinprojection "
          "isospin");
      break;
    }
    custom_nucleus.push_back(nucleon);
  }
  return custom_nucleus;
}
 
}  // namespace smash