20 : sig_(param.gaussian_sigma),
21 rc_(conf.take({
"Gaussian_Cutoff"}, 2.2)),
22 rr_(conf.take({
"Spatial_Averaging_Radius"}, 1.86)),
23 rp_(conf.take({
"Momentum_Averaging_Radius"}, 0.08)),
24 ntest_(param.testparticles),
25 n_ensembles_(param.n_ensembles) {
26 if (ntest_ * n_ensembles_ < 20) {
28 "Phase-space density calculation in Pauli blocking will not work "
29 "reasonably.\nEither use testparticles or ensembles, or both.\n"
30 "The recommended testparticles * ensembles is at least 20.");
33 if (rc_ < rr_ || rr_ < 0.0 || rp_ < 0) {
35 "Please choose reasonable parameters for Pauli blocking:\n"
36 "All radii have to be positive and Gaussian_Cutoff should\n"
37 "be larger than Spatial_Averaging_Radius");
40 init_weights_analytical();
46 const std::vector<Particles> &ensembles,
48 const ParticleList &disregard)
const {
54 for (
const Particles &particles : ensembles) {
57 if (part.pdgcode() != pdg) {
61 const double pdist_sqr = (part.momentum().threevec() -
p).sqr();
62 if (pdist_sqr >
rp_ *
rp_) {
65 const double rdist_sqr = (part.position().threevec() - r).sqr();
71 bool to_disregard =
false;
72 for (
const auto &disregard_part : disregard) {
73 if (part.id() == disregard_part.id()) {
83 const size_t i = std::floor(i_real);
84 const double rest = i_real - i;
94 const double pi = M_PI;
95 const double sqrt2 = std::sqrt(2.);
96 const double sqrt_2pi = std::sqrt(2. * pi);
98 const double phase_volume =
108 const double d_pos = (
rr_ +
rc_) /
static_cast<double>(
weights_.size());
110 for (
size_t k = 0; k <
weights_.size(); k++) {
112 const double rj = d_pos * k;
115 const double A =
rr_ / sqrt2 /
sig_;
116 integral = sqrt_2pi *
sig_ * std::erf(A) - 2 *
rr_ * std::exp(-A * A);
118 }
else if (
rc_ > rj +
rr_) {
119 const double A = (rj +
rr_) / sqrt2 /
sig_;
120 const double B = (rj -
rr_) / sqrt2 /
sig_;
121 integral =
sig_ / rj * (std::exp(-A * A) - std::exp(-B * B)) +
122 0.5 * sqrt_2pi * (std::erf(A) - std::erf(B));
125 const double A =
rc_ / sqrt2 /
sig_;
126 const double B = (rj -
rr_) / sqrt2 /
sig_;
129 (0.5 * std::exp(-A * A) * C -
sig_ *
sig_ * std::exp(-B * B)) / rj +
130 0.5 * sqrt_2pi *
sig_ * (std::erf(A) - std::erf(B));
133 integral *= 2 * pi / std::pow(2 * pi *
sig_ *
sig_, 1.5);
134 weights_[k] = integral / norm / phase_volume;
Interface to the SMASH configuration files.
ParticleData contains the dynamic information of a certain particle.
The Particles class abstracts the storage and manipulation of particles.
void init_weights_analytical()
Analytical calculation of weights.
double rr_
Radius of averaging in coordinate space, fm.
int n_ensembles_
Number of ensembles.
double sig_
Standard deviation of the gaussian used for smearing.
~PauliBlocker()
Destructor.
int ntest_
Testparticles number.
double rc_
Radius, after which gaussians (used for averaging) are cut, fm.
double phasespace_dens(const ThreeVector &r, const ThreeVector &p, const std::vector< Particles > &ensembles, const PdgCode pdg, const ParticleList &disregard) const
Calculate phase-space density of a particle species at the point (r,p).
PauliBlocker(Configuration conf, const ExperimentParameters ¶meters)
PauliBlocker constructor.
double rp_
Radius of averaging in momentum space, GeV.
std::array< double, 30 > weights_
Weights: tabulated results of numerical integration.
PdgCode stores a Particle Data Group Particle Numbering Scheme particle type number.
The ThreeVector class represents a physical three-vector with the components .
Collection of useful constants that are known at compile time.
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
#define likely(x)
Tell the branch predictor that this expression is likely true.
static constexpr int LPauliBlocking
constexpr double hbarc
GeV <-> fm conversion factor.
constexpr double really_small
Numerical error tolerance.
Helper structure for Experiment.