diff --git a/#.gitlab-ci.yml# b/#.gitlab-ci.yml#
new file mode 100644
index 0000000000000000000000000000000000000000..fb8f6a95b188b27e5e77d969a4c244a8d3788dde
--- /dev/null
+++ b/#.gitlab-ci.yml#
@@ -0,0 +1,37 @@
+image: ubuntu:bionic
+
+variables:
+ GIT_SSL_NO_VERIFY: "1"
+
+before_script:
+ - apt-get update --yes
+ - apt-get install --yes cmake libboost-dev libeigen3-dev python3 gfortran
+
+build:
+ stage: build
+ tags:
+ - run1
+ script:
+ - mkdir build
+ - cd build
+ - cmake ..
+ - cmake --build .
+ - ctest -V
+
+# code_quality:
+# image: docker:stable
+# variables:
+# DOCKER_DRIVER: overlay2
+# allow_failure: true
+# services:
+# - docker:stable-dind
+# script:
+# - export SP_VERSION=$(echo "$CI_SERVER_VERSION" | sed 's/^\([0-9]*\)\.\([0-9]*\).*/\1-\2-stable/')
+# - docker run
+# --env SOURCE_CODE="$PWD"
+# --volume "$PWD":/code
+# --volume /var/run/docker.sock:/var/run/docker.sock
+# "registry.gitlab.com/gitlab-org/security-products/codequality:$SP_VERSION" /code
+# artifacts:
+# reports:
+# codequality: gl-code-quality-report.json
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 567dd5ce531516ee2446f4b3443fc1daa9f30257..d2d4df9b12cb9bc92a2251bc99abda96b8507a0a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -28,7 +28,7 @@ if(EXISTS "${CMAKE_SOURCE_DIR}/.git")
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.")
+ message(STATUS "Setting build type to '${default_build_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
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index e221b28fbc3b5d977cd6ae8ac193708e6e149b7a..60e436a4e7574f233d642bd6673a635a6db25823 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -119,8 +119,10 @@ public:
template <typename Particle>
LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
+ cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
+ cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
const Code pid = p.GetPID();
- if (isEmParticle(pid) || isInvisible(pid)) {
+ if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
return 0_m;
} else {
@@ -134,40 +136,26 @@ public:
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
const Code pid = p.GetPID();
+ EProcessReturn ret = EProcessReturn::eOk;
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += p.GetEnergy();
fEmCount += 1;
p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += p.GetEnergy();
fInvCount += 1;
p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy();
p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
}
- // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
- // cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
- // cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
- // const Code pid = p.GetPID();
- // if( isEmParticle( pid ) ){
- // cout << "removing em. particle..." << endl;
- // fEmEnergy += p.GetEnergy();
- // fEmCount += 1;
- // p.Delete();
- // return EProcessReturn::eParticleAbsorbed;
- // }
- // if ( isInvisible( pid ) ){
- // cout << "removing inv. particle..." << endl;
- // fInvEnergy += p.GetEnergy();
- // fInvCount += 1;
- // p.Delete();
- // return EProcessReturn::eParticleAbsorbed;
- // }
- return EProcessReturn::eOk;
+ return ret;
}
void Init() {
@@ -182,10 +170,11 @@ public:
void ShowResults() {
cout << " ******************************" << endl
<< " ParticleCut: " << endl
- << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
- << " no. of em. particles injected: " << fEmCount << endl
- << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
- << " no. of inv. particles injected: " << fInvCount << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
<< " ******************************" << endl;
}
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 23954fef70ffe62140074010a1207624b1cfcaa2..24fd807d2723343be6a827f0783fa8fd5fdd0383 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -42,13 +42,23 @@ namespace corsika::process {
}
}
+ void setUnstable(const corsika::particles::Code pCode) const {
+ int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
+ }
+
+ void setStable(const corsika::particles::Code pCode) const {
+ int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ }
+
void setAllStable() const {
+ // name? also makes EM particles stable
- // using namespace corsika::io;
using std::cout;
using std::endl;
- // name? also makes EM particles stable
+ std::cout << "Decay: setting all particles stable.." << std::endl;
// loop over all particles in sibyll
// should be changed to loop over human readable list
@@ -58,11 +68,7 @@ namespace corsika::process {
const int sibCode = (int)p;
// skip unknown and antiparticles
if (sibCode < 1) continue;
- const corsika::particles::Code pid =
- ConvertFromSibyll(static_cast<SibyllCode>(sibCode));
- std::cout << "Sibyll: Decay: setting " << pid << " stable" << std::endl;
s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
- std::cout << "decay table value: " << s_csydec_.idb[sibCode - 1] << std::endl;
}
}
@@ -124,6 +130,7 @@ namespace corsika::process {
corsika::particles::GetLifetime(p.GetPID());
cout << "Decay: code: " << p.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
+ cout << "Decay: MinStep: energy: " << E / 1_GeV << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
// cout << "Decay: MinStep: density: " << density << endl;
// return as column density
@@ -139,20 +146,29 @@ namespace corsika::process {
template <typename Particle, typename Stack>
void DoDecay(Particle& p, Stack& s) const {
+ using corsika::geometry::Point;
+ using namespace corsika::units::si;
SibStack ss;
ss.Clear();
// copy particle to sibyll stack
auto pin = ss.NewParticle();
- pin.SetPID(process::sibyll::ConvertToSibyllRaw(p.GetPID()));
+ const corsika::particles::Code pCode = p.GetPID();
+ pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
pin.SetEnergy(p.GetEnergy());
pin.SetMomentum(p.GetMomentum());
+ // remember position
+ Point decayPoint = p.GetPosition();
+ TimeType t0 = p.GetTime();
// remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable
- setHadronsUnstable();
+ // setHadronsUnstable();
+ setUnstable(pCode);
// call sibyll decay
std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
decsib_();
+ // reset to stable
+ setStable(pCode);
// print output
int print_unit = 6;
sib_list_(print_unit);
@@ -169,6 +185,8 @@ namespace corsika::process {
pnew.SetEnergy(psib.GetEnergy());
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
pnew.SetMomentum(psib.GetMomentum());
+ pnew.SetPosition(decayPoint);
+ pnew.SetTime(t0);
}
// empty sibyll stack
ss.Clear();
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index a82dac38034c26a6373fe6ff07a4a278dcdb400e..fbaf7c5f519391af653da149d7c9aa3d173d6b0e 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -28,12 +28,6 @@ namespace corsika::process::sibyll {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm");
- // test random number generator
- std::cout << "Interaction: "
- << " 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
sibyll_ini_();
}
@@ -56,9 +50,8 @@ namespace corsika::process::sibyll {
// 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);
+ const int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
+ const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/*
the target should be defined by the Environment,
@@ -67,7 +60,7 @@ namespace corsika::process::sibyll {
*/
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
- int kTarget = 16;
+ const int kTarget = corsika::particles::Oxygen::GetNucleusA();
hep::EnergyType Etot =
p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
@@ -77,8 +70,8 @@ namespace corsika::process::sibyll {
Ptot += p.GetMomentum();
Ptot += pTarget;
// calculate cm. energy
- hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
- double Ecm = sqs / 1_GeV;
+ const hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
+ const double Ecm = sqs / 1_GeV;
std::cout << "Interaction: "
<< "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
@@ -149,8 +142,8 @@ namespace corsika::process::sibyll {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
- Point pOrig(rootCS, coordinates);
+ Point pOrig = p.GetPosition();
+ TimeType tOrig = p.GetTime();
/*
the target should be defined by the Environment,
@@ -160,8 +153,9 @@ namespace corsika::process::sibyll {
here we need: GetTargetMassNumber() or GetTargetPID()??
GetTargetMomentum() (zero in EAS)
*/
- // FOR NOW: set target to proton
- int kTarget = 1; // env.GetTargetParticle().GetPID();
+ // FOR NOW: set target to oxygen
+ const int kTarget = corsika::particles::Oxygen::
+ GetNucleusA(); // env.GetTargetParticle().GetPID();
cout << "defining target momentum.." << endl;
// FOR NOW: target is always at rest
@@ -171,6 +165,8 @@ namespace corsika::process::sibyll {
<< pTarget.GetComponents() / 1_GeV * constants::c << endl;
cout << "beam momentum (GeV/c): "
<< p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
+ cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
+ cout << "time: " << tOrig << endl;
// get energy of particle from stack
/*
@@ -184,7 +180,7 @@ namespace corsika::process::sibyll {
EnergyType E = p.GetEnergy();
EnergyType Etot = E + Etarget;
// total momentum
- MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+ MomentumVector Ptot = p.GetMomentum();
// invariant mass, i.e. cm. energy
EnergyType Ecm =
sqrt(Etot * Etot - Ptot.squaredNorm()); // sqrt( 2. * E * 0.93827_GeV );
@@ -208,8 +204,8 @@ namespace corsika::process::sibyll {
<< " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
- << " DoInteraction: dropping particle.." << std::endl;
- p.Delete();
+ << " DoInteraction: should have dropped particle.." << std::endl;
+ // p.Delete(); delete later... different process
} else {
// Sibyll does not know about units..
double sqs = Ecm / 1_GeV;
@@ -265,6 +261,8 @@ namespace corsika::process::sibyll {
corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]};
pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
+ pnew.SetPosition(pOrig);
+ pnew.SetTime(tOrig);
Ptot_final += pnew.GetMomentum();
}
// cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
index 5d2d1e69644e812246b731aa4bf276d090d79910..0a0ecc1133ce2fce469e93747368ab9224568449 100644
--- a/Processes/Sibyll/sibyll2.3c.h
+++ b/Processes/Sibyll/sibyll2.3c.h
@@ -59,7 +59,7 @@ extern struct {
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
// sibyll main subroutine
-void sibyll_(int&, int&, double&);
+void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
@@ -79,8 +79,9 @@ 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&);
+void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&);
+void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
+ double&);
double s_rndm_(int&);
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 3a12c052fcd833adf7254b139bac311c8a977893..3e9bfe831d1dfbe3d3d16fb1f67347dcaf0e51fa 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -74,12 +74,9 @@ namespace corsika::stack {
return GetStackData().GetTime(GetIndex());
}
- //#warning this does not really work, nor makes sense:
corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>
GetDirection() const {
- auto P = GetMomentum();
- return P / P.norm() * 1e10 *
- (corsika::units::si::meter / corsika::units::si::second);
+ return GetMomentum() / GetEnergy() * corsika::units::constants::c;
}
};