IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 1601 deletions
Processes/AnalyticProcessors/ImplDecay.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/DecayProcess.h>
#include <corsika/process/analytic_processors/ExecTimeImpl.h>
namespace corsika::process {
namespace analytic_processors {
namespace detail {
template <typename T>
class ExecTimeImpl;
}
/// Base for Decay Implementation
template <class T, bool TCheck>
class Decay;
/// Specialisation if class is not DecayProcess
template <class T>
class Decay<T, false> {};
/// Specialisation if class is a DecayProcess
template <class T>
class Decay<T, true> : public detail::ExecTimeImpl<T> {
private:
public:
template <typename Particle>
EProcessReturn DoDecay(Particle& p) {
this->start();
auto r = detail::ExecTimeImpl<T>::DoDecay(p);
this->stop();
return r;
}
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle& p) {
this->start();
auto r = T::GetLifetime(p);
this->stop();
return r;
}
};
} // namespace analytic_processors
} // namespace corsika::process
\ No newline at end of file
Processes/AnalyticProcessors/ImplInteraction.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/analytic_processors/ExecTimeImpl.h>
namespace corsika::process {
namespace analytic_processors {
namespace detail {
template <typename T>
class ExecTimeImpl;
}
template <class T, bool TCheck>
class Interaction;
template <class T>
class Interaction<T, false> {};
template <class T>
class Interaction<T, true> : public detail::ExecTimeImpl<T> {
private:
public:
template <typename Particle>
EProcessReturn DoInteraction(Particle& p) {
this->start();
auto r = detail::ExecTimeImpl<T>::DoInteraction(p);
this->stop();
return r;
}
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle& p) {
this->start();
auto r = T::GetInteractionLength(p);
this->stop();
return r;
}
};
} // namespace analytic_processors
} // namespace corsika::process
\ No newline at end of file
Processes/AnalyticProcessors/ImplSecondaries.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/analytic_processors/ExecTimeImpl.h>
namespace corsika::process {
namespace analytic_processors {
namespace detail {
template <typename T>
class ExecTimeImpl;
}
template <class T, bool TCheck>
class Secondaries;
template <class T>
class Secondaries<T, false> {};
template <class T>
class Secondaries<T, true> : public detail::ExecTimeImpl<T> {
private:
public:
template <typename Secondaries>
inline EProcessReturn DoSecondaries(Secondaries& sec) {
this->start();
auto r = detail::ExecTimeImpl<T>::DoSecondaries(sec);
this->stop();
return r;
}
};
} // namespace analytic_processors
} // namespace corsika::process
\ No newline at end of file
Processes/CONEXSourceCut/CMakeLists.txt deleted 100644 → 0
View file @ eabed243
set (
MODEL_SOURCES
CONEXSourceCut.cc
)
set (
MODEL_HEADERS
CONEXSourceCut.h
CONEX_f.h
)
set (
MODEL_NAMESPACE
corsika/process/conex_source_cut
)
add_library (ProcessCONEXSourceCut STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessCONEXSourceCut ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessCONEXSourceCut
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessCONEXSourceCut
CORSIKAunits
CORSIKAparticles
CORSIKAprocesssequence
CORSIKAsetup
C8::ext::conex
ProcessSibyll
ProcessUrQMD
)
target_include_directories (
ProcessCONEXSourceCut
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessCONEXSourceCut
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testCONEXSourceCut SOURCES
testCONEXSourceCut.cc
${MODEL_HEADERS}
)
#
target_compile_definitions (
testCONEXSourceCut
PRIVATE
REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}"
)
#
target_link_libraries (
testCONEXSourceCut
ProcessCONEXSourceCut
CORSIKAunits
CORSIKAstackinterface
CORSIKAprocesssequence
CORSIKAsetup
CORSIKAgeometry
CORSIKAenvironment
CORSIKAtesting
)
Processes/CONEXSourceCut/CONEXSourceCut.cc deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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/logging/Logging.h>
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/conex_source_cut/CONEX_f.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalConstants.h>
#include <algorithm>
#include <fstream>
#include <iomanip>
#include <utility>
using namespace corsika::process::conex_source_cut;
using namespace corsika::units::si;
using namespace corsika::particles;
using namespace corsika::setup;
corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
corsika::setup::StackView& vS) {
auto p = vS.begin();
while (p != vS.end()) {
Code const pid = p.GetPID();
if (addParticle(pid, p.GetEnergy(), p.GetMass(), p.GetPosition(),
p.GetMomentum().normalized(), p.GetTime())) {
p.Delete();
}
++p;
}
return corsika::process::EProcessReturn::eOk;
}
bool CONEXSourceCut::addParticle(particles::Code pid, HEPEnergyType energy,
HEPEnergyType mass, geometry::Point const& position,
geometry::Vector<dimensionless_d> const& direction,
TimeType t) {
auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
[=](auto const& p) { return pid == p.first; });
if (it == egs_em_codes_.cend()) { return false; }
// EM particle
auto const egs_pid = it->second;
std::cout << "position conexObs: " << position.GetCoordinates(conexObservationCS_)
<< std::endl;
auto const coords = position.GetCoordinates(conexObservationCS_) / 1_m;
double const x = coords[0].magnitude();
double const y = coords[1].magnitude();
double const altitude = ((position - center_).norm() - conex::earthRadius) / 1_m;
auto const d = position - showerCore_;
// distance from core to particle projected along shower axis
double const slantDistance = -d.dot(showerAxis_.GetDirection()) / 1_m;
// lateral coordinates in CONEX shower frame
auto const dShowerPlane = d - d.parallelProjectionOnto(showerAxis_.GetDirection());
double const lateralX = dShowerPlane.dot(x_sf_) / 1_m;
double const lateralY = dShowerPlane.dot(y_sf_) / 1_m;
double const slantX = showerAxis_.projectedX(position) * (1_cm * 1_cm / 1_g);
double const time = (t * units::constants::c - groundDist_) / 1_m;
// fill u,v,w momentum direction in EGS frame
double const u = direction.dot(y_sf_).magnitude();
double const v = direction.dot(x_sf_).magnitude();
double const w = direction.dot(showerAxis_.GetDirection()).magnitude();
double const weight = 1; // NEEDS TO BE CHANGED WHEN WE HAVE WEIGHTS!
// generation, TO BE CHANGED WHEN WE HAVE THAT INFORMATION AVAILABLE
int const latchin = 1;
double const E = energy / 1_GeV;
double const m = mass / 1_GeV;
std::cout << "CONEXSourceCut: removing " << egs_pid << " " << std::scientific << energy
<< " GeV" << std::endl;
std::cout << "#### parameters to cegs4_() ####" << std::endl;
std::cout << "egs_pid = " << egs_pid << std::endl;
std::cout << "E = " << E << std::endl;
std::cout << "m = " << m << std::endl;
std::cout << "x = " << x << std::endl;
std::cout << "y = " << y << std::endl;
std::cout << "altitude = " << altitude << std::endl;
std::cout << "slantDistance = " << slantDistance << std::endl;
std::cout << "lateralX = " << lateralX << std::endl;
std::cout << "lateralY = " << lateralY << std::endl;
std::cout << "slantX = " << slantX << std::endl;
std::cout << "time = " << time << std::endl;
std::cout << "u = " << u << std::endl;
std::cout << "v = " << v << std::endl;
std::cout << "w = " << w << std::endl;
conex::cxoptl_.dptl[10 - 1] = egs_pid;
conex::cxoptl_.dptl[4 - 1] = E;
conex::cxoptl_.dptl[5 - 1] = m;
conex::cxoptl_.dptl[6 - 1] = x;
conex::cxoptl_.dptl[7 - 1] = y;
conex::cxoptl_.dptl[8 - 1] = altitude;
conex::cxoptl_.dptl[9 - 1] = time;
conex::cxoptl_.dptl[11 - 1] = weight;
conex::cxoptl_.dptl[12 - 1] = latchin;
conex::cxoptl_.dptl[13 - 1] = slantX;
conex::cxoptl_.dptl[14 - 1] = lateralX;
conex::cxoptl_.dptl[15 - 1] = lateralY;
conex::cxoptl_.dptl[16 - 1] = slantDistance;
conex::cxoptl_.dptl[2 - 1] = u;
conex::cxoptl_.dptl[1 - 1] = v;
conex::cxoptl_.dptl[3 - 1] = w;
int n = 1, i = 1;
conex::cegs4_(n, i);
return true;
}
void CONEXSourceCut::SolveCE() {
conex::conexcascade_();
int nX = conex::get_number_of_depth_bins_(); // make sure this works!
int icut = 1;
int icutg = 2;
int icute = 3;
int icutm = 2;
int icuth = 3;
int iSec = 0;
const int maxX = nX;
auto X = std::make_unique<float[]>(maxX);
auto H = std::make_unique<float[]>(maxX);
auto D = std::make_unique<float[]>(maxX);
auto N = std::make_unique<float[]>(maxX);
auto dEdX = std::make_unique<float[]>(maxX);
auto Mu = std::make_unique<float[]>(maxX);
auto dMu = std::make_unique<float[]>(maxX);
auto Gamma = std::make_unique<float[]>(maxX);
auto Electrons = std::make_unique<float[]>(maxX);
auto Hadrons = std::make_unique<float[]>(maxX);
float EGround[3], fitpars[13];
conex::get_shower_data_(icut, iSec, nX, X[0], N[0], fitpars[0], H[0], D[0]);
conex::get_shower_edep_(icut, nX, dEdX[0], EGround[0]);
conex::get_shower_muon_(icutm, nX, Mu[0], dMu[0]);
conex::get_shower_gamma_(icutg, nX, Gamma[0]);
conex::get_shower_electron_(icute, nX, Electrons[0]);
conex::get_shower_hadron_(icuth, nX, Hadrons[0]);
std::ofstream file{"conex_output.txt"};
file << fmt::format("#{:>8} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13}\n", "X",
"N", "dEdX", "Mu", "dMu", "Gamma", "El", "Had");
for (int i = 0; i < nX; ++i) {
file << fmt::format(
" {:>8.2f} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3}\n",
X[i], N[i], dEdX[i], Mu[i], dMu[i], Gamma[i], Electrons[i], Hadrons[i]);
}
file.close();
std::ofstream fitout{"conex_fit.txt"};
fitout << fitpars[1 - 1] << " # log10(eprima/eV)" << std::endl;
fitout << fitpars[2 - 1] << " # theta" << std::endl;
fitout << fitpars[3 - 1] << " # X1 (first interaction)" << std::endl;
fitout << fitpars[4 - 1] << " # Nmax" << std::endl;
fitout << fitpars[5 - 1] << " # X0" << std::endl;
fitout << fitpars[6 - 1] << " # P1" << std::endl;
fitout << fitpars[7 - 1] << " # P2" << std::endl;
fitout << fitpars[8 - 1] << " # P3" << std::endl;
fitout << fitpars[9 - 1] << " # chi^2 / sqrt(Nmax)" << std::endl;
fitout << fitpars[10 - 1] << " # Xmax" << std::endl;
fitout << fitpars[11 - 1] << " # phi" << std::endl;
fitout << fitpars[12 - 1] << " # inelasticity 1st int." << std::endl;
fitout << fitpars[13 - 1] << " # ???" << std::endl;
fitout.close();
}
CONEXSourceCut::CONEXSourceCut(geometry::Point center,
environment::ShowerAxis const& showerAxis,
units::si::LengthType groundDist,
units::si::LengthType injectionHeight,
units::si::HEPEnergyType primaryEnergy,
particles::PDGCode primaryPDG)
: center_{center}
, showerAxis_{showerAxis}
, groundDist_{groundDist}
, showerCore_{showerAxis_.GetStart() + showerAxis_.GetDirection() * groundDist_}
, conexObservationCS_{std::invoke([&]() {
auto const& c8cs = center.GetCoordinateSystem();
auto const translation = showerCore_ - center;
auto const intermediateCS = c8cs.translate(translation.GetComponents(c8cs));
auto const intermediateCS2 = intermediateCS.RotateToZ(translation);
std::cout << "translation C8/CONEX obs: " << translation.GetComponents()
<< std::endl;
auto const transform = geometry::CoordinateSystem::GetTransformation(
intermediateCS2, c8cs); // either this way or vice versa... TODO: test this!
std::cout << transform.matrix() << std::endl << std::endl;
std::cout
<< geometry::CoordinateSystem::GetTransformation(intermediateCS, c8cs).matrix()
<< std::endl
<< std::endl;
std::cout << geometry::CoordinateSystem::GetTransformation(intermediateCS2,
intermediateCS)
.matrix()
<< std::endl;
return geometry::CoordinateSystem(c8cs, transform);
})}
, x_sf_{std::invoke([&]() {
geometry::Vector<length_d> const a{conexObservationCS_, 0._m, 0._m, 1._m};
auto b = a.cross(showerAxis_.GetDirection());
auto const lengthB = b.norm();
if (lengthB < 1e-10_m) {
b = geometry::Vector<length_d>{conexObservationCS_, 1_m, 0_m, 0_m};
}
return b.normalized();
})}
, y_sf_{showerAxis_.GetDirection().cross(x_sf_)} {
std::cout << "x_sf (conexObservationCS): " << x_sf_.GetComponents(conexObservationCS_)
<< std::endl;
std::cout << "x_sf (C8): " << x_sf_.GetComponents(center.GetCoordinateSystem())
<< std::endl;
std::cout << "y_sf (conexObservationCS): " << y_sf_.GetComponents(conexObservationCS_)
<< std::endl;
std::cout << "y_sf (C8): " << y_sf_.GetComponents(center.GetCoordinateSystem())
<< std::endl;
std::cout << "showerAxisDirection (conexObservationCS): "
<< showerAxis_.GetDirection().GetComponents(conexObservationCS_) << std::endl;
std::cout << "showerAxisDirection (C8): "
<< showerAxis_.GetDirection().GetComponents(center.GetCoordinateSystem())
<< std::endl;
std::cout << "showerCore (conexObservationCS): "
<< showerCore_.GetCoordinates(conexObservationCS_) << std::endl;
std::cout << "showerCore (C8): "
<< showerCore_.GetCoordinates(center.GetCoordinateSystem()) << std::endl;
int randomSeeds[3] = {1234, 0, 0}; // will be overwritten later??
int heModel = eSibyll23;
int nShower = 1; // large to avoid final stats.
int maxDetail = 0;
#ifdef CONEX_EXTENSIONS
int particleListMode = 0;
#endif
std::string configPath = CONEX_CONFIG_PATH;
conex::initconex_(nShower, randomSeeds, heModel, maxDetail,
#ifdef CONEX_EXTENSIONS
particleListMode,
#endif
configPath.c_str(), configPath.size());
double eprima = primaryEnergy / 1_GeV;
// set phi, theta
geometry::Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}};
auto const c = showerAxis_.GetDirection().dot(ez) / 1_m;
double theta = std::acos(c) * 180 / M_PI;
auto const showerAxisConex =
showerAxis_.GetDirection().GetComponents(conexObservationCS_);
double phi = std::atan2(-showerAxisConex.GetY().magnitude(),
showerAxisConex.GetX().magnitude()) *
180 / M_PI;
std::cout << "theta (deg) = " << theta << "; phi (deg) = " << phi << std::endl;
int ipart = static_cast<int>(primaryPDG);
auto rng = corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
double dimpact = 0.; // valid only if shower core is fixed on the observation plane; for
// skimming showers an offset is needed like in CONEX
std::array<int, 3> ioseed{static_cast<int>(rng()), static_cast<int>(rng()),
static_cast<int>(rng())};
double xminp = injectionHeight / 1_m;
conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
}
Processes/CONEXSourceCut/CONEXSourceCut.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/environment/ShowerAxis.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/process/conex_source_cut/CONEX_f.h>
namespace conex {
corsika::units::si::LengthType constexpr earthRadius{6371315 *
corsika::units::si::meter};
} // namespace conex
namespace corsika::process {
namespace conex_source_cut {
class CONEXSourceCut : public process::SecondariesProcess<CONEXSourceCut> {
public:
CONEXSourceCut(geometry::Point center, environment::ShowerAxis const& showerAxis,
units::si::LengthType groundDist,
units::si::LengthType injectionHeight,
units::si::HEPEnergyType primaryEnergy, particles::PDGCode pdg);
corsika::process::EProcessReturn DoSecondaries(corsika::setup::StackView&);
void SolveCE();
bool addParticle(particles::Code pid, units::si::HEPEnergyType energy,
units::si::HEPEnergyType mass, geometry::Point const& position,
geometry::Vector<units::si::dimensionless_d> const& direction,
units::si::TimeType t);
auto const& GetObserverCS() const { return conexObservationCS_; }
private:
//! CONEX e.m. particle codes
static std::array<std::pair<particles::Code, int>, 3> constexpr egs_em_codes_{
{{particles::Code::Gamma, 0},
{particles::Code::Electron, -1},
{particles::Code::Positron, -1}}};
geometry::Point const center_; //!< center of CONEX Earth
environment::ShowerAxis const& showerAxis_;
units::si::LengthType groundDist_; //!< length from injection point to shower core
geometry::Point const showerCore_; //!< shower core
geometry::CoordinateSystem const conexObservationCS_; //!< CONEX observation frame
geometry::Vector<units::si::dimensionless_d> const x_sf_,
y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction
};
} // namespace conex_source_cut
} // namespace corsika::process
Processes/CONEXSourceCut/conex_fit_REF.txt deleted 100644 → 0
View file @ eabed243
15 # log10(eprima/eV)
60 # theta
1e+30 # X1 (first interaction)
1093.9 # Nmax
475.585 # X0
50.3301 # P1
-0.00584492 # P2
3.55546e-06 # P3
4.49378e-05 # chi^2 / sqrt(Nmax)
880.549 # Xmax
-0 # phi
-1 # inelasticity 1st int.
100000 # ???
Processes/CONEXSourceCut/conex_output_REF.txt deleted 100644 → 0
View file @ eabed243
# X N dEdX Mu dMu Gamma El Had
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
10.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
20.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
30.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
40.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
50.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
60.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
70.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
80.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
90.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
100.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
110.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
120.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
130.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
140.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
150.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
160.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
170.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
180.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
190.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
200.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
210.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
220.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
230.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
240.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
250.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
260.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
270.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
280.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
290.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
300.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
310.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
320.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
330.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
340.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
350.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
360.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
370.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
380.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
390.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
400.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
410.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
420.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
430.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
440.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
450.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
460.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
470.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
480.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
490.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
500.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
510.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
520.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
530.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
540.00 0.00 0.00062 0.00 0.00 1.00 0.00 0.00
550.00 0.512 0.00246 0.000784 0.000802 1.69 0.272 0.00518
560.00 1.45 0.00598 0.00326 0.00256 4.23 0.787 0.00985
570.00 3.23 0.0123 0.00719 0.00411 9.26 1.78 0.0145
580.00 6.30 0.0227 0.0126 0.00569 18.0 3.51 0.0195
590.00 11.2 0.0385 0.0194 0.00734 32.0 6.27 0.0248
600.00 18.4 0.0613 0.0278 0.00906 53.3 10.4 0.0304
610.00 28.6 0.0925 0.0376 0.0108 84.1 16.2 0.0364
620.00 42.3 0.134 0.0489 0.0125 127.0 24.1 0.0425
630.00 60.1 0.186 0.0615 0.0142 184.0 34.5 0.0489
640.00 82.5 0.251 0.0753 0.0158 258.0 47.5 0.0553
650.00 110.0 0.328 0.0902 0.0173 350.0 63.4 0.0618
660.00 142.0 0.416 0.106 0.0186 464.0 82.5 0.0681
670.00 180.0 0.518 0.123 0.0198 600.0 105.0 0.0744
680.00 223.0 0.632 0.14 0.0208 759.0 130.0 0.0803
690.00 270.0 0.757 0.157 0.0216 942.0 158.0 0.086
700.00 322.0 0.892 0.175 0.0223 1.15e+03 189.0 0.0913
710.00 378.0 1.04 0.193 0.0227 1.38e+03 222.0 0.0962
720.00 437.0 1.19 0.21 0.023 1.62e+03 257.0 0.101
730.00 497.0 1.34 0.227 0.023 1.89e+03 294.0 0.104
740.00 559.0 1.50 0.244 0.0229 2.17e+03 331.0 0.108
750.00 621.0 1.65 0.26 0.0227 2.46e+03 369.0 0.111
760.00 683.0 1.81 0.275 0.0223 2.76e+03 406.0 0.113
770.00 742.0 1.96 0.289 0.0218 3.06e+03 443.0 0.114
780.00 799.0 2.10 0.303 0.0212 3.36e+03 478.0 0.116
790.00 853.0 2.23 0.315 0.0205 3.66e+03 511.0 0.116
800.00 902.0 2.35 0.327 0.0197 3.94e+03 542.0 0.116
810.00 947.0 2.45 0.337 0.0189 4.21e+03 570.0 0.116
820.00 986.0 2.55 0.347 0.018 4.47e+03 595.0 0.115
830.00 1.02e+03 2.63 0.355 0.0171 4.70e+03 616.0 0.113
840.00 1.05e+03 2.69 0.362 0.0162 4.92e+03 634.0 0.112
850.00 1.07e+03 2.74 0.369 0.0153 5.10e+03 648.0 0.109
860.00 1.08e+03 2.77 0.374 0.0144 5.27e+03 659.0 0.107
870.00 1.09e+03 2.79 0.379 0.0135 5.40e+03 666.0 0.104
880.00 1.10e+03 2.79 0.382 0.0126 5.50e+03 669.0 0.101
890.00 1.09e+03 2.78 0.385 0.0117 5.58e+03 669.0 0.0983
900.00 1.08e+03 2.76 0.387 0.0109 5.63e+03 665.0 0.095
910.00 1.07e+03 2.72 0.388 0.0101 5.64e+03 658.0 0.0916
920.00 1.05e+03 2.67 0.388 0.00933 5.64e+03 649.0 0.0881
930.00 1.03e+03 2.62 0.388 0.00861 5.60e+03 636.0 0.0845
940.00 1.01e+03 2.55 0.388 0.00793 5.54e+03 622.0 0.0809
950.00 978.0 2.48 0.387 0.00729 5.46e+03 605.0 0.0773
960.00 947.0 2.39 0.385 0.00669 5.36e+03 587.0 0.0737
970.00 913.0 2.31 0.383 0.00613 5.24e+03 567.0 0.0702
980.00 878.0 2.22 0.38 0.00561 5.11e+03 546.0 0.0667
990.00 841.0 2.12 0.378 0.00512 4.96e+03 524.0 0.0633
1000.00 803.0 2.03 0.375 0.00467 4.80e+03 501.0 0.0599
1010.00 765.0 1.93 0.371 0.00426 4.63e+03 478.0 0.0566
1020.00 726.0 1.83 0.368 0.00387 4.45e+03 455.0 0.0535
1030.00 688.0 1.73 0.364 0.00352 4.26e+03 431.0 0.0504
1040.00 650.0 1.64 0.36 0.0032 4.07e+03 408.0 0.0475
1050.00 612.0 1.54 0.356 0.0029 3.88e+03 385.0 0.0447
1060.00 575.0 1.45 0.352 0.00263 3.69e+03 362.0 0.042
1070.00 539.0 1.36 0.348 0.00239 3.50e+03 340.0 0.0394
1080.00 504.0 1.27 0.344 0.00216 3.31e+03 319.0 0.037
1090.00 470.0 1.19 0.34 0.00196 3.12e+03 298.0 0.0346
1100.00 438.0 1.10 0.336 0.00177 2.94e+03 278.0 0.0324
1110.00 407.0 1.03 0.331 0.0016 2.76e+03 258.0 0.0303
1120.00 377.0 0.952 0.327 0.00145 2.58e+03 240.0 0.0283
1130.00 349.0 0.881 0.323 0.00131 2.41e+03 222.0 0.0265
1140.00 323.0 0.814 0.319 0.00119 2.25e+03 206.0 0.0247
1150.00 298.0 0.751 0.314 0.00108 2.10e+03 190.0 0.023
1160.00 274.0 0.692 0.31 0.000973 1.95e+03 175.0 0.0215
1170.00 252.0 0.636 0.306 0.000881 1.81e+03 161.0 0.02
1180.00 231.0 0.584 0.302 0.000798 1.67e+03 148.0 0.0186
1190.00 212.0 0.535 0.298 0.000722 1.55e+03 136.0 0.0173
1200.00 194.0 0.49 0.294 0.000655 1.43e+03 125.0 0.0161
1210.00 177.0 0.447 0.29 0.000593 1.32e+03 114.0 0.015
1220.00 161.0 0.408 0.286 0.000538 1.21e+03 104.0 0.014
1230.00 147.0 0.372 0.283 0.000489 1.11e+03 94.8 0.013
1240.00 134.0 0.338 0.279 0.000444 1.02e+03 86.3 0.0121
1250.00 121.0 0.307 0.275 0.000403 933.0 78.5 0.0112
1260.00 110.0 0.279 0.272 0.000367 853.0 71.3 0.0104
1270.00 99.8 0.253 0.268 0.000335 779.0 64.7 0.0097
1280.00 90.3 0.229 0.265 0.000305 710.0 58.6 0.00901
1290.00 81.6 0.207 0.262 0.000278 647.0 53.0 0.00838
1300.00 73.7 0.187 0.258 0.000254 588.0 47.9 0.00779
1310.00 66.5 0.169 0.255 0.000232 534.0 43.2 0.00724
1320.00 59.9 0.152 0.252 0.000212 485.0 39.0 0.00673
1330.00 53.9 0.137 0.249 0.000194 439.0 35.1 0.00626
1340.00 48.5 0.123 0.246 0.000178 397.0 31.6 0.00582
1350.00 43.5 0.111 0.243 0.000163 359.0 28.4 0.00541
1360.00 39.1 0.0995 0.24 0.00015 324.0 25.5 0.00503
1370.00 35.0 0.0892 0.237 0.000138 293.0 22.9 0.00469
1380.00 31.4 0.08 0.234 0.000127 264.0 20.5 0.00436
1390.00 28.1 0.0716 0.231 0.000117 237.0 18.4 0.00406
1400.00 25.1 0.0641 0.229 0.000108 213.0 16.4 0.00379
1410.00 22.4 0.0573 0.226 9.93e-05 192.0 14.7 0.00353
1420.00 20.0 0.0512 0.223 9.17e-05 172.0 13.1 0.00329
1430.00 17.9 0.0457 0.221 8.48e-05 154.0 11.7 0.00307
1440.00 16.0 0.0408 0.218 7.85e-05 138.0 10.4 0.00287
1450.00 14.2 0.0364 0.216 7.28e-05 124.0 9.29 0.00268
1460.00 12.7 0.0325 0.214 6.75e-05 111.0 8.27 0.0025
1470.00 11.3 0.0289 0.211 6.27e-05 98.9 7.36 0.00234
1480.00 10.0 0.0258 0.209 5.83e-05 88.3 6.54 0.00219
1490.00 8.93 0.0229 0.207 5.42e-05 78.7 5.81 0.00205
1500.00 7.94 0.0204 0.204 5.05e-05 70.2 5.16 0.00192
1510.00 7.05 0.0182 0.202 4.70e-05 62.5 4.58 0.0018
1520.00 6.27 0.0161 0.2 4.39e-05 55.7 4.06 0.00169
1530.00 5.57 0.0144 0.198 4.10e-05 49.5 3.60 0.00159
1540.00 4.95 0.0128 0.196 3.83e-05 44.0 3.18 0.0015
1550.00 4.40 0.0114 0.194 3.58e-05 39.1 2.82 0.00141
1560.00 3.91 0.0101 0.192 3.35e-05 34.7 2.49 0.00133
1570.00 3.47 0.00899 0.19 3.14e-05 30.8 2.20 0.00125
1580.00 3.09 0.008 0.188 2.94e-05 27.3 1.95 0.00118
1590.00 2.75 0.00712 0.186 2.76e-05 24.2 1.72 0.00112
1600.00 2.44 0.00635 0.185 2.59e-05 21.4 1.52 0.00106
1610.00 2.18 0.00566 0.183 2.44e-05 19.0 1.34 0.001
1620.00 1.94 0.00505 0.181 2.29e-05 16.8 1.19 0.000948
1630.00 1.73 0.00451 0.179 2.16e-05 14.9 1.05 0.0009
1640.00 1.55 0.00404 0.177 2.04e-05 13.1 0.923 0.000855
1650.00 1.38 0.00362 0.176 1.92e-05 11.6 0.814 0.000814
1660.00 1.24 0.00324 0.174 1.81e-05 10.3 0.718 0.000775
1670.00 1.11 0.00292 0.173 1.71e-05 9.07 0.634 0.00074
1680.00 0.998 0.00263 0.171 1.62e-05 8.01 0.559 0.000706
1690.00 0.9 0.00237 0.169 1.53e-05 7.08 0.494 0.000675
1700.00 0.812 0.00214 0.168 1.45e-05 6.25 0.436 0.000647
1710.00 0.736 0.00195 0.166 1.38e-05 5.52 0.385 0.00062
1720.00 0.668 0.00177 0.165 1.31e-05 4.87 0.34 0.000595
1730.00 0.608 0.00161 0.163 1.24e-05 4.30 0.301 0.000572
1740.00 0.556 0.00148 0.162 1.18e-05 3.80 0.266 0.00055
1750.00 0.509 0.00136 0.161 1.12e-05 3.36 0.236 0.00053
1760.00 0.469 0.00125 0.159 1.07e-05 2.97 0.209 0.000512
1770.00 0.433 0.00116 0.158 1.02e-05 2.63 0.185 0.000494
1780.00 0.401 0.00108 0.156 9.71e-06 2.33 0.165 0.000478
1790.00 0.373 0.001 0.155 9.27e-06 2.06 0.147 0.000463
1800.00 0.349 0.000939 0.154 8.86e-06 1.83 0.131 0.000448
1810.00 0.327 0.000882 0.152 8.48e-06 1.62 0.117 0.000435
1820.00 0.308 0.000831 0.151 8.12e-06 1.44 0.105 0.000423
1830.00 0.291 0.000787 0.15 7.78e-06 1.28 0.0941 0.000411
1840.00 0.276 0.000747 0.149 7.46e-06 1.14 0.0847 0.0004
1850.00 0.262 0.000712 0.147 7.16e-06 1.02 0.0764 0.00039
1860.00 0.251 0.000681 0.146 6.88e-06 0.914 0.0692 0.000381
1870.00 0.24 0.000653 0.145 6.62e-06 0.82 0.0629 0.000372
1880.00 0.231 0.000629 0.144 6.37e-06 0.737 0.0573 0.000364
1890.00 0.222 0.000607 0.143 6.14e-06 0.665 0.0524 0.000356
1900.00 0.215 0.000587 0.142 5.92e-06 0.601 0.0481 0.000349
1910.00 0.208 0.000569 0.141 5.72e-06 0.545 0.0443 0.000342
1920.00 0.202 0.000553 0.139 5.52e-06 0.495 0.041 0.000335
1930.00 0.197 0.000539 0.138 5.34e-06 0.452 0.0381 0.000329
1940.00 0.192 0.000526 0.137 5.17e-06 0.414 0.0355 0.000324
1950.00 0.188 0.000514 0.136 5.01e-06 0.381 0.0332 0.000318
1960.00 0.184 0.000503 0.135 4.85e-06 0.351 0.0312 0.000313
1970.00 0.18 0.000493 0.134 4.71e-06 0.325 0.0295 0.000309
1980.00 0.177 0.000484 0.133 4.57e-06 0.303 0.0279 0.000304
1990.00 0.174 0.000476 0.132 4.44e-06 0.283 0.0265 0.0003
2000.00 0.171 0.000469 0.131 4.32e-06 0.265 0.0252 0.000296
2010.00 0.168 0.000461 0.13 4.20e-06 0.249 0.0241 0.000292
2020.00 0.166 0.000455 0.129 4.10e-06 0.235 0.0232 0.000289
2030.00 0.164 0.000449 0.128 3.99e-06 0.223 0.0223 0.000285
2040.00 0.162 0.000443 0.127 3.89e-06 0.212 0.0215 0.000282
2050.00 0.16 0.000437 0.127 3.80e-06 0.202 0.0208 0.000279
2060.00 0.158 0.00 0.126 3.71e-06 0.193 0.0201 0.000276
Processes/CONEXSourceCut/testCONEXSourceCut.cc deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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/setup/SetupEnvironment.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/MediumPropertyModel.h>
#include <corsika/environment/UniformMagneticField.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::environment;
using namespace corsika::geometry;
using namespace corsika::units::si;
const std::string refDataDir = std::string(REFDATADIR); // from cmake
template <typename T>
using MExtraEnvirnoment =
environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
TEST_CASE("CONEXSourceCut") {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
feenableexcept(FE_INVALID);
// setup environment, geometry
setup::Environment env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
auto builder = environment::make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface,
MExtraEnvirnoment>::create(center, conex::earthRadius,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 50_mT, 0_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env);
const HEPEnergyType E0 = 1_PeV;
double thetaDeg = 60.;
auto const thetaRad = thetaDeg / 180. * M_PI;
auto const observationHeight = 1.4_km + conex::earthRadius;
auto const injectionHeight = 112.75_km + conex::earthRadius;
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
units::static_pow<2>(injectionHeight));
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
showerCore +
Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
environment::ShowerAxis const showerAxis{injectionPos,
(showerCore - injectionPos) * 1.02, env};
// need to initialize Sibyll, done in constructor:
process::sibyll::Interaction sibyll;
[[maybe_unused]] process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
corsika::process::conex_source_cut::CONEXSourceCut conex(
center, showerAxis, t, injectionHeight, E0,
particles::GetPDG(particles::Code::Proton));
HEPEnergyType const Eem{1_PeV};
auto const momentum = showerAxis.GetDirection() * Eem;
auto const emPosition = showerCore + showerAxis.GetDirection() * (-20_km);
std::cout << "position injection: "
<< injectionPos.GetCoordinates(conex.GetObserverCS()) << " "
<< injectionPos.GetCoordinates(rootCS) << std::endl;
std::cout << "position core: " << showerCore.GetCoordinates(conex.GetObserverCS())
<< " " << showerCore.GetCoordinates(rootCS) << std::endl;
std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " "
<< emPosition.GetCoordinates(rootCS) << std::endl;
conex.addParticle(particles::Code::Proton, Eem, 0_eV, emPosition, momentum.normalized(),
0_s);
// supperimpose a photon
auto const momentumPhoton = showerAxis.GetDirection() * 1_TeV;
conex.addParticle(particles::Code::Gamma, 1_TeV, 0_eV, emPosition,
momentumPhoton.normalized(), 0_s);
conex.SolveCE();
}
#include <algorithm>
#include <iterator>
#include <string>
#include <fstream>
TEST_CASE("ConexOutput", "[output validation]") {
auto file = GENERATE(as<std::string>{}, "conex_fit", "conex_output");
SECTION(std::string("check saved data, ") + file + ".txt") {
// compare to binary reference data
std::ifstream file1(file + ".txt");
std::ifstream file1ref(refDataDir + "/" + file + "_REF.txt");
std::istreambuf_iterator<char> begin1(file1);
std::istreambuf_iterator<char> begin1ref(file1ref);
std::istreambuf_iterator<char> end;
while (begin1 != end && begin1ref != end) {
CHECK(*begin1 == *begin1ref);
++begin1;
++begin1ref;
}
CHECK(begin1 == end);
CHECK(begin1ref == end);
file1.close();
file1ref.close();
}
}
Processes/ExampleProcessors/CMakeLists.txt deleted 100644 → 0
View file @ eabed243
set (
ExampleProcessors_HEADERS
DummyBoundaryCrossingProcess.h
DummyContinuousProcess.h
DummyDecayProcess.h
DummyInteractionProcess.h
DummySecondariesProcess.h
)
set (
ExampleProcessors_NAMESPACE
corsika/process/example_processors
)
add_library (ExampleProcessors INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ExampleProcessors ${ExampleProcessors_NAMESPACE} ${ExampleProcessors_HEADERS})
target_include_directories (
ExampleProcessors
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
FILES ${ExampleProcessors_HEADERS}
DESTINATION include/${ExampleProcessors_NAMESPACE}
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST (TestDummy)
target_link_libraries (
TestDummy
ExampleProcessors
CORSIKAunits
CORSIKAstackinterface
CORSIKAprocesssequence
CORSIKAsetup
CORSIKAgeometry
CORSIKAenvironment
CORSIKAtesting
)
Processes/ExampleProcessors/DummyBoundaryCrossingProcess.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/BoundaryCrossingProcess.h>
#include <chrono>
#include <thread>
namespace corsika::process {
namespace example_processors {
template <int ISleep>
class DummyBoundaryCrossingProcess
: public BoundaryCrossingProcess<DummyBoundaryCrossingProcess<ISleep>> {
private:
protected:
public:
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle&, VTNType const&, VTNType const&) {
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return EProcessReturn::eOk;
}
};
} // namespace example_processors
} // namespace corsika::process
\ No newline at end of file
Processes/ExampleProcessors/DummyContinuousProcess.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/ContinuousProcess.h>
#include <chrono>
#include <thread>
namespace corsika::process {
namespace example_processors {
template <int ISleep>
class DummyContinuousProcess
: public ContinuousProcess<DummyContinuousProcess<ISleep>> {
private:
public:
template <typename Particle, typename Track>
inline EProcessReturn DoContinuous(Particle&, Track const&) const {
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return process::EProcessReturn::eOk;
}
template <typename Particle, typename Track>
inline units::si::LengthType MaxStepLength(Particle const&, Track const&) const {
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return units::si::meter * std::numeric_limits<double>::infinity();
}
std::string name() { return "DummyContinuousProcess"; }
};
} // namespace example_processors
} // namespace corsika::process
\ No newline at end of file
Processes/ExampleProcessors/DummyDecayProcess.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/DecayProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <chrono>
#include <thread>
namespace corsika::process {
namespace example_processors {
template <int ISleep>
class DummyDecayProcess : public DecayProcess<DummyDecayProcess<ISleep>> {
private:
public:
template <typename Particle>
EProcessReturn DoDecay(Particle&) {
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return process::EProcessReturn::eOk;
}
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle&) {
using namespace corsika::units::si;
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return std::numeric_limits<double>::infinity() * 1_s;
}
};
} // namespace example_processors
} // namespace corsika::process
\ No newline at end of file
Processes/ExampleProcessors/DummyInteractionProcess.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/InteractionProcess.h>
#include <chrono>
#include <thread>
namespace corsika::process {
namespace example_processors {
template <int ISleep>
class DummyInteractionProcess
: public InteractionProcess<DummyInteractionProcess<ISleep> > {
private:
public:
template <typename Particle>
EProcessReturn DoInteraction(Particle&) {
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return process::EProcessReturn::eOk;
}
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle&) {
using namespace corsika::units::si;
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return std::numeric_limits<double>::infinity() * (1_g / 1_cm / 1_cm);
}
};
} // namespace example_processors
} // namespace corsika::process
\ No newline at end of file
Processes/ExampleProcessors/DummySecondariesProcess.h deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/SecondariesProcess.h>
#include <chrono>
#include <thread>
namespace corsika::process {
namespace example_processors {
template <int ISleep>
class DummySecondariesProcess
: public SecondariesProcess<DummySecondariesProcess<ISleep>> {
private:
public:
template <typename TSecondaries>
inline EProcessReturn DoSecondaries(TSecondaries&) {
std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
return process::EProcessReturn::eOk;
}
};
} // namespace example_processors
} // namespace corsika::process
\ No newline at end of file
Processes/ExampleProcessors/TestDummy.cc deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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/example_processors/DummyBoundaryCrossingProcess.h>
#include <corsika/process/example_processors/DummyContinuousProcess.h>
#include <corsika/process/example_processors/DummyDecayProcess.h>
#include <corsika/process/example_processors/DummyInteractionProcess.h>
#include <corsika/process/example_processors/DummySecondariesProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <catch2/catch.hpp>
#include <chrono>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::process::example_processors;
using namespace corsika::units::si;
TEST_CASE("Dummy Processes") {
DummyBoundaryCrossingProcess<1000> dbc;
DummyContinuousProcess<1000> dc;
DummyDecayProcess<1000> dd;
DummyInteractionProcess<1000> di;
DummySecondariesProcess<1000> dse;
int tmp = 0;
SECTION("BoundaryCrossing") {
auto start = std::chrono::steady_clock::now();
REQUIRE(dbc.DoBoundaryCrossing(tmp, 0, 0) == EProcessReturn::eOk);
auto end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
}
SECTION("Continuous") {
auto start = std::chrono::steady_clock::now();
REQUIRE(dc.DoContinuous(tmp, nullptr) == EProcessReturn::eOk);
auto end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
start = std::chrono::steady_clock::now();
REQUIRE(dc.MaxStepLength(nullptr, nullptr) ==
units::si::meter * std::numeric_limits<double>::infinity());
end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
}
SECTION("Decay") {
auto start = std::chrono::steady_clock::now();
REQUIRE(dd.DoDecay(tmp) == EProcessReturn::eOk);
auto end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
start = std::chrono::steady_clock::now();
REQUIRE(dd.GetLifetime(tmp) ==
units::si::second * std::numeric_limits<double>::infinity());
end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
}
SECTION("Interaction") {
auto start = std::chrono::steady_clock::now();
REQUIRE(di.DoInteraction(tmp) == EProcessReturn::eOk);
auto end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
start = std::chrono::steady_clock::now();
REQUIRE(di.GetInteractionLength(tmp) ==
(units::si::gram / 1_cm / 1_cm) * std::numeric_limits<double>::infinity());
end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
}
SECTION("Secondaries") {
auto start = std::chrono::steady_clock::now();
REQUIRE(dse.DoSecondaries(tmp) == EProcessReturn::eOk);
auto end = std::chrono::steady_clock::now();
REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
Approx(1000).margin(1));
}
}
Processes/InteractionCounter/CMakeLists.txt deleted 100644 → 0
View file @ eabed243
set (
MODEL_HEADERS
InteractionCounter.hpp
InteractionHistogram.hpp
)
set (
MODEL_SOURCES
InteractionHistogram.cc
)
set (
MODEL_NAMESPACE
corsika/process/interaction_counter
)
add_library (ProcessInteractionCounter STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessInteractionCounter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessInteractionCounter
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessInteractionCounter
CORSIKAunits
CORSIKAutilities
CORSIKAprocesssequence
CORSIKAthirdparty
C8::ext::boost
)
target_include_directories (
ProcessInteractionCounter
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testInteractionCounter SOURCES
testInteractionCounter.cc
${MODEL_HEADERS}
)
#
target_compile_definitions (
testInteractionCounter
PRIVATE
REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}"
)
#
target_link_libraries (
testInteractionCounter
ProcessInteractionCounter
CORSIKAsetup
CORSIKAunits
CORSIKAenvironment
CORSIKAtesting
)
Processes/InteractionCounter/InteractionCounter.hpp deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/process/interaction_counter/InteractionHistogram.hpp>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/setup/SetupStack.h>
namespace corsika::process::interaction_counter {
/*!
* Wrapper around an InteractionProcess that fills histograms of the number
* of calls to DoInteraction() binned in projectile energy (both in
* lab and center-of-mass frame) and species
*/
template <class TCountedProcess>
class InteractionCounter
: public InteractionProcess<InteractionCounter<TCountedProcess>> {
TCountedProcess& process_;
InteractionHistogram histogram_;
public:
InteractionCounter(TCountedProcess& process)
: process_(process) {}
template <typename TSecondaryView>
auto DoInteraction(TSecondaryView& view) {
auto const projectile = view.GetProjectile();
auto const massNumber = projectile.GetNode()
->GetModelProperties()
.GetNuclearComposition()
.GetAverageMassNumber();
auto const massTarget = massNumber * units::constants::nucleonMass;
if (auto const projectile_id = projectile.GetPID();
projectile_id == particles::Code::Nucleus) {
auto const A = projectile.GetNuclearA();
auto const Z = projectile.GetNuclearZ();
histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget, A, Z);
} else {
histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget);
}
return process_.DoInteraction(view);
}
template <typename TParticle>
auto GetInteractionLength(TParticle const& particle) const {
return process_.GetInteractionLength(particle);
}
InteractionHistogram const& GetHistogram() const { return histogram_; }
};
} // namespace corsika::process::interaction_counter
Processes/InteractionCounter/InteractionHistogram.cc deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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/interaction_counter/InteractionHistogram.hpp>
#include <fstream>
#include <string>
using namespace corsika::process::interaction_counter;
InteractionHistogram::InteractionHistogram()
: inthist_cms_{detail::hist_factory(num_bins_cms, lower_edge_cms, upper_edge_cms)}
, inthist_lab_{detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab)} {}
void InteractionHistogram::fill(particles::Code projectile_id,
units::si::HEPEnergyType lab_energy,
units::si::HEPEnergyType mass_target, int A, int Z) {
using namespace units::si;
auto constexpr inv_eV = 1 / 1_eV;
if (projectile_id == particles::Code::Nucleus) {
auto const sqrtS =
sqrt(A * A * (units::constants::nucleonMass * units::constants::nucleonMass) +
mass_target * mass_target + 2 * lab_energy * mass_target);
int32_t const pdg = 1'000'000'000l + Z * 10'000l + A * 10l;
inthist_lab_(pdg, lab_energy * inv_eV);
inthist_cms_(pdg, sqrtS * inv_eV);
} else {
auto const projectile_mass = particles::GetMass(projectile_id);
auto const sqrtS = sqrt(projectile_mass * projectile_mass +
mass_target * mass_target + 2 * lab_energy * mass_target);
inthist_cms_(static_cast<int>(particles::GetPDG(projectile_id)), sqrtS * inv_eV);
inthist_lab_(static_cast<int>(particles::GetPDG(projectile_id)), lab_energy * inv_eV);
}
}
void InteractionHistogram::saveLab(std::string const& filename,
utl::SaveMode mode) const {
corsika::utl::save_hist(inthist_lab_, filename, mode);
}
void InteractionHistogram::saveCMS(std::string const& filename,
utl::SaveMode mode) const {
corsika::utl::save_hist(inthist_cms_, filename, mode);
}
InteractionHistogram& InteractionHistogram::operator+=(
InteractionHistogram const& other) {
inthist_lab_ += other.inthist_lab_;
inthist_cms_ += other.inthist_cms_;
return *this;
}
InteractionHistogram InteractionHistogram::operator+(InteractionHistogram other) const {
other.inthist_lab_ += inthist_lab_;
other.inthist_cms_ += inthist_cms_;
return other;
}
Processes/InteractionCounter/testInteractionCounter.cc deleted 100644 → 0
View file @ eabed243
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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/interaction_counter/InteractionCounter.hpp>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/setup/SetupStack.h>
#include <catch2/catch.hpp>
#include <numeric>
using namespace corsika;
using namespace corsika::process::interaction_counter;
using namespace corsika::units;
using namespace corsika::units::si;
const std::string refDataDir = std::string(REFDATADIR); // from cmake
struct DummyProcess {
template <typename TParticle>
GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
return 100_g / 1_cm / 1_cm;
}
template <typename TParticle>
auto DoInteraction([[maybe_unused]] TParticle& projectile) {
return nullptr;
}
};
TEST_CASE("InteractionCounter", "[process]") {
logging::SetLevel(logging::level::debug);
DummyProcess d;
InteractionCounter countedProcess(d);
SECTION("GetInteractionLength") {
REQUIRE(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
}
auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
[[maybe_unused]] auto& env_dummy = env;
SECTION("DoInteraction nucleus") {
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
Z, 105_TeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
auto const ret = countedProcess.DoInteraction(*secViewPtr);
REQUIRE(ret == nullptr);
auto const& h = countedProcess.GetHistogram().labHist();
REQUIRE(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
auto const& h2 = countedProcess.GetHistogram().CMSHist();
REQUIRE(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
utl::SaveMode::overwrite);
countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz",
utl::SaveMode::overwrite);
}
SECTION("DoInteraction Lambda") {
auto constexpr code = particles::Code::Lambda0;
auto [stackPtr, secViewPtr] =
setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
auto const ret = countedProcess.DoInteraction(*secViewPtr);
REQUIRE(ret == nullptr);
auto const& h = countedProcess.GetHistogram().labHist();
REQUIRE(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
auto const& h2 = countedProcess.GetHistogram().CMSHist();
REQUIRE(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
}
}
#include <algorithm>
#include <iterator>
#include <string>
#include <fstream>
TEST_CASE("InteractionCounterOutput", "[output validation]") {
auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
"testInteractionCounter_file2");
SECTION(std::string("check saved data, ") + file + ".npz") {
std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
// compare to binary reference data
std::ifstream file1(file + ".npz");
std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");
std::istreambuf_iterator<char> begin1(file1);
std::istreambuf_iterator<char> begin1ref(file1ref);
std::istreambuf_iterator<char> end;
while (begin1 != end && begin1ref != end) {
CHECK(*begin1 == *begin1ref);
++begin1;
++begin1ref;
}
CHECK(begin1 == end);
CHECK(begin1ref == end);
file1.close();
file1ref.close();
}
}

data privacy - impressum