diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index f91ed4e43a7a3a7a50faf8c91c4505f616390b6b..a42cf04f1349c948e4902b76fa0b2bb5ed079756 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -1,3 +1,4 @@
+
#include <corsika/process/proposal/Interaction.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -9,7 +10,8 @@
#include <memory>
#include <limits>
#include <random>
-
+#include <tuple>
+#include <corsika/utl/COMBoost.h>
namespace corsika::process::proposal {
@@ -169,7 +171,7 @@ corsika::units::si::GrammageType Interaction<SetupEnvironment>::GetInteractionLe
double rndiMin = corsika_particle_to_utility_map[vP.GetPID()]->Calculate(particle_energy, 0, 0);
if (rndi >= rndiMin || rndiMin <= 0) {
- return 1 * 1_g / 1_cm / 1_cm;
+ return std::numeric_limits<double>::max() * 1_g / 1_cm / 1_cm;
}
double final_energy = corsika_particle_to_utility_map[vP.GetPID()]->GetUpperLimit(particle_energy, rndi) ;
return corsika_particle_to_displacement_map[vP.GetPID()]->Calculate(
@@ -190,6 +192,8 @@ process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(SetupPartic
double rnd3 = distribution(fRNG);
double particle_energy = vP.GetEnergy() / 1_GeV * 1000;
+ corsika::geometry::Point pOrig = vP.GetPosition();
+ corsika::units::si::TimeType tOrig = vP.GetTime();
corsika::stack::MomentumVector Ptot = vP.GetMomentum();
corsika::geometry::QuantityVector q_momentum = Ptot.GetComponents();
PROPOSAL::Vector3D initial_direction(q_momentum[0] / 1_GeV, q_momentum[1] / 1_GeV, q_momentum[2] / 1_GeV);
@@ -198,31 +202,46 @@ process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(SetupPartic
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<CrossSection*> cross_sections = utility_.GetCrosssections();
- // for (unsigned int i = 0; i < cross_sections.size(); i++) {
- // if (cross_sections[i]->GetTypeId() == energy_loss.second)
- // {
- // deflection_angles = cross_sections[i]->StochasticDeflection(particle_energy, energy_loss.first);
- // // calculate produced particles, get an empty list if no particles
- // // are produced in interaction
- // products = cross_sections[i]->CalculateProducedParticles(particle_energy, energy_loss.first, particle_.GetDirection());
- // }
- // }
-
- // auto pnew = vP.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- // geometry::Point, units::si::TimeType>
- // {
- // process::sibyll::ConvertFromSibyll(psib.GetPID()),
- // Plab.GetTimeLikeComponent(),
- // Plab.GetSpaceLikeComponents(),
- // pOrig,
- // tOrig
- // }
- // );
-
+ 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
- // Point pOrig = vP.GetPosition();
- // TimeType tOrig = vP.GetTime();
// // add to corsika stack
// auto pnew = vP.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index a574e1eb147084cc3caed0e6bf2bf6f1436ebcc3..3defff2a05ea53461fd0862fc61725dcf820da1d 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -41,6 +41,33 @@ namespace corsika::process::proposal {
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 = {
+ {particles::Code::Gamma, PROPOSAL::GammaDef::Get()},
+ {particles::Code::Electron, PROPOSAL::EMinusDef::Get()},
+ {particles::Code::Positron, PROPOSAL::EPlusDef::Get()},
+ {particles::Code::MuMinus, PROPOSAL::MuMinusDef::Get()},
+ {particles::Code::MuPlus, PROPOSAL::MuPlusDef::Get()},
+ {particles::Code::TauPlus, PROPOSAL::TauPlusDef::Get()},
+ {particles::Code::TauMinus, PROPOSAL::TauMinusDef::Get()},
+ };
+
+ // fTrackedParticles[proposal_particle.GetName()] return particle::Code
+
+ bool IsTracked(particles::Code pcode) {
+ for (auto i : convert_corsika_code_to_proposal_particle_def) if (i.first==pcode) return true;
+ return false;
+ };
+
public:
Interaction(TEnvironment const& env, ParticleCut const& cut);