IAP GITLAB

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

secondary add still not yet work

parent e708331b
No related branches found
No related tags found
1 merge request!245Include proposal process rebase
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/process/proposal/Interaction.h> #include <corsika/process/proposal/Interaction.h>
#include <corsika/setup/SetupStack.h> #include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h> #include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/units/PhysicalUnits.h> #include <corsika/units/PhysicalUnits.h>
#include <set> #include <corsika/utl/COMBoost.h>
#include <limits>
#include <map> #include <map>
#include <memory> #include <memory>
#include <limits>
#include <random> #include <random>
#include <set>
#include <tuple> #include <tuple>
#include <corsika/utl/COMBoost.h>
using Component_PROPOSAL = PROPOSAL::Components::Component;
namespace corsika::process::proposal { namespace corsika::process::proposal {
using namespace corsika::setup; using namespace corsika::setup;
using namespace corsika::environment; using namespace corsika::environment;
using namespace corsika::units::si; using namespace corsika::units::si;
using namespace corsika; using namespace corsika;
using SetupParticle = setup::Stack::StackIterator; using SetupParticle = corsika::setup::Stack::StackIterator;
template <> template <>
Interaction<SetupEnvironment>::Interaction(SetupEnvironment const& env, ParticleCut const& cut) std::unordered_map<particles::Code, PROPOSAL::ParticleDef>
: fEnvironment(env), Interaction<SetupEnvironment>::particle_map{
pCut(cut) {particles::Code::Gamma, PROPOSAL::GammaDef()},
{ {particles::Code::Electron, PROPOSAL::EMinusDef()},
using namespace corsika::particles; {particles::Code::Positron, PROPOSAL::EPlusDef()},
using namespace units::si; {particles::Code::MuMinus, PROPOSAL::MuMinusDef()},
{particles::Code::MuPlus, PROPOSAL::MuPlusDef()},
auto& universe = *(fEnvironment.GetUniverse()); {particles::Code::TauPlus, PROPOSAL::TauPlusDef()},
{particles::Code::TauMinus, PROPOSAL::TauMinusDef()},
auto const allCompositionsInUniverse = std::invoke([&]() {
std::vector<NuclearComposition> allCompositionsInUniverse;
auto collectCompositions = [&](auto& vtn) {
if (vtn.HasModelProperties()) {
auto const& comp =
vtn.GetModelProperties().GetNuclearComposition();
allCompositionsInUniverse.push_back(comp);
}
}; };
universe.walk(collectCompositions);
return allCompositionsInUniverse;
});
std::map<const NuclearComposition*, PROPOSAL::Medium> composition_medium_map;
for (const NuclearComposition& ncarg : allCompositionsInUniverse) {
std::vector<particles::Code> pcode { ncarg.GetComponents() };
std::vector<float> fractions { ncarg.GetFractions() };
std::vector<std::shared_ptr<PROPOSAL::Components::Component>> comp_vec;
for(unsigned int i=0; i< pcode.size(); i++ ) { template <>
comp_vec.push_back( Interaction<SetupEnvironment>::Interaction(SetupEnvironment const& env,
std::make_shared<PROPOSAL::Components::Component>( CORSIKA_ParticleCut const& e_cut)
PROPOSAL::Components::Component( : fEnvironment(env)
GetName(pcode[i]), , cut(make_shared<const PROPOSAL::EnergyCutSettings>(e_cut.GetCutEnergy() / 1_GeV,
GetNucleusZ(pcode[i]), 0, false)) {
GetNucleusA(pcode[i]),
fractions[i]) auto all_compositions = std::vector<NuclearComposition>();
) fEnvironment.GetUniverse()->walk([&](auto& vtn) {
); if (vtn.HasModelProperties())
} all_compositions.push_back(vtn.GetModelProperties().GetNuclearComposition());
});
// use ionization constants from air, until CORSIKA implements them for (auto& ncarg : all_compositions) {
PROPOSAL::Air tmp_air; auto comp_vec = std::vector<Component_PROPOSAL>();
const PROPOSAL::Medium air( auto frac_iter = ncarg.GetFractions().cbegin();
tmp_air.GetName(), for (auto& pcode : ncarg.GetComponents()) {
1., // density correction set to 1 because it cancels out comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
tmp_air.GetI(), *frac_iter);
tmp_air.GetC(), ++frac_iter;
tmp_air.GetA(), }
tmp_air.GetM(), medium_map[&ncarg] = PROPOSAL::Medium(
tmp_air.GetX0(), "Modified Air", 1., PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
tmp_air.GetX1(), PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
tmp_air.GetD0(), PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
1.0,// massDensity set to 1, because it cancels out
comp_vec
);
composition_medium_map[&ncarg] = air;
} }
const double e_cut = 500.;
const double v_cut = -1.;
PROPOSAL::EnergyCutSettings loss_cut(e_cut,v_cut); // energy_loss, relatic_energy loss (negativ means unused)
std::string bremsstrahlung_param_str ="BremsKelnerKokoulinPetrukhin";
std::string photonuclear_param_str = "PhotoAbramowiczLevinLevyMaor97";
std::string epairproduction_param_str = "EpairKelnerKokoulinPetrukhin";
std::string ionization_param_str = "IonizBetheBlochRossi"; // IonizBergerSeltzerBhabha IonizBergerSeltzerMollerg
std::string shadowing_str = "ShadowButkevichMikhailov";
std::string annihilation_param_str = "None"; // Heitler
std::string compton_param_str = "None"; // KleinNishina
std::string muPair_param_str = "None"; // KelnerKokoulinPetrukhin
std::string photoPair_param_str = "None"; // Tsai
std::string weak_param_str = "None"; // CooperSarkarMertsch
PROPOSAL::Utility::Definition utility_def;
utility_def.brems_def.lpm_effect = false;
utility_def.epair_def.lpm_effect = false;
utility_def.photo_def.hard_component = false;
utility_def.photo_def.shadow = PROPOSAL::PhotonuclearFactory::Get().GetShadowEnumFromString(shadowing_str);
utility_def.brems_def.parametrization = PROPOSAL::BremsstrahlungFactory::Get().GetEnumFromString(bremsstrahlung_param_str);
utility_def.epair_def.parametrization = PROPOSAL::EpairProductionFactory::Get().GetEnumFromString(epairproduction_param_str);
utility_def.ioniz_def.parametrization = PROPOSAL::IonizationFactory::Get().GetEnumFromString(ionization_param_str);
utility_def.photo_def.parametrization = PROPOSAL::PhotonuclearFactory::Get().GetEnumFromString(photonuclear_param_str);
utility_def.annihilation_def.parametrization = PROPOSAL::AnnihilationFactory::Get().GetEnumFromString(annihilation_param_str);
utility_def.compton_def.parametrization = PROPOSAL::ComptonFactory::Get().GetEnumFromString(compton_param_str);
utility_def.mupair_def.parametrization = PROPOSAL::MupairProductionFactory::Get().GetEnumFromString(muPair_param_str);
utility_def.photopair_def.parametrization = PROPOSAL::PhotoPairFactory::Get().GetEnumFromString(photoPair_param_str);
utility_def.weak_def.parametrization = PROPOSAL::WeakInteractionFactory::Get().GetEnumFromString(weak_param_str);
PROPOSAL::InterpolationDef interpolation_def;
interpolation_def.path_to_tables = "~/.local/share/PROPOSAL/tables";
interpolation_def.path_to_tables_readonly = "~/.local/share/PROPOSAL/tables";
interpolation_def.order_of_interpolation = 5; // order of interpolation
interpolation_def.max_node_energy = 1e14; // upper energy bound for Interpolation (MeV)
interpolation_def.nodes_cross_section = 100; // number of interpolation in cross section
interpolation_def.nodes_propagate = 1000; // number of interpolation in propagate
interpolation_def.do_binary_tables = true;
interpolation_def.just_use_readonly_path = false;
// std::map<particles::Code, double> corsika_particle_to_utility_map;
std::cout << "go PROPOSAL !!!" << std::endl;
for(auto const& medium:composition_medium_map){
corsika_particle_to_utility_map[Code::MuMinus] = std::make_unique<PROPOSAL::UtilityInterpolantInteraction>(
PROPOSAL::UtilityInterpolantInteraction(PROPOSAL::Utility(PROPOSAL::MuMinusDef::Get(), medium.second, loss_cut, utility_def, interpolation_def),
interpolation_def)
);
corsika_particle_to_displacement_map[Code::MuMinus] = std::make_unique<PROPOSAL::UtilityInterpolantDisplacement>(
PROPOSAL::UtilityInterpolantDisplacement(PROPOSAL::Utility(PROPOSAL::MuMinusDef::Get(), medium.second, loss_cut, utility_def, interpolation_def),
interpolation_def)
);
}
std::cout << "PROPOSAL done." << std::endl; std::cout << "PROPOSAL done." << std::endl;
}
}
template <>
// Interaction::~Interaction {} template <>
auto Interaction<SetupEnvironment>::GetCalculator(SetupParticle const& vP) {
// template <> auto& mediumComposition = vP.GetNode()->GetModelProperties().GetNuclearComposition();
// template <> auto calc_it = calculators.find(&mediumComposition);
// corsika::process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(SetupParticle const& vP) { if (calc_it != calculators.end()) return calc_it;
auto const& p_def = particle_map[vP.GetPID()];
// // double rndi = dist(rng); auto& medium = medium_map[&mediumComposition];
// // double rndiMin = corsika_particle_to_utility_map[pcode]->Calculate(initial_energy, p_low, 0); auto cross = PROPOSAL::GetStdCrossSections(p_def, medium, cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
// // if (rndd >= rnddMin || rnddMin <= 0) { auto [insert_it, success] = calculators.insert(
// // // return particle_.GetLow(); {&mediumComposition,
// // return process::EProcessReturn::eOk; make_tuple(PROPOSAL::SecondariesCalculator(inter_types, p_def, medium),
// // } else { PROPOSAL::make_interaction(cross, true),
// // return corsika_particle_to_utility_map[pcode]->GetUpperLimit(initial_energy, rndi) ; PROPOSAL::make_displacement(cross, true))});
// // } if (success) return insert_it;
throw std::logic_error("insert failed.");
// return process::EProcessReturn::eOk; }
// }
template <>
template <> template <>
template <> corsika::process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(
corsika::units::si::GrammageType Interaction<SetupEnvironment>::GetInteractionLength(SetupParticle const& vP) { SetupParticle const& vP) {
std::uniform_real_distribution<double> distribution(0., 1.); auto calc = GetCalculator(vP); // [CrossSections]
double rndi = distribution(fRNG); std::uniform_real_distribution<double> distr(0., 1.);
double particle_energy = vP.GetEnergy() / 1_GeV * 1000; auto [type, comp_ptr, v] = std::get<INTERACTION>(calc->second)
double rndiMin = corsika_particle_to_utility_map[vP.GetPID()]->Calculate(particle_energy, 0, 0); ->TypeInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
auto rnd = std::vector<double>();
if (rndi >= rndiMin || rndiMin <= 0) { for (auto i = 0; i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type);
return std::numeric_limits<double>::max() * 1_g / 1_cm / 1_cm; ++i)
} rnd.push_back(distr(fRNG));
double final_energy = corsika_particle_to_utility_map[vP.GetPID()]->GetUpperLimit(particle_energy, rndi) ; double primary_energy = vP.GetEnergy() / 1_GeV;
return corsika_particle_to_displacement_map[vP.GetPID()]->Calculate( auto point =
particle_energy, PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm, vP.GetPosition().GetY() / 1_cm,
final_energy, vP.GetPosition().GetZ() / 1_cm);
std::numeric_limits<double>::max(), auto p = vP.GetMomentum().GetComponents();
PROPOSAL::Vector3D(), auto direction = PROPOSAL::Vector3D(p[0] / 1_GeV, p[1] / 1_GeV, p[2] / 1_GeV);
PROPOSAL::Vector3D() auto loss =
) * 1_g / 1_cm / 1_cm; 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);
template <> for (auto& s : sec) {
template <> auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(SetupParticle& vP) { auto vec = corsika::geometry::QuantityVector(
std::uniform_real_distribution<double> distribution(0., 1.); std::get<PROPOSAL::Loss::DIRECTION>(s).GetX() * energy,
double rnd1 = distribution(fRNG); std::get<PROPOSAL::Loss::DIRECTION>(s).GetY() * energy,
double rnd2 = distribution(fRNG); std::get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * energy);
double rnd3 = distribution(fRNG); auto momentum = corsika::stack::MomentumVector(
double particle_energy = vP.GetEnergy() / 1_GeV * 1000; corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(),
corsika::geometry::Point pOrig = vP.GetPosition(); vec);
corsika::units::si::TimeType tOrig = vP.GetTime(); // particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
corsika::stack::MomentumVector Ptot = vP.GetMomentum(); // geometry::Point, units::si::TimeType
corsika::geometry::QuantityVector q_momentum = Ptot.GetComponents(); auto sec = make_tuple(get<PROPOSAL::Loss::TYPE>(s), energy, momentum,
PROPOSAL::Vector3D initial_direction(q_momentum[0] / 1_GeV, q_momentum[1] / 1_GeV, q_momentum[2] / 1_GeV); vP.GetPosition(), vP.GetTime());
initial_direction.normalise(); vP.AddSecondary(sec);
initial_direction.CalculateSphericalCoordinates();
std::pair<double, PROPOSAL::DynamicData::Type> energy_loss;
energy_loss = corsika_particle_to_utility_map[vP.GetPID()]->GetUtility().StochasticLoss(particle_energy, rnd1, rnd2, rnd3);
std::vector<PROPOSAL::CrossSection*> cross_sections = corsika_particle_to_utility_map[vP.GetPID()]->GetUtility().GetCrosssections();
std::pair<std::vector<PROPOSAL::Particle*>, bool> products;
for (unsigned int i = 0; i < cross_sections.size(); i++) {
if (cross_sections[i]->GetTypeId() == energy_loss.second)
{
// calculate produced particles, get an empty list if no particles
// are produced in interaction
products = cross_sections[i]->CalculateProducedParticles(particle_energy, energy_loss.first, initial_direction);
break;
}
}
double particle_energy_after_loss = particle_energy - energy_loss.first;
if (particle_energy_after_loss == 0.)
{
PROPOSAL::Vector3D particle_direction = initial_direction * particle_energy_after_loss;
corsika::geometry::QuantityVector q_vec(
particle_direction.GetX() * 1_MeV,
particle_direction.GetY() * 1_MeV,
particle_direction.GetZ() * 1_MeV);
// vp.GetCoordinateSystem();
const corsika::geometry::CoordinateSystem& rootCS = corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
corsika::stack::MomentumVector p_particle(rootCS, q_vec);
// std::cout << p_particle.GetTimeLikeComponent() << " and " << particle_energy_after_loss << std::endl;
units::si::HEPEnergyType const converted_particle_energy = particle_energy * 1_GeV;
HEPEnergyType const eCoM = particle_energy * 1_GeV ;
//auto const Plab = boost.fromCoM(FourVector(eCoM, p_particle));
auto pnew = vP.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>
{
//convert_proposal_particle_name_to_corsika_code[PROPOSAL::MuMinusDef::Get().name],
vP.GetPID(),
eCoM,
p_particle,
pOrig,
tOrig
}
);
} }
// // position and time of interaction, not used in PROPOSAL return process::EProcessReturn::eOk;
}
// // add to corsika stack
// auto pnew = vP.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector, template <>
// geometry::Point, units::si::TimeType> template <>
// { corsika::units::si::GrammageType Interaction<SetupEnvironment>::GetInteractionLength(
// process::sibyll::ConvertFromSibyll(psib.GetPID()), SetupParticle const& vP) {
// energy_lost, auto calc = GetCalculator(vP); // [CrossSections]
// Plab.GetTimeLikeComponent(), std::uniform_real_distribution<double> distr(0., 1.);
// Plab.GetSpaceLikeComponents(), auto energy = get<INTERACTION>(calc->second)
// pOrig, ->EnergyInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
// tOrig return get<DISPLACEMENT>(calc->second)
// } ->SolveTrackIntegral(vP.GetEnergy() / 1_GeV, energy) *
// ); 1_g / 1_cm / 1_cm;
}
return process::EProcessReturn::eOk; } // namespace corsika::process::proposal
}
} // corsika::process::proposal
...@@ -11,76 +11,80 @@ ...@@ -11,76 +11,80 @@
#ifndef _corsika_process_proposal_interaction_h_ #ifndef _corsika_process_proposal_interaction_h_
#define _corsika_process_proposalythia_interaction_h_ #define _corsika_process_proposalythia_interaction_h_
#include "PROPOSAL/PROPOSAL.h" #include <corsika/environment/Environment.h>
#include <corsika/particles/ParticleProperties.h> #include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h> #include <corsika/process/InteractionProcess.h>
#include <corsika/environment/Environment.h>
#include <corsika/process/particle_cut/ParticleCut.h> #include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/random/UniformRealDistribution.h>
#include <corsika/random/RNGManager.h> #include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
#include <random> #include <random>
#include <unordered_map>
#include "PROPOSAL/PROPOSAL.h"
using namespace corsika::environment; using namespace corsika::environment;
using namespace corsika::process::particle_cut;
using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
namespace corsika::process::proposal { namespace corsika::process::proposal {
/* static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particle_map{ */
/* {particles::Code::Gamma, PROPOSAL::GammaDef()}, */
/* {particles::Code::Electron, PROPOSAL::EMinusDef()}, */
/* {particles::Code::Positron, PROPOSAL::EPlusDef()}, */
/* {particles::Code::MuMinus, PROPOSAL::MuMinusDef()}, */
/* {particles::Code::MuPlus, PROPOSAL::MuPlusDef()}, */
/* {particles::Code::TauPlus, PROPOSAL::TauPlusDef()}, */
/* {particles::Code::TauMinus, PROPOSAL::TauMinusDef()}, */
/* }; */
template <class TEnvironment> template <class TEnvironment>
class Interaction : public corsika::process::InteractionProcess<Interaction<TEnvironment>> { class Interaction
: public corsika::process::InteractionProcess<Interaction<TEnvironment>> {
private: private:
TEnvironment const& fEnvironment; TEnvironment const& fEnvironment;
ParticleCut const& pCut; shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
// Initializing of uniform_real_distribution class
// double min{0.0};
// double max{1.0};
// std::uniform_real_distribution<double> rnd_uniform(min,max);
std::map<particles::Code, std::unique_ptr<PROPOSAL::UtilityInterpolantInteraction>> corsika_particle_to_utility_map;
std::map<particles::Code, std::unique_ptr<PROPOSAL::UtilityInterpolantDisplacement>> corsika_particle_to_displacement_map;
corsika::random::RNG& fRNG = corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
std::map<std::string, particles::Code> const convert_proposal_particle_name_to_corsika_code = {
{"Gamma", particles::Code::Gamma},
{"EMinus", particles::Code::Electron},
{"EPlus", particles::Code::Positron},
{"MuMinu", particles::Code::MuMinus},
{"MuPlus", particles::Code::MuPlus},
{"TauPlus", particles::Code::TauPlus},
{"TauMinus",particles::Code::TauMinus},
};
std::map<particles::Code, PROPOSAL::ParticleDef> const convert_corsika_code_to_proposal_particle_def = { static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particle_map;
{particles::Code::Gamma, PROPOSAL::GammaDef::Get()}, std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> medium_map;
{particles::Code::Electron, PROPOSAL::EMinusDef::Get()},
{particles::Code::Positron, PROPOSAL::EPlusDef::Get()}, enum { SECONDARIES, INTERACTION, DISPLACEMENT };
{particles::Code::MuMinus, PROPOSAL::MuMinusDef::Get()}, /* std::uniform_real_distribution<> rnd_uniform(0., 1.); */
{particles::Code::MuPlus, PROPOSAL::MuPlusDef::Get()},
{particles::Code::TauPlus, PROPOSAL::TauPlusDef::Get()}, /* std::map<particles::Code, std::unique_ptr<PROPOSAL::UtilityInterpolantInteraction>>
{particles::Code::TauMinus, PROPOSAL::TauMinusDef::Get()}, * corsika_particle_to_utility_map; */
}; /* std::map<particles::Code,
* std::unique_ptr<PROPOSAL::UtilityInterpolantDisplacement>>
* corsika_particle_to_displacement_map; */
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
// fTrackedParticles[proposal_particle.GetName()] return particle::Code // fTrackedParticles[proposal_particle.GetName()] return particle::Code
bool IsTracked(particles::Code pcode) { auto IsTracked(particles::Code pcode) const noexcept {
for (auto i : convert_corsika_code_to_proposal_particle_def) if (i.first==pcode) return true; auto search = particle_map.find(pcode);
return false; if (search != particle_map.end()) return true;
return false;
}; };
using calculator_t =
tuple<PROPOSAL::SecondariesCalculator, unique_ptr<PROPOSAL::Interaction>,
unique_ptr<PROPOSAL::Displacement>>;
std::unordered_map<const NuclearComposition*, calculator_t> calculators;
public: template <typename Particle>
Interaction(TEnvironment const& env, ParticleCut const& cut); auto GetCalculator(Particle&);
// ~Interaction; public:
Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
template <typename Particle> template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&); corsika::process::EProcessReturn DoInteraction(Particle&);
template <typename TParticle> template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle& p); corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
}; };
} // namespace corsika::process } // namespace corsika::process::proposal
#endif #endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment