Version: SMASH-3.1
pauliblocking.cc
Go to the documentation of this file.
1 /*
2  *
3  * Copyright (c) 2015-2020,2022
4  * SMASH Team
5  *
6  * GNU General Public License (GPLv3 or later)
7  *
8  */
9 
10 #include "smash/pauliblocking.h"
11 
12 #include "smash/constants.h"
13 #include "smash/logging.h"
14 
15 namespace smash {
16 static constexpr int LPauliBlocking = LogArea::PauliBlocking::id;
17 
19  const ExperimentParameters &param)
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) {
27  logg[LPauliBlocking].warn(
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.");
31  }
32 
33  if (rc_ < rr_ || rr_ < 0.0 || rp_ < 0) {
34  logg[LPauliBlocking].error(
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");
38  }
39 
40  init_weights_analytical();
41 }
42 
44 
46  const std::vector<Particles> &ensembles,
47  const PdgCode pdg,
48  const ParticleList &disregard) const {
49  double f = 0.0;
50 
51  /* TODO(oliiny): looping over all particles is inefficient,
52  * I need only particles within rp_ radius in momentum and
53  * within rr_+rc_ in coordinate space. Some search algorithm might help. */
54  for (const Particles &particles : ensembles) {
55  for (const ParticleData &part : particles) {
56  // Only consider identical particles
57  if (part.pdgcode() != pdg) {
58  continue;
59  }
60  // Only consider momenta in sphere of radius rp_ with center at p
61  const double pdist_sqr = (part.momentum().threevec() - p).sqr();
62  if (pdist_sqr > rp_ * rp_) {
63  continue;
64  }
65  const double rdist_sqr = (part.position().threevec() - r).sqr();
66  // Only consider coordinates in sphere of radius rr_+rc_ with center at r
67  if (rdist_sqr >= (rr_ + rc_) * (rr_ + rc_)) {
68  continue;
69  }
70  // Do not count particles that should be disregarded.
71  bool to_disregard = false;
72  for (const auto &disregard_part : disregard) {
73  if (part.id() == disregard_part.id()) {
74  to_disregard = true;
75  }
76  }
77  if (to_disregard) {
78  continue;
79  }
80  // 1st order interpolation using tabulated values
81  const double i_real =
82  std::sqrt(rdist_sqr) / (rr_ + rc_) * weights_.size();
83  const size_t i = std::floor(i_real);
84  const double rest = i_real - i;
85  if (likely(i + 1 < weights_.size())) {
86  f += weights_[i] * rest + weights_[i + 1] * (1. - rest);
87  }
88  } // loop over particles in one ensemble
89  } // loop over ensembles
90  return f / ntest_ / n_ensembles_;
91 }
92 
94  const double pi = M_PI;
95  const double sqrt2 = std::sqrt(2.);
96  const double sqrt_2pi = std::sqrt(2. * pi);
97  // Volume of the phase-space area; Factor 2 stands for spin.
98  const double phase_volume =
99  2 * (4. / 3. * pi * rr_ * rr_ * rr_) * (4. / 3. * pi * rp_ * rp_ * rp_) /
100  ((2 * pi * hbarc) * (2 * pi * hbarc) * (2 * pi * hbarc));
101  // Analytical expression for integral in denominator
102  const double norm =
103  std::erf(rc_ / sqrt2 / sig_) -
104  rc_ * 2 / sqrt_2pi / sig_ * std::exp(-0.5 * rc_ * rc_ / sig_ / sig_);
105 
106  double integral;
107  // Step of the table for tabulated integral
108  const double d_pos = (rr_ + rc_) / static_cast<double>(weights_.size());
109 
110  for (size_t k = 0; k < weights_.size(); k++) {
111  // rdist = 0 ... rc_ (gauss cut) + rr_ (position cut)
112  const double rj = d_pos * k;
113  if (rj < really_small) {
114  // Assuming rc_ > rr_
115  const double A = rr_ / sqrt2 / sig_;
116  integral = sqrt_2pi * sig_ * std::erf(A) - 2 * rr_ * std::exp(-A * A);
117  integral *= sig_ * sig_;
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));
123  integral *= sig_ * sig_ * sig_;
124  } else {
125  const double A = rc_ / sqrt2 / sig_;
126  const double B = (rj - rr_) / sqrt2 / sig_;
127  const double C = (rc_ - rj) * (rc_ - rj) - rr_ * rr_ + 2 * sig_ * sig_;
128  integral =
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));
131  integral *= sig_ * sig_;
132  }
133  integral *= 2 * pi / std::pow(2 * pi * sig_ * sig_, 1.5);
134  weights_[k] = integral / norm / phase_volume;
135  logg[LPauliBlocking].debug("Analytical weights[", k, "] = ", weights_[k]);
136  }
137 }
138 
139 } // namespace smash
Interface to the SMASH configuration files.
ParticleData contains the dynamic information of a certain particle.
Definition: particledata.h:58
The Particles class abstracts the storage and manipulation of particles.
Definition: particles.h:33
void init_weights_analytical()
Analytical calculation of weights.
double rr_
Radius of averaging in coordinate space, fm.
Definition: pauliblocking.h:85
int n_ensembles_
Number of ensembles.
Definition: pauliblocking.h:94
double sig_
Standard deviation of the gaussian used for smearing.
Definition: pauliblocking.h:79
~PauliBlocker()
Destructor.
int ntest_
Testparticles number.
Definition: pauliblocking.h:91
double rc_
Radius, after which gaussians (used for averaging) are cut, fm.
Definition: pauliblocking.h:82
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 &parameters)
PauliBlocker constructor.
double rp_
Radius of averaging in momentum space, GeV.
Definition: pauliblocking.h:88
std::array< double, 30 > weights_
Weights: tabulated results of numerical integration.
Definition: pauliblocking.h:97
PdgCode stores a Particle Data Group Particle Numbering Scheme particle type number.
Definition: pdgcode.h:111
The ThreeVector class represents a physical three-vector with the components .
Definition: threevector.h:31
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.
Definition: logging.cc:39
#define likely(x)
Tell the branch predictor that this expression is likely true.
Definition: macros.h:14
constexpr int p
Proton.
Definition: action.h:24
static constexpr int LPauliBlocking
Definition: action.cc:27
constexpr double hbarc
GeV <-> fm conversion factor.
Definition: constants.h:25
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:37
Helper structure for Experiment.