17 const ExperimentParameters ¶m)
18 : sig_(param.gaussian_sigma),
19 rc_(conf.take({
"Gaussian_Cutoff"}, 2.2)),
20 rr_(conf.take({
"Spatial_Averaging_Radius"}, 1.86)),
21 rp_(conf.take({
"Momentum_Averaging_Radius"}, 0.08)),
22 ntest_(param.testparticles) {
36 const auto &log = logger<LogArea::PauliBlocking>();
40 "Phase-space density calculation in Pauli blocking" 41 " will not work reasonably for a small number of testparticles." 42 " The recommended number of testparticles is 20.");
45 if (rc_ < rr_ || rr_ < 0.0 || rp_ < 0) {
47 "Please choose reasonable parameters for Pauli blocking:" 48 "All radii have to be positive and Gaussian_Cutoff should" 49 "be larger than Spatial_Averaging_Radius");
52 init_weights_analytical();
55 PauliBlocker::~PauliBlocker() {}
60 const ParticleList &disregard)
const {
66 for (
const auto &part : particles) {
68 if (part.pdgcode() != pdg) {
72 const double pdist_sqr = (part.momentum().threevec() -
p).sqr();
73 if (pdist_sqr > rp_ * rp_) {
76 const double rdist_sqr = (part.position().threevec() - r).sqr();
78 if (rdist_sqr >= (rr_ + rc_) * (rr_ + rc_)) {
82 bool to_disregard =
false;
83 for (
const auto &disregard_part : disregard) {
84 if (part.id() == disregard_part.id()) {
92 const double i_real = std::sqrt(rdist_sqr) / (rr_ + rc_) * weights_.size();
93 const size_t i = std::floor(i_real);
94 const double rest = i_real - i;
95 if (
likely(i + 1 < weights_.size())) {
96 f += weights_[i] * rest + weights_[i + 1] * (1. - rest);
102 void PauliBlocker::init_weights_analytical() {
103 const auto &log = logger<LogArea::PauliBlocking>();
105 const double pi = M_PI;
106 const double sqrt2 = std::sqrt(2.);
107 const double sqrt_2pi = std::sqrt(2. * pi);
109 const double phase_volume =
110 2 * (4. / 3. * pi * rr_ * rr_ * rr_) * (4. / 3. * pi * rp_ * rp_ * rp_) /
114 std::erf(rc_ / sqrt2 / sig_) -
115 rc_ * 2 / sqrt_2pi / sig_ * std::exp(-0.5 * rc_ * rc_ / sig_ / sig_);
119 const double d_pos = (rr_ + rc_) / static_cast<double>(weights_.size());
121 for (
size_t k = 0; k < weights_.size(); k++) {
123 const double rj = d_pos * k;
126 const double A = rr_ / sqrt2 / sig_;
127 integral = sqrt_2pi * sig_ * std::erf(A) - 2 * rr_ * std::exp(-A * A);
128 integral *= sig_ * sig_;
129 }
else if (rc_ > rj + rr_) {
130 const double A = (rj + rr_) / sqrt2 / sig_;
131 const double B = (rj - rr_) / sqrt2 / sig_;
132 integral = sig_ / rj * (std::exp(-A * A) - std::exp(-B * B)) +
133 0.5 * sqrt_2pi * (std::erf(A) - std::erf(B));
134 integral *= sig_ * sig_ * sig_;
136 const double A = rc_ / sqrt2 / sig_;
137 const double B = (rj - rr_) / sqrt2 / sig_;
138 const double C = (rc_ - rj) * (rc_ - rj) - rr_ * rr_ + 2 * sig_ * sig_;
140 (0.5 * std::exp(-A * A) * C - sig_ * sig_ * std::exp(-B * B)) / rj +
141 0.5 * sqrt_2pi * sig_ * (std::erf(A) - std::erf(B));
142 integral *= sig_ * sig_;
144 integral *= 2 * pi / std::pow(2 * pi * sig_ * sig_, 1.5);
145 weights_[k] = integral / norm / phase_volume;
146 log.debug(
"Analytical weights[", k,
"] = ", weights_[k]);
The ThreeVector class represents a physical three-vector with the components .
constexpr double really_small
Numerical error tolerance.
#define likely(x)
Tell the branch predictor that this expression is likely true.
Collection of useful constants that are known at compile time.
PauliBlocker(Configuration conf, const ExperimentParameters ¶meters)
PauliBlocker constructor.
constexpr double hbarc
GeV <-> fm conversion factor.
PdgCode stores a Particle Data Group Particle Numbering Scheme particle type number.
The Particles class abstracts the storage and manipulation of particles.