#include <pauliblocking.h>

A class that stores parameters needed for Pauli blocking, tabulates necessary integrals and computes phase-space density.

Pauli blocking is the way to go from the classical Boltzmann equation to the quantum one, effectively reducing cross-sections by \(1 - f(r,p) \) factors, where \(f(r,p)\) is a phase-space density at coordinate and momentum of a final-state fermion. Effective reduction of cross-section is done via random rejection of reaction with probability \( 1 - f\). More details can be found in Gaitanos:2010fd [22], section III B. Our implementation mainly follows this article (and therefore GiBUU, see http://gibuu.hepforge.org).

Definition at line 38 of file pauliblocking.h.

Public Member Functions
	PauliBlocker (Configuration conf, const ExperimentParameters &parameters)
	PauliBlocker constructor. More...

	~PauliBlocker ()
	Destructor. More...

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). More...

Private Member Functions
void	init_weights ()
	Tabulate integrals for weights. More...

void	init_weights_analytical ()
	Analytical calculation of weights. More...

Private Attributes
double	sig_
	Standard deviation of the gaussian used for smearing. More...

double	rc_
	Radius, after which gaussians (used for averaging) are cut, fm. More...

double	rr_
	Radius of averaging in coordinate space, fm. More...

double	rp_
	Radius of averaging in momentum space, GeV. More...

int	ntest_
	Testparticles number. More...

int	n_ensembles_
	Number of ensembles. More...

std::array< double, 30 >	weights_
	Weights: tabulated results of numerical integration. More...

Constructor & Destructor Documentation

◆ PauliBlocker()

smash::PauliBlocker::PauliBlocker	(	Configuration	conf,
		const ExperimentParameters &	parameters
	)

PauliBlocker constructor.

Gets parameters from configuration and experiment. Tabulates necessary integrals.

Parameters

[in]	conf	Configurations from config.yaml.
[in]	parameters	Parameters given by Experiment.

Returns: The constructed object.

Definition at line 18 of file pauliblocking.cc.

     : sig_(param.gaussian_sigma),
       rc_(conf.take({"Gaussian_Cutoff"}, 2.2)),
       rr_(conf.take({"Spatial_Averaging_Radius"}, 1.86)),
       rp_(conf.take({"Momentum_Averaging_Radius"}, 0.08)),
       ntest_(param.testparticles),
       n_ensembles_(param.n_ensembles) {
   if (ntest_ * n_ensembles_ < 20) {
     logg[LPauliBlocking].warn(
         "Phase-space density calculation in Pauli blocking will not work "
         "reasonably.\nEither use testparticles or ensembles, or both.\n"
         "The recommended testparticles * ensembles is at least 20.");
   }
  
   if (rc_ < rr_ || rr_ < 0.0 || rp_ < 0) {
     logg[LPauliBlocking].error(
         "Please choose reasonable parameters for Pauli blocking:\n"
         "All radii have to be positive and Gaussian_Cutoff should\n"
         "be larger than Spatial_Averaging_Radius");
   }
  
   init_weights_analytical();
 }

◆ ~PauliBlocker()

smash::PauliBlocker::~PauliBlocker ( )

Destructor.

Definition at line 43 of file pauliblocking.cc.

43 {}

Member Function Documentation

◆ phasespace_dens()

double smash::PauliBlocker::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).

Parameters

[in]	r	Position vector of the particle.
[in]	p	Momentum vector of the particle.
[in]	ensembles	Current list of particles in all ensembles.
[in]	pdg	PDG number of species for which density to be calculated.
[in]	disregard	Do not count particles that should be disregarded. This is intended to avoid counting incoming particles when the phase-space density for outgoing ones is estimated.

Returns: Phase-space density

Definition at line 45 of file pauliblocking.cc.

                                                                           {
   double f = 0.0;
  
   /* TODO(oliiny): looping over all particles is inefficient,
    * I need only particles within rp_ radius in momentum and
    * within rr_+rc_ in coordinate space. Some search algorithm might help. */
   for (const Particles &particles : ensembles) {
     for (const ParticleData &part : particles) {
       // Only consider identical particles
       if (part.pdgcode() != pdg) {
         continue;
       }
       // Only consider momenta in sphere of radius rp_ with center at p
       const double pdist_sqr = (part.momentum().threevec() - p).sqr();
       if (pdist_sqr > rp_ * rp_) {
         continue;
       }
       const double rdist_sqr = (part.position().threevec() - r).sqr();
       // Only consider coordinates in sphere of radius rr_+rc_ with center at r
       if (rdist_sqr >= (rr_ + rc_) * (rr_ + rc_)) {
         continue;
       }
       // Do not count particles that should be disregarded.
       bool to_disregard = false;
       for (const auto &disregard_part : disregard) {
         if (part.id() == disregard_part.id()) {
           to_disregard = true;
         }
       }
       if (to_disregard) {
         continue;
       }
       // 1st order interpolation using tabulated values
       const double i_real =
           std::sqrt(rdist_sqr) / (rr_ + rc_) * weights_.size();
       const size_t i = std::floor(i_real);
       const double rest = i_real - i;
       if (likely(i + 1 < weights_.size())) {
         f += weights_[i] * rest + weights_[i + 1] * (1. - rest);
       }
     }  // loop over particles in one ensemble
   }    // loop over ensembles
   return f / ntest_ / n_ensembles_;
 }

