diff --git a/CMakeLists.txt b/CMakeLists.txt
index 524f0ff9c3064967fa5310ad2c10d529eff4fd23..e423548404e3247bffbe4fce3deab8bd89f19e48 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,6 +7,9 @@ project (
LANGUAGES CXX
)
+# as long as there still are modules using it:
+enable_language (Fortran)
+
# ignore many irrelevant Up-to-date messages during install
set (CMAKE_INSTALL_MESSAGE LAZY)
@@ -21,7 +24,21 @@ set (CMAKE_CXX_STANDARD 17)
enable_testing ()
set (CTEST_OUTPUT_ON_FAILURE 1)
-ENABLE_LANGUAGE (Fortran)
+# Set a default build type if none was specified
+set(default_build_type "Release")
+if(EXISTS "${CMAKE_SOURCE_DIR}/.git")
+ set(default_build_type "Debug")
+endif()
+
+if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
+ message(STATUS "Setting build type to '${default_bluild_type}' as no other was specified.")
+ set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
+ STRING "Choose the type of build." FORCE)
+ # Set the possible values of build type for cmake-gui
+ set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
+ "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
+endif()
+
# unit testing coverage, does not work yet
#include (CodeCoverage)
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 35b9908c94473ee9137ba4413e362c0509ea856d..2a66a3470584d14733e8a4651238e13b95756879 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -31,6 +31,7 @@ add_executable (staticsequence_example staticsequence_example.cc)
target_link_libraries (staticsequence_example
CORSIKAprocesssequence
CORSIKAunits
+ CORSIKAgeometry
CORSIKAlogging)
install (TARGETS staticsequence_example DESTINATION share/examples)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 7b4ab0ad4e344ce3ccae06b2555247a14e47f4f3..08e8cd21e8f48a2d18d427a58c1c067bfdcd3b49 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -12,14 +12,15 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/stack_inspector/StackInspector.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/random/RNGManager.h>
-#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/cascade/SibStack.h>
+#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
@@ -29,9 +30,10 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
+using namespace corsika::setup;
-#include <typeinfo>
#include <iostream>
+#include <typeinfo>
using namespace std;
static int fCount = 0;
@@ -41,23 +43,25 @@ public:
ProcessSplit() {}
template <typename Particle>
- double MinStepLength(Particle& p) const {
+ double MinStepLength(Particle& p, setup::Trajectory&) const {
const Code corsikaBeamId = p.GetPID();
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
+
int kBeam = process::sibyll::GetSibyllXSCode( corsikaBeamId );
bool kInteraction = process::sibyll::CanInteract( corsikaBeamId );
/*
the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
+ ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
+
int kTarget = 16;
double beamEnergy = p.GetEnergy() / 1_GeV;
#warning boost to cm. still missing, sibyll cross section input is cm. energy!
@@ -98,25 +102,27 @@ public:
/*
what are the units of the output? slant depth or 3space length?
-
+
*/
- std::cout << "ProcessSplit: " << "next step (g/cm2): " << next_step << std::endl;
+ std::cout << "ProcessSplit: "
+ << "next step (g/cm2): " << next_step << std::endl;
return next_step;
}
- template <typename Particle, typename Trajectory, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
// corsika::utls::ignore(p);
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
+
cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
if( process::sibyll::CanInteract( p.GetPID() ) ){
cout << "defining coordinates" << endl;
// coordinate system, get global frame of reference
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point pOrig(rootCS, coordinates);
@@ -253,84 +259,87 @@ public:
//cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
}
}else
+
p.Delete();
}
-
-
- void Init()
- {
+
+ void Init() {
fCount = 0;
// define reference frame? --> defines boosts between corsika stack and model stack.
-
+
// initialize random numbers for sibyll
// FOR NOW USE SIBYLL INTERNAL !!!
// rnd_ini_();
-
- corsika::random::RNGManager & rmng = corsika::random::RNGManager::GetInstance();;
+
+ corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+ ;
const std::string str_name = "s_rndm";
rmng.RegisterRandomStream(str_name);
-
+
// // corsika::random::RNG srng;
// auto srng = rmng.GetRandomStream("s_rndm");
// test random number generator
- std::cout << "ProcessSplit: " << " test sequence of random numbers." << std::endl;
+ std::cout << "ProcessSplit: "
+ << " test sequence of random numbers." << std::endl;
int a = 0;
- for(int i=0; i<8; ++i)
- std::cout << i << " " << s_rndm_(a) << std::endl;
-
- //initialize Sibyll
+ for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
+
+ // initialize Sibyll
sibyll_ini_();
// set particles stable / unstable
// use stack to loop over particles
setup::Stack ds;
- ds.NewParticle().SetPID( Code::Proton );
- ds.NewParticle().SetPID( Code::Neutron );
- ds.NewParticle().SetPID( Code::PiPlus );
- ds.NewParticle().SetPID( Code::PiMinus );
- ds.NewParticle().SetPID( Code::KPlus );
- ds.NewParticle().SetPID( Code::KMinus );
- ds.NewParticle().SetPID( Code::K0Long );
- ds.NewParticle().SetPID( Code::K0Short );
-
- for( auto &p: ds){
- int s_id = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+ ds.NewParticle().SetPID(Code::Proton);
+ ds.NewParticle().SetPID(Code::Neutron);
+ ds.NewParticle().SetPID(Code::PiPlus);
+ ds.NewParticle().SetPID(Code::PiMinus);
+ ds.NewParticle().SetPID(Code::KPlus);
+ ds.NewParticle().SetPID(Code::KMinus);
+ ds.NewParticle().SetPID(Code::K0Long);
+ ds.NewParticle().SetPID(Code::K0Short);
+
+ for (auto& p : ds) {
+ int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
// set particle stable by setting table value negative
- cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")" << " stable in Sibyll .." << endl;
- s_csydec_.idb[ s_id ] = -s_csydec_.idb[ s_id-1 ];
+ cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")"
+ << " stable in Sibyll .." << endl;
+ s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1];
p.Delete();
}
}
-
+
int GetCount() { return fCount; }
private:
};
-double s_rndm_(int &)
-{
- static corsika::random::RNG& rmng = corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");;
- return rmng()/(double)rmng.max();
+double s_rndm_(int&) {
+ static corsika::random::RNG& rmng =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ ;
+ return rmng() / (double)rmng.max();
}
-
-int main(){
+int main() {
// coordinate system, get global frame of reference
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point pOrig(rootCS, coordinates);
- stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+ // stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+ tracking_line::TrackingLine<setup::Stack> tracking;
+ stack_inspector::StackInspector<setup::Stack> p0(true);
+
ProcessSplit p1;
const auto sequence = p0 + p1;
setup::Stack stack;
-
- corsika::cascade::Cascade EAS(sequence, stack);
+ corsika::cascade::Cascade EAS(tracking, sequence, stack);
stack.Clear();
auto particle = stack.NewParticle();
@@ -341,10 +350,10 @@ int main(){
0. * 1_GeV / si::constants::c,
P0);
particle.SetEnergy(E0);
+
particle.SetMomentum(plab);
particle.SetPID( Code::Proton );
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
-
}
diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index 048347aee5be19df25571cfc96752f49e1a0dad7..5bacd1489ab9bff7add254ad8edc70ccada9ea43 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/Documentation/Examples/geometry_example.cc
@@ -9,8 +9,8 @@
* the license.
*/
-#include <corsika/geometry/CoordinateSystem.h>
#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
@@ -19,12 +19,14 @@
#include <iostream>
#include <typeinfo>
+using namespace corsika;
using namespace corsika::geometry;
using namespace corsika::units::si;
int main() {
// define the root coordinate system
- auto const root = CoordinateSystem::CreateRootCS();
+ geometry::CoordinateSystem& root =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCS();
// another CS defined by a translation relative to the root CS
CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});
diff --git a/Documentation/Examples/helix_example.cc b/Documentation/Examples/helix_example.cc
index 2797a30a1e9eb61d54cb51b4a3040e4dea80ca80..c51afce291181a4629b4fceaf776b0bfa9fc40cd 100644
--- a/Documentation/Examples/helix_example.cc
+++ b/Documentation/Examples/helix_example.cc
@@ -9,20 +9,22 @@
* the license.
*/
-#include <corsika/geometry/CoordinateSystem.h>
#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <array>
#include <cstdlib>
#include <iostream>
+using namespace corsika;
using namespace corsika::geometry;
using namespace corsika::units::si;
int main() {
- auto const root = CoordinateSystem::CreateRootCS();
+ geometry::CoordinateSystem& root =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCS();
Point const r0(root, {0_m, 0_m, 0_m});
auto const omegaC = 2 * M_PI * 1_Hz;
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index 8049983f767a5ea28544f83cda9bffc84a35d95a..07cc7561e0a8222e26d9e5b7eeba8ced68f6f63b 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/Documentation/Examples/staticsequence_example.cc
@@ -15,6 +15,11 @@
#include <corsika/process/ProcessSequence.h>
+#include <corsika/setup/SetupTrajectory.h> // TODO: try to break this dependence later
+using corsika::setup::Trajectory;
+#include <corsika/units/PhysicalUnits.h> // dito
+using namespace corsika::units::si;
+
using namespace std;
using namespace corsika::process;
@@ -72,7 +77,6 @@ struct DummyData {
double p[10];
};
struct DummyStack {};
-struct DummyTrajectory {};
void modular() {
@@ -84,8 +88,8 @@ void modular() {
const auto sequence = m1 + m2 + m3 + m4;
DummyData p;
- DummyTrajectory t;
DummyStack s;
+ Trajectory t;
const int n = 100000000;
for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
index 8681fae994b1ff91cb8bb1d69f98f5a23da1e50c..ad48313530b2df7dfcbda79cf8dcaf2791a0bf18 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -15,6 +15,7 @@ set (
#set (
# CORSIKAcascade_SOURCES
+# TrackingStep.cc
# Cascade.cc
# )
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 0e920cc4bb6e2fa1eb5e25378848ab57c8d8a0a7..e46960dc9c3d8e8c33254c2751a601ab72843a18 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -9,19 +9,21 @@
* the license.
*/
-#ifndef _include_Cascade_h_
-#define _include_Cascade_h_
+#ifndef _include_corsika_cascade_Cascade_h_
+#define _include_corsika_cascade_Cascade_h_
#include <corsika/process/ProcessReturn.h>
#include <corsika/units/PhysicalUnits.h>
+#include <type_traits>
+
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika::units::si;
namespace corsika::cascade {
- template <typename ProcessList, typename Stack> //, typename Trajectory>
+ template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
typedef typename Stack::ParticleType Particle;
@@ -29,32 +31,26 @@ namespace corsika::cascade {
Cascade() = delete;
public:
- Cascade(ProcessList& pl, Stack& stack)
- : fProcesseList(pl)
- , fStack(stack) {}
+ Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
+ : fTracking(tr)
+ , fProcesseList(pl)
+ , fStack(stack) {
+ // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoDiscrete)>::value,
+ //"ProcessList has not function DoDiscrete.");
+ // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoContinuous)>::value,
+ // "ProcessList has not function DoContinuous.");
+ }
void Init() {
- fStack.Init();
+ fTracking.Init();
fProcesseList.Init();
+ fStack.Init();
}
void Run() {
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
Particle& pNext = *fStack.GetNextParticle();
- /*
- EnergyType Emin;
- typename Stack::StackIterator pNext(fStack, 0);
- bool first = true;
- for (typename Stack::StackIterator ip = fStack.begin(); ip != fStack.end();
- ++ip) {
- if (first || ip.GetEnergy() < Emin) {
- first = false;
- pNext = ip;
- Emin = pMin.GetEnergy();
- }
- }
- */
Step(pNext);
}
// do cascade equations, which can put new particles on Stack,
@@ -62,27 +58,25 @@ namespace corsika::cascade {
// DoCascadeEquations(); //
}
}
-
+
void Step(Particle& particle) {
- [[maybe_unused]] double nextStep = fProcesseList.MinStepLength(particle);
- // corsika::utls::ignore(nextStep);
- auto const root = corsika::geometry::CoordinateSystem::CreateRootCS();
- corsika::geometry::Trajectory<corsika::geometry::Line>
- trajectory( // trajectory is not yet used. this is a dummy.
- corsika::geometry::Line(corsika::geometry::Point(root, {0_m, 0_m, 0_m}),
- corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>(
- root, 0 * 1_m / second, 0 * 1_m / second, 1 * 1_m / second)),
- 0_s, 1_s);
+ corsika::setup::Trajectory step = fTracking.GetTrack(particle);
+ fProcesseList.MinStepLength(particle, step);
+
+ /// here the particle is actually moved along the trajectory to new position:
+ std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
+
corsika::process::EProcessReturn status =
- fProcesseList.DoContinuous(particle, trajectory, fStack);
+ fProcesseList.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
- fStack.Delete(particle);
+ fStack.Delete(particle); // TODO: check if this is really needed
} else {
fProcesseList.DoDiscrete(particle, fStack);
}
}
private:
+ Tracking& fTracking;
ProcessList& fProcesseList;
Stack& fStack;
};
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index c10ea7d8d97e0959da6fe911fc87347c9324ec6d..57a3f9c16b099c6d00e71b405223f6c3f7eb41c5 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -1,12 +1,12 @@
#ifndef _include_sibstack_h_
#define _include_sibstack_h_
-#include <vector>
#include <string>
+#include <vector>
-#include <corsika/stack/Stack.h>
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/stack/Stack.h>
using namespace std;
using namespace corsika::stack;
@@ -14,10 +14,10 @@ using namespace corsika::units;
using namespace corsika::geometry;
class SibStackData {
-
- public:
+
+public:
void Init();
-
+
void Clear() { s_plist_.np = 0; }
int GetSize() const { return s_plist_.np; }
@@ -40,7 +40,7 @@ class SibStackData {
super_stupid::MomentumVector GetMomentum(const int i) const
{
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
super_stupid::MomentumVector v1(rootCS,components);
return v1;
@@ -50,18 +50,20 @@ class SibStackData {
s_plist_.llist[i1] = s_plist_.llist[i2];
s_plist_.p[3][i1] = s_plist_.p[3][i2];
}
-
- protected:
+
+protected:
void IncrementSize() { s_plist_.np++; }
- void DecrementSize() { if ( s_plist_.np>0) { s_plist_.np--; } }
+ void DecrementSize() {
+ if (s_plist_.np > 0) { s_plist_.np--; }
+ }
};
-
-template<typename StackIteratorInterface>
+template <typename StackIteratorInterface>
class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetStackData;
using ParticleBase<StackIteratorInterface>::GetIndex;
- public:
+
+public:
void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
@@ -70,7 +72,6 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
};
-
typedef Stack<SibStackData, ParticleInterface> SibStack;
#endif
diff --git a/Framework/Cascade/sibyll2.3c.h b/Framework/Cascade/sibyll2.3c.h
index e66f1f74426ede7caf3bbf19ff60097ae6c2aa3c..5d2d1e69644e812246b731aa4bf276d090d79910 100644
--- a/Framework/Cascade/sibyll2.3c.h
+++ b/Framework/Cascade/sibyll2.3c.h
@@ -3,94 +3,88 @@
//----------------------------------------------
// C++ interface for the SIBYLL event generator
//----------------------------------------------
-//wrapper
+// wrapper
+extern "C" {
+typedef char s_name[6];
+// SIBYLL particle stack (FORTRAN COMMON)
+// variables are: np : numer of particles on stack
+// p : 4momentum + mass of particles on stack
+// llist : id of particles on stack
+extern struct {
+ double p[5][8000];
+ int llist[8000];
+ int np;
+} s_plist_;
-extern"C"{
+extern struct {
+ double cbr[223 + 16 + 12 + 8];
+ int kdec[1338 + 6 * (16 + 12 + 8)];
+ int lbarp[99];
+ int idb[99];
+} s_csydec_;
- typedef char s_name[6];
+// additional particle stack for the mother particles of unstable particles
+// stable particles have entry zero
+extern struct { int llist1[8000]; } s_plist1_;
- // SIBYLL particle stack (FORTRAN COMMON)
- // variables are: np : numer of particles on stack
- // p : 4momentum + mass of particles on stack
- // llist : id of particles on stack
- extern struct {
- double p[5][8000];
- int llist[8000];
- int np;
- }s_plist_;
+// tables with particle properties
+// charge, strangeness and baryon number
+extern struct {
+ int ichp[99];
+ int istr[99];
+ int ibar[99];
+} s_chp_;
+// tables with particle properties
+// mass and mass squared
+extern struct {
+ double am[99];
+ double am2[99];
+} s_mass1_;
- extern struct {
- double cbr[223+16+12+8];
- int kdec[1338+6*(16+12+8)];
- int lbarp[99];
- int idb[99];
- }s_csydec_;
+// table with particle names
+extern struct { char namp[6][99]; } s_cnam_;
-
- // additional particle stack for the mother particles of unstable particles
- // stable particles have entry zero
- extern struct {
- int llist1[8000];
- }s_plist1_;
+// debug info
+extern struct {
+ int ncall;
+ int ndebug;
+ int lun;
+} s_debug_;
+// lund random generator setup
+// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
- // tables with particle properties
- // charge, strangeness and baryon number
- extern struct {
- int ichp[99];
- int istr[99];
- int ibar[99];
- }s_chp_;
+// sibyll main subroutine
+void sibyll_(int&, int&, double&);
- // tables with particle properties
- // mass and mass squared
- extern struct {
- double am[99];
- double am2[99];
- }s_mass1_;
+// subroutine to initiate sibyll
+void sibyll_ini_();
- // table with particle names
- extern struct {char namp[6][99];}s_cnam_;
+// subroutine to SET DECAYS
+void dec_ini_();
- // debug info
- extern struct {int ncall; int ndebug; int lun; }s_debug_;
+// subroutine to initiate random number generator
+// void rnd_ini_();
- // lund random generator setup
- //extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
+// print event
+void sib_list_(int&);
+// decay routine
+void decsib_();
- // sibyll main subroutine
- void sibyll_(int&,int&,double&);
-
- // subroutine to initiate sibyll
- void sibyll_ini_();
+// interaction length
+// double fpni_(double&, int&);
- // subroutine to SET DECAYS
- void dec_ini_();
+void sib_sigma_hnuc_(int&, int&, double&, double&, double&);
+void sib_sigma_hp_(int&, double&, double&, double&, double&, double*, double&, double&);
- // subroutine to initiate random number generator
- //void rnd_ini_();
-
- // print event
- void sib_list_(int&);
+double s_rndm_(int&);
- // decay routine
- void decsib_();
-
- // interaction length
- //double fpni_(double&, int&);
-
- void sib_sigma_hnuc_(int&,int&,double&,double&,double&);
- void sib_sigma_hp_(int&,double&,double&,double&,double&,double*,double&,double&);
-
- double s_rndm_(int&);
-
- // phojet random generator setup
- void pho_rndin_(int&,int&,int&,int&);
+// phojet random generator setup
+void pho_rndin_(int&, int&, int&, int&);
}
#endif
-
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index edb825bc78e70961920dcba98c28b1de1d02d337..ad4daa2e4df0a554727f1a0584cc88a3af3cff72 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -10,11 +10,19 @@
*/
#include <corsika/cascade/Cascade.h>
+
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/stack_inspector/StackInspector.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/stack/super_stupid/SuperStupidStack.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
@@ -23,6 +31,7 @@
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
+using namespace corsika::geometry;
#include <iostream>
using namespace std;
@@ -34,13 +43,10 @@ public:
ProcessSplit() {}
template <typename Particle>
- double MinStepLength(Particle&) const {
- return 0;
- }
+ void MinStepLength(Particle&, setup::Trajectory&) const {}
- template <typename Particle, typename Trajectory, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
- // corsika::utls::ignore(p);
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
return EProcessReturn::eOk;
}
@@ -52,30 +58,42 @@ public:
fCount++;
} else {
p.SetEnergy(E / 2);
- s.NewParticle().SetEnergy(E / 2);
+ auto pnew = s.NewParticle();
+ // s.Copy(p, pnew);
+ pnew.SetEnergy(E / 2);
+ pnew.SetPosition(p.GetPosition());
+ pnew.SetMomentum(p.GetMomentum());
}
}
void Init() { fCount = 0; }
- int GetCount() { return fCount; }
+ int GetCount() const { return fCount; }
private:
};
+
TEST_CASE("Cascade", "[Cascade]") {
- stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+ tracking_line::TrackingLine<setup::Stack> tracking;
+
+ stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
setup::Stack stack;
- corsika::cascade::Cascade EAS(sequence, stack);
+ corsika::cascade::Cascade EAS(tracking, sequence, stack);
+
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
particle.SetEnergy(E0);
+ particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
+ particle.SetMomentum(
+ corsika::stack::super_stupid::MomentumVector(rootCS, {0*newton*second, 0*newton*second, -1*newton*second}));
EAS.Init();
EAS.Run();
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index 580066ae0fbe9e9f6c79c3980dd6a17bfc9c1fd9..2b0c8d3af13a5ba00b5fbe6674febb0b21461f96 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -11,11 +11,12 @@ set (
Line.h
Sphere.h
CoordinateSystem.h
+ RootCoordinateSystem.h
Helix.h
BaseVector.h
QuantityVector.h
Trajectory.h
- BaseTrajectory.h
+ # BaseTrajectory.h
)
set (
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index 089a3f80c96a8c7c2e9b2ebceeae4c6dc5bb2c8e..3a06dfbd8bc65b06a83a58184088d5bec42e8b03 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -21,6 +21,8 @@ typedef Eigen::Translation<double, 3> EigenTranslation;
namespace corsika::geometry {
+ class RootCoordinateSystem;
+
using corsika::units::si::length_d;
class CoordinateSystem {
@@ -30,14 +32,17 @@ namespace corsika::geometry {
CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
: reference(&reference)
, transf(transf) {}
-
+
CoordinateSystem()
: // for creating the root CS
transf(EigenTransform::Identity()) {}
+ protected:
+ static auto CreateCS() { return CoordinateSystem(); }
+ friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
+ /// creat ONE unique root CS
+
public:
- static auto CreateRootCS() { return CoordinateSystem(); }
-
static EigenTransform GetTransformation(CoordinateSystem const& c1,
CoordinateSystem const& c2);
@@ -55,9 +60,9 @@ namespace corsika::geometry {
auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
if (axis.eVector.isZero()) {
- throw std::string("null-vector given as axis parameter");
+ throw std::string("null-vector given as axis parameter");
}
-
+
EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
return CoordinateSystem(*this, rotation);
@@ -66,9 +71,9 @@ namespace corsika::geometry {
auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
QuantityVector<phys::units::length_d> axis, double angle) {
if (axis.eVector.isZero()) {
- throw std::string("null-vector given as axis parameter");
+ throw std::string("null-vector given as axis parameter");
}
-
+
EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
EigenTranslation(translation.eVector)};
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index b647688d45c012d5188dd5a8a6147eabdca75ffc..8eb215820097fae23eef0e535631384d31db01c4 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -31,7 +31,7 @@ namespace corsika::geometry {
*/
class Helix {
-
+
using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
Point const r0;
@@ -58,8 +58,8 @@ namespace corsika::geometry {
auto GetRadius() const { return radius; }
- LengthType DistanceBetween(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const {
+ LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
+ corsika::units::si::TimeType t2) const {
return (vPar + vPerp).norm() * (t2 - t1);
}
};
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index 829c16044a274785f838f4975f78af4c42354744..e093a0055fb82ceaeab2c656227b2c005f15958e 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -19,7 +19,7 @@
namespace corsika::geometry {
class Line {
-
+
using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
Point const r0;
@@ -30,13 +30,11 @@ namespace corsika::geometry {
: r0(pR0)
, v0(pV0) {}
- Point GetPosition(corsika::units::si::TimeType t) const {
- return r0 + v0 * t;
- }
-
- LengthType DistanceBetween(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const {
- assert(t2 >= t1);
+ Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
+
+ LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
+ corsika::units::si::TimeType t2) const {
+ // assert(t2 >= t1);
return v0.norm() * (t2 - t1);
}
};
diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h
index 26de8510662aefc6961c536541483d28c136c3ce..c9c88d7bf599e3f41d33bd92dfe746706f5e5339 100644
--- a/Framework/Geometry/Point.h
+++ b/Framework/Geometry/Point.h
@@ -34,8 +34,11 @@ namespace corsika::geometry {
Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
: BaseVector<phys::units::length_d>(cs, {x, y, z}) {}
+ // TODO: this should be private or protected, we don NOT want to expose numbers
+ // without reference to outside:
auto GetCoordinates() const { return BaseVector<phys::units::length_d>::qVector; }
+ /// this always returns a QuantityVector as triple
auto GetCoordinates(CoordinateSystem const& pCS) const {
if (&pCS == BaseVector<phys::units::length_d>::cs) {
return BaseVector<phys::units::length_d>::qVector;
diff --git a/Framework/Geometry/RootCoordinateSystem.h b/Framework/Geometry/RootCoordinateSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..4c3133becd9c59e4b50a9730b893a933657c5518
--- /dev/null
+++ b/Framework/Geometry/RootCoordinateSystem.h
@@ -0,0 +1,33 @@
+#ifndef _include_corsika_geometry_rootcoordinatesystem_h_
+#define _include_corsika_geometry_rootcoordinatesystem_h_
+
+#include <corsika/utl/Singleton.h>
+
+#include <corsika/geometry/CoordinateSystem.h>
+
+/*!
+ * This is the only way to get a root-coordinate system, and it is a
+ * singleton. All other CoordinateSystems must be relative to the
+ * RootCoordinateSystem
+ */
+
+namespace corsika::geometry {
+
+ class RootCoordinateSystem : public corsika::utl::Singleton<RootCoordinateSystem> {
+
+ friend class corsika::utl::Singleton<RootCoordinateSystem>;
+
+ protected:
+ RootCoordinateSystem() {}
+
+ public:
+ corsika::geometry::CoordinateSystem& GetRootCS() { return fRootCS; }
+ const corsika::geometry::CoordinateSystem& GetRootCS() const { return fRootCS; }
+
+ private:
+ corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT
+ };
+
+} // namespace corsika::geometry
+
+#endif
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 97adfdd3bba6a4c449316b37fec428d77060c8db..686799756504138e349aacfa428b325166ef1e23 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -12,37 +12,39 @@
#ifndef _include_TRAJECTORY_H
#define _include_TRAJECTORY_H
-#include <corsika/geometry/BaseTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
+using corsika::units::si::LengthType;
+using corsika::units::si::TimeType;
+
namespace corsika::geometry {
template <typename T>
- class Trajectory : public BaseTrajectory {
+ class Trajectory : public T {
- T fTraj;
+ corsika::units::si::TimeType fTimeLength;
public:
- Trajectory(T const& theT, corsika::units::si::TimeType pTStart,
- corsika::units::si::TimeType pTEnd)
- //: T(theT), fTStart(pTStart), fTEnd(pTEnd) {}
- : BaseTrajectory(pTStart, pTEnd)
- , fTraj(theT) {}
+ using T::GetDistanceBetween;
+ using T::GetPosition;
+
+ Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
+ : T(theT)
+ , fTimeLength(timeLength) {}
- Point GetPosition(corsika::units::si::TimeType t) const {
+ /*Point GetPosition(corsika::units::si::TimeType t) const {
return fTraj.GetPosition(t + fTStart);
- }
+ }*/
- Point GetPosition(double u) const {
- return GetPosition(fTEnd * u + fTStart * (1 - u));
- }
-
- LengthType GetDistance(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const {
- return fTraj.DistanceBetween(t1, t2);
- }
-
+ Point GetPosition(const double u) const { return T::GetPosition(fTimeLength * u); }
+
+ TimeType GetDuration() const { return fTimeLength; }
+ LengthType GetDistance(const corsika::units::si::TimeType t) const {
+ assert(t > fTimeLength);
+ assert(t >= 0 * corsika::units::si::second);
+ return T::DistanceBetween(0, t);
+ }
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 742938ffec44f3e65054313634f8c246c33f40e2..9e3be3310bbaf3ed6c9ae22bb6faa67701f20687 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -17,6 +17,7 @@
#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
@@ -28,7 +29,7 @@ using namespace corsika::units::si;
double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
.isApprox(EigenTransform::Identity()));
@@ -44,10 +45,11 @@ TEST_CASE("transformations between CoordinateSystems") {
REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() ==
Approx(0).margin(absMargin));
+ /*
SECTION("unconnected CoordinateSystems") {
CoordinateSystem rootCS2 = CoordinateSystem::CreateRootCS();
REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2));
- }
+ }*/
SECTION("translations") {
QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
@@ -126,7 +128,7 @@ TEST_CASE("transformations between CoordinateSystems") {
}
TEST_CASE("Sphere") {
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
@@ -137,7 +139,7 @@ TEST_CASE("Sphere") {
}
TEST_CASE("Trajectories") {
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
Point r0(rootCS, {0_m, 0_m, 0_m});
SECTION("Line") {
@@ -145,16 +147,16 @@ TEST_CASE("Trajectories") {
{1_m / second, 0_m / second, 0_m / second});
Line const line(r0, v0);
- CHECK((line.GetPosition(2_s).GetCoordinates() -
- QuantityVector<length_d>(2_m, 0_m, 0_m))
- .norm()
- .magnitude() == Approx(0).margin(absMargin));
+ CHECK(
+ (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(2_m, 0_m, 0_m))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
- Trajectory<Line> base(line, 0_s, 1_s);
+ Trajectory<Line> base(line, 1_s);
CHECK(line.GetPosition(2_s).GetCoordinates() ==
base.GetPosition(2_s).GetCoordinates());
- CHECK(base.GetDistance(1_s, 2_s) / 1_m == Approx(1));
+ CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
}
SECTION("Helix") {
@@ -178,10 +180,10 @@ TEST_CASE("Trajectories") {
.norm()
.magnitude() == Approx(0).margin(absMargin));
- Trajectory<Helix> const base(helix, 0_s, 1_s);
+ Trajectory<Helix> const base(helix, 1_s);
CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
base.GetPosition(1234_s).GetCoordinates());
- CHECK(base.GetDistance(1_s, 2_s) / 1_m == Approx(5));
+ CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(5));
}
}
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index c24567821d938b68fa27e9f2a7be29eec33bb295..c1f5a45862faffc92aba606e8fc1abde94220d89 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -42,6 +42,8 @@ add_executable (
target_link_libraries (
testProcessSequence
+ CORSIKAsetup
+ CORSIKAgeometry
CORSIKAprocesssequence
CORSIKAthirdparty # for catch2
)
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 4188fb0aa56e9ff8735560cde5b9121a91dce22f..a07e3ef84c91e6e4e54743c4a4e87c09d933bf2f 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -17,8 +17,14 @@
#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/ProcessReturn.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <variant>
+
//#include <type_traits> // still needed ?
+using corsika::setup::Trajectory;
+
namespace corsika::process {
/* namespace detail { */
@@ -143,7 +149,7 @@ namespace corsika::process {
// example for a trait-based call:
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
- template <typename Particle, typename Trajectory, typename Stack>
+ template <typename Particle, typename Stack>
inline EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
@@ -155,9 +161,10 @@ namespace corsika::process {
return ret;
}
- template <typename D>
- inline double MinStepLength(D& d) const {
- return std::min(A.MinStepLength(d), B.MinStepLength(d));
+ template <typename Particle>
+ inline void MinStepLength(Particle& p, Trajectory& step) const {
+ A.MinStepLength(p, step);
+ B.MinStepLength(p, step);
}
/*
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index 2d4b46b15f11a8d9b396bf24e0073bdc2cdbb147..9757298cdc2269a90fd875cdedfa68a16f19ab5a 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -19,17 +19,24 @@
#include <corsika/process/ProcessSequence.h>
+#include <corsika/setup/SetupTrajectory.h> // TODO: maybe try to break this dependency later!
+using corsika::setup::Trajectory;
+#include <corsika/units/PhysicalUnits.h>
+using namespace corsika::units::si;
+
using namespace std;
using namespace corsika::process;
+static const int nData = 10;
+
class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
public:
ContinuousProcess1() {}
void Init() { cout << "ContinuousProcess1::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
cout << "ContinuousProcess1::DoContinuous" << endl;
- for (int i = 0; i < 10; ++i) d.p[i] += 0.933;
+ for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
}
@@ -44,10 +51,10 @@ class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
public:
ContinuousProcess2() {}
void Init() { cout << "ContinuousProcess2::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
cout << "ContinuousProcess2::DoContinuous" << endl;
- for (int i = 0; i < 20; ++i) d.p[i] += 0.933;
+ for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
}
@@ -62,9 +69,9 @@ class Process1 : public DiscreteProcess<Process1> {
public:
Process1() {}
void Init() { cout << "Process1::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
- for (int i = 0; i < 10; ++i) d.p[i] += 1 + i;
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+ for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
@@ -78,9 +85,9 @@ class Process2 : public DiscreteProcess<Process2> {
public:
Process2() {}
void Init() { cout << "Process2::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
- for (int i = 0; i < 10; ++i) d.p[i] *= 0.7;
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+ for (int i = 0; i < nData; ++i) d.p[i] *= 0.7;
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
@@ -94,9 +101,9 @@ class Process3 : public DiscreteProcess<Process3> {
public:
Process3() {}
void Init() { cout << "Process3::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
- for (int i = 0; i < 10; ++i) d.p[i] += 0.933;
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+ for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
@@ -110,9 +117,9 @@ class Process4 : public BaseProcess<Process4> {
public:
Process4() {}
void Init() { cout << "Process4::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
- for (int i = 0; i < 10; ++i) d.p[i] /= 1.2;
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+ for (int i = 0; i < nData; ++i) d.p[i] /= 1.2;
return EProcessReturn::eOk;
}
// inline double MinStepLength(D& d) {
@@ -120,10 +127,9 @@ public:
};
struct DummyData {
- double p[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
struct DummyStack {};
-struct DummyTrajectory {};
TEST_CASE("Cascade", "[Cascade]") {
@@ -142,12 +148,19 @@ TEST_CASE("Cascade", "[Cascade]") {
const auto sequence2 = cp1 + m2 + m3 + cp2;
DummyData p;
- DummyTrajectory t;
DummyStack s;
cout << "-->init" << endl;
sequence2.Init();
cout << "-->docont" << endl;
+
+ // auto const root = corsika::geometry::CoordinateSystem::CreateRootCS();
+ // corsika::geometry::Point pos(root, {0_m, 0_m, 0_m});
+ // corsika::geometry::Vector<SpeedType::dimension_type> vec(root,
+ // {1_m/1_s,0_m/1_s,0_m/1_s}); corsika::geometry::Line traj(pos, vec);
+ Trajectory
+ t; //(corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns));
+
sequence2.DoContinuous(p, t, s);
cout << "-->dodisc" << endl;
sequence2.DoDiscrete(p, s);
@@ -155,11 +168,11 @@ TEST_CASE("Cascade", "[Cascade]") {
sequence.Init();
- const int n = 100;
- cout << "Running loop with n=" << n << endl;
- for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
-
- for (int i = 0; i < 10; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
+ const int nLoop = 5;
+ cout << "Running loop with n=" << nLoop << endl;
+ for (int i = 0; i < nLoop; ++i) { sequence.DoContinuous(p, t, s); }
+ for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
+ cout << "done" << endl;
}
SECTION("sectionThree") {}
diff --git a/Framework/Random/RNGManager.h b/Framework/Random/RNGManager.h
index c9fab9ef70e45dfc737c5cd2e4ed5dcad3146212..c440e7aaae0beb3ae745562fa4949eff5ff504ee 100644
--- a/Framework/Random/RNGManager.h
+++ b/Framework/Random/RNGManager.h
@@ -30,13 +30,13 @@ namespace corsika::random {
class RNGManager : public corsika::utl::Singleton<RNGManager> {
friend class corsika::utl::Singleton<RNGManager>;
-
+
std::map<std::string, RNG> rngs;
- std::map<std::string, std::seed_seq> seeds;
+ std::map<std::string, std::seed_seq> seeds;
protected:
RNGManager() {}
-
+
public:
/*!
* This function is to be called by a module requiring a random-number
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 65a98e482f7727af4176ded2c2b84339e6fcdf6e..94f98ac2ccf6a0f9c96f06e7b4429b9c022e1e19 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -56,8 +56,7 @@ namespace corsika::units::si::constants {
// barn moved to PhysicalUnits
// constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
-
-
+
// etc.
} // namespace corsika::units::si::constants
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index f1a1b845ab854284ed6c348f0c0df40b55fea4b7..e97003eee0a0e649021e20c420776afd6e586e59 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -16,7 +16,7 @@ namespace corsika::units::si {
using namespace phys::units;
using namespace phys::units::literals;
// namespace literals = phys::units::literals;
-
+
/// defining momentum you suckers
/// dimensions, i.e. composition in base SI dimensions
using momentum_d = phys::units::dimensions<1, 1, -1>;
@@ -35,14 +35,12 @@ namespace corsika::units::si {
using ElectricChargeType =
phys::units::quantity<phys::units::electric_charge_d, double>;
using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
- using MassType = phys::units::quantity<phys::units::mass_d, double>;
+ using MassType = phys::units::quantity<phys::units::mass_d, double>;
using MomentumType = phys::units::quantity<momentum_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
} // end namespace corsika::units::si
-
-
/**
* @file PhysicalUnits
*
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 712146c6b657fff07000ac891c71d29f632990ad..5785cc4aef1b285c3ac0676e692bbeaecdfa8d17 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -32,7 +32,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
LengthType l1 = 10_nm;
l1 = l1;
-
+
LengthType arr0[5];
arr0[0] = 5_m;
@@ -58,7 +58,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
REQUIRE(1_K / 1_zK == Approx(1e21));
REQUIRE(1_K / 1_yK == Approx(1e24));
REQUIRE(1_barn / 1_mbarn == Approx(1e3));
-
+
REQUIRE(1_A / 1_hA == Approx(1e-2));
REQUIRE(1_m / 1_km == Approx(1e-3));
REQUIRE(1_m / 1_Mm == Approx(1e-6));
@@ -82,8 +82,8 @@ TEST_CASE("PhysicalUnits", "[Units]") {
const MassType m = 1_kg;
const SpeedType v = 1_m / 1_s;
- REQUIRE( m*v == 1_newton_second);
-
+ REQUIRE(m * v == 1_newton_second);
+
const double lgE = log10(E2 / 1_GeV);
REQUIRE(lgE == Approx(log10(40.)));
@@ -95,7 +95,6 @@ TEST_CASE("PhysicalUnits", "[Units]") {
const LengthType farAway = std::numeric_limits<double>::infinity() * meter;
REQUIRE(farAway > 100000_m);
- REQUIRE_FALSE(farAway < 1e19*meter);
-
+ REQUIRE_FALSE(farAway < 1e19 * meter);
}
}
diff --git a/Framework/Utilities/Bit.h b/Framework/Utilities/Bit.h
index 4c21449d710d04aa21cbb6082c4a1c7615796e97..67473b37a131e2baff58cf9ca5f4b0d744a84a18 100644
--- a/Framework/Utilities/Bit.h
+++ b/Framework/Utilities/Bit.h
@@ -14,27 +14,27 @@
// #include <utl/AugerException.h>
-
namespace corsika::utl {
namespace Bit {
- template<typename T>
+ template <typename T>
class Array {
public:
- Array(T& target) : fTarget(target) { }
+ Array(T& target)
+ : fTarget(target) {}
class Bit {
public:
- Bit(T& target, T mask) : fTarget(target), fMask(mask) { }
+ Bit(T& target, T mask)
+ : fTarget(target)
+ , fMask(mask) {}
operator bool() const { return fTarget & fMask; }
bool operator~() const { return !bool(*this); }
- Bit&
- operator=(const bool value)
- {
+ Bit& operator=(const bool value) {
if (value)
fTarget |= fMask;
else
@@ -49,41 +49,38 @@ namespace corsika::utl {
T fMask;
};
- Bit operator[](unsigned int position)
- { return Bit(fTarget, T(1) << position); }
+ Bit operator[](unsigned int position) { return Bit(fTarget, T(1) << position); }
- Bit
- At(unsigned int position)
- {
- if (position >= 8*sizeof(T))
- //throw std::exceptionOutOfBoundException("Running out of bits.");
- throw std::exception("Running out of bits.");
+ Bit At(unsigned int position) {
+ if (position >= 8 * sizeof(T))
+ // throw std::exceptionOutOfBoundException("Running out of bits.");
+ throw std::exception("Running out of bits.");
return (*this)[position];
}
- template<typename M>
- Array& Mask(const M mask, const bool value)
- { Bit(fTarget, mask) = value; return *this; }
+ template <typename M>
+ Array& Mask(const M mask, const bool value) {
+ Bit(fTarget, mask) = value;
+ return *this;
+ }
- template<typename M>
- T Get(const M mask) { return fTarget & T(mask); }
+ template <typename M>
+ T Get(const M mask) {
+ return fTarget & T(mask);
+ }
private:
T& fTarget;
};
- }
+ } // namespace Bit
// helper
- template<typename T>
- inline
- Bit::Array<T>
- AsBitArray(T& target)
- {
+ template <typename T>
+ inline Bit::Array<T> AsBitArray(T& target) {
return Bit::Array<T>(target);
}
-}
-
+} // namespace corsika::utl
#endif
diff --git a/Framework/Utilities/Dummy.h b/Framework/Utilities/Dummy.h
index 09260f7eecc293d334aed170b614bad360e15056..dac2021b49df401d221c5160f825909440468e22 100644
--- a/Framework/Utilities/Dummy.h
+++ b/Framework/Utilities/Dummy.h
@@ -4,7 +4,6 @@
namespace corsika::utl {
// void....
-
}
#endif
diff --git a/Framework/Utilities/Singleton.h b/Framework/Utilities/Singleton.h
index c718065a19d96ef9fb17320139ec16dba17a4cc1..e06d974cc3c7de8003991227ded2cd6df0b0ab53 100644
--- a/Framework/Utilities/Singleton.h
+++ b/Framework/Utilities/Singleton.h
@@ -34,25 +34,23 @@ namespace corsika::utl {
* \ingroup stl
*/
- template<typename T>
+ template <typename T>
class Singleton {
public:
- static
- T&
- GetInstance()
-# ifdef __MAKECINT__
- ;
-# else
+ static T& GetInstance()
+#ifdef __MAKECINT__
+ ;
+#else
{
static T instance;
return instance;
}
-# endif
+#endif
protected:
// derived class can call ctor and dtor
- Singleton() { }
- ~Singleton() { }
+ Singleton() {}
+ ~Singleton() {}
private:
// no one should do copies
@@ -63,22 +61,18 @@ namespace corsika::utl {
#else
/// classical Gamma singleton
- template<typename T>
+ template <typename T>
class Singleton {
public:
- static
- T&
- GetInstance()
- {
- if (!fgInstance)
- fgInstance = new T;
+ static T& GetInstance() {
+ if (!fgInstance) fgInstance = new T;
return *fgInstance;
}
protected:
// derived class can call ctor and dtor
- Singleton() { }
- ~Singleton() { }
+ Singleton() {}
+ ~Singleton() {}
private:
// no one should do copies
@@ -90,7 +84,6 @@ namespace corsika::utl {
#endif
-
/**
* \class LeakingSingleton Singleton.h utl/Singleton.h
*
@@ -123,21 +116,18 @@ namespace corsika::utl {
* \ingroup stl
*/
- template<class T>
+ template <class T>
class LeakingSingleton {
public:
- static
- T&
- GetInstance()
- {
+ static T& GetInstance() {
static T* const instance = new T;
return *instance;
}
protected:
// derived class can call ctor and dtor
- LeakingSingleton() { }
- ~LeakingSingleton() { }
+ LeakingSingleton() {}
+ ~LeakingSingleton() {}
private:
// no one should do copies
@@ -145,7 +135,6 @@ namespace corsika::utl {
LeakingSingleton& operator=(const LeakingSingleton&);
};
-}
-
+} // namespace corsika::utl
#endif
diff --git a/Framework/Utilities/Test.h b/Framework/Utilities/Test.h
index 8135219a0d71e5913721e4e904e367aac3970771..41c7cecd7a9856e2c5b68f845d22a25f64a3cf0a 100644
--- a/Framework/Utilities/Test.h
+++ b/Framework/Utilities/Test.h
@@ -13,73 +13,73 @@
\ingroup testing
*/
-#include <cmath>
+#include <utl/Triple.h>
#include <boost/format.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_comparison.hpp>
#include <boost/tuple/tuple_io.hpp>
-#include <utl/Triple.h>
-
+#include <cmath>
namespace utl {
//! Predicate used in STL for searching for whitespace
struct IsSpace {
- bool operator()(const char x) const
- { return x == ' ' || x == '\r' || x == '\n' || x == '\t'; }
+ bool operator()(const char x) const {
+ return x == ' ' || x == '\r' || x == '\n' || x == '\t';
+ }
};
-
/// Predicate for equality
class Equal {
public:
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return lhs == rhs; }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return lhs == rhs;
+ }
static const char* Name() { return "equal"; }
};
-
/// Predicate for less
class Less {
public:
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return lhs < rhs; }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return lhs < rhs;
+ }
static const char* Name() { return "less"; }
};
-
/// Predicate for less or equal
class LessOrEqual {
public:
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return lhs <= rhs; }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return lhs <= rhs;
+ }
static const char* Name() { return "less or equal"; }
};
-
/// Predicate for greater
class Greater {
public:
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return lhs > rhs; }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return lhs > rhs;
+ }
static const char* Name() { return "greater"; }
};
-
/// Predicate for greater or equal
class GreaterOrEqual {
public:
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return lhs >= rhs; }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return lhs >= rhs;
+ }
static const char* Name() { return "greater or equal"; }
};
@@ -90,41 +90,41 @@ namespace utl {
*/
class CloseTo {
public:
- CloseTo(const double eps = 1e-6) : fEpsilon(eps) { }
+ CloseTo(const double eps = 1e-6)
+ : fEpsilon(eps) {}
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return IsCloseTo(lhs, rhs); }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return IsCloseTo(lhs, rhs);
+ }
- boost::format Name() const
- { return boost::format("close (@%g) to") % fEpsilon; }
+ boost::format Name() const { return boost::format("close (@%g) to") % fEpsilon; }
protected:
- template<typename T>
- bool IsCloseAbs(const T& lhs, const T& rhs) const
- { return std::abs(double(lhs) - double(rhs)) < fEpsilon; }
-
- bool IsCloseAbs(const utl::Triple& lhs, const utl::Triple& rhs) const
- { return std::sqrt(TupleDist2(lhs, rhs)) < fEpsilon; }
-
- template<typename T>
- bool IsCloseRel(const T& lhs, const T& rhs) const
- { return 2*std::abs(double(lhs) - double(rhs)) / (std::abs(double(lhs)) + std::abs(double(rhs))) < fEpsilon; }
-
- bool
- IsCloseRel(const utl::Triple& lhs, const utl::Triple& rhs)
- const
- {
- return (2*sqrt(TupleDist2(lhs, rhs)) /
- (sqrt(TupleDist2(lhs, utl::Triple(0,0,0))) +
- sqrt(TupleDist2(rhs, utl::Triple(0,0,0))))) < fEpsilon;
- }
-
- template<typename T>
- bool
- IsCloseTo(const T& lhs, const T& rhs)
- const
- {
+ template <typename T>
+ bool IsCloseAbs(const T& lhs, const T& rhs) const {
+ return std::abs(double(lhs) - double(rhs)) < fEpsilon;
+ }
+
+ bool IsCloseAbs(const utl::Triple& lhs, const utl::Triple& rhs) const {
+ return std::sqrt(TupleDist2(lhs, rhs)) < fEpsilon;
+ }
+
+ template <typename T>
+ bool IsCloseRel(const T& lhs, const T& rhs) const {
+ return 2 * std::abs(double(lhs) - double(rhs)) /
+ (std::abs(double(lhs)) + std::abs(double(rhs))) <
+ fEpsilon;
+ }
+
+ bool IsCloseRel(const utl::Triple& lhs, const utl::Triple& rhs) const {
+ return (2 * sqrt(TupleDist2(lhs, rhs)) /
+ (sqrt(TupleDist2(lhs, utl::Triple(0, 0, 0))) +
+ sqrt(TupleDist2(rhs, utl::Triple(0, 0, 0))))) < fEpsilon;
+ }
+
+ template <typename T>
+ bool IsCloseTo(const T& lhs, const T& rhs) const {
if (IsCloseAbs(lhs, rhs))
return true;
else
@@ -132,107 +132,106 @@ namespace utl {
}
// tuple distance
- template<typename Head, typename Tail>
- static
- double
- TupleDist2(const boost::tuples::cons<Head, Tail>& lhs,
- const boost::tuples::cons<Head, Tail>& rhs)
- {
+ template <typename Head, typename Tail>
+ static double TupleDist2(const boost::tuples::cons<Head, Tail>& lhs,
+ const boost::tuples::cons<Head, Tail>& rhs) {
const double t = lhs.get_head() - rhs.get_head();
- return t*t + TupleDist2(lhs.get_tail(), rhs.get_tail());
+ return t * t + TupleDist2(lhs.get_tail(), rhs.get_tail());
}
static double TupleDist2(const boost::tuples::null_type& /*lhs*/,
- const boost::tuples::null_type& /*rhs*/)
- { return 0; }
+ const boost::tuples::null_type& /*rhs*/) {
+ return 0;
+ }
double fEpsilon;
};
-
class CloseAbs : public CloseTo {
public:
- CloseAbs(const double eps = 1e-6) : CloseTo(eps) { }
+ CloseAbs(const double eps = 1e-6)
+ : CloseTo(eps) {}
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return IsCloseAbs(lhs, rhs); }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return IsCloseAbs(lhs, rhs);
+ }
- boost::format Name() const
- { return boost::format("absolutely close (@%g) to") % fEpsilon; }
+ boost::format Name() const {
+ return boost::format("absolutely close (@%g) to") % fEpsilon;
+ }
};
-
class CloseRel : public CloseTo {
public:
- CloseRel(const double eps = 1e-6) : CloseTo(eps) { }
+ CloseRel(const double eps = 1e-6)
+ : CloseTo(eps) {}
- template<typename T>
- bool operator()(const T& lhs, const T& rhs) const
- { return IsCloseRel(lhs, rhs); }
+ template <typename T>
+ bool operator()(const T& lhs, const T& rhs) const {
+ return IsCloseRel(lhs, rhs);
+ }
- boost::format Name() const
- { return boost::format("relatively close (@%g) to") % fEpsilon; }
+ boost::format Name() const {
+ return boost::format("relatively close (@%g) to") % fEpsilon;
+ }
};
-
- template<typename Predicate>
+ template <typename Predicate>
class Not : public Predicate {
public:
- Not() : Predicate() { }
-
- Not(const double eps) : Predicate(eps) { }
-
- template<typename T>
- bool operator()(const T& x) const
- { return !Predicate::operator()(x); }
-
- template<typename T, typename U>
- bool operator()(const T& x, const U& y) const
- { return !Predicate::operator()(x, y); }
-
- template<typename T, typename U, typename W>
- bool operator()(const T& x, const U& y, const W& z) const
- { return !Predicate::operator()(x, y, z); }
-
- static boost::format Name()
- { return boost::format("not-%s") % Predicate().Name(); }
- };
+ Not()
+ : Predicate() {}
+
+ Not(const double eps)
+ : Predicate(eps) {}
+
+ template <typename T>
+ bool operator()(const T& x) const {
+ return !Predicate::operator()(x);
+ }
+ template <typename T, typename U>
+ bool operator()(const T& x, const U& y) const {
+ return !Predicate::operator()(x, y);
+ }
+
+ template <typename T, typename U, typename W>
+ bool operator()(const T& x, const U& y, const W& z) const {
+ return !Predicate::operator()(x, y, z);
+ }
- inline
- utl::Triple
- Diff(const utl::Triple& lhs, const utl::Triple& rhs)
- {
- return utl::Triple(lhs.get<0>() - rhs.get<0>(),
- lhs.get<1>() - rhs.get<1>(),
+ static boost::format Name() { return boost::format("not-%s") % Predicate().Name(); }
+ };
+
+ inline utl::Triple Diff(const utl::Triple& lhs, const utl::Triple& rhs) {
+ return utl::Triple(lhs.get<0>() - rhs.get<0>(), lhs.get<1>() - rhs.get<1>(),
lhs.get<2>() - rhs.get<2>());
}
-
/// Test condition by evaluating a predicate
/** If the predicate evaluates to false, we print a failure message
with the values of the left- and right-hand side and the name of
the predicate. This information is normally not available when
using the CPPUNIT_ASSERT macro.
*/
- template<class Predicate, typename T>
- inline bool Test(const Predicate& pred, const T& lhs, const T& rhs)
- { return pred(lhs, rhs); }
-
+ template <class Predicate, typename T>
+ inline bool Test(const Predicate& pred, const T& lhs, const T& rhs) {
+ return pred(lhs, rhs);
+ }
/// Main test function
- template<class Predicate, typename T>
- inline bool Test(const T& lhs, const T& rhs)
- { return Test(Predicate(), lhs, rhs); }
-
+ template <class Predicate, typename T>
+ inline bool Test(const T& lhs, const T& rhs) {
+ return Test(Predicate(), lhs, rhs);
+ }
/// Test function for predicates that take an option
- template<class Predicate, typename T, typename U>
- inline bool Test(const T& lhs, const T& rhs, const U& eps)
- { return Test(Predicate(eps), lhs, rhs); }
-
-}
+ template <class Predicate, typename T, typename U>
+ inline bool Test(const T& lhs, const T& rhs, const U& eps) {
+ return Test(Predicate(eps), lhs, rhs);
+ }
+} // namespace utl
#endif
diff --git a/Framework/Utilities/testBit.cc b/Framework/Utilities/testBit.cc
index dc4f454368d80f4388baad7b5b9e772e2b319bf5..a4d36d1bcae35a23c3baaa0fa96a0c05821c8fa5 100644
--- a/Framework/Utilities/testBit.cc
+++ b/Framework/Utilities/testBit.cc
@@ -10,18 +10,17 @@
*/
#include <corsika/utl/Test.h>
+#include <cppunit/extensions/HelperMacros.h>
#include <tst/Verify.h>
#include <utl/Bit.h>
-#include <cppunit/extensions/HelperMacros.h>
+#include <bitset>
#include <cstdio>
#include <iostream>
-#include <bitset>
using namespace tst;
using namespace utl;
using namespace std;
-
/**
\ingroup testing
*/
@@ -34,14 +33,12 @@ class TestBit : public CppUnit::TestFixture {
CPPUNIT_TEST_SUITE_END();
public:
- void setUp() { }
+ void setUp() {}
- void tearDown() { }
+ void tearDown() {}
- void
- TestGet()
- {
- const int size = sizeof(int)*8;
+ void TestGet() {
+ const int size = sizeof(int) * 8;
const int bc2 = 12345;
int b2 = bc2;
bitset<size> b1(bc2);
@@ -59,23 +56,18 @@ public:
CPPUNIT_ASSERT(Verify<Equal>(out1.str(), out3.str()));
}
- void
- TestSet()
- {
- const int size = sizeof(int)*8;
+ void TestSet() {
+ const int size = sizeof(int) * 8;
const int number = 12345;
bitset<size> b1(number);
int b2 = 11111;
- for (int i = 0; i < size; ++i)
- AsBitArray(b2)[i] = b1[i];
+ for (int i = 0; i < size; ++i) AsBitArray(b2)[i] = b1[i];
CPPUNIT_ASSERT(Verify<Equal>(b2, number));
}
- void
- TestMask()
- {
+ void TestMask() {
const int n = (1 << 18) | (1 << 5);
int m = 0;
@@ -83,19 +75,16 @@ public:
AsBitArray(m)[5] = true;
CPPUNIT_ASSERT(Verify<Equal>(n, m));
- for (unsigned int i = 0; i < 8*sizeof(int); ++i)
- AsBitArray(m)[i] = 0;
+ for (unsigned int i = 0; i < 8 * sizeof(int); ++i) AsBitArray(m)[i] = 0;
CPPUNIT_ASSERT(Verify<Equal>(m, 0));
m = 1;
AsBitArray(m).Mask(n, true);
- CPPUNIT_ASSERT(Verify<Equal>(m, n+1));
+ CPPUNIT_ASSERT(Verify<Equal>(m, n + 1));
AsBitArray(m).Mask(n, false);
CPPUNIT_ASSERT(Verify<Equal>(m, 1));
}
-
};
-
CPPUNIT_TEST_SUITE_REGISTRATION(TestBit);
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index c8f8557b074a16b0d152f84bc7e96273b10d063b..41a8277d9c736f387d68b255ff2af7415133c10c 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -2,9 +2,11 @@
add_subdirectory (NullModel)
add_subdirectory (Sibyll)
add_subdirectory (StackInspector)
+add_subdirectory (TrackingLine)
#add_custom_target(CORSIKAprocesses)
add_library (CORSIKAprocesses INTERFACE)
add_dependencies(CORSIKAprocesses ProcessNullModel)
add_dependencies(CORSIKAprocesses ProcessSibyll)
add_dependencies(CORSIKAprocesses ProcessStackInspector)
+add_dependencies(CORSIKAprocesses ProcessTrackingLine)
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index 9d7227ce97b324ba33977ce351c476cb1a2fe0c9..35bfbc617df904ba0e859fbf451ab1356092f207 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -35,21 +35,23 @@ namespace corsika::process::sibyll {
SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) {
//~ assert(handledBySibyll(pCode));
- return static_cast<SibyllCode>(corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
+ return static_cast<SibyllCode>(
+ corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
}
-
+
corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) {
return sibyll2corsika[static_cast<SibyllCodeIntType>(pCode) - minSibyll];
}
- int ConvertToSibyllRaw(corsika::particles::Code pCode){
- return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>( corsika::process::sibyll::ConvertToSibyll( pCode ) );
+ int ConvertToSibyllRaw(corsika::particles::Code pCode) {
+ return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>(
+ corsika::process::sibyll::ConvertToSibyll(pCode));
}
int GetSibyllXSCode(corsika::particles::Code pCode) {
return corsika2sibyllXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
}
-
+
} // namespace corsika::process::sibyll
#endif
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index 77c1f5dd8216ad48ebf5c4aa476ef57695bb7fea..a593a1f2ce5a5c3a281c758e2404c72347286ec0 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -159,8 +159,6 @@ if __name__ == "__main__":
pythia_db = load_pythiadb(sys.argv[1])
read_sibyll_codes(sys.argv[2], pythia_db)
- print (str(pythia_db))
-
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
print(generate_sibyll_enum(pythia_db), file=f)
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 12917bfd91cf85d9f3f310256386f70bcfce483d..67c950e5585f831f8a95ae4f0e00f59df95c8919 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -50,7 +50,7 @@ TEST_CASE("Sibyll", "[processes]") {
REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::Electron::GetCode()));
REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::SigmaC0::GetCode()));
}
-
+
SECTION("cross-section type") {
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::Electron) == 0);
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 5f8471177ae2dc351a440693937a3ee3e3b756da..b0acefc4a86c72acba373feafc82b7434aa2b9a0 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -9,11 +9,14 @@
* the license.
*/
+#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/logging/Logger.h>
+#include <corsika/setup/SetupTrajectory.h>
+
#include <iostream>
using namespace std;
@@ -21,18 +24,16 @@ using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
-template <typename Stack, typename Trajectory>
-StackInspector<Stack, Trajectory>::StackInspector(const bool aReport)
+template <typename Stack>
+StackInspector<Stack>::StackInspector(const bool aReport)
: fReport(aReport) {}
-template <typename Stack, typename Trajectory>
-StackInspector<Stack, Trajectory>::~StackInspector() {}
-
-template <typename Stack, typename Trajectory>
-process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle&,
- Trajectory&,
- Stack& s) const {
+template <typename Stack>
+StackInspector<Stack>::~StackInspector() {}
+template <typename Stack>
+process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
+ Stack& s) const {
static int countStep = 0;
if (!fReport) return EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
@@ -41,24 +42,27 @@ process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle
for (auto& iterP : s) {
EnergyType E = iterP.GetEnergy();
Etot += E;
- std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()<< ", en=" << iterP.GetEnergy() / 1_GeV << " | " << std::endl;
+ geometry::CoordinateSystem& rootCS =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCS(); // for printout
+ auto pos = iterP.GetPosition().GetCoordinates(rootCS);
+ cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
+ << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
+ << " pos=" << pos << endl;
}
countStep++;
- // cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl;
- cout << "call: " << countStep << " stack size: " << s.GetSize() << " energy: " << Etot / 1_GeV << endl;
+ cout << "StackInspector: nStep=" << countStep << " stackSize=" << s.GetSize()
+ << " Estack=" << Etot / 1_GeV << " GeV" << endl;
return EProcessReturn::eOk;
}
-template <typename Stack, typename Trajectory>
-double StackInspector<Stack, Trajectory>::MinStepLength(Particle&) const {
- return 0;
+template <typename Stack>
+void StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
+ // return 0;
}
-template <typename Stack, typename Trajectory>
-void StackInspector<Stack, Trajectory>::Init() {}
+template <typename Stack>
+void StackInspector<Stack>::Init() {}
#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-template class corsika::process::stack_inspector::StackInspector<setup::Stack,
- setup::Trajectory>;
+template class process::stack_inspector::StackInspector<setup::Stack>;
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 118c06d36a562dc8ed703ebca0eb6ffb3489a5ec..e40517499d7c07b984c5d9b30caac469269cf649 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -14,13 +14,16 @@
#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+
namespace corsika::process {
namespace stack_inspector {
- template <typename Stack, typename Trajectory>
+ template <typename Stack>
class StackInspector
- : public corsika::process::ContinuousProcess<StackInspector<Stack, Trajectory>> {
+ : public corsika::process::ContinuousProcess<StackInspector<Stack>> {
typedef typename Stack::ParticleType Particle;
@@ -31,10 +34,10 @@ namespace corsika::process {
void Init();
// template <typename Particle, typename Trajectory, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack& s) const;
+ EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
- // template <typename Particle>
- double MinStepLength(Particle&) const;
+ // template <typename Particle>
+ void MinStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2911c9b787787f43ba837e64fd67c89921c699b0
--- /dev/null
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -0,0 +1,38 @@
+
+set (
+ MODEL_HEADERS
+ TrackingLine.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/tracking_line
+ )
+
+add_library (ProcessTrackingLine INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLine ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+
+target_include_directories (
+ ProcessTrackingLine
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
+
+
+# # --------------------
+# # code unit testing
+# add_executable (testStackInspector testStackInspector.cc)
+
+# target_link_libraries (
+# testStackInspector
+# CORSIKAunits
+# CORSIKAthirdparty # for catch2
+# )
+
+# add_test (NAME testStackInspector COMMAND testStackInspector)
+
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
new file mode 100644
index 0000000000000000000000000000000000000000..602d5473269c8e15561967fe9adc8b16ff4dfa0c
--- /dev/null
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -0,0 +1,35 @@
+#ifndef _include_corsika_processes_TrackinLine_h_
+#define _include_corsika_processes_TrackinLine_h_
+
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/Point.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using namespace corsika;
+
+namespace corsika::process {
+
+ namespace tracking_line {
+
+ template <typename Stack>
+ class TrackingLine { // Newton-step, naja.. not yet
+ typedef typename Stack::ParticleType Particle;
+
+ public:
+ void Init() {}
+ setup::Trajectory GetTrack(Particle& p) {
+ geometry::Vector<SpeedType::dimension_type> v = p.GetDirection();
+ geometry::Line traj(p.GetPosition(), v);
+ return geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
+ }
+ };
+
+ } // namespace stack_inspector
+
+} // namespace corsika::process
+
+#endif
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index d7f5c482d207ea04cfae1365ed7f11af13541f7f..f1dfd497a6bec31eeb9e4d84967a8b80f380ebb0 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -12,12 +12,53 @@
#ifndef _corsika_setup_setuptrajectory_h_
#define _corsika_setup_setuptrajectory_h_
+#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Trajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <variant>
+
namespace corsika::setup {
- typedef corsika::geometry::Trajectory<corsika::geometry::Line> Trajectory;
-}
+ using corsika::geometry::Helix;
+ using corsika::geometry::Line;
+
+ /// definition of Trajectory base class, to be used in tracking and cascades
+ typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
+ corsika::geometry::Trajectory<Helix>>
+ Trajectory;
+
+ /// helper visitor to modify Particle by moving along Trajectory
+ template <typename Particle>
+ class ParticleUpdate {
+
+ Particle& fP;
+
+ public:
+ ParticleUpdate(Particle& p)
+ : fP(p) {}
+ void operator()(std::monostate const&) {}
+
+ template <typename T>
+ void operator()(T const& trajectory) {
+ fP.SetPosition(trajectory.GetPosition(1));
+ }
+ };
+
+ /// helper visitor to modify Particle by moving along Trajectory
+ class GetDuration {
+ public:
+ corsika::units::si::TimeType operator()(std::monostate const&) {
+ return 0 * corsika::units::si::second;
+ }
+ template <typename T>
+ corsika::units::si::TimeType operator()(T const& trajectory) {
+ return trajectory.GetDuration();
+ }
+ };
+
+} // namespace corsika::setup
#endif
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 2b03c54431ca71335b026c5b892a704eb9ebf418..88ec5d58baa4d6d0da208628a41dbdf92e28e569 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -16,53 +16,66 @@
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/geometry/CoordinateSystem.h> // remove
#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h> // remove
#include <corsika/geometry/Vector.h>
-#include <vector>
#include <algorithm>
+#include <vector>
+
+using namespace corsika;
namespace corsika::stack {
namespace super_stupid {
+ using corsika::geometry::Point;
+ using corsika::geometry::Vector;
using corsika::particles::Code;
+ using corsika::units::si::energy_d;
using corsika::units::si::EnergyType;
- using corsika::units::si::TimeType;
- using corsika::units::si::second;
- using corsika::units::si::meter;
using corsika::units::si::joule;
- using corsika::units::si::newton_second;
- using corsika::units::si::energy_d;
+ using corsika::units::si::meter;
using corsika::units::si::momentum_d;
- using corsika::geometry::Point;
- using corsika::geometry::Vector;
+ using corsika::units::si::newton_second;
+ using corsika::units::si::second;
+ using corsika::units::si::SpeedType;
+ using corsika::units::si::TimeType;
- typedef Vector<momentum_d> MomentumVector;
+ typedef Vector<momentum_d> MomentumVector;
/**
* Example of a particle object on the stack.
*/
template <typename StackIteratorInterface>
- class ParticleInterface : public ParticleBase<StackIteratorInterface> {
-
+ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
+
using ParticleBase<StackIteratorInterface>::GetStackData;
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
void SetPID(const Code id) { GetStackData().SetPID(GetIndex(), id); }
void SetEnergy(const EnergyType& e) { GetStackData().SetEnergy(GetIndex(), e); }
- void SetMomentum(const MomentumVector& v) { GetStackData().SetMomentum(GetIndex(), v); }
+ void SetMomentum(const MomentumVector& v) {
+ GetStackData().SetMomentum(GetIndex(), v);
+ }
void SetPosition(const Point& v) { GetStackData().SetPosition(GetIndex(), v); }
void SetTime(const TimeType& v) { GetStackData().SetTime(GetIndex(), v); }
Code GetPID() const { return GetStackData().GetPID(GetIndex()); }
EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
- MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
+ MomentumVector GetMomentum() const {
+ return GetStackData().GetMomentum(GetIndex());
+ }
Point GetPosition() const { return GetStackData().GetPosition(GetIndex()); }
TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); }
+
+#warning this does not really work, nor make sense:
+ Vector<SpeedType::dimension_type> GetDirection() const {
+ auto P = GetMomentum();
+ return P / P.norm() * 1e10 * (units::si::meter / units::si::second);
+ }
};
/**
@@ -82,13 +95,13 @@ namespace corsika::stack {
int GetSize() const { return fDataPID.size(); }
int GetCapacity() const { return fDataPID.size(); }
-
+
void SetPID(const int i, const Code id) { fDataPID[i] = id; }
void SetEnergy(const int i, const EnergyType e) { fDataE[i] = e; }
void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; }
void SetPosition(const int i, const Point& v) { fPosition[i] = v; }
void SetTime(const int i, const TimeType& v) { fTime[i] = v; }
-
+
Code GetPID(const int i) const { return fDataPID[i]; }
EnergyType GetEnergy(const int i) const { return fDataE[i]; }
MomentumVector GetMomentum(const int i) const { return fMomentum[i]; }
@@ -101,41 +114,41 @@ namespace corsika::stack {
void Copy(const int i1, const int i2) {
fDataPID[i2] = fDataPID[i1];
fDataE[i2] = fDataE[i1];
- fMomentum[i2] = fMomentum[i1];
- fPosition[i2] = fPosition[i1];
- fTime[i2] = fTime[i1];
+ fMomentum[i2] = fMomentum[i1];
+ fPosition[i2] = fPosition[i1];
+ fTime[i2] = fTime[i1];
}
/**
* Function to copy particle at location i2 in stack to i1
*/
void Swap(const int i1, const int i2) {
- std::swap(fDataPID[i2], fDataPID[i1]);
- std::swap(fDataE[i2], fDataE[i1]);
- std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!!
- std::swap(fPosition[i2], fPosition[i1]);
- std::swap(fTime[i2], fTime[i1]);
+ std::swap(fDataPID[i2], fDataPID[i1]);
+ std::swap(fDataE[i2], fDataE[i1]);
+ std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!!
+ std::swap(fPosition[i2], fPosition[i1]);
+ std::swap(fTime[i2], fTime[i1]);
}
protected:
void IncrementSize() {
fDataPID.push_back(Code::Unknown);
fDataE.push_back(0 * joule);
-#warning this here makes no sense: see issue #48
- auto const dummyCS = corsika::geometry::CoordinateSystem::CreateRootCS();
- fMomentum.push_back(MomentumVector(dummyCS,
- {0 * newton_second, 0 * newton_second, 0 * newton_second}));
- fPosition.push_back(Point(dummyCS,
- {0 * meter, 0 * meter, 0 * meter}));
- fTime.push_back(0 * second);
+ //#TODO this here makes no sense: see issue #48
+ geometry::CoordinateSystem& dummyCS =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+ fMomentum.push_back(MomentumVector(
+ dummyCS, {0 * newton_second, 0 * newton_second, 0 * newton_second}));
+ fPosition.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
+ fTime.push_back(0 * second);
}
void DecrementSize() {
if (fDataE.size() > 0) {
- fDataPID.pop_back();
+ fDataPID.pop_back();
fDataE.pop_back();
- fMomentum.pop_back();
- fPosition.pop_back();
- fTime.pop_back();
+ fMomentum.pop_back();
+ fPosition.pop_back();
+ fTime.pop_back();
}
}
@@ -144,7 +157,7 @@ namespace corsika::stack {
std::vector<Code> fDataPID;
std::vector<EnergyType> fDataE;
- std::vector<MomentumVector> fMomentum;
+ std::vector<MomentumVector> fMomentum;
std::vector<Point> fPosition;
std::vector<TimeType> fTime;
diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc
index bc0814b2e91c2096494dd43faef8a17a9503efdd..28cb891bdc8dfead5e3b04e8113d46f7bc127108 100644
--- a/Stack/SuperStupidStack/testSuperStupidStack.cc
+++ b/Stack/SuperStupidStack/testSuperStupidStack.cc
@@ -9,6 +9,7 @@
* the license.
*/
+#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/stack/super_stupid/SuperStupidStack.h>
#include <corsika/units/PhysicalUnits.h>
@@ -22,46 +23,44 @@ using namespace corsika::units::si;
using namespace corsika;
using namespace corsika::stack::super_stupid;
-
#include <iostream>
using namespace std;
-
TEST_CASE("SuperStupidStack", "[stack]") {
SECTION("read+write") {
SuperStupidStack s;
auto p = s.NewParticle();
- p.SetPID(corsika::particles::Code::Electron);
+ p.SetPID(particles::Code::Electron);
p.SetEnergy(1.5_GeV);
- auto const dummyCS = corsika::geometry::CoordinateSystem::CreateRootCS();
- p.SetMomentum(MomentumVector(dummyCS, {1 * newton_second, 1 * newton_second, 1 * newton_second}));
+ geometry::CoordinateSystem& dummyCS =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+ p.SetMomentum(MomentumVector(
+ dummyCS, {1 * newton_second, 1 * newton_second, 1 * newton_second}));
p.SetPosition(Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));
p.SetTime(100_s);
-
+
// read
REQUIRE(s.GetSize() == 1);
auto pout = s.GetNextParticle();
- REQUIRE(pout.GetPID() == corsika::particles::Code::Electron);
+ REQUIRE(pout.GetPID() == particles::Code::Electron);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
#warning Fix the next two lines:
- //REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 * joule}));
- //REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));
+ // REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 *
+ // joule})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter, 1 *
+ // meter}));
REQUIRE(pout.GetTime() == 100_s);
}
-
+
SECTION("write+delete") {
SuperStupidStack s;
- for (int i=0; i<99; ++i)
- s.NewParticle();
+ for (int i = 0; i < 99; ++i) s.NewParticle();
REQUIRE(s.GetSize() == 99);
- for (int i=0; i<99; ++i)
- s.GetNextParticle().Delete();
-
+ for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete();
REQUIRE(s.GetSize() == 0);
}