9 #ifndef SRC_INCLUDE_DENSITY_H_ 10 #define SRC_INCLUDE_DENSITY_H_ 73 const double tmp = two_sigma_sqr * M_PI;
74 return tmp * std::sqrt(tmp);
95 const double x = rcut_in_sigma / std::sqrt(2.0);
96 return -2.0 / std::sqrt(M_PI) * x * std::exp(-x * x) + std::erf(x);
114 :
sig_(par.gaussian_sigma),
115 r_cut_(par.gauss_cutoff_in_sigma * par.gaussian_sigma),
116 ntest_(par.testparticles) {
118 const double two_sig_sqr = 2 *
sig_ *
sig_;
121 const double corr_factor =
213 std::tuple<double, ThreeVector, ThreeVector, ThreeVector>
rho_eckart(
217 std::tuple<double, ThreeVector, ThreeVector, ThreeVector>
rho_eckart(
264 if (FactorTimesSf > 0.0) {
265 jmu_pos_ += PartFourVelocity * FactorTimesSf;
267 jmu_neg_ += PartFourVelocity * FactorTimesSf;
286 for (
int k = 1; k <= 3; k++) {
287 djmu_dx_[k] += factor * PartFourVelocity * sf_grad[k - 1];
289 factor * PartFourVelocity * sf_grad[k - 1] * part.
velocity()[k - 1];
310 double density(
const double norm_factor = 1.0) {
325 j_rot *= norm_factor;
337 for (
int i = 1; i < 4; i++) {
338 rho_grad[i - 1] =
djmu_dx_[i].x0() * norm_factor;
350 return djmu_dx_[0].threevec() * norm_factor;
385 template <
typename T>
389 const bool compute_gradient =
false) {
391 if (lat ==
nullptr || lat->
when_update() != update) {
396 for (
const auto &part : particles) {
397 const double dens_factor =
density_factor(part.type(), dens_type);
402 const double m =
p.abs();
404 const auto &log = logger<LogArea::Density>();
405 log.warn(
"Gaussian smearing is undefined for momentum ",
p);
408 const double m_inv = 1.0 / m;
410 const ThreeVector pos = part.position().threevec();
412 pos, par.
r_cut(), [&](T &node,
int ix,
int iy,
int iz) {
417 node.add_particle(part, sf.first * norm_factor * dens_factor);
419 if (compute_gradient) {
420 node.add_particle_for_derivatives(part, dens_factor,
421 sf.second * norm_factor);
429 #endif // SRC_INCLUDE_DENSITY_H_ A class to pre-calculate and store parameters relevant for density calculation.
std::array< FourVector, 4 > djmu_dx_
A class for time-efficient (time-memory trade-off) calculation of density on the lattice.
The ThreeVector class represents a physical three-vector with the components .
DensityParameters(const ExperimentParameters &par)
Constructor of DensityParameters.
constexpr double really_small
Numerical error tolerance.
double norm_factor_sf() const
FourVector jmu_neg_
Four-current density of the negatively charged particle.
const int ntest_
Testparticle number.
void add_particle(const ParticleData &part, double FactorTimesSf)
Adds particle to 4-current: .
const double sig_
Gaussian smearing width [fm].
void update_lattice(RectangularLattice< T > *lat, const LatticeUpdate update, const DensityType dens_type, const DensityParameters &par, const Particles &particles, const bool compute_gradient=false)
Updates the contents on the lattice.
std::tuple< double, ThreeVector, ThreeVector, ThreeVector > rho_eckart(const ThreeVector &r, const ParticleList &plist, const DensityParameters &par, DensityType dens_type, bool compute_gradient)
Calculates Eckart rest frame density and optionally the gradient of the density in an arbitary frame...
FourVector jmu_net() const
LatticeUpdate
Enumerator option for lattice updates.
double two_sig_sqr_inv() const
void iterate_in_radius(const ThreeVector &point, const double r_cut, F &&func)
Iterates only nodes, whose cell centers lie not further than r_cut in x, y, z directions from the giv...
LatticeUpdate when_update() const
void set_x1(double x)
set first component
double two_sig_sqr_inv_
[fm ]
A container class to hold all the arrays on the lattice and access them.
void reset()
Sets all values on lattice to zeros.
DensityOnLattice()
Default constructor.
double gauss_cutoff_in_sigma
Distance at which gaussian is cut, i.e. set to zero, IN SIGMA (not fm)
ThreeVector rot_j(const double norm_factor=1.0)
Compute curl of the current on the local lattice.
double norm_factor_sf_
Normalization for smearing factor.
double smearing_factor_rcut_correction(const double rcut_in_sigma)
Gaussians used for smearing are cut at radius for calculation speed-up.
double density_factor(const ParticleType &type, DensityType dens_type)
Get the factor that determines how much a particle contributes to the density type that is computed...
ThreeVector dj_dt(const double norm_factor=1.0)
Compute time derivative of the current density on the local lattice.
void set_x3(double z)
set third component
ThreeVector cell_center(int ix, int iy, int iz) const
Find the coordinates of a given cell.
double smearing_factor_norm(const double two_sigma_sqr)
Norm of the Gaussian smearing function.
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)
Implements gaussian smearing for any quantity.
const double r_cut_
Cut-off radius [fm].
ThreeVector velocity() const
Get the velocity 3-vector.
#define unlikely(x)
Tell the branch predictor that this expression is likely false.
FourVector jmu_pos_
Four-current density of the positively charged particle.
The Particles class abstracts the storage and manipulation of particles.
DensityType
Allows to choose which kind of density to calculate.
std::ostream & operator<<(std::ostream &out, const ActionPtr &action)
Convenience: dereferences the ActionPtr to Action.
ThreeVector grad_rho(const double norm_factor=1.0)
Compute gradient of the density on the local lattice.
The FourVector class holds relevant values in Minkowski spacetime with (+, −, −, −) metric signature.
double density(const double norm_factor=1.0)
Compute the net Eckart density on the local lattice.
double abs() const
calculate the lorentz invariant absolute value
RectangularLattice< DensityOnLattice > DensityLattice
Conveniency typedef for lattice of density.
Helper structure for Experiment.
ParticleData contains the dynamic information of a certain particle.
void set_x2(double y)
set second component
double r_cut_sqr_
Squared cut-off radius [fm ].
void add_particle_for_derivatives(const ParticleData &part, double factor, ThreeVector sf_grad)
Adds particle to the time and spatial derivatives of the 4-current.