diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc index 523c9400a67d94fd8b9a4152d69ffa3eb1daa0de..5914d529606bf264da310ad46eed8717ddbe6e7e 100644 --- a/Processes/Proposal/ContinuousProcess.cc +++ b/Processes/Proposal/ContinuousProcess.cc @@ -36,22 +36,22 @@ namespace corsika::process::proposal { template <> ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env, CORSIKA_ParticleCut const& _cut) - : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetCutEnergy() / 1_GeV, 1, + : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1, false)) { - auto all_compositions = std::vector<NuclearComposition>(); + auto all_compositions = std::vector<const NuclearComposition*>(); _env.GetUniverse()->walk([&](auto& vtn) { if (vtn.HasModelProperties()) - all_compositions.push_back(vtn.GetModelProperties().GetNuclearComposition()); + all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition()); }); for (auto& ncarg : all_compositions) { auto comp_vec = std::vector<PROPOSAL::Components::Component>(); - auto frac_iter = ncarg.GetFractions().cbegin(); - for (auto& pcode : ncarg.GetComponents()) { + auto frac_iter = ncarg->GetFractions().cbegin(); + for (auto& pcode : ncarg->GetComponents()) { comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode), *frac_iter); ++frac_iter; } - media[&ncarg] = PROPOSAL::Medium( + media[ncarg] = PROPOSAL::Medium( "Modified Air", 1., PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(), PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(), PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec); @@ -62,14 +62,16 @@ namespace corsika::process::proposal { GrammageType const vDX) { auto calc_ptr = GetCalculator(vP); auto upper_energy = calc_ptr->second->UpperLimitTrackIntegral( - vP.GetEnergy() / 1_GeV, vDX / 1_g * 1_cm * 1_cm); - return upper_energy * 1_GeV; + vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm); + std::cout << "upper_energy: " << upper_energy << "MeV" << std::endl; + return upper_energy * 1_MeV; } + void ContinuousProcess::Init() {} template <> EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP, SetupTrack const& vT) { if (vP.GetChargeNumber() == 0) return process::EProcessReturn::eOk; - + std::cout << "DoContinuous..." << std::endl; GrammageType const dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength()); HEPEnergyType dE = TotalEnergyLoss(vP, dX); diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h index ef0782c70f344a090b0fdf868056e362a467fe97..e123ee938f5a0c6a5e20b10489687fa419bccc93 100644 --- a/Processes/Proposal/ContinuousProcess.h +++ b/Processes/Proposal/ContinuousProcess.h @@ -42,18 +42,45 @@ namespace corsika::process::proposal { bool CanInteract(particles::Code pcode) const noexcept; - unordered_map<const NuclearComposition*, unique_ptr<PROPOSAL::Displacement>> calc; - + 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; + } + }; + + unordered_map<std::pair<const NuclearComposition*, particles::Code>, unique_ptr<PROPOSAL::Displacement>, interaction_hash> calc; + + auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) { + auto medium = media.at(&comp); + if (corsika_code == particles::Code::Gamma) { + auto cross = GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true); + auto [insert_it, success] = + calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)}); + return insert_it; + } + if (corsika_code == particles::Code::Electron) { + auto cross = GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true); + auto [insert_it, success] = + calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)}); + return insert_it; + } + if (corsika_code == particles::Code::Positron) { + auto cross = GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true); + auto [insert_it, success] = + calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)}); + return insert_it; + } + } + template <typename Particle> auto GetCalculator(Particle& vP) { auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition(); - auto calc_it = calc.find(&comp); + auto calc_it = calc.find(std::make_pair(&comp, vP.GetPID())); if (calc_it != calc.end()) return calc_it; - auto cross = - PROPOSAL::GetStdCrossSections(particles[vP.GetPID()], media[&comp], cut, true); - auto [insert_it, success] = - calc.insert({&comp, PROPOSAL::make_displacement(cross, true)}); - return insert_it; + return BuildCalculator(vP.GetPID(), comp); } units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&, @@ -63,6 +90,8 @@ namespace corsika::process::proposal { template <typename TEnvironment> ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut); + void Init(); + template <typename Particle, typename Track> EProcessReturn DoContinuous(Particle&, Track const&) ;