IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 1f9e5d4c authored by Maximilian Sackel's avatar Maximilian Sackel Committed by Ralf Ulrich
Browse files

stochastic losses are now possible for electrons and gammas

parent 7e3df7d7
No related branches found
No related tags found
1 merge request!245Include proposal process rebase
Hide whitespace changes
Inline Side-by-side
Processes/Proposal/Interaction.cc
+ 22
− 36
View file @ 1f9e5d4c
......@@ -30,17 +30,15 @@ namespace corsika::process::proposal {
};
bool Interaction::CanInteract(particles::Code pcode) const noexcept {
auto search = particles.find(pcode);
if (search != particles.end()) return true;
if (particles.find(pcode) != particles.end()) return true;
return false;
}
template <>
Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut const& _cut)
: cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
false)) {
std::cout << "corsika set cut to " << _cut.GetECut() / 1_GeV << " [GeV]" << std::endl;
std::cout << "proposal set ecut to " << cut->GetEcut() << " [MeV]" << std::endl;
false))
, fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
auto all_compositions = std::vector<const NuclearComposition*>();
_env.GetUniverse()->walk([&](auto& vtn) {
if (vtn.HasModelProperties()) {
......@@ -62,8 +60,6 @@ namespace corsika::process::proposal {
}
}
void Interaction::Init() {}
template <>
corsika::process::EProcessReturn Interaction::DoInteraction(
setup::StackView::StackIterator& vP) {
......@@ -71,39 +67,34 @@ namespace corsika::process::proposal {
auto calc = GetCalculator(vP); // [CrossSections]
std::uniform_real_distribution<double> distr(0., 1.);
auto [type, comp_ptr, v] =
std::get<INTERACTION>(calc->second)
get<INTERACTION>(calc->second)
->TypeInteraction(vP.GetEnergy() / 1_MeV, distr(fRNG));
std::cout << "InteractionType: "<< static_cast<int>(type) << std::endl;
auto rnd = std::vector<double>();
for (size_t i = 0;
i < std::get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type); ++i)
rnd.push_back(distr(fRNG));
double primary_energy = vP.GetEnergy() / 1_MeV;
auto rnd =
vector<double>(get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type));
for (auto& it : rnd) it = distr(fRNG);
auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
vP.GetPosition().GetY() / 1_cm,
vP.GetPosition().GetZ() / 1_cm);
auto d = vP.GetDirection().GetComponents();
auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(),
d.GetY().magnitude(),
auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
d.GetZ().magnitude());
auto loss =
make_tuple(static_cast<int>(type), point, direction, v * primary_energy, 0.);
auto sec = std::get<SECONDARIES>(calc->second)
->CalculateSecondaries(primary_energy, loss, *comp_ptr, rnd);
auto loss = make_tuple(static_cast<int>(type), point, direction,
v * vP.GetEnergy() / 1_MeV, 0.);
auto sec = get<SECONDARIES>(calc->second)
->CalculateSecondaries(vP.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd);
for (auto& s : sec) {
auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
auto vec = corsika::geometry::QuantityVector(
std::get<PROPOSAL::Loss::DIRECTION>(s).GetX() * energy,
std::get<PROPOSAL::Loss::DIRECTION>(s).GetY() * energy,
std::get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * energy);
auto momentum = corsika::stack::MomentumVector(
get<PROPOSAL::Loss::DIRECTION>(s).GetX() * E,
get<PROPOSAL::Loss::DIRECTION>(s).GetY() * E,
get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * E);
auto p = corsika::stack::MomentumVector(
corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(),
vec);
particles::Code sec_code = corsika::particles::ConvertFromPDG(
auto sec_code = corsika::particles::ConvertFromPDG(
static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
vP.AddSecondary(
make_tuple(sec_code, energy, momentum, vP.GetPosition(), vP.GetTime()));
vP.AddSecondary(make_tuple(sec_code, E, p, vP.GetPosition(), vP.GetTime()));
}
}
return process::EProcessReturn::eOk;
......@@ -113,14 +104,9 @@ namespace corsika::process::proposal {
corsika::units::si::GrammageType Interaction::GetInteractionLength(
setup::Stack::StackIterator const& vP) {
if (CanInteract(vP.GetPID())) {
auto calc = GetCalculator(vP); // [CrossSections]
std::uniform_real_distribution<double> distr(0., 1.);
auto rnd = distr(fRNG);
auto energy = get<INTERACTION>(calc->second)
->EnergyInteraction(vP.GetEnergy() / 1_MeV, rnd);
return get<DISPLACEMENT>(calc->second)
->SolveTrackIntegral(vP.GetEnergy() / 1_MeV, energy) *
1_g / 1_cm / 1_cm;
auto calc = GetCalculator(vP);
return get<INTERACTION>(calc->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
(1_cm * 1_cm);
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
......
Processes/Proposal/Interaction.h
+ 43
− 104
View file @ 1f9e5d4c
......@@ -23,6 +23,8 @@
using namespace corsika::environment;
using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
using std::make_pair;
using std::make_tuple;
namespace corsika::process::proposal {
......@@ -33,130 +35,67 @@ namespace corsika::process::proposal {
static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("proposal");
corsika::random::RNG& fRNG;
bool CanInteract(particles::Code pcode) const noexcept;
using calculator_t =
tuple<unique_ptr<PROPOSAL::SecondariesCalculator>, unique_ptr<PROPOSAL::Interaction>,
unique_ptr<PROPOSAL::Displacement>>;
using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
unique_ptr<PROPOSAL::Interaction>>;
using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
struct interaction_hash {
size_t operator()(const std::pair<const NuclearComposition*, particles::Code>& p) const
{
auto hash1 = std::hash<const NuclearComposition*>{}(p.first);
auto hash2 = std::hash<particles::Code>{}(p.second);
return hash1 ^ hash2;
}
size_t operator()(const calc_key_t& p) const {
return std::hash<const NuclearComposition*>{}(p.first) ^
std::hash<particles::Code>{}(p.second);
}
};
std::unordered_map<std::pair<const NuclearComposition*, particles::Code>, calculator_t, interaction_hash> calculators;
std::unordered_map<calc_key_t, calculator_t, interaction_hash> calculators;
template <typename Particle>
auto BuildCalculator(particles::Code corsika_code, Particle p_def,
NuclearComposition const& comp) {
auto cross = GetStdCrossSections(p_def, media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{make_pair(&comp, corsika_code),
make_tuple(PROPOSAL::make_secondaries(inter_types, p_def, media.at(&comp)),
PROPOSAL::make_interaction(cross, true))});
return insert_it;
}
auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
auto medium = media.at(&comp);
if (corsika_code == particles::Code::Gamma) {
std::cout << "Build gamma tables" << std::endl;
auto cross =
GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
inter_types, PROPOSAL::GammaDef(), media.at(&comp)),
PROPOSAL::make_interaction(cross, true),
PROPOSAL::make_displacement(cross, true))});
return insert_it;
}
if (corsika_code == particles::Code::Electron) {
std::cout << "Build electron tables" << std::endl;
auto cross =
GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
inter_types, PROPOSAL::EMinusDef(), media.at(&comp)),
PROPOSAL::make_interaction(cross, true),
PROPOSAL::make_displacement(cross, true))});
return insert_it;
}
if (corsika_code == particles::Code::Positron) {
std::cout << "Build positron tables" << std::endl;
auto cross =
GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
inter_types, PROPOSAL::EPlusDef(), media.at(&comp)),
PROPOSAL::make_interaction(cross, true),
PROPOSAL::make_displacement(cross, true))});
return insert_it;
}
if (corsika_code == particles::Code::MuMinus) {
std::cout << "Build MuMinus tables" << std::endl;
auto cross =
GetStdCrossSections(PROPOSAL::MuMinusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
inter_types, PROPOSAL::MuMinusDef(), media.at(&comp)),
PROPOSAL::make_interaction(cross, true),
PROPOSAL::make_displacement(cross, true))});
return insert_it;
}
if (corsika_code == particles::Code::MuPlus) {
std::cout << "Build MuPlus tables" << std::endl;
auto cross =
GetStdCrossSections(PROPOSAL::MuPlusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
inter_types, PROPOSAL::MuPlusDef(), media.at(&comp)),
PROPOSAL::make_interaction(cross, true),
PROPOSAL::make_displacement(cross, true))});
return insert_it;
}
if (corsika_code == particles::Code::TauMinus) {
std::cout << "Build TauMinus tables" << std::endl;
auto cross =
GetStdCrossSections(PROPOSAL::TauMinusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
inter_types, PROPOSAL::TauMinusDef(), media.at(&comp)),
PROPOSAL::make_interaction(cross, true),
PROPOSAL::make_displacement(cross, true))});
return insert_it;
}
if (corsika_code == particles::Code::TauPlus) {
std::cout << "Build TauPlus tables" << std::endl;
auto cross =
GetStdCrossSections(PROPOSAL::TauPlusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
{std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
inter_types, PROPOSAL::TauPlusDef(), media.at(&comp)),
PROPOSAL::make_interaction(cross, true),
PROPOSAL::make_displacement(cross, true))});
return insert_it;
}
if (corsika_code == particles::Code::Gamma)
return BuildCalculator(particles::Code::Gamma, PROPOSAL::GammaDef(), comp);
if (corsika_code == particles::Code::Electron)
return BuildCalculator(particles::Code::Electron, PROPOSAL::EMinusDef(), comp);
if (corsika_code == particles::Code::Positron)
return BuildCalculator(particles::Code::Positron, PROPOSAL::EPlusDef(), comp);
if (corsika_code == particles::Code::MuMinus)
return BuildCalculator(particles::Code::MuMinus, PROPOSAL::MuMinusDef(), comp);
if (corsika_code == particles::Code::MuPlus)
return BuildCalculator(particles::Code::MuPlus, PROPOSAL::MuPlusDef(), comp);
if (corsika_code == particles::Code::TauMinus)
return BuildCalculator(particles::Code::TauMinus, PROPOSAL::TauMinusDef(), comp);
if (corsika_code == particles::Code::TauPlus)
return BuildCalculator(particles::Code::TauPlus, PROPOSAL::TauPlusDef(), comp);
throw std::runtime_error("PROPOSAL could not find corresponding builder");
} // namespace corsika::process::proposal
}
enum { SECONDARIES, INTERACTION, DISPLACEMENT };
enum { SECONDARIES, INTERACTION };
template <typename Particle>
auto GetCalculator(Particle& vP) {
auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
auto calc_it = calculators.find(std::make_pair(&comp, vP.GetPID()));
if (calc_it != calculators.end()) return calc_it;
return BuildCalculator(vP.GetPID(), comp);
auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
auto calc_it = calculators.find(make_pair(&comp, vP.GetPID()));
if (calc_it != calculators.end()) return calc_it;
return BuildCalculator(vP.GetPID(), comp);
}
public:
template <typename TEnvironment>
Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
void Init();
void Init() {};
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum