density.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2013-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
#ifndef SRC_INCLUDE_DENSITY_H_
#define SRC_INCLUDE_DENSITY_H_

#include <iostream>
#include <typeinfo>
#include <utility>
#include <vector>

#include "energymomentumtensor.h"
#include "experimentparameters.h"
#include "forwarddeclarations.h"
#include "fourvector.h"
#include "lattice.h"
#include "particledata.h"
#include "particles.h"
#include "pdgcode.h"
#include "threevector.h"

namespace smash {

enum class DensityType {
  None = 0,
  Hadron = 1,
  Baryon = 2,
  BaryonicIsospin = 3,
  Pion = 4,
  Isospin3_tot = 5,
  Charge = 6,
  Strangeness = 7,
};

std::ostream &operator<<(std::ostream &os, DensityType dt);

double density_factor(const ParticleType &type, DensityType dens_type);

inline double smearing_factor_norm(const double two_sigma_sqr) {
  const double tmp = two_sigma_sqr * M_PI;
  return tmp * std::sqrt(tmp);
}

inline double smearing_factor_rcut_correction(const double rcut_in_sigma) {
  const double x = rcut_in_sigma / std::sqrt(2.0);
  return -2.0 / std::sqrt(M_PI) * x * std::exp(-x * x) + std::erf(x);
}

class DensityParameters {
 public:
  DensityParameters(const ExperimentParameters &par)  // NOLINT
      : sig_(par.gaussian_sigma),
        r_cut_(par.gauss_cutoff_in_sigma * par.gaussian_sigma),
        ntest_(par.testparticles) {
    r_cut_sqr_ = r_cut_ * r_cut_;
    const double two_sig_sqr = 2 * sig_ * sig_;
    two_sig_sqr_inv_ = 1. / two_sig_sqr;
    const double norm = smearing_factor_norm(two_sig_sqr);
    const double corr_factor =
        smearing_factor_rcut_correction(par.gauss_cutoff_in_sigma);
    norm_factor_sf_ = 1. / (norm * ntest_ * corr_factor);
  }
  int ntest() const { return ntest_; }
  double r_cut() const { return r_cut_; }
  double r_cut_sqr() const { return r_cut_sqr_; }
  double two_sig_sqr_inv() const { return two_sig_sqr_inv_; }
  double norm_factor_sf() const { return norm_factor_sf_; }

 private:
  const double sig_;
  const double r_cut_;
  double r_cut_sqr_;
  double two_sig_sqr_inv_;
  double norm_factor_sf_;
  const int ntest_;
};

std::pair<double, ThreeVector> unnormalized_smearing_factor(
    const ThreeVector &r, const FourVector &p, const double m_inv,
    const DensityParameters &dens_par, const bool compute_gradient = false);

std::tuple<double, FourVector, ThreeVector, ThreeVector, ThreeVector>
current_eckart(const ThreeVector &r, const ParticleList &plist,
               const DensityParameters &par, DensityType dens_type,
               bool compute_gradient, bool smearing);
std::tuple<double, FourVector, ThreeVector, ThreeVector, ThreeVector>
current_eckart(const ThreeVector &r, const Particles &plist,
               const DensityParameters &par, DensityType dens_type,
               bool compute_gradient, bool smearing);

class DensityOnLattice {
 public:
  DensityOnLattice()
      : jmu_pos_(FourVector()),
        jmu_neg_(FourVector()),
        djmu_dx_({FourVector(), FourVector(), FourVector(), FourVector()}) {}

  void add_particle(const ParticleData &part, double FactorTimesSf) {
    const FourVector PartFourVelocity = FourVector(1.0, part.velocity());
    if (FactorTimesSf > 0.0) {
      jmu_pos_ += PartFourVelocity * FactorTimesSf;
    } else {
      jmu_neg_ += PartFourVelocity * FactorTimesSf;
    }
  }

  void add_particle_for_derivatives(const ParticleData &part, double factor,
                                    ThreeVector sf_grad) {
    const FourVector PartFourVelocity = FourVector(1.0, part.velocity());
    for (int k = 1; k <= 3; k++) {
      djmu_dx_[k] += factor * PartFourVelocity * sf_grad[k - 1];
      djmu_dx_[0] -=
          factor * PartFourVelocity * sf_grad[k - 1] * part.velocity()[k - 1];
    }
  }

  double density(const double norm_factor = 1.0) {
    return (jmu_pos_.abs() - jmu_neg_.abs()) * norm_factor;
  }

  ThreeVector rot_j(const double norm_factor = 1.0) {
    ThreeVector j_rot = ThreeVector();
    j_rot.set_x1(djmu_dx_[2].x3() - djmu_dx_[3].x2());
    j_rot.set_x2(djmu_dx_[3].x1() - djmu_dx_[1].x3());
    j_rot.set_x3(djmu_dx_[1].x2() - djmu_dx_[2].x1());
    j_rot *= norm_factor;
    return j_rot;
  }

  ThreeVector grad_rho(const double norm_factor = 1.0) {
    ThreeVector rho_grad = ThreeVector();
    for (int i = 1; i < 4; i++) {
      rho_grad[i - 1] = djmu_dx_[i].x0() * norm_factor;
    }
    return rho_grad;
  }

  ThreeVector dj_dt(const double norm_factor = 1.0) {
    return djmu_dx_[0].threevec() * norm_factor;
  }

  FourVector jmu_net() const { return jmu_pos_ + jmu_neg_; }

 private:
  FourVector jmu_pos_;
  FourVector jmu_neg_;
  std::array<FourVector, 4> djmu_dx_;
};

typedef RectangularLattice<DensityOnLattice> DensityLattice;

template <typename T>
void update_lattice(RectangularLattice<T> *lat, const LatticeUpdate update,
                    const DensityType dens_type, const DensityParameters &par,
                    const Particles &particles,
                    const bool compute_gradient = false) {
  // Do not proceed if lattice does not exists/update not required
  if (lat == nullptr || lat->when_update() != update) {
    return;
  }
  lat->reset();
  const double norm_factor = par.norm_factor_sf();
  for (const auto &part : particles) {
    const double dens_factor = density_factor(part.type(), dens_type);
    if (std::abs(dens_factor) < really_small) {
      continue;
    }
    const FourVector p = part.momentum();
    const double m = p.abs();
    if (unlikely(m < really_small)) {
      const auto &log = logger<LogArea::Density>();
      log.warn("Gaussian smearing is undefined for momentum ", p);
      continue;
    }
    const double m_inv = 1.0 / m;

    const ThreeVector pos = part.position().threevec();
    lat->iterate_in_radius(
        pos, par.r_cut(), [&](T &node, int ix, int iy, int iz) {
          const ThreeVector r = lat->cell_center(ix, iy, iz);
          const auto sf = unnormalized_smearing_factor(pos - r, p, m_inv, par,
                                                       compute_gradient);
          if (sf.first * norm_factor > really_small) {
            node.add_particle(part, sf.first * norm_factor * dens_factor);
          }
          if (compute_gradient) {
            node.add_particle_for_derivatives(part, dens_factor,
                                              sf.second * norm_factor);
          }
        });
  }
}

}  // namespace smash

#endif  // SRC_INCLUDE_DENSITY_H_