◆ init_weights()

void smash::PauliBlocker::init_weights ( )

private

Tabulate integrals for weights.

◆ init_weights_analytical()

void smash::PauliBlocker::init_weights_analytical ( )

private

Analytical calculation of weights.

Definition at line 93 of file pauliblocking.cc.

                                            {
   const double pi = M_PI;
   const double sqrt2 = std::sqrt(2.);
   const double sqrt_2pi = std::sqrt(2. * pi);
   // Volume of the phase-space area; Factor 2 stands for spin.
   const double phase_volume =
       2 * (4. / 3. * pi * rr_ * rr_ * rr_) * (4. / 3. * pi * rp_ * rp_ * rp_) /
       ((2 * pi * hbarc) * (2 * pi * hbarc) * (2 * pi * hbarc));
   // Analytical expression for integral in denominator
   const double norm =
       std::erf(rc_ / sqrt2 / sig_) -
       rc_ * 2 / sqrt_2pi / sig_ * std::exp(-0.5 * rc_ * rc_ / sig_ / sig_);
  
   double integral;
   // Step of the table for tabulated integral
   const double d_pos = (rr_ + rc_) / static_cast<double>(weights_.size());
  
   for (size_t k = 0; k < weights_.size(); k++) {
     // rdist = 0 ... rc_ (gauss cut) + rr_ (position cut)
     const double rj = d_pos * k;
     if (rj < really_small) {
       // Assuming rc_ > rr_
       const double A = rr_ / sqrt2 / sig_;
       integral = sqrt_2pi * sig_ * std::erf(A) - 2 * rr_ * std::exp(-A * A);
       integral *= sig_ * sig_;
     } else if (rc_ > rj + rr_) {
       const double A = (rj + rr_) / sqrt2 / sig_;
       const double B = (rj - rr_) / sqrt2 / sig_;
       integral = sig_ / rj * (std::exp(-A * A) - std::exp(-B * B)) +
                  0.5 * sqrt_2pi * (std::erf(A) - std::erf(B));
       integral *= sig_ * sig_ * sig_;
     } else {
       const double A = rc_ / sqrt2 / sig_;
       const double B = (rj - rr_) / sqrt2 / sig_;
       const double C = (rc_ - rj) * (rc_ - rj) - rr_ * rr_ + 2 * sig_ * sig_;
       integral =
           (0.5 * std::exp(-A * A) * C - sig_ * sig_ * std::exp(-B * B)) / rj +
           0.5 * sqrt_2pi * sig_ * (std::erf(A) - std::erf(B));
       integral *= sig_ * sig_;
     }
     integral *= 2 * pi / std::pow(2 * pi * sig_ * sig_, 1.5);
     weights_[k] = integral / norm / phase_volume;
     logg[LPauliBlocking].debug("Analytical weights[", k, "] = ", weights_[k]);
   }
 }

Member Data Documentation

◆ sig_

double smash::PauliBlocker::sig_

private

Standard deviation of the gaussian used for smearing.

Definition at line 79 of file pauliblocking.h.

◆ rc_

double smash::PauliBlocker::rc_

private

Radius, after which gaussians (used for averaging) are cut, fm.

Definition at line 82 of file pauliblocking.h.

◆ rr_

double smash::PauliBlocker::rr_

private

Radius of averaging in coordinate space, fm.

Definition at line 85 of file pauliblocking.h.

◆ rp_

double smash::PauliBlocker::rp_

private

Radius of averaging in momentum space, GeV.

Definition at line 88 of file pauliblocking.h.

◆ ntest_

int smash::PauliBlocker::ntest_

private

Testparticles number.

Definition at line 91 of file pauliblocking.h.

◆ n_ensembles_

int smash::PauliBlocker::n_ensembles_

private

Number of ensembles.

Definition at line 94 of file pauliblocking.h.

◆ weights_

std::array<double, 30> smash::PauliBlocker::weights_

private

Weights: tabulated results of numerical integration.

Definition at line 97 of file pauliblocking.h.

The documentation for this class was generated from the following files:

src/include/smash/pauliblocking.h
src/pauliblocking.cc

Public Member Functions

Private Member Functions

Private Attributes

Constructor & Destructor Documentation

◆ PauliBlocker()

◆ ~PauliBlocker()

Member Function Documentation

◆ phasespace_dens()

◆ init_weights()

◆ init_weights_analytical()

Member Data Documentation

◆ sig_

◆ rc_

◆ rr_

◆ rp_

◆ ntest_

◆ n_ensembles_

◆ weights_