Initalize interpolation tables

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

@@ -27,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

@@ -5,16 +5,16 @@ 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})
+ADD_LIBRARY (ProcessProposal STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessProposal ${MODEL_NAMESPACE} ${MODEL_HEADERS})
 
-SET_TARGET_PROPERTIES ( ProcessPROPOSAL PROPERTIES VERSION ${PROJECT_VERSION}
+SET_TARGET_PROPERTIES ( ProcessProposal PROPERTIES VERSION ${PROJECT_VERSION}
     SOVERSION 1
     #  PUBLIC_HEADER "${MODEL_HEADERS}"
   )
 
 TARGET_LINK_LIBRARIES (
-  ProcessPROPOSAL
+  ProcessProposal
   CORSIKAenvironment
   ${PROPOSAL_LIBRARY}
   )

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

@@ -8,23 +8,26 @@
  * the license.
  */
 
-#ifndef _corsika_process_proposal_interaction_h_
-#define _corsika_process_proposalythia_interaction_h_
+// #ifndef _corsika_process_proposal_interaction_h_
+// #define _corsika_process_proposalythia_interaction_h_
 
-#include <PROPOSAL/PROPOSAL.h>
+// #include <PROPOSAL/PROPOSAL.h>
+// #include <corsika/process/ContinuousProcess.h>
 
 
-namespace corsika::process::proposal {
+// namespace corsika::process::proposal {
 
-  class Interaction : public corsika::process::ContinuousProcess<Continuous> {
-  private:
+//   class ContinuousProcess : public corsika::process::ContinuousProcess<Continuous> {
+//   private:
 
 
-  public:
-    template <typename Particle, typename Track>
-    EProcessReturn DoContinuous(Particle&, Track const&) const;
+//   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
+//     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 <set>
+#include <map>
+#include <memory>
+
+
+namespace corsika::process::proposal {
+
+using namespace corsika::setup;
+using namespace corsika::environment;
+using namespace corsika;
+
+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::map<particles::Code, std::unique_ptr<const PROPOSAL::UtilityInterpolantInteraction>> 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<const PROPOSAL::UtilityInterpolantInteraction>(
+			PROPOSAL::UtilityInterpolantInteraction(PROPOSAL::Utility(PROPOSAL::MuMinusDef::Get(), medium.second, loss_cut, utility_def, interpolation_def),
+													interpolation_def)
+		);
+	}
+    std::cout << "PROPOSAL done." << std::endl;
+}
+
+// Interaction::~Interaction {}
+
+// template <>
+// corsika::process::EProcessReturn Interaction::DoInteraction(Particle&) {}
+
+// template <>
+// corsika::units::si::GrammageType Interaction::GetInteractionLength(TParticle& p) {}
+
+} // corsika::process::proposal

Processes/Proposal/Interaction.h

@@ -11,21 +11,36 @@
 #ifndef _corsika_process_proposal_interaction_h_
 #define _corsika_process_proposalythia_interaction_h_
 
-#include <PROPOSAL/PROPOSAL.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>
 
+using namespace corsika::environment;
+using namespace corsika::process::particle_cut;
 
 namespace corsika::process::proposal {
 
-  class Interaction : public corsika::process::InteractionProcess<Interaction> {
-  private:
+  template <class TEnvironment>
+  class Interaction : public corsika::process::InteractionProcess<Interaction<TEnvironment>> {
 
+  private:
+    TEnvironment const& fEnvironment;
+    ParticleCut const& pCut;
 
   public:
-    template <typename Particle>
-    EProcessReturn DoInteraction(Particle&);
+    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);
 
-    template <typename TParticle>
-    corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
+  };
 
-  }
 } // namespace corsika::process
+#endif

Processes/Proposal/Interfaces.cc

-#include <iostream>
-#include <corsika/process/Proposal/Proposal.h>
-
-Interface::Interface() { std::cout<< "Hello World!" << std::endl; }

Processes/Proposal/Interfaces.h

-class Interface {
-    Interface();
-};

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	//