From 7fbb4a9bb549ce0290c6293be085a3fa29f290d0 Mon Sep 17 00:00:00 2001
From: Jean-Marco Alameddine <jean-marco.alameddine@udo.edu>
Date: Mon, 22 May 2023 17:05:20 +0200
Subject: [PATCH] start refactoring variable names in InteractionModel.inl
---
.../modules/proposal/InteractionModel.inl | 97 ++++++++++---------
1 file changed, 52 insertions(+), 45 deletions(-)
diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index 7a34c1365..3f531f2bd 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -64,19 +64,19 @@ namespace corsika::proposal {
if (canInteract(projectileId)) {
// get or build corresponding calculators
- auto c = getCalculator(projectile, calc_);
+ auto calculator = getCalculator(projectile, calc_);
// get the rates of the interaction types for every component.
std::uniform_real_distribution<double> distr(0., 1.);
- // sample a interaction-type, loss and component
+ // sample an interaction-type, loss and component
auto rates =
- std::get<eINTERACTION>(c->second)->Rates(projectile.getEnergy() / 1_MeV);
- auto [type, target_hash, v] = std::get<eINTERACTION>(c->second)->SampleLoss(
+ std::get<eINTERACTION>(calculator->second)->Rates(projectile.getEnergy() / 1_MeV);
+ auto [interaction_type, interaction_target_hash, v] = std::get<eINTERACTION>(calculator->second)->SampleLoss(
projectile.getEnergy() / 1_MeV, rates, distr(RNG_));
// TODO: This should become obsolete as soon #482 is fixed
- if (type == PROPOSAL::InteractionType::Undefined) {
+ if (interaction_type == PROPOSAL::InteractionType::Undefined) {
CORSIKA_LOG_WARN(
"PROPOSAL: No particle interaction possible. "
"Put initial particle back on stack.");
@@ -86,57 +86,64 @@ namespace corsika::proposal {
return ProcessReturn::Ok;
}
- // Read how much random numbers are required to calculate the secondaries.
- // Calculate the secondaries and deploy them on the corsika stack.
- auto rnd = std::vector<double>(
- std::get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type));
- for (auto& it : rnd) it = distr(RNG_);
- Point const& place = projectile.getPosition();
- CoordinateSystemPtr const& labCS = place.getCoordinateSystem();
-
- auto point = PROPOSAL::Cartesian3D(
- place.getX(labCS) / 1_cm, place.getY(labCS) / 1_cm, place.getZ(labCS) / 1_cm);
- auto projectile_dir = projectile.getDirection();
- auto d = projectile_dir.getComponents(labCS);
- auto direction = PROPOSAL::Cartesian3D(d.getX().magnitude(), d.getY().magnitude(),
- d.getZ().magnitude());
-
- auto loss = PROPOSAL::StochasticLoss(
- static_cast<int>(type), v * projectile.getEnergy() / 1_MeV, point, direction,
+ // Read how much random numbers are required by proposal to sample the secondary particles
+ auto random_numbers_secondaries = std::vector<double>(
+ std::get<eSECONDARIES>(calculator->second)->RequiredRandomNumbers(interaction_type));
+ for (auto& it : random_numbers_secondaries) it = distr(RNG_);
+
+ // calculate position of interaction
+ Point const& interaction_position = projectile.getPosition();
+ CoordinateSystemPtr const& labCS = interaction_position.getCoordinateSystem();
+ auto interaction_position_pp = PROPOSAL::Cartesian3D(
+ interaction_position.getX(labCS) / 1_cm,
+ interaction_position.getY(labCS) / 1_cm,
+ interaction_position.getZ(labCS) / 1_cm);
+
+ // calculate direction of initial projection
+ auto projectile_direction = projectile.getDirection();
+ auto projectile_direction_components = projectile_direction.getComponents(labCS);
+ auto projectile_direction_pp = PROPOSAL::Cartesian3D(
+ projectile_direction_components.getX().magnitude(),
+ projectile_direction_components.getY().magnitude(),
+ projectile_direction_components.getZ().magnitude());
+
+ // calculate secondary particles from interaction
+ auto energy_loss_interaction = PROPOSAL::StochasticLoss(
+ static_cast<int>(interaction_type), v * projectile.getEnergy() / 1_MeV, interaction_position_pp, projectile_direction_pp,
projectile.getTime() / 1_s, 0., projectile.getEnergy() / 1_MeV);
- PROPOSAL::Component target;
- if (type != PROPOSAL::InteractionType::Ioniz)
- target = PROPOSAL::Component::GetComponentForHash(target_hash);
-
- auto sec =
- std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
- for (auto& s : sec) {
- auto E = s.energy * 1_MeV;
- auto vecProposal = s.direction;
- auto dir = DirectionVector(
- labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
-
- if (static_cast<PROPOSAL::ParticleType>(s.type) ==
+ PROPOSAL::Component interaction_target;
+ if (interaction_type != PROPOSAL::InteractionType::Ioniz)
+ interaction_target = PROPOSAL::Component::GetComponentForHash(interaction_target_hash);
+ auto secondary_particles =
+ std::get<eSECONDARIES>(calculator->second)->CalculateSecondaries(energy_loss_interaction, interaction_target, random_numbers_secondaries);
+
+ // deploy sampled secondary particles on stack
+ for (auto& secondary_particle : secondary_particles) {
+ auto energy = secondary_particle.energy * 1_MeV;
+ auto direction_pp = secondary_particle.direction;
+ auto direction = DirectionVector(
+ labCS, {direction_pp.GetX(), direction_pp.GetY(), direction_pp.GetZ()});
+
+ if (static_cast<PROPOSAL::ParticleType>(secondary_particle.type) ==
PROPOSAL::ParticleType::Hadron) {
- FourMomentum const photonP4(E, E * dir);
+ FourMomentum const photonP4(energy, energy * direction);
// for PROPOSAL media target.GetAtomicNum() returns the atomic number in units
// of atomic mass, so Nitrogen is 14.0067. When using media in CORSIKA the same
// field is filled with the number of nucleons (ie. 14 for Nitrogen) To be sure
// we explicitly cast to int
- auto const A = int(target.GetAtomicNum());
- auto const Z = int(target.GetNucCharge());
+ auto const A = int(interaction_target.GetAtomicNum());
+ auto const Z = int(interaction_target.GetNucCharge());
Code const targetId = get_nucleus_code(A, Z);
CORSIKA_LOGGER_DEBUG(
logger_,
"photo-hadronic interaction of projectile={} with target={}! Energy={} GeV",
- projectileId, targetId, E / 1_GeV);
+ projectileId, targetId, energy / 1_GeV);
this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
} else {
- auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
+ auto sec_code = convert_from_PDG(static_cast<PDGCode>(secondary_particle.type));
// use mass provided by PROPOSAL to ensure correct conversion to kinetic energy
- auto massProposal =
- PROPOSAL::ParticleDef::GetParticleDefForType(s.type).mass * 1_MeV;
- view.addSecondary(std::make_tuple(sec_code, E - massProposal, dir));
+ auto mass = PROPOSAL::ParticleDef::GetParticleDefForType(secondary_particle.type).mass * 1_MeV;
+ view.addSecondary(std::make_tuple(sec_code, energy - mass, direction));
}
}
}
@@ -161,8 +168,8 @@ namespace corsika::proposal {
// ==============================================
if (canInteract(projectileId)) {
- auto c = getCalculator(projectile, calc_);
- return meanMass / (std::get<eINTERACTION>(c->second)->MeanFreePath(
+ auto calculator = getCalculator(projectile, calc_);
+ return meanMass / (std::get<eINTERACTION>(calculator->second)->MeanFreePath(
projectileP4.getTimeLikeComponent() / 1_MeV) *
1_g / (1_cm * 1_cm));
}
--
GitLab