IAP GITLAB

Commit 4159e9bc authored by Ralf Ulrich's avatar Ralf Ulrich Committed by Felix Riehn

Resolve "Start process interface for Pythia8 (decays)"

parent cd8e728b
......@@ -72,6 +72,8 @@ if (Boost_FOUND)
include_directories(${Boost_INCLUDE_DIRS})
endif (Boost_FOUND)
find_package (Pythia8) # optional
find_package (Eigen3 REQUIRED)
#find_package (HDF5) # not yet needed
......
#
# (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.
#
#################################################
#
# takes a list of input files and prepends a path
#
......@@ -12,6 +22,7 @@ function (CORSIKA_PREPEND_PATH return prefix)
endfunction (CORSIKA_PREPEND_PATH)
#################################################
#
# use: CORSIKA_COPY_HEADERS_TO_NAMESPACE theLib theNamesapce header1.h header2.h ...
#
......@@ -58,6 +69,7 @@ endfunction (CORSIKA_COPY_HEADERS_TO_NAMESPACE)
#################################################
#
# use: CORSIKA_ADD_FILES_ABSOLUTE varname
#
......@@ -81,6 +93,7 @@ endmacro(CORSIKA_ADD_FILES_ABSOLUTE)
#################################################
#
# central macro to activate unit tests in cmake
#
......
#
# (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.
#
#################################################
#
# run pythia8-config and interpret result
#
function (_PYTHIA8_CONFIG_ option variable type doc)
execute_process(COMMAND ${PYTHIA8_CONFIG} ${option}
OUTPUT_VARIABLE _local_out_
RESULT_VARIABLE _local_res_)
string(REGEX REPLACE "\n$" "" _local_out_ "${_local_out_}")
if (NOT ${_local_res_} EQUAL 0)
message ("Error in running ${PYTHIA8_CONFIG} ${option}")
else ()
set(${variable} "${_local_out_}" CACHE ${type} ${doc})
endif ()
endfunction (_PYTHIA8_CONFIG_)
#################################################
#
# Searched PYTHIA8 on system. Expect pythia8-config in PATH, or typical installation location
#
# This module defines
# HAVE_PYTHIA8
# PYTHIA8_INCLUDE_DIR where to locate Pythia.h file
# PYTHIA8_LIBRARY where to find the libpythia8 library
# PYTHIA8_LIBRARIES (not cached) the libraries to link against to use Pythia8
# PYTHIA8_VERSION version of Pythia8 if found
#
set(_SEARCH_PYTHIA8_
${PROJECT_BINARY_DIR}/ThirdParty/pythia8-install
${PYTHIA8}
$ENV{PYTHIA8}
${PYTHIA8DIR}
$ENV{PYTHIA8DIR}
${PYTHIA8_DIR}
$ENV{PYTHIA8_DIR}
${PYTHIA8_ROOT}
$ENV{PYTHIA8_ROOT}
/opt/pythia8)
find_program(PYTHIA8_CONFIG
NAME pythia8-config
PATHS ${_SEARCH_PYTHIA8_}
PATH_SUFFIXES "/bin"
DOC "The location of the pythia8-config script")
if (PYTHIA8_CONFIG)
set(HAVE_PYTHIA8 1 CACHE BOOL "presence of pythia8, found via pythia8-config")
_PYTHIA8_CONFIG_("--prefix" PYTHIA8_PREFIX PATH "location of pythia8 installation")
_PYTHIA8_CONFIG_("--includedir" PYTHIA8_INCLUDE_DIR PATH "pythia8 include directory")
_PYTHIA8_CONFIG_("--libs" PYTHIA8_LIBRARY STRING "the pythia8 libs")
_PYTHIA8_CONFIG_("--datadir" PYTHIA8_DATA_DIR PATH "the pythia8 data dir")
_PYTHIA8_CONFIG_("--cxxflags" PYTHIA8_CXXFLAGS STRING "the pythia8 cxxflags")
endif ()
# standard cmake infrastructure:
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(Pythia8 DEFAULT_MSG PYTHIA8_PREFIX PYTHIA8_INCLUDE_DIR PYTHIA8_LIBRARY)
mark_as_advanced(PYTHIA8_DATA_DIR PYTHIA8_CXXFLAGS PYTHIA8_INCLUDE_DIR PYTHIA8_LIBRARY)
......@@ -18,14 +18,17 @@ CORSIKA_ADD_TEST (logger_example)
add_executable (stack_example stack_example.cc)
target_compile_options(stack_example PRIVATE -g) # do not skip asserts
target_link_libraries (stack_example SuperStupidStack CORSIKAunits CORSIKAlogging)
target_link_libraries (stack_example SuperStupidStack CORSIKAunits
CORSIKAlogging)
CORSIKA_ADD_TEST (stack_example)
add_executable (cascade_example cascade_example.cc)
target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
CORSIKArandom
ProcessSibyll
target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia
CORSIKAcascade
ProcessStackInspector
ProcessTrackWriter
......
......@@ -28,6 +28,8 @@
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
......@@ -251,10 +253,17 @@ int main() {
tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll(env);
process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
process::sibyll::Decay decay;
process::sibyll::Decay decay(trackedHadrons);
//random::RNGManager::GetInstance().RegisterRandomStream("pythia");
//process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
......@@ -274,7 +283,7 @@ int main() {
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Nucleus;
const int nuclA = 56;
const int nuclA = 4;
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) * particles::Neutron::GetMass();
......
......@@ -47,10 +47,10 @@ namespace corsika::geometry {
void LimitEndTo(corsika::units::si::LengthType limit) {
fTimeLength = T::TimeFromArclength(limit);
}
auto NormalizedDirection() const {
static_assert(std::is_same_v<T, corsika::geometry::Line>);
return T::GetV0().normalized();
static_assert(std::is_same_v<T, corsika::geometry::Line>);
return T::GetV0().normalized();
}
};
......
......@@ -177,8 +177,10 @@ TEST_CASE("Trajectories") {
CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
CHECK((base.NormalizedDirection().GetComponents(rootCS) - QuantityVector<dimensionless_d>{1, 0, 0}).eVector.norm() == Approx(0).margin(absMargin));
CHECK((base.NormalizedDirection().GetComponents(rootCS) -
QuantityVector<dimensionless_d>{1, 0, 0})
.eVector.norm() == Approx(0).margin(absMargin));
}
SECTION("Helix") {
......
add_subdirectory (NullModel)
add_subdirectory (Sibyll)
if (PYTHIA8_FOUND)
add_subdirectory (Pythia)
endif (PYTHIA8_FOUND)
add_subdirectory (StackInspector)
add_subdirectory (TrackWriter)
add_subdirectory (TrackingLine)
......@@ -10,5 +13,8 @@ add_subdirectory (HadronicElasticModel)
add_library (CORSIKAprocesses INTERFACE)
add_dependencies(CORSIKAprocesses ProcessNullModel)
add_dependencies(CORSIKAprocesses ProcessSibyll)
if (PYTHIA8_FOUND)
add_dependencies(CORSIKAprocesses ProcessPythia)
endif (PYTHIA8_FOUND)
add_dependencies(CORSIKAprocesses ProcessStackInspector)
add_dependencies(CORSIKAprocesses ProcessTrackingLine)
set(Python_ADDITIONAL_VERSIONS 3)
find_package(PythonInterp 3 REQUIRED)
set (
MODEL_SOURCES
Decay.cc
Random.cc
)
set (
MODEL_HEADERS
Decay.h
Random.h
)
set (
MODEL_NAMESPACE
corsika/process/pythia
)
add_library (ProcessPythia STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessPythia ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessPythia
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessPythia
CORSIKAparticles
CORSIKAutilities
CORSIKAunits
CORSIKAthirdparty
CORSIKAgeometry
CORSIKAenvironment
${PYTHIA8_LIBRARY}
)
target_include_directories (
ProcessPythia
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
# do it again, since we want -isystem for pythia exteranls
target_include_directories (
ProcessPythia
SYSTEM PUBLIC
${PYTHIA8_INCLUDE_DIR}
)
target_include_directories (
ProcessPythia
PRIVATE
${PYTHIA8_INCLUDE_DIR}
)
install (
TARGETS ProcessPythia
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
)
# --------------------
# code unit testing
add_executable (testPythia
testPythia.cc
${MODEL_HEADERS}
)
target_link_libraries (
testPythia
ProcessPythia
CORSIKAsetup
CORSIKArandom
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
${PYTHIA8_LIBRARY}
)
CORSIKA_ADD_TEST(testPythia)
/*
* (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/process/pythia/Decay.h>
#include <corsika/process/pythia/Random.h>
#include <Pythia8/Pythia.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using std::cout;
using std::endl;
using std::tuple;
using std::vector;
using namespace corsika;
using namespace corsika::setup;
using Particle = Stack::StackIterator; // ParticleType;
using Track = Trajectory;
namespace corsika::process::pythia {
Decay::Decay(vector<particles::Code> pParticles)
: fTrackedParticles(pParticles) {}
Decay::~Decay() { cout << "Pythia::Decay n=" << fCount << endl; }
void Decay::Init() {
Decay::SetParticleListStable(fTrackedParticles);
// set random number generator in pythia
Pythia8::RndmEngine* rndm = new corsika::process::pythia::Random();
fPythia.setRndmEnginePtr( rndm );
fPythia.readString("Next:numberShowInfo = 0");
fPythia.readString("Next:numberShowProcess = 0");
fPythia.readString("Next:numberShowEvent = 0");
fPythia.readString("Print:quiet = on");
fPythia.readString("ProcessLevel:all = off");
fPythia.readString("ProcessLevel:resonanceDecays = off");
fPythia.particleData.readString("59:m0 = 101.00");
fPythia.init();
}
void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
for (auto p : particleList)
Decay::SetStable( p );
}
void Decay::SetUnstable(const particles::Code pCode) {
cout << "Pythia::Decay: setting " << pCode << " unstable.." << endl;
fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , true);
}
void Decay::SetStable(const particles::Code pCode) {
cout << "Pythia::Decay: setting " << pCode << " stable.." << endl;
fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , false);
}
template <>
units::si::TimeType Decay::GetLifetime(Particle const& p) {
using namespace units::si;
HEPEnergyType E = p.GetEnergy();
HEPMassType m = particles::GetMass(p.GetPID());
const double gamma = E / m;
const TimeType t0 = particles::GetLifetime(p.GetPID());
auto const lifetime = gamma * t0;
return lifetime;
}
template <>
void Decay::DoDecay(Particle& p, Stack&) {
using geometry::Point;
using namespace units;
using namespace units::si;
auto const decayPoint = p.GetPosition();
auto const t0 = p.GetTime();
// coordinate system, get global frame of reference
geometry::CoordinateSystem& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
fCount++;
// pythia stack
Pythia8::Event& event = fPythia.event;
event.reset();
// set particle unstable
Decay::SetUnstable( p.GetPID() );
// input particle PDG
auto const pdgCode = static_cast<int>( particles::GetPDG( p.GetPID() ));
auto const pcomp = p.GetMomentum().GetComponents();
double px = pcomp[0] / 1_GeV ;
double py = pcomp[1] / 1_GeV ;
double pz = pcomp[2] / 1_GeV ;
double en = p.GetEnergy() / 1_GeV;
double m = particles::GetMass( p.GetPID() ) / 1_GeV;
// add particle to pythia stack
event.append( pdgCode, 1, 0, 0, px, py, pz, en, m);
if (!fPythia.next())
cout << "Pythia::Decay: decay failed!" << endl;
else
cout << "Pythia::Decay: particles after decay: " << event.size() << endl;
// list final state
event.list();
// loop over final state
for (int i = 0; i < event.size(); ++i)
if (event[i].isFinal()) {
auto const pyId = particles::ConvertFromPDG(static_cast<particles::PDGCode>(event[i].id()));
HEPEnergyType pyEn = event[i].e() * 1_GeV;
MomentumVector pyP(rootCS, { event[i].px() * 1_GeV,
event[i].py() * 1_GeV,
event[i].pz() * 1_GeV});
cout << "particle: id=" << pyId << " momentum=" << pyP.GetComponents() / 1_GeV
<< " energy=" << pyEn << endl;
p.AddSecondary( tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
pyId, pyEn, pyP, decayPoint, t0} );
}
// set particle stable
Decay::SetStable( p.GetPID() );
// remove original particle from corsika stack
p.Delete();
// if (fCount>10) throw std::runtime_error("stop here");
}
} // namespace corsika::process::pythia
/*
* (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_corsika_process_pythia_decay_h_
#define _include_corsika_process_pythia_decay_h_
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/DecayProcess.h>
#include <Pythia8/Pythia.h>
namespace corsika::process {
namespace pythia {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
class Decay : public corsika::process::DecayProcess<Decay> {
const std::vector<particles::Code> fTrackedParticles;
int fCount = 0;
public:
Decay(std::vector<corsika::particles::Code> );
~Decay();
void Init();
void SetParticleListStable(const std::vector<particles::Code>);
void SetUnstable(const corsika::particles::Code );
void SetStable(const corsika::particles::Code );
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle const& p);
template <typename Particle, typename Stack>
void DoDecay(Particle& p, Stack&);
private:
Pythia8::Pythia fPythia;
};
} // namespace pythia
} // namespace corsika::process
#endif
/*
* (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/process/pythia/Random.h>
namespace corsika::process::pythia {
double Random::flat(){
std::uniform_real_distribution<double> dist;
return dist(fRNG);
}
} // namespace corsika::process::pythia
/*
* (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_corsika_process_pythia_random_h_
#define _include_corsika_process_pythia_random_h_
#include <corsika/random/RNGManager.h>
#include <Pythia8/Pythia.h>
namespace corsika::process {
namespace pythia {
class Random : public Pythia8::RndmEngine {
double flat();
private:
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
};
} // namespace pythia
} // namespace corsika::process
#endif
/*
* (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/process/pythia/Decay.h>
#include <Pythia8/Pythia.h>
#include <corsika/random/RNGManager.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
TEST_CASE("Pythia", "[processes]") {
SECTION("linking pythia") {
using namespace Pythia8;
using std::cout;
using std::endl;
// Generator. Process selection. LHC initialization. Histogram.
Pythia pythia;
pythia.readString("Next:numberShowInfo = 0");
pythia.readString("Next:numberShowProcess = 0");
pythia.readString("Next:numberShowEvent = 0");
pythia.readString("ProcessLevel:all = off");
pythia.init();
Event& event = pythia.event;
event.reset();
pythia.particleData.mayDecay(321, true);
double pz = 100.;
double m = 0.49368;
event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
if (!pythia.next())
cout << "decay failed!" << endl;
else
cout << "particles after decay: " << event.size() << endl;
event.list();
// loop over final state
for (int i = 0; i < pythia.event.size(); ++i)
if (pythia.event[i].isFinal()) {
cout << "particle: id=" << pythia.event[i].id() << endl;
}
}