diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 1592372c2c6a346def1c866508841097c4bbd1c2..93cb0509a93aedeb38864cfbc91e98fa5d988383 100644 --- a/Documentation/Examples/CMakeLists.txt +++ b/Documentation/Examples/CMakeLists.txt @@ -19,6 +19,7 @@ target_link_libraries (boundary_example CORSIKAlogging CORSIKArandom ProcessSibyll + ProcessProposal CORSIKAcascade ProcessTrackWriter ProcessParticleCut diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 93902e146ae7b223644eb8a0e34fb73cf93ebcab..8ff257cc5f29da29d05ec2bf2bc0cea6de7de50b 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -23,6 +23,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> @@ -140,6 +142,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{showerAxis}; diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index dc7682842ac6139edd90a9fc84937c806c3cf98c..c405641bbcbbae6bc614cb8f6f377465c41a700f 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -36,6 +36,7 @@ add_subdirectory (SwitchProcess) 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) diff --git a/Processes/Proposal/CMakeLists.txt b/Processes/Proposal/CMakeLists.txt index 737cdc20ce5b2c53a1f6cf6952377514bd611f89..970abe29dc0ed296479e2ba20df2f134f0c0d007 100644 --- a/Processes/Proposal/CMakeLists.txt +++ b/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} ) diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc new file mode 100644 index 0000000000000000000000000000000000000000..4293977794951c7be14862ea39744d277d0b7d3a --- /dev/null +++ b/Processes/Proposal/ContinuousProcess.cc @@ -0,0 +1,14 @@ +// #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 diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h index 91d35174672983fa5ed78c80e18a447d32bd3c38..9ed747dac682af11b9e55ce4cc184090c9653407 100644 --- a/Processes/Proposal/ContinuousProcess.h +++ b/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 diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc new file mode 100644 index 0000000000000000000000000000000000000000..44a17868444bd7a70a466462a05683e98a910561 --- /dev/null +++ b/Processes/Proposal/Interaction.cc @@ -0,0 +1,142 @@ +#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 diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h index 4a3f533a1721a02cb72cf70ce416a2687f912b8d..ed3fd2e57e0f0af804c4f5327c60d2e41ce7c5d4 100644 --- a/Processes/Proposal/Interaction.h +++ b/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 diff --git a/Processes/Proposal/Interfaces.cc b/Processes/Proposal/Interfaces.cc deleted file mode 100644 index d6f47e28ec65c33b19e3efaac5108637f2fbb919..0000000000000000000000000000000000000000 --- a/Processes/Proposal/Interfaces.cc +++ /dev/null @@ -1,4 +0,0 @@ -#include <iostream> -#include <corsika/process/Proposal/Proposal.h> - -Interface::Interface() { std::cout<< "Hello World!" << std::endl; } diff --git a/Processes/Proposal/Interfaces.h b/Processes/Proposal/Interfaces.h deleted file mode 100644 index 7dfd75b374335c9e6be82030c961621c3cda5a61..0000000000000000000000000000000000000000 --- a/Processes/Proposal/Interfaces.h +++ /dev/null @@ -1,3 +0,0 @@ -class Interface { - Interface(); -}; diff --git a/tags b/tags new file mode 100644 index 0000000000000000000000000000000000000000..d01d39bc1217a3e782aafdbd41b0cc83b1770d8f --- /dev/null +++ b/tags @@ -0,0 +1,6 @@ +!_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 //