boxmodus.cc

Go to the documentation of this file.

/*
 *    Copyright (c) 2013-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <list>
#include <map>
#include <utility>
#include <vector>

#include "smash/algorithms.h"
#include "smash/angles.h"
#include "smash/boxmodus.h"
#include "smash/configuration.h"
#include "smash/constants.h"
#include "smash/cxx14compat.h"
#include "smash/distributions.h"
#include "smash/experimentparameters.h"
#include "smash/hadgas_eos.h"
#include "smash/logging.h"
#include "smash/macros.h"
#include "smash/outputinterface.h"
#include "smash/particles.h"
#include "smash/processbranch.h"
#include "smash/quantumnumbers.h"
#include "smash/random.h"
#include "smash/threevector.h"
#include "smash/wallcrossingaction.h"

namespace smash {

/* console output on startup of box specific parameters */
std::ostream &operator<<(std::ostream &out, const BoxModus &m) {
  out << "-- Box Modus:\nSize of the box: (" << m.length_ << " fm)³\n";
  if (m.use_thermal_) {
    out << "Thermal multiplicities "
        << "(T = " << m.temperature_ << " GeV, muB = " << m.mub_
        << " GeV, muS = " << m.mus_ << " GeV)\n";
  } else {
    for (const auto &p : m.init_multipl_) {
      ParticleTypePtr ptype = &ParticleType::find(p.first);
      out << ptype->name() << " initial multiplicity " << p.second << '\n';
    }
  }
  if (m.initial_condition_ == BoxInitialCondition::PeakedMomenta) {
    out << "All initial momenta = 3T = " << 3 * m.temperature_ << " GeV\n";
  } else {
    out << "Boltzmann momentum distribution with T = " << m.temperature_
        << " GeV.\n";
  }
  if (m.insert_jet_) {
    ParticleTypePtr ptype = &ParticleType::find(m.jet_pdg_);
    out << "Adding a " << ptype->name() << " as a jet in the middle "
        << "of the box with " << m.jet_mom_ << " GeV initial momentum.\n";
  }
  return out;
}

BoxModus::BoxModus(Configuration modus_config, const ExperimentParameters &)
    : initial_condition_(modus_config.take({"Box", "Initial_Condition"})),
      length_(modus_config.take({"Box", "Length"})),
      temperature_(modus_config.take({"Box", "Temperature"})),
      start_time_(modus_config.take({"Box", "Start_Time"}, 0.)),
      use_thermal_(
          modus_config.take({"Box", "Use_Thermal_Multiplicities"}, false)),
      mub_(modus_config.take({"Box", "Baryon_Chemical_Potential"}, 0.)),
      mus_(modus_config.take({"Box", "Strange_Chemical_Potential"}, 0.)),
      account_for_resonance_widths_(
          modus_config.take({"Box", "Account_Resonance_Widths"}, true)),
      init_multipl_(use_thermal_
                        ? std::map<PdgCode, int>()
                        : modus_config.take({"Box", "Init_Multiplicities"})
                              .convert_for(init_multipl_)),
      insert_jet_(modus_config.has_value({"Box", "Jet", "Jet_PDG"})),
      jet_pdg_(insert_jet_ ? modus_config.take({"Box", "Jet", "Jet_PDG"})
                                 .convert_for(jet_pdg_)
                           : pdg::p),  // dummy default; never used
      jet_mom_(modus_config.take({"Box", "Jet", "Jet_Momentum"}, 20.)) {}

