Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add LPM effect for photopair, fixes for BremsLPM #336

Merged
merged 6 commits into from
Jan 5, 2023
Merged

Add LPM effect for photopair, fixes for BremsLPM #336

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
453f145
Change definition of s1 such that it is applicable to all lepton masses
asandrock Jan 4, 2023
afadcfe
Add LPM effect for PhotoPairProduction
asandrock Jan 4, 2023
87f12c3
Fix ifstream in UnitTests to correctly read in ecut==INF using std::stod
Jean1995 Jan 4, 2023
b9be763
Merge branch 'fix_ifstream_for_infinite_ecuts' into add_photopairlpm_…
Jean1995 Jan 4, 2023
0216967
Add LPM UnitTests to PhotoPair
Jean1995 Jan 4, 2023
ecf9628
Add density_corretion tests to PhotoPair_TEST
Jean1995 Jan 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/config_docu.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -446,6 +446,7 @@ For there parametrizations, the parametrization of the shadowing effect can be c
| ----------------- | ------- | ------- | -------------------------- |
| `parametrization` | String | `-` | Parametrization to be used |
| `multiplier` | Number | `1.0` | Multiplier to scale the cross section with a coefficient |
| `lpm` | Boolean | `true` | Enabling or disabling the reduction of the cross section at high energies by the Landau-Pomeranchuk-Migdal effect. |

Creation of an electron-positron pair by an ingoing photon. Available `photopair` parametrizations are:

