Documentation/Examples/CMakeLists.txt

@@ -22,6 +22,7 @@ target_link_libraries (cascade_example
   CORSIKAlogging
   CORSIKArandom
   ProcessSibyll
+  ProcessProposal
   CORSIKAcascade
   ProcessEnergyLoss
   ProcessStackInspector

Documentation/Examples/cascade_example.cc

@@ -24,6 +24,8 @@
 
 #include <corsika/geometry/Sphere.h>
 
+#include <corsika/process/proposal/Interaction.h>
+
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/NuclearInteraction.h>
@@ -139,6 +141,7 @@ int main() {
   process::sibyll::Decay decay;
   // cascade with only HE model ==> HE cut
   process::particle_cut::ParticleCut cut(80_GeV);
+  process::proposal::Interaction proposal(env, cut);
 
   process::track_writer::TrackWriter trackWriter("tracks.dat");
   process::energy_loss::EnergyLoss eLoss(

Processes/CMakeLists.txt

 # general
 add_subdirectory (NullModel)
 # tracking
-add_subdirectory (TrackingLine) 
+add_subdirectory (TrackingLine)
 # hadron interaction models
 add_subdirectory (Sibyll)
 if (PYTHIA8_FOUND)
@@ -10,6 +10,7 @@ endif (PYTHIA8_FOUND)
 add_subdirectory (HadronicElasticModel)
 add_subdirectory (UrQMD)
 add_subdirectory (SwitchProcess)
+add_subdirectory (Proposal)
 
 # continuous physics
 add_subdirectory (EnergyLoss)
@@ -26,6 +27,7 @@ add_subdirectory (ParticleCut)
 add_library (CORSIKAprocesses INTERFACE)
 add_dependencies(CORSIKAprocesses ProcessNullModel)
 add_dependencies(CORSIKAprocesses ProcessSibyll)
+add_dependencies(CORSIKAprocesses ProcessProposal)
 if (PYTHIA8_FOUND)
   add_dependencies(CORSIKAprocesses ProcessPythia)
 endif (PYTHIA8_FOUND)

Processes/Proposal/CMakeLists.txt

+# find_package(PROPOSAL REQUIRED ${LIB_INCLUDE})
+
+
+FILE (GLOB MODEL_SOURCES *.cc)
+FILE (GLOB MODEL_HEADERS *.h)
+SET (MODEL_NAMESPACE corsika/process/proposal)
+
+ADD_LIBRARY (ProcessProposal STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessProposal ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+SET_TARGET_PROPERTIES ( ProcessProposal PROPERTIES VERSION ${PROJECT_VERSION}
+    SOVERSION 1
+    #  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+TARGET_LINK_LIBRARIES (
+  ProcessProposal
+  CORSIKAenvironment
+  ${PROPOSAL_LIBRARY}
+  )
+
+TARGET_INCLUDE_DIRECTORIES (
+  ProcessProposal
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+TARGET_INCLUDE_DIRECTORIES (
+  ProcessProposal
+  SYSTEM
+  PUBLIC ${PROPOSAL_INCLUDE_DIR}
+  INTERFACE ${PROPOSAL_INCLUDE_DIR}
+  )
+
+INSTALL (
+  TARGETS ProcessProposal
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )

Processes/Proposal/ContinuousProcess.cc

+// #include <PROPOSAL/PROPOSAL.h>
+// #include <corsika/process/proposal/ContinuousProcess.h>
+// // #include <corsika/process/sibyll/Interaction.h>
+
+
+// namespace corsika::process::proposal {
+
+// template <>
+// EProcessReturn ContinuousProcess::DoContinuous(Particle&, Track const&) const {};
+
+// template <>
+// units::si::LengthType ContinuousProcess::MaxStepLength(Particle const& p, Track const& track) const {};
+
+// } // namespace corsika::process

Processes/Proposal/ContinuousProcess.h

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+// #ifndef _corsika_process_proposal_interaction_h_
+// #define _corsika_process_proposalythia_interaction_h_
+
+// #include <PROPOSAL/PROPOSAL.h>
+// #include <corsika/process/ContinuousProcess.h>
+
+
+// namespace corsika::process::proposal {
+
+//   class ContinuousProcess : public corsika::process::ContinuousProcess<Continuous> {
+//   private:
+
+
+//   public:
+//     template <typename Particle, typename Track>
+//     EProcessReturn DoContinuous(Particle&, Track const&) const;
+
+//     template <typename Particle, typename Track>
+//     units::si::LengthType MaxStepLength(Particle const& p, Track const& track) const;
+//   }
+// } // namespace corsika::process
+
+// #endif

Processes/Proposal/Interaction.cc

+
+#include <corsika/process/proposal/Interaction.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <set>
+#include <map>
+#include <memory>
+#include <limits>
+#include <random>
+#include <tuple>
+#include <corsika/utl/COMBoost.h>
+
+namespace corsika::process::proposal {
+
+using namespace corsika::setup;
+using namespace corsika::environment;
+using namespace corsika::units::si;
+using namespace corsika;
+
+using SetupParticle = setup::Stack::StackIterator;
+
+template <>
+Interaction<SetupEnvironment>::Interaction(SetupEnvironment const& env, ParticleCut const& cut)
+      : fEnvironment(env),
+        pCut(cut)
+{
+    using namespace corsika::particles;
+    using namespace units::si;
+
+    auto& universe = *(fEnvironment.GetUniverse());
+
+    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++ ) {
+            comp_vec.push_back(
+				std::make_shared<PROPOSAL::Components::Component>(
+				  PROPOSAL::Components::Component(
+					GetName(pcode[i]),
+					GetNucleusZ(pcode[i]),
+					GetNucleusA(pcode[i]),
+					fractions[i])
+				)
+			);
+        }
+
+        // use ionization constants from air, until CORSIKA implements them
+		PROPOSAL::Air tmp_air;
+        const PROPOSAL::Medium air(
+            tmp_air.GetName(),
+            1., // density correction set to 1 because it cancels out
+            tmp_air.GetI(),
+            tmp_air.GetC(),
+            tmp_air.GetA(),
+            tmp_air.GetM(),
+            tmp_air.GetX0(),
+            tmp_air.GetX1(),
+            tmp_air.GetD0(),
+            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;
+
+}
+
+// Interaction::~Interaction {}
+
+// template <>
+// template <>
+// corsika::process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(SetupParticle const& vP) {
+
+//     // double rndi = dist(rng);
+//     // double rndiMin = corsika_particle_to_utility_map[pcode]->Calculate(initial_energy, p_low, 0);
+
+//     // if (rndd >= rnddMin || rnddMin <= 0) {
+//     //     // return particle_.GetLow();
+//     //     return process::EProcessReturn::eOk;
+//     // } else {
+//     //     return corsika_particle_to_utility_map[pcode]->GetUpperLimit(initial_energy, rndi) ;
+//     // }
+
+//     return process::EProcessReturn::eOk;
+// }
+
+template <>
+template <>
+corsika::units::si::GrammageType Interaction<SetupEnvironment>::GetInteractionLength(SetupParticle const& vP) {
+        std::uniform_real_distribution<double> distribution(0., 1.);
+        double rndi = distribution(fRNG);
+        double particle_energy = vP.GetEnergy()  / 1_GeV * 1000;
+        double rndiMin = corsika_particle_to_utility_map[vP.GetPID()]->Calculate(particle_energy, 0, 0);
+
+        if (rndi >= rndiMin || rndiMin <= 0) {
+            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(
+                    particle_energy,
+                    final_energy,
+                    std::numeric_limits<double>::max(),
+                    PROPOSAL::Vector3D(),
+                    PROPOSAL::Vector3D()
+                ) * 1_g / 1_cm / 1_cm;
+}
+
+template <>
+template <>
+process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(SetupParticle& vP) {
+    std::uniform_real_distribution<double> distribution(0., 1.);
+    double rnd1 = distribution(fRNG);
+    double rnd2 = distribution(fRNG);
+    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);
+	initial_direction.normalise();
+	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
+
+	// // add to corsika stack
+	// auto pnew = vP.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+	// 							geometry::Point, units::si::TimeType>
+	// 							{
+	// 							  process::sibyll::ConvertFromSibyll(psib.GetPID()),
+	// 							  energy_lost,
+                  				  // Plab.GetTimeLikeComponent(),
+	// 							  Plab.GetSpaceLikeComponents(),
+	// 							  pOrig,
+	// 							  tOrig
+	// 							}
+	// 			);
+
+	return process::EProcessReturn::eOk;
+}
+
+} // corsika::process::proposal

Processes/Proposal/Interaction.h

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _corsika_process_proposal_interaction_h_
+#define _corsika_process_proposalythia_interaction_h_
+
+#include "PROPOSAL/PROPOSAL.h"
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <corsika/random/RNGManager.h>
+#include <random>
+
+using namespace corsika::environment;
+using namespace corsika::process::particle_cut;
+
+namespace corsika::process::proposal {
+
+  template <class TEnvironment>
+  class Interaction : public corsika::process::InteractionProcess<Interaction<TEnvironment>> {
+
+  private:
+    TEnvironment const& fEnvironment;
+    ParticleCut const& pCut;
+    // 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 = {
+        {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);
+
+    // ~Interaction;
+
+    template <typename Particle>
+    corsika::process::EProcessReturn DoInteraction(Particle&);
+
+    template <typename TParticle>
+    corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
+
+  };
+
+} // namespace corsika::process
+#endif

tags

+!_TAG_FILE_FORMAT	2	/extended format; --format=1 will not append ;" to lines/
+!_TAG_FILE_SORTED	1	/0=unsorted, 1=sorted, 2=foldcase/
+!_TAG_PROGRAM_AUTHOR	Darren Hiebert	/dhiebert@users.sourceforge.net/
+!_TAG_PROGRAM_NAME	Exuberant Ctags	//
+!_TAG_PROGRAM_URL	http://ctags.sourceforge.net	/official site/
+!_TAG_PROGRAM_VERSION	5.8	//