double BoxModus::initial_conditions(Particles *particles,
                                    const ExperimentParameters &parameters) {
  const auto &log = logger<LogArea::Box>();
  double momentum_radial = 0, mass;
  Angles phitheta;
  FourVector momentum_total(0, 0, 0, 0);
  auto uniform_length = random::make_uniform_distribution(0.0, this->length_);
  const double T = this->temperature_;
  /* Create NUMBER OF PARTICLES according to configuration, or thermal case */
  if (use_thermal_) {
    const double V = length_ * length_ * length_;
    if (average_multipl_.empty()) {
      for (const ParticleType &ptype : ParticleType::list_all()) {
        if (HadronGasEos::is_eos_particle(ptype)) {
          const double n = HadronGasEos::partial_density(
              ptype, T, mub_, mus_, account_for_resonance_widths_);
          average_multipl_[ptype.pdgcode()] = n * V * parameters.testparticles;
        }
      }
    }
    double nb_init = 0.0, ns_init = 0.0;
    for (const auto &mult : average_multipl_) {
      const int thermal_mult_int = random::poisson(mult.second);
      particles->create(thermal_mult_int, mult.first);
      nb_init += mult.second * mult.first.baryon_number();
      ns_init += mult.second * mult.first.strangeness();
      log.debug(mult.first, " initial multiplicity ", thermal_mult_int);
    }
    log.info("Initial hadron gas baryon density ", nb_init);
    log.info("Initial hadron gas strange density ", ns_init);
  } else {
    for (const auto &p : init_multipl_) {
      particles->create(p.second * parameters.testparticles, p.first);
      log.debug("Particle ", p.first, " initial multiplicity ", p.second);
    }
  }

  for (ParticleData &data : *particles) {
    /* Set MOMENTUM SPACE distribution */
    if (this->initial_condition_ == BoxInitialCondition::PeakedMomenta) {
      /* initial thermal momentum is the average 3T */
      momentum_radial = 3.0 * T;
      mass = data.pole_mass();
    } else {
      /* thermal momentum according Maxwell-Boltzmann distribution */
      mass = (!account_for_resonance_widths_)
                 ? data.type().mass()
                 : HadronGasEos::sample_mass_thermal(data.type(), 1.0 / T);
      momentum_radial = sample_momenta_from_thermal(T, mass);
    }
    phitheta.distribute_isotropically();
    log.debug(data.type().name(), "(id ", data.id(), ") radial momentum ",
              momentum_radial, ", direction", phitheta);
    data.set_4momentum(mass, phitheta.threevec() * momentum_radial);
    momentum_total += data.momentum();

    /* Set COORDINATE SPACE distribution */
    ThreeVector pos{uniform_length(), uniform_length(), uniform_length()};
    data.set_4position(FourVector(start_time_, pos));
    data.set_formation_time(start_time_);
  }

  /* Make total 3-momentum 0 */
  for (ParticleData &data : *particles) {
    data.set_4momentum(data.momentum().abs(),
                       data.momentum().threevec() -
                           momentum_total.threevec() / particles->size());
  }

  /* Add a single highly energetic particle in the center of the box (jet) */
  if (insert_jet_) {
    auto &jet_particle = particles->create(jet_pdg_);
    jet_particle.set_formation_time(start_time_);
    jet_particle.set_4position(FourVector(start_time_, 0., 0., 0.));
    jet_particle.set_4momentum(ParticleType::find(jet_pdg_).mass(),
                               ThreeVector(jet_mom_, 0., 0.));
  }

  /* Recalculate total momentum */
  momentum_total = FourVector(0, 0, 0, 0);
  for (ParticleData &data : *particles) {
    momentum_total += data.momentum();
    /* IC: debug checks */
    log.debug() << data;
  }
  /* allows to check energy conservation */
  log.debug() << "Initial total 4-momentum [GeV]: " << momentum_total;
  return start_time_;
}

int BoxModus::impose_boundary_conditions(Particles *particles,
                                         const OutputsList &output_list) {
  const auto &log = logger<LogArea::Box>();
  int wraps = 0;

  for (ParticleData &data : *particles) {
    FourVector position = data.position();
    bool wall_hit = enforce_periodic_boundaries(position.begin() + 1,
                                                position.end(), length_);
    if (wall_hit) {
      const ParticleData incoming_particle(data);
      data.set_4position(position);
      ++wraps;
      ActionPtr action =
          make_unique<WallcrossingAction>(incoming_particle, data);
      for (const auto &output : output_list) {
        if (!output->is_dilepton_output() && !output->is_photon_output()) {
          output->at_interaction(*action, 0.);
        }
      }
    }
  }
  log.debug("Moved ", wraps, " particles back into the box.");
  return wraps;
}

}  // namespace smash