Expand Down
3 changes: 2 additions & 1 deletion src/PROPOSAL/PROPOSAL/crosssection/Factories/PhotoPairProductionFactory.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,5 +12,6 @@ namespace PROPOSAL {

namespace PROPOSAL {
cross_ptr make_photopairproduction(const ParticleDef&, const Medium&, bool,
const nlohmann::json&);
const nlohmann::json&,
double density_correction = 1.0);
}
2 changes: 1 addition & 1 deletion src/PROPOSAL/PROPOSAL/crosssection/ParticleDefaultCrossSectionList.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ template<>
template <typename CrossVec, typename P, typename M>
void DefaultCrossSections<GammaDef>::Append(CrossVec& cross_vec, P p, M m,std::shared_ptr<const EnergyCutSettings> cut, bool interpolate)
{
auto photopair = std::make_tuple(crosssection::PhotoPairKochMotz{}, p, m, nullptr, interpolate);
auto photopair = std::make_tuple(crosssection::PhotoPairKochMotz{false}, p, m, nullptr, interpolate);
auto compton = std::make_tuple(crosssection::ComptonKleinNishina{}, p, m, cut, interpolate);
auto photoproduction = std::make_tuple(crosssection::PhotoproductionRhode{}, p, m, nullptr, interpolate);
auto photoeffect = std::make_tuple(crosssection::PhotoeffectSauter{}, p, m, nullptr, interpolate);
Expand Down
31 changes: 28 additions & 3 deletions src/PROPOSAL/PROPOSAL/crosssection/parametrization/PhotoPairProduction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,9 +34,14 @@

namespace PROPOSAL {
namespace crosssection {
class PhotoPairLPM;

class PhotoPairProduction : public Parametrization<Component> {
protected:
std::shared_ptr<PhotoPairLPM> lpm_;
double density_correction_; // correction to standard medium density for LPM
public:
PhotoPairProduction() = default;
PhotoPairProduction();
virtual ~PhotoPairProduction() = default;

double GetLowerEnergyLim(const ParticleDef&) const noexcept final;
Expand All @@ -49,7 +54,9 @@ namespace crosssection {
};

struct PhotoPairTsai : public PhotoPairProduction {
PhotoPairTsai() { hash_combine(hash, std::string("tsai")); }
PhotoPairTsai(bool lpm = false);
PhotoPairTsai(bool lpm, const ParticleDef&, const Medium&,
double density_correction = 1.0);
using base_param_t = PhotoPairProduction;
std::unique_ptr<Parametrization<Component>> clone() const final;

Expand All @@ -58,7 +65,9 @@ namespace crosssection {
};

struct PhotoPairKochMotz : public PhotoPairProduction {
PhotoPairKochMotz();
PhotoPairKochMotz(bool lpm = false);
PhotoPairKochMotz(bool lpm, const ParticleDef&, const Medium&,
double density_correction = 1.0);
using base_param_t = PhotoPairProduction;
std::unique_ptr<Parametrization<Component>> clone() const final;

Expand Down Expand Up @@ -86,5 +95,21 @@ namespace crosssection {
template <> struct ParametrizationId<PhotoPairKochMotz> {
static constexpr size_t value = 1000000013;
};

// LPM effect object
class PhotoPairLPM {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we really need this constructor with mol_density, mass_density, sum_charge, eLPM? Some variables are not really intuitive.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

eLPM is a computationally expensive calculation that needs to be executed only once, so we should save it.
mol_density_ and sum_charge_ are medium variables. Since those are the only two variables we need for the calcluations, we save them to avoid saving the whole medium in the PhotoPairLPM object. So I believe these three values are useful to save. Their names are just adopted from what they are called in the Medium class (although I agree that the variable names are not always well-chosen, but for the sake of consistency, I would use the same names here).

mass_density_ is only used in the constructor, so I would agree that we don't need to save them as a private member of the object. But in my opinion, this is the only change I would make to the constructor.

size_t hash;
double mol_density_;
double mass_density_;
double sum_charge_;
double eLpm_;

public:
PhotoPairLPM(const ParticleDef&, const Medium&, const PhotoPairProduction&);
double suppression_factor(double energy, double x, const Component&,
double density_correction = 1.0) const;
size_t GetHash() const noexcept { return hash; }
};

} // namespace crosssection
} // namespace PROPOSAL
3 changes: 2 additions & 1 deletion src/PROPOSAL/detail/PROPOSAL/Propagator.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -485,7 +485,8 @@ Propagator::CreateCrossSectionList(const ParticleDef& p_def,
p_def, medium, config["photoproduction"]));
if (config.contains("photopair"))
cross.emplace_back(make_photopairproduction(
p_def, medium, interpolate, config["photopair"]));
p_def, medium, interpolate, config["photopair"],
density_correction));
if (config.contains("weak"))
cross.emplace_back(
make_weakinteraction(p_def, medium, interpolate, config["weak"]));
Expand Down
15 changes: 10 additions & 5 deletions src/PROPOSAL/detail/PROPOSAL/crosssection/Factories/PhotoPairProductionFactory.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,15 @@
#include <nlohmann/json.hpp>

using namespace PROPOSAL;
using photopair_func_ptr = cross_ptr (*)(const ParticleDef&, const Medium&, bool);
using photopair_func_ptr = cross_ptr (*)(const ParticleDef&, const Medium&,
bool, bool, double);

template <typename Param>
cross_ptr create_photopairproduction(const ParticleDef& p_def, const Medium& medium, bool interpol)
cross_ptr create_photopairproduction(
const ParticleDef& p_def, const Medium& medium, bool interpol,
bool lpm, double density_correction)
{
auto param = Param();
auto param = Param(lpm, p_def, medium, density_correction);
return make_crosssection(param, p_def, medium, nullptr, interpol);
}

Expand All @@ -24,16 +27,18 @@ std::map<std::string, photopair_func_ptr, Helper::case_insensitive_comp> photopa

namespace PROPOSAL {
cross_ptr make_photopairproduction(const ParticleDef &p_def, const Medium &medium,
bool interpol, const nlohmann::json &config) {
bool interpol, const nlohmann::json &config,
double density_correction) {
if (!config.contains("parametrization"))
throw std::logic_error("No parametrization passed for photopairproduction");
std::string param_name = config["parametrization"];
bool lpm = config.value("lpm", true);

auto it = photopair_map.find(param_name);
if (it == photopair_map.end())
throw std::logic_error("Unknown parametrization for photopairproduction");

auto cross = it->second(p_def, medium, interpol);
auto cross = it->second(p_def, medium, interpol, lpm, density_correction);

double multiplier = config.value("multiplier", 1.0);
if (multiplier != 1.0)
Expand Down
8 changes: 5 additions & 3 deletions src/PROPOSAL/detail/PROPOSAL/crosssection/parametrization/Bremsstrahlung.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -595,7 +595,7 @@ crosssection::BremsLPM::BremsLPM(const ParticleDef& p_def, const Medium& medium,
}
sum = sum * mass_density_;
eLpm_ = ALPHA * mass_;
eLpm_ *= eLpm_ / (4. * PI * ME * RE * sum);
eLpm_ *= 2* eLpm_ / (PI * ME * RE * sum);
}

double crosssection::BremsLPM::suppression_factor(
Expand All @@ -610,8 +610,10 @@ double crosssection::BremsLPM::suppression_factor(
double Z3 = std::pow(comp.GetNucCharge(), -1. / 3);

double Dn = 1.54 * std::pow(comp.GetAtomicNum(), 0.27);
s1 = ME * Dn / (mass_ * Z3 * comp.GetLogConstant());
s1 = s1 * s1 * SQRT2; // TODO: is SQRT2 correct here, this factor is not in the paper?
double aux_nuc = mass_ / MMU * Dn;
s1 = ME / (mass_ * Z3 * comp.GetLogConstant());
s1 = s1 * s1 * (1 + aux_nuc * aux_nuc) * SQRT2;
// SQRT2 appears together with s1 everywhere in Stanev et al.

// Calc xi(s') from Stanev, Vankow, Streitmatter, Ellsworth, Bowen
// Phys. Rev. D 25 (1982), 1291
Expand Down