diff --git a/CMakeLists.txt b/CMakeLists.txt index 14d59a7fdebd00fc3e45c50ef493677c669d595a..216bb098c81608a19999f7f975b8a99cb16a4aa5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,6 +45,13 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O3 -g") # -O2 would not trade speed for size, neit # clang produces a lot of unecessary warnings without this: add_compile_options("$<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>>:-Wno-nonportable-include-path>") +# COAST - interface +if (WITH_COAST) + message(STATUS "Compiling CORSIKA8 for the use with COAST/corsika7.") + add_compile_options("-fPIC") +endif() + + # unit testing coverage, does not work yet #include (CodeCoverage) ##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*') @@ -79,3 +86,6 @@ add_subdirectory (Processes) add_subdirectory (Documentation) add_subdirectory (Main) add_subdirectory (Tools) +if (WITH_COAST) + add_subdirectory (COAST) +endif() diff --git a/COAST/CMakeLists.txt b/COAST/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..ae85a9176aee3e55efb226d3fae4df0ae5a7e212 --- /dev/null +++ b/COAST/CMakeLists.txt @@ -0,0 +1,61 @@ + +set ( + COAST_HEADERS + COASTProcess.h + COASTStack.h + ParticleConversion.h + ) + +set ( + COAST_SOURCES + COASTUserLib.cc + COASTProcess.cc + ParticleConversion.cc + ) + +set ( + COAST_NAMESPACE + corsika/coast + ) + +add_library (COAST SHARED ${COAST_SOURCES}) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (COAST ${COAST_NAMESPACE} ${COAST_HEADERS}) + +set_target_properties ( + COAST + PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION 1 +# PUBLIC_HEADER "${MODEL_HEADERS}" + ) + +target_link_libraries ( + COAST + PUBLIC + CORSIKAgeometry + CORSIKAunits + CORSIKAparticles + CORSIKAgeometry + CORSIKAsetup + # SuperStupidStack + ) + +target_include_directories ( + COAST + PUBLIC + $ENV{COAST_DIR}/include + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) + +install ( + TARGETS COAST + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib + PUBLIC_HEADER DESTINATION include/${COAST_NAMESPACE} + ) + +#install ( +# FILES ${COAST_HEADERS} +# DESTINATION include/${COAST_NAMESPACE} +# ) diff --git a/COAST/COASTProcess.cc b/COAST/COASTProcess.cc new file mode 100644 index 0000000000000000000000000000000000000000..d362bc044a619e6b327c1f5dfc4d2b022ec4e19e --- /dev/null +++ b/COAST/COASTProcess.cc @@ -0,0 +1,44 @@ +/** + * (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. + */ + +#include <corsika/coast/COASTProcess.h> +#include <corsika/coast/COASTStack.h> + +#include <iostream> + +using namespace std; + +namespace corsika::coast { + + /** + the init function is called during the start of corsika. + */ + void COASTProcess::Init() { + cout << "************* Greetings from CORSIKA8 *****************" << endl; + } + + /** + the docontinous function is called for each tracking step in + corsika. Take care: you cannot modify the particle, the track or + the stack from here (docontinuous) inside corisika7. In corsika8 + you will be able to do that. + */ + corsika::process::EProcessReturn COASTProcess::DoContinuous(Particle& p, Track& t, + Stack&) { + using namespace corsika::units::si; + auto const start = t.GetPosition(0).GetCoordinates(); + auto const delta = t.GetPosition(1).GetCoordinates() - start; + auto const name = corsika::particles::GetName(p.GetPID()); + cout << "CORSIKA8: particle=" << name << ", pos=" << start + << " track-l=" << delta.norm() << ", track-t=" << t.GetDuration() << endl; + return corsika::process::EProcessReturn::eOk; + } + +} // namespace corsika::coast diff --git a/COAST/COASTProcess.h b/COAST/COASTProcess.h new file mode 100644 index 0000000000000000000000000000000000000000..19d1db2c6ff2fc407ad91fd84b6f6e0d7144c805 --- /dev/null +++ b/COAST/COASTProcess.h @@ -0,0 +1,35 @@ +#ifndef _include_corsika_coast_coastprocess_h_ +#define _include_corsika_coast_coastprocess_h_ + +#include <corsika/particles/ParticleProperties.h> +#include <corsika/process/ContinuousProcess.h> +#include <corsika/setup/SetupTrajectory.h> +#include <corsika/units/PhysicalUnits.h> + +#include <corsika/coast/COASTStack.h> + +#include <limits> + +typedef corsika::coast::COASTStack Stack; +typedef corsika::coast::COASTStack::ParticleType Particle; +typedef corsika::geometry::Trajectory<corsika::geometry::Line> Track; + +namespace corsika::coast { + + class COASTProcess : public corsika::process::ContinuousProcess<COASTProcess> { + + public: + void Init(); + + corsika::process::EProcessReturn DoContinuous(Particle&, Track&, Stack&); + + corsika::units::si::LengthType MaxStepLength(Particle&, Track&) { + return corsika::units::si::meter * std::numeric_limits<double>::infinity(); + } + + private: + }; + +} // namespace corsika::coast + +#endif diff --git a/COAST/COASTStack.h b/COAST/COASTStack.h new file mode 100644 index 0000000000000000000000000000000000000000..d5ac8c738869f3680faaebed8117c9ccff9cfb26 --- /dev/null +++ b/COAST/COASTStack.h @@ -0,0 +1,193 @@ + +/** + * (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 _include_coaststack_h_ +#define _include_coaststack_h_ + +#include <corsika/coast/ParticleConversion.h> +#include <corsika/particles/ParticleProperties.h> +#include <corsika/stack/ParticleBase.h> +#include <corsika/stack/Stack.h> +#include <corsika/units/PhysicalUnits.h> + +#include <corsika/geometry/Point.h> +#include <corsika/geometry/RootCoordinateSystem.h> // remove +#include <corsika/geometry/Vector.h> + +#include <crs/CParticle.h> +#include <crs/CorsikaTypes.h> + +#include <algorithm> +#include <vector> + +namespace corsika::coast { + + typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector; + + /** + * Example of a particle object on the stack. + */ + + template <typename StackIteratorInterface> + class ParticleInterface : public corsika::stack::ParticleBase<StackIteratorInterface> { + + using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData; + using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex; + + public: + corsika::particles::Code GetPID() const { return GetStackData().GetPID(GetIndex()); } + + corsika::units::si::HEPEnergyType GetEnergy() const { + return GetStackData().GetEnergy(GetIndex()); + } + + MomentumVector GetMomentum(const corsika::geometry::CoordinateSystem& cs) const { + using namespace corsika::units::si; + const HEPEnergyType mass = corsika::particles::GetMass(GetPID()); + auto P = sqrt((GetEnergy() - mass) * (GetEnergy() + mass)); + auto p = GetStackData().GetDirection(cs, GetIndex()); + return p * P; + } + + corsika::geometry::Point GetPosition( + const corsika::geometry::CoordinateSystem& cs) const { + return GetStackData().GetPosition(cs, GetIndex()); + } + + corsika::geometry::Vector<corsika::units::si::speed_d> GetVelocity( + const corsika::geometry::CoordinateSystem& cs) const { + return GetStackData().GetVelocity(cs, GetIndex()); + } + + corsika::units::si::TimeType GetTime() const { + return GetStackData().GetTime(GetIndex()); + } + + corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection( + const corsika::geometry::CoordinateSystem& cs) const { + return GetStackData().GetDirection(cs, GetIndex()); + } + + corsika::units::si::TimeType GetTimeInterval() const { + return GetStackData().GetTimeInterval(); + } + }; + + /** + * + * Memory implementation of the most simple (stupid) particle stack object. + */ + + class COASTStackImpl { + + const crs::CParticle* fParticle1 = 0; + const crs::CParticle* fParticle2 = 0; + + public: + void Init() {} + void SetParticle(const crs::CParticle* v1, const crs::CParticle* v2) { + fParticle1 = v1; + fParticle2 = v2; + } + void Clear() {} + + // there is one particle only + int GetSize() const { return 1; } + int GetCapacity() const { return 1; } + + // you cannot modify the particle: + // there are no Set... function defined + + // readout particle data, there is just one particle! + /* + double x; + double y; + double z; + double depth; + double time; + double energy; + double weight; + int particleId; + int hadronicGeneration; + */ + corsika::particles::Code GetPID(const int) const { + return ConvertFromCoast(static_cast<CoastCode>(fParticle1->particleId)); + } + corsika::units::si::HEPEnergyType GetEnergy(const int) const { + using namespace corsika::units::si; + return fParticle1->energy * 1_GeV; + } + corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection( + const corsika::geometry::CoordinateSystem& cs, const int) const { + using namespace corsika::units::si; + corsika::geometry::Point p1( + cs, {fParticle1->x * 1_cm, fParticle1->y * 1_cm, fParticle1->y * 1_cm}); + corsika::geometry::Point p2( + cs, {fParticle2->x * 1_cm, fParticle2->y * 1_cm, fParticle2->y * 1_cm}); + const corsika::geometry::Vector D = p2 - p1; + const auto magD = D.norm(); + const corsika::geometry::Vector dir = D / magD; + return dir; + } + corsika::geometry::Vector<corsika::units::si::speed_d> GetVelocity( + const corsika::geometry::CoordinateSystem& cs, const int) const { + using namespace corsika::units::si; + corsika::geometry::Vector<corsika::units::si::dimensionless_d> dir = + GetDirection(cs, 0); + corsika::geometry::Point p1( + cs, {fParticle1->x * 1_cm, fParticle1->y * 1_cm, fParticle1->y * 1_cm}); + corsika::geometry::Point p2( + cs, {fParticle2->x * 1_cm, fParticle2->y * 1_cm, fParticle2->y * 1_cm}); + const corsika::geometry::Vector D = p2 - p1; + const LengthType magD = D.norm(); + const TimeType deltaT = (fParticle2->time - fParticle1->time) * 1_ns; + return dir * magD / deltaT; + } + corsika::geometry::Point GetPosition(const corsika::geometry::CoordinateSystem& cs, + const int) const { + using namespace corsika::units::si; + return corsika::geometry::Point( + cs, {fParticle1->x * 1_cm, fParticle1->y * 1_cm, fParticle1->z * 1_cm}); + } + corsika::units::si::TimeType GetTime(const int) const { + using namespace corsika::units::si; + return fParticle1->time * 1_ns; + } + + corsika::units::si::TimeType GetTimeInterval() const { + using namespace corsika::units::si; + return (fParticle2->time - fParticle1->time) * 1_ns; + } + + /** + * We do not allow copying! + */ + void Copy(const int, const int) {} + + /** + * We do not allow swapping particles, there is just one... + */ + void Swap(const int, const int) {} + + protected: + // size cannot be increased + void IncrementSize() {} + + // size cannot be decremented + void DecrementSize() {} + + }; // end class COASTStackImpl + + typedef corsika::stack::Stack<COASTStackImpl, ParticleInterface> COASTStack; + +} // namespace corsika::coast + +#endif diff --git a/COAST/COASTUserLib.cc b/COAST/COASTUserLib.cc new file mode 100644 index 0000000000000000000000000000000000000000..ab4c3d987c81f65dcb41e1b6a6e61f3764e15ae0 --- /dev/null +++ b/COAST/COASTUserLib.cc @@ -0,0 +1,102 @@ +#include <interface/CorsikaInterface.h> + +#include <corsika/coast/COASTProcess.h> +#include <corsika/coast/COASTStack.h> +#include <corsika/geometry/CoordinateSystem.h> +#include <corsika/geometry/Line.h> +#include <corsika/geometry/Point.h> +#include <corsika/geometry/Trajectory.h> +#include <corsika/geometry/Vector.h> +#include <corsika/units/PhysicalUnits.h> + +#include <crs/CInteraction.h> +#include <crs/CParticle.h> +#include <crs/CorsikaTypes.h> +#include <crs/TSubBlock.h> + +#include <iostream> +#include <sstream> + +using namespace std; +using namespace corsika; +using namespace corsika::units::si; + +coast::COASTStack gCOASTStack; +corsika::coast::COASTProcess gCorsikaProcess; + +/* + Data is one CORSIKA data-block constining of 21 SubBlocks. + A SubBlock can be: + - thinned mode: 39 (Particles) * 8 (ENTRIES) * 4 (BYTES) + - not-thinned mode: 39 (Particles) * 7 (ENTRIES) * 4 (BYTES) +*/ +extern "C" void wrida_([[maybe_unused]] const CREAL* Data) { + // crs::CParticleFortranPtr p; + // const bool isF = prminfo_(p); +} + +extern "C" void inida_([[maybe_unused]] const char* filename, + [[maybe_unused]] const int& thinning, + [[maybe_unused]] const int& /*curved*/, + [[maybe_unused]] const int& /*slant*/, + [[maybe_unused]] const int& /*stackinput*/, + [[maybe_unused]] const int& /*preshower*/, + [[maybe_unused]] int str_length) { + gCorsikaProcess.Init(); +} + +extern "C" void cloda_() { + // crs::CParticleFortranPtr pptr; + // const bool isF = prminfo_(pptr); + // gCorsikaProcess.Close(); +} + +void interaction_([[maybe_unused]] const crs::CInteraction& interaction) { + /* + all interactions in the shower are available in this function ! + the information availabel in the CInteraction class are: + double x; + double y; + double z; + double etot; // lab energy + double sigma; // cross-section of process + double kela; // elasticity + int projId; // projectile + int targetId; // target + double time; + */ +} + +extern "C" void track_([[maybe_unused]] const crs::CParticle& pre, + [[maybe_unused]] const crs::CParticle& post) { + /* + all particles in the shower are available in this function ! + The pre and post objecte are the two endpoints for one single track + in the shower, where the information available in CParticle is: + double x; + double y; + double z; + double depth; + double time; + double energy; + double weight; + int particleId; + int hadronicGeneration; + */ + gCOASTStack.SetParticle(&pre, &post); + auto particle = gCOASTStack.GetNextParticle(); + + geometry::CoordinateSystem& rootCS = + geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); + + geometry::Line const line(particle.GetPosition(rootCS), particle.GetVelocity(rootCS)); + const TimeType time = particle.GetTimeInterval(); + geometry::Trajectory<geometry::Line> track(line, time); + gCorsikaProcess.DoContinuous(particle, track, gCOASTStack); +} + +extern "C" void tabularizedatmosphere_([[maybe_unused]] const int& nPoints, + [[maybe_unused]] const double* height, + [[maybe_unused]] const double* refractiveIndex) { + // for special use only but should be defined because it is delcared in CORSIKA.F +} diff --git a/COAST/ParticleConversion.cc b/COAST/ParticleConversion.cc new file mode 100644 index 0000000000000000000000000000000000000000..9b1fd3f0e6d75daf45023d2cfb6b86560ce66880 --- /dev/null +++ b/COAST/ParticleConversion.cc @@ -0,0 +1,36 @@ + +/** + * (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. + */ + +#include <corsika/coast/ParticleConversion.h> +#include <corsika/particles/ParticleProperties.h> + +#include <exception> +#include <iostream> +#include <sstream> + +using namespace std; + +namespace corsika::coast { + + corsika::particles::Code ConvertFromCoast(CoastCode pCode) { + if (coast2corsika.count(pCode) == 0) { + ostringstream err; + err << "corsika::coast::ConvertFromCoast CoastCode does not exists=" + << static_cast<CoastCodeIntType>(pCode) << endl; + cout << err.str() << endl; + throw std::runtime_error(err.str()); + } + return coast2corsika.find(pCode)->second; + } + // process::sibyll::Sibyll2Corsika = { + // {PID::E_MINUS, InternalParticleCode::Electron}, + //}; +} // namespace corsika::coast diff --git a/COAST/ParticleConversion.h b/COAST/ParticleConversion.h new file mode 100644 index 0000000000000000000000000000000000000000..d9e87b8e38071051e527cf70640a10aa57c5caf1 --- /dev/null +++ b/COAST/ParticleConversion.h @@ -0,0 +1,181 @@ + +/** + * (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 _include_processes_sibyll_particles_h_ +#define _include_processes_sibyll_particles_h_ + +#include <corsika/particles/ParticleProperties.h> + +#include <bitset2/bitset2.hpp> + +#include <map> + +namespace corsika::coast { + + enum class CoastCode : int16_t { + Gamma = 1, + Positron = 2, + Electron = 3, + MuonBar = 5, + Muon = 6, + Pi0 = 7, + PiP = 8, + PiM = 9, + Klong = 10, + KP = 11, + KM = 12, + Neutron = 13, + Proton = 14, + ProtonBar = 15, + Kshort = 16, + Eta = 17, + Lambda = 18, + SigmaPlus = 19, + Sigma0 = 20, + SigmaMinus = 21, + Xi0 = 22, + XiMinus = 23, + OmegaMinus = 24, + NeutronBar = 25, + LambdaBar = 26, + SigmaMinusBar = 27, + Sigma0Bar = 28, + SigmaPlusBar = 29, + Xi0Bar = 30, + XiPlusBar = 31, + OmegaPlusBar = 32, + + EtaPrime = 48, + Phi = 49, + omega = 50, + Rho0 = 51, + RhoPlus = 52, + RhoMinus = 53, + DeltaPlusPlus = 54, + DeltaPlus = 55, + Delta0 = 56, + DeltaMinus = 57, + DeltaMinusMinusBar = 58, + DeltaMinusBar = 59, + Delta0Bar = 60, + DeltaPlusBar = 61, + KStar0 = 62, + KStarPlus = 63, + KStarMinus = 64, + KStar0Bar = 65, + + NeutrinoE = 66, + NeutrinoEBar = 67, + NeutrinoMu = 68, + NeutrinoMuBar = 69, + + Code71 = 71, + Code72 = 72, + Code73 = 73, + Code74 = 74, + Code75 = 75, + Code76 = 76, + + Code85 = 85, + Code86 = 86, + + Code95 = 95, + Code96 = 96, + + Helium = 402, + Carbon = 612, + Nitrogen = 714, + Oxygen = 816, + Iron = 5628, + }; + using CoastCodeIntType = std::underlying_type<CoastCode>::type; + + const std::map<corsika::coast::CoastCode, corsika::particles::Code> coast2corsika = { + {CoastCode::Gamma, corsika::particles::Code::Gamma}, + {CoastCode::Positron, corsika::particles::Code::Positron}, + {CoastCode::Electron, corsika::particles::Code::Electron}, + //{CoastCode:: ,corsika::particles::Code::Unknown}, // 4 + {CoastCode::MuonBar, corsika::particles::Code::MuPlus}, + {CoastCode::Muon, corsika::particles::Code::MuMinus}, + {CoastCode::Pi0, corsika::particles::Code::Pi0}, + {CoastCode::PiP, corsika::particles::Code::PiPlus}, + {CoastCode::PiM, corsika::particles::Code::PiMinus}, + {CoastCode::Klong, corsika::particles::Code::K0Long}, // 10 + {CoastCode::KP, corsika::particles::Code::KPlus}, + {CoastCode::KM, corsika::particles::Code::KMinus}, + {CoastCode::Neutron, corsika::particles::Code::Neutron}, + {CoastCode::Proton, corsika::particles::Code::Proton}, // 14 + {CoastCode::ProtonBar, corsika::particles::Code::AntiProton}, + {CoastCode::Kshort, corsika::particles::Code::K0Short}, + {CoastCode::Eta, corsika::particles::Code::Eta}, // 17 + {CoastCode::Lambda, corsika::particles::Code::Lambda0}, + {CoastCode::SigmaPlus, corsika::particles::Code::SigmaPlus}, + {CoastCode::Sigma0, corsika::particles::Code::Sigma0}, // 20 + {CoastCode::SigmaMinus, corsika::particles::Code::SigmaMinus}, + {CoastCode::Xi0, corsika::particles::Code::Xi0}, + {CoastCode::XiMinus, corsika::particles::Code::XiMinus}, + {CoastCode::OmegaMinus, corsika::particles::Code::OmegaMinus}, + {CoastCode::NeutronBar, corsika::particles::Code::AntiNeutron}, // 25 + {CoastCode::LambdaBar, corsika::particles::Code::Lambda0Bar}, + {CoastCode::SigmaMinusBar, corsika::particles::Code::SigmaMinusBar}, + {CoastCode::Sigma0Bar, corsika::particles::Code::Sigma0Bar}, + {CoastCode::SigmaPlusBar, corsika::particles::Code::SigmaPlusBar}, + {CoastCode::Xi0Bar, corsika::particles::Code::Xi0Bar}, + {CoastCode::XiPlusBar, corsika::particles::Code::XiPlusBar}, + {CoastCode::OmegaPlusBar, corsika::particles::Code::OmegaPlusBar}, // 32 + //{CoastCode:: ,corsika::particles::Code::Unknown}, // eta-prime + //{CoastCode:: ,corsika::particles::Code::Unknown}, // PHI + //{CoastCode:: ,corsika::particles::Code::Unknown}, // omega + {CoastCode::Rho0, corsika::particles::Code::Rho0}, // 51 + {CoastCode::RhoPlus, corsika::particles::Code::RhoPlus}, + {CoastCode::RhoMinus, corsika::particles::Code::RhoMinus}, + {CoastCode::DeltaPlusPlus, corsika::particles::Code::DeltaPlusPlus}, + {CoastCode::DeltaPlus, corsika::particles::Code::DeltaPlus}, + {CoastCode::Delta0, corsika::particles::Code::Delta0}, // 56 + //{CoastCode:: ,corsika::particles::Code::Unknown}, // DeltaMinus}, + {CoastCode::DeltaMinusMinusBar, corsika::particles::Code::DeltaMinusMinusBar}, + {CoastCode::DeltaMinusBar, corsika::particles::Code::DeltaMinusBar}, + {CoastCode::Delta0Bar, corsika::particles::Code::Delta0Bar}, + //{CoastCode:: ,corsika::particles::Code::Unknown}, // DeltaPlusBar + {CoastCode::KStar0, corsika::particles::Code::KStar0}, // 62 + {CoastCode::KStarPlus, corsika::particles::Code::KStarPlus}, + {CoastCode::KStarMinus, corsika::particles::Code::KStarMinus}, + {CoastCode::KStar0Bar, corsika::particles::Code::KStar0Bar}, + {CoastCode::NeutrinoE, corsika::particles::Code::NuE}, + {CoastCode::NeutrinoEBar, corsika::particles::Code::NuEBar}, + {CoastCode::NeutrinoMu, corsika::particles::Code::NuMu}, + {CoastCode::NeutrinoMuBar, corsika::particles::Code::NuMuBar}, // 69 + {CoastCode::Code71, corsika::particles::Code::Unknown}, + {CoastCode::Code72, corsika::particles::Code::Unknown}, + {CoastCode::Code73, corsika::particles::Code::Unknown}, + {CoastCode::Code74, corsika::particles::Code::Unknown}, + {CoastCode::Code75, corsika::particles::Code::Unknown}, + {CoastCode::Code76, corsika::particles::Code::Unknown}, + {CoastCode::Code85, corsika::particles::Code::Unknown}, + {CoastCode::Code86, corsika::particles::Code::Unknown}, + {CoastCode::Code95, corsika::particles::Code::Unknown}, + {CoastCode::Code96, corsika::particles::Code::Unknown}, + + {CoastCode::Helium, corsika::particles::Code::Helium}, // 402 + {CoastCode::Carbon, corsika::particles::Code::Carbon}, + {CoastCode::Nitrogen, corsika::particles::Code::Nitrogen}, + {CoastCode::Oxygen, corsika::particles::Code::Oxygen}, + {CoastCode::Iron, corsika::particles::Code::Iron}, // 5628 + + //{CoastCode::, corsika::particles::Code::} + + }; + + corsika::particles::Code ConvertFromCoast(CoastCode pCode); + +} // namespace corsika::coast + +#endif diff --git a/COAST/README.md b/COAST/README.md new file mode 100644 index 0000000000000000000000000000000000000000..215d5f10c3e772fedab9d81ff03e8bb325722f30 --- /dev/null +++ b/COAST/README.md @@ -0,0 +1,58 @@ +This is based on corsika7/trunk/coast/CoastOptions/example with an +additional interface to CORSIKA8. + +This is an example for a "COAST user library" using CORSIKA8 +technology. It explains the steps, how to use the COAST_USER_LIB +option of CORSIKA together with CORSIKA8/COAST + + +Step 1: +------- + +Configure CORSIKA8 with 'cmake -DWITH_COAST=1' and other options you +prefer. You have to define COAST_DIR environment variable to the +location of your existing CORSIKA7 installation. Check that 'ls +$COAST_DIR/include/interface' contains the file +'CorsikaInterface.h'. If this is not the case, this is not a valid or +proper installation of CORSIKA. First compile CORSIKA, e.g. with +ROOTOUT option, or CoReas to get COAST installed! + + +Step 2: +------- + +Compile CORSIKA8, edit COAST/COASTProcess.cc to modify/add your +physics module. + +There should be a libCOAST.so as a result! This is what you need. + + +Step 3: +------- + +Create COAST_USER_LIB environment variable to point at your current +directory, where now the 'libCOAST.so' library is located. +Add the path in $COAST_USER_LIB to your LD_LIBRARY_PATH environment +variable. + + +Step 4: +------- + +Go back to your CORSIKA directory and re-start 'coconut'. The option +COAST_USER_LIB will now be visible. Please select it, copile CORSIKA +and start the executable. In the generated console output you will the +statements from your COAST library during Init and Close of the +simulation. +Add any kind of code now working on the CParticle or +CInteraction class to start using the full power of COAST. + +Note: the default COASTProcess just prints out tons of ASCII. This is +not very useful, don't run a full MC with this... + + +Step 5: +------- + +enjoy... + diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 8545e3ab1648191920dda0f8973033675b796e2b..f18ff72efcf8db5c95cc370b8624229d65fccdc0 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -251,7 +251,8 @@ int main() { // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); - const HEPEnergyType E0 = 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later) + const HEPEnergyType E0 = + 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later) double theta = 0.; double phi = 0.; { diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 33914fe527128826155a449452f0c5af0f8e8ac4..cd22f884089092dfdbc373d538168c346b6e35fd 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -73,9 +73,10 @@ class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> { int fCount = 0; int fCalls = 0; HEPEnergyType fEcrit; - + public: - ProcessSplit(HEPEnergyType e) : fEcrit(e) {} + ProcessSplit(HEPEnergyType e) + : fEcrit(e) {} template <typename Particle, typename T> LengthType MaxStepLength(Particle&, T&) const { @@ -83,7 +84,7 @@ public: } template <typename Particle, typename T, typename Stack> - EProcessReturn DoContinuous(Particle& p, T&, Stack& s) { + EProcessReturn DoContinuous(Particle& p, T&, Stack& s) { fCalls++; HEPEnergyType E = p.GetEnergy(); if (E < fEcrit) { @@ -102,7 +103,10 @@ public: return EProcessReturn::eOk; } - void Init() { fCount = 0; fCalls =0; } + void Init() { + fCount = 0; + fCalls = 0; + } int GetCount() const { return fCount; } int GetCalls() const { return fCalls; } @@ -140,8 +144,8 @@ TEST_CASE("Cascade", "[Cascade]") { EAS.Init(); EAS.Run(); - CHECK( p1.GetCount() == 2048 ); - CHECK( p1.GetCalls() == 4095 ); + CHECK(p1.GetCount() == 2048); + CHECK(p1.GetCalls() == 4095); /* SECTION("sectionTwo") { diff --git a/Framework/Geometry/testFourVector.cc b/Framework/Geometry/testFourVector.cc index 301fb97e93d71dccc79fe0eda844376b7ecc133d..6d3d758cbf9a76be337ce8a6927d11f704604eb3 100644 --- a/Framework/Geometry/testFourVector.cc +++ b/Framework/Geometry/testFourVector.cc @@ -63,15 +63,15 @@ TEST_CASE("four vectors") { CHECK(p0.IsSpacelike()); CHECK(!p0.IsTimelike()); - //CHECK(!p0.IsPhotonlike()); + // CHECK(!p0.IsPhotonlike()); CHECK(!p1.IsSpacelike()); CHECK(p1.IsTimelike()); - //CHECK(!p1.IsPhotonlike()); + // CHECK(!p1.IsPhotonlike()); CHECK(!p2.IsSpacelike()); CHECK(!p2.IsTimelike()); - //CHECK(p2.IsPhotonlike()); + // CHECK(p2.IsPhotonlike()); } /* diff --git a/Framework/Particles/NuclearData.xml b/Framework/Particles/NuclearData.xml index d6d2a46dee4c0ac3beb53e8ec5ca5db94bf20680..c707f3a2716a6fbbe61e15df98fcdaaf0e3e71bb 100644 --- a/Framework/Particles/NuclearData.xml +++ b/Framework/Particles/NuclearData.xml @@ -24,4 +24,7 @@ <particle id="1000180400" name="argon" A="40" Z="18" > </particle> +<particle id="1000280560" name="iron" A="56" Z="28" > +</particle> + </chapter> diff --git a/Framework/Random/testRandom.cc b/Framework/Random/testRandom.cc index 97923eb09e19a7e6bc1cd8a4c8b6971907af9528..a8f8595512815574fd3abb0b9c29274a37ac22d4 100644 --- a/Framework/Random/testRandom.cc +++ b/Framework/Random/testRandom.cc @@ -12,9 +12,9 @@ // cpp file #include <catch2/catch.hpp> +#include <corsika/random/ExponentialDistribution.h> #include <corsika/random/RNGManager.h> #include <corsika/random/UniformRealDistribution.h> -#include <corsika/random/ExponentialDistribution.h> #include <corsika/units/PhysicalUnits.h> #include <iostream> #include <limits> @@ -93,17 +93,17 @@ TEST_CASE("ExponentialDistribution") { auto const beta = 15_m; corsika::random::ExponentialDistribution dist(beta); - + SECTION("mean") { std::remove_const<decltype(beta)>::type mean = beta * 0; - + int constexpr N = 1'000'000; - + for (int i{0}; i < N; ++i) { decltype(beta) x = dist(rng); mean += x / N; } - + CHECK(mean / beta == Approx(1).margin(1e-2)); } } diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc index e7c877c58cc550146b8ce0e44db17be50d8d75b4..a9bd4115362de0685dd8bbcfd58cda9ff3064a2b 100644 --- a/Framework/Utilities/COMBoost.cc +++ b/Framework/Utilities/COMBoost.cc @@ -17,7 +17,8 @@ using namespace corsika::utl; using namespace corsika::units::si; template <typename FourVector> -COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const HEPMassType massTarget) +COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, + const HEPMassType massTarget) : fRotation(Eigen::Matrix3d::Identity()) , fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) { // calculate matrix for rotating pProjectile to z-axis first @@ -44,8 +45,7 @@ COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const HEPMassType } // calculate boost - double const beta = - pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget); + double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget); /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */ double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta)); @@ -83,8 +83,8 @@ FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const { (p.GetSpaceLikeComponents().GetComponents().eVector(2) * (1 / 1_GeV).magnitude()); std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << " GeV, " - << " pcm=" << p.GetSpaceLikeComponents().GetComponents().squaredNorm() / 1_GeV << " GeV" - << std::endl; + << " pcm=" << p.GetSpaceLikeComponents().GetComponents().squaredNorm() / 1_GeV + << " GeV" << std::endl; auto const boostedZ = fInverseBoost * com; auto const E_lab = boostedZ(0) * 1_GeV; diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h index 46e808c71aeb57d384fe09fc65d8ce898904e28b..e48de10ba1ea87f76a3a5d715495061c35c2efeb 100644 --- a/Framework/Utilities/COMBoost.h +++ b/Framework/Utilities/COMBoost.h @@ -31,7 +31,8 @@ namespace corsika::utl { public: //! construct a COMBoost given four-vector of prjectile and mass of target - COMBoost(const FourVector& Pprojectile, const corsika::units::si::HEPEnergyType massTarget); + COMBoost(const FourVector& Pprojectile, + const corsika::units::si::HEPEnergyType massTarget); //! transforms a 4-momentum from lab frame to the center-of-mass frame FourVector toCoM(const FourVector& p) const; diff --git a/Framework/Utilities/CorsikaFenv.h b/Framework/Utilities/CorsikaFenv.h index e1977d6d2ec5178ae518538e2ef07dc0d74ac6ea..d8aa00b3f0987973c71af606d2d5dbd93304dfe4 100644 --- a/Framework/Utilities/CorsikaFenv.h +++ b/Framework/Utilities/CorsikaFenv.h @@ -27,11 +27,8 @@ */ extern "C" { - int - feenableexcept(int excepts); - int - fedisableexcept(int excepts); - +int feenableexcept(int excepts); +int fedisableexcept(int excepts); } -#endif //CORSIKA_CORSIKAFENV_H +#endif // CORSIKA_CORSIKAFENV_H diff --git a/Framework/Utilities/CorsikaFenvFallback.cc b/Framework/Utilities/CorsikaFenvFallback.cc index ad495c4d3b0bc7396490fbcededc1d44c767806d..f3e2f00f26840838365b45fa7426930a8b818e6e 100644 --- a/Framework/Utilities/CorsikaFenvFallback.cc +++ b/Framework/Utilities/CorsikaFenvFallback.cc @@ -24,14 +24,7 @@ extern "C" { #warning No enabling/disabling of floating point exceptions - platform needs better implementation - int feenableexcept(int excepts) - { - return -1; - } - - int fedisableexcept(int excepts) - { - return -1; - } +int feenableexcept(int excepts) { return -1; } +int fedisableexcept(int excepts) { return -1; } } diff --git a/Framework/Utilities/CorsikaFenvOSX.cc b/Framework/Utilities/CorsikaFenvOSX.cc index 224e04e733da6cf760e91fe5e94873afe7848898..1028c4b7e0988bb7ced605f71990c9c9f344fbbc 100644 --- a/Framework/Utilities/CorsikaFenvOSX.cc +++ b/Framework/Utilities/CorsikaFenvOSX.cc @@ -14,47 +14,42 @@ #include <cfenv> // Implementation for OS X on intel X64_86 -// code from https://stackoverflow.com/questions/37819235/how-do-you-enable-floating-point-exceptions-for-clang-in-os-x -// based on http://www-personal.umich.edu/~williams/archive/computation/fe-handling-example.c +// code from +// https://stackoverflow.com/questions/37819235/how-do-you-enable-floating-point-exceptions-for-clang-in-os-x +// based on +// http://www-personal.umich.edu/~williams/archive/computation/fe-handling-example.c extern "C" { - int - feenableexcept(int excepts) { - static fenv_t fenv; - int new_excepts = excepts & FE_ALL_EXCEPT; - // previous masks - int old_excepts; +int feenableexcept(int excepts) { + static fenv_t fenv; + int new_excepts = excepts & FE_ALL_EXCEPT; + // previous masks + int old_excepts; - if (fegetenv(&fenv)) { - return -1; - } - old_excepts = fenv.__control & FE_ALL_EXCEPT; + if (fegetenv(&fenv)) { return -1; } + old_excepts = fenv.__control & FE_ALL_EXCEPT; - // unmask - fenv.__control &= ~new_excepts; - fenv.__mxcsr &= ~(new_excepts << 7); + // unmask + fenv.__control &= ~new_excepts; + fenv.__mxcsr &= ~(new_excepts << 7); - return fesetenv(&fenv) ? -1 : old_excepts; - } - - int - fedisableexcept(int excepts) { - static fenv_t fenv; - int new_excepts = excepts & FE_ALL_EXCEPT; - // all previous masks - int old_excepts; + return fesetenv(&fenv) ? -1 : old_excepts; +} - if (fegetenv(&fenv)) { - return -1; - } - old_excepts = fenv.__control & FE_ALL_EXCEPT; +int fedisableexcept(int excepts) { + static fenv_t fenv; + int new_excepts = excepts & FE_ALL_EXCEPT; + // all previous masks + int old_excepts; - // mask - fenv.__control |= new_excepts; - fenv.__mxcsr |= new_excepts << 7; + if (fegetenv(&fenv)) { return -1; } + old_excepts = fenv.__control & FE_ALL_EXCEPT; - return fesetenv(&fenv) ? -1 : old_excepts; - } + // mask + fenv.__control |= new_excepts; + fenv.__mxcsr |= new_excepts << 7; + return fesetenv(&fenv) ? -1 : old_excepts; +} } diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc index 712c4996c89a986ac3449c5b4d24ccbaf38e316c..996da626138ebdd0ec7b989672f6161638729173 100644 --- a/Framework/Utilities/testCOMBoost.cc +++ b/Framework/Utilities/testCOMBoost.cc @@ -36,7 +36,7 @@ TEST_CASE("boosts") { auto energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) { return sqrt(m * m + p.squaredNorm()); }; - + // helper function for mandelstam-s auto s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) { return E * E - p.squaredNorm(); @@ -50,7 +50,7 @@ TEST_CASE("boosts") { /* General tests check the interface and basic operation */ - + SECTION("General tests") { // define projectile kinematics in lab frame @@ -124,7 +124,7 @@ TEST_CASE("boosts") { /* special case: projectile with arbitrary direction */ - + SECTION("Test boost along tilted axis") { const HEPMomentumType P0 = 1_PeV; @@ -161,7 +161,7 @@ TEST_CASE("boosts") { /* test the ultra-high energy behaviour: E=ZeV */ - + SECTION("High energy") { // define projectile kinematics in lab frame HEPMassType const projectileMass = 1_GeV; diff --git a/Framework/Utilities/testCorsikaFenv.cc b/Framework/Utilities/testCorsikaFenv.cc index 9cae946918047166c347e73b907a5ea562a61479..6f76930c86c7d3ee22f5a155700af35de7959098 100644 --- a/Framework/Utilities/testCorsikaFenv.cc +++ b/Framework/Utilities/testCorsikaFenv.cc @@ -11,19 +11,14 @@ #include <corsika/utl/CorsikaFenv.h> #include <cmath> -#include <iostream> #include <csignal> +#include <iostream> extern "C" { - static void - handle_fpe(int /*signo*/ ) { - exit(0); - } +static void handle_fpe(int /*signo*/) { exit(0); } } - -int -main() { +int main() { feenableexcept(FE_ALL_EXCEPT); signal(SIGFPE, handle_fpe); diff --git a/Framework/Utilities/try_feenableexcept.cc b/Framework/Utilities/try_feenableexcept.cc index 270307b407a85881bcbc273c67c44e8d4e2dad60..1ee8ecc1eeb67c267578b74388d3054d12c16872 100644 --- a/Framework/Utilities/try_feenableexcept.cc +++ b/Framework/Utilities/try_feenableexcept.cc @@ -16,9 +16,7 @@ #include <cfenv> -int -main() -{ +int main() { feenableexcept(FE_ALL_EXCEPT); return 0; } \ No newline at end of file