diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9d33bf5e4f9fa0a6785c4fa50b73f1033aaaf80e..ad27460006921a1e4c0a5e7334510dedc0349d3b 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -16,7 +16,7 @@ build:
- cd build
- cmake ..
- cmake --build .
- - ctest -V
+ - ctest
# code_quality:
# image: docker:stable
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index cc37b2131566ad1ad39afcba3311f66532298154..232e7d002945f41a613266c0738015a2dd34717e 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -45,24 +45,24 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::hep;
-static EnergyType fEnergy = 0. * 1_GeV;
+class ProcessCut : public corsika::process::ContinuousProcess<ProcessCut> {
-// FOR NOW: global static variables for ParticleCut process
-// this is just wrong...
-static EnergyType fEmEnergy;
-static int fEmCount;
+ EnergyType fECut;
-static EnergyType fInvEnergy;
-static int fInvCount;
+ mutable EnergyType fEnergy = 0_GeV;
+ mutable EnergyType fEmEnergy = 0_GeV;
+ mutable int fEmCount = 0;
+ mutable EnergyType fInvEnergy = 0_GeV;
+ mutable int fInvCount = 0;
-class ProcessEMCut : public corsika::process::ContinuousProcess<ProcessEMCut> {
public:
- ProcessEMCut() {}
+ ProcessCut(const EnergyType v)
+ : fECut(v) {}
template <typename Particle>
bool isBelowEnergyCut(Particle& p) const {
// FOR NOW: center-of-mass energy hard coded
const EnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
- if (p.GetEnergy() < 50_GeV || Ecm < 10_GeV)
+ if (p.GetEnergy() < fECut || Ecm < 10_GeV)
return true;
else
return false;
@@ -134,25 +134,27 @@ public:
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
- cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
const Code pid = p.GetPID();
+ EnergyType energy = p.GetEnergy();
+ cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
+ << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
EProcessReturn ret = EProcessReturn::eOk;
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
- fEmEnergy += p.GetEnergy();
+ fEmEnergy += energy;
fEmCount += 1;
- p.Delete();
+ // p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
- fInvEnergy += p.GetEnergy();
+ fInvEnergy += energy;
fInvCount += 1;
- p.Delete();
+ // p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
- fEnergy += p.GetEnergy();
- p.Delete();
+ fEnergy += energy;
+ // p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
}
return ret;
@@ -178,20 +180,21 @@ public:
<< " ******************************" << endl;
}
- EnergyType GetInvEnergy() { return fInvEnergy; }
-
- EnergyType GetCutEnergy() { return fEnergy; }
-
- EnergyType GetEmEnergy() { return fEmEnergy; }
-
-private:
+ EnergyType GetInvEnergy() const { return fInvEnergy; }
+ EnergyType GetCutEnergy() const { return fEnergy; }
+ EnergyType GetEmEnergy() const { return fEmEnergy; }
};
+//
+// The example main program for a particle cascade
+//
int main() {
+ // initialize random number sequence(s)
corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
- corsika::environment::Environment env; // dummy environment
+ // setup environment, geometry
+ corsika::environment::Environment env;
auto& universe = *(env.GetUniverse());
auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
@@ -201,47 +204,47 @@ int main() {
using MyHomogeneousModel =
corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
- 1_g / (1_m * 1_m * 1_m),
+ 1_kg / (1_m * 1_m * 1_m),
corsika::environment::NuclearComposition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+ std::vector<corsika::particles::Code>{corsika::particles::Code::Oxygen},
std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
- CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ // setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
-
corsika::process::sibyll::Interaction sibyll;
corsika::process::sibyll::Decay decay;
- ProcessEMCut cut;
+ ProcessCut cut(8_GeV);
+
+ // assemble all processes into an ordered process list
const auto sequence = /*p0 <<*/ sibyll << decay << cut;
+ // setup particle stack, and add primary particle
setup::Stack stack;
+ stack.Clear();
+ const hep::EnergyType E0 = 1_TeV;
+ {
+ auto particle = stack.NewParticle();
+ particle.SetPID(Code::Proton);
+ hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
+ auto plab = stack::super_stupid::MomentumVector(rootCS, 0_GeV, 0_GeV, -P0);
+ particle.SetEnergy(E0);
+ particle.SetMomentum(plab);
+ particle.SetTime(0_ns);
+ Point p(rootCS, 0_m, 0_m, 10_km);
+ particle.SetPosition(p);
+ }
+ // define air shower object, run simulation
corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
-
- stack.Clear();
- auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV;
- hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
- auto plab = stack::super_stupid::MomentumVector(rootCS, 0. * 1_GeV, 0. * 1_GeV, P0);
- particle.SetEnergy(E0);
- particle.SetMomentum(plab);
- particle.SetPID(Code::Proton);
- particle.SetTime(0_ns);
- Point p(rootCS, 0_m, 0_m, 0_m);
- particle.SetPosition(p);
EAS.Init();
EAS.Run();
- cout << "Result: E0="
- << E0 / 1_GeV
- //<< "GeV, particles below energy threshold =" << p1.GetCount()
- << endl;
- cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
- << std::endl;
+
+ cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
cout << "total energy (GeV): "
<< (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy()) / 1_GeV << endl;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index c2815c3ed4578804cf3619b2c496298fc6f29e64..275f9702e7e5cacc7c7448b01aae14b9e1a19744 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -71,9 +71,8 @@ namespace corsika::cascade {
fProcessSequence.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length in grammage
- auto constexpr grammageConversion = 1_g / (1_m * 1_m);
- std::exponential_distribution expDist(1 / (grammageConversion * total_inv_lambda));
- GrammageType const next_interact = grammageConversion * expDist(fRNG);
+ std::exponential_distribution expDist(total_inv_lambda * (1_g / (1_m * 1_m)));
+ GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
@@ -100,7 +99,6 @@ namespace corsika::cascade {
<< ", next_decay=" << next_decay << std::endl;
// convert next_decay from time to length [m]
- // Environment::GetDistance(step, next_decay);
LengthType const distance_decay = next_decay * particle.GetMomentum().norm() /
particle.GetEnergy() *
corsika::units::constants::c;
@@ -109,10 +107,11 @@ namespace corsika::cascade {
auto const min_distance =
std::min({distance_interact, distance_decay, distance_max});
+ std::cout << " move particle by : " << min_distance << std::endl;
+
// here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.PositionFromArclength(min_distance));
-
// .... also update time, momentum, direction, ...
// apply all continuous processes on particle + track
@@ -120,39 +119,40 @@ namespace corsika::cascade {
fProcessSequence.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
- // fStack.Delete(particle); // TODO: check if this is really needed
- } else {
-
- std::cout << "sth. happening before geometric limit ?"
- << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
-
- if (min_distance < distance_max) { // interaction to happen within geometric limit
- // check weather decay or interaction limits this step
-
- if (min_distance == distance_interact) {
- std::cout << "collide" << std::endl;
-
- InverseGrammageType const actual_inv_length =
- fProcessSequence.GetTotalInverseInteractionLength(particle, step);
-
- corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
- actual_inv_length);
- const auto sample_process = uniDist(fRNG);
- InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
- fProcessSequence.SelectInteraction(particle, fStack, sample_process,
- inv_lambda_count);
- } else {
- std::cout << "decay" << std::endl;
- InverseTimeType const actual_decay_time =
- fProcessSequence.GetTotalInverseLifetime(particle);
-
- corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
- actual_decay_time);
- const auto sample_process = uniDist(fRNG);
- InverseTimeType inv_decay_count = 0 / second;
- fProcessSequence.SelectDecay(particle, fStack, sample_process,
- inv_decay_count);
- }
+ std::cout << "Cascade: delete absorbed particle " << particle.GetPID() << " "
+ << particle.GetEnergy() / 1_GeV << "GeV" << std::endl;
+ particle.Delete();
+ return;
+ }
+
+ std::cout << "sth. happening before geometric limit ?"
+ << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
+
+ if (min_distance < distance_max) { // interaction to happen within geometric limit
+ // check weather decay or interaction limits this step
+
+ if (min_distance == distance_interact) {
+ std::cout << "collide" << std::endl;
+
+ InverseGrammageType const actual_inv_length =
+ fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+
+ corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
+ actual_inv_length);
+ const auto sample_process = uniDist(fRNG);
+ InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
+ fProcessSequence.SelectInteraction(particle, fStack, sample_process,
+ inv_lambda_count);
+ } else {
+ std::cout << "decay" << std::endl;
+ InverseTimeType const actual_decay_time =
+ fProcessSequence.GetTotalInverseLifetime(particle);
+
+ corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
+ actual_decay_time);
+ const auto sample_process = uniDist(fRNG);
+ InverseTimeType inv_decay_count = 0 / second;
+ fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count);
}
}
}
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 80ebf4b21b1999754a2ccf2d244e4dfd7e5e0b4f..081325f96038dc8c733f4599603e0ee97a22d04f 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -80,7 +80,7 @@ public:
}
template <typename Particle, typename T, typename Stack>
- void DoContinuous(Particle& p, T&, Stack& s) const {
+ EProcessReturn DoContinuous(Particle& p, T&, Stack& s) const {
EnergyType E = p.GetEnergy();
if (E < 85_MeV) {
p.Delete();
@@ -95,6 +95,7 @@ public:
pnew.SetPosition(p.GetPosition());
pnew.SetMomentum(p.GetMomentum());
}
+ return EProcessReturn::eOk;
}
void Init() { fCount = 0; }
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index 19effade79572ef87d0091f9e86f416dbfdac30e..dd9572ac94cbb8b24352df1c77a22227679ea1c2 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -35,7 +35,7 @@ namespace corsika::process {
// -> enforce derived to implement DoContinuous...
template <typename Particle, typename Track, typename Stack>
EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
-
+
// -> enforce derived to implement MaxStepLength...
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const;
diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h
index ae4869ea84ba228e1b59b2fe89778f2e932efc25..215994b42fb74daaba3b382513d94e50a45d430b 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -20,13 +20,25 @@ namespace corsika::process {
that can be accumulated easily with "|="
*/
- enum class EProcessReturn {
- eOk = 1,
- eParticleAbsorbed = 2,
- eInteracted = 3,
- eDecayed = 4,
+ enum class EProcessReturn : int {
+ eOk = (1 << 0),
+ eParticleAbsorbed = (1 << 2),
+ eInteracted = (1 << 3),
+ eDecayed = (1 << 4),
};
+ inline EProcessReturn operator|(EProcessReturn a, EProcessReturn b) {
+ return static_cast<EProcessReturn>(static_cast<int>(a) | static_cast<int>(b));
+ }
+
+ inline EProcessReturn& operator|=(EProcessReturn& a, EProcessReturn b) {
+ return a = a | b;
+ }
+
+ inline bool operator==(EProcessReturn a, EProcessReturn b) {
+ return (static_cast<int>(a) & static_cast<int>(b)) != 0;
+ }
+
} // namespace corsika::process
#endif
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index beeeb469f99106b8b18f636e075ee0abbae11c10..a75ad6adc3d9aeae6159028f0d905683ff8d88d0 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -60,20 +60,21 @@ namespace corsika::process {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
- A.DoContinuous(p, t, s);
+ ret |= A.DoContinuous(p, t, s);
}
if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
is_process_sequence<T2>::value) {
- B.DoContinuous(p, t, s);
+ ret |= B.DoContinuous(p, t, s);
}
return ret;
}
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const {
- corsika::units::si::LengthType max_length = // if no other process in the sequence implements it
+ corsika::units::si::LengthType
+ max_length = // if no other process in the sequence implements it
std::numeric_limits<double>::infinity() * corsika::units::si::meter;
-
+
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
corsika::units::si::LengthType const len = A.MaxStepLength(p, track);
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 1ad2d6428ef870aeaccd3fd55400d8e03a5ab2bf..fe4e284acd7f039147acebb39512390bcef1431c 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -12,7 +12,6 @@
on. This breaks ADL (argument-dependent lookup). Here we "fix" this:
*/
namespace phys::units {
- // using namespace phys::units::io;
using phys::units::io::operator<<;
} // namespace phys::units
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index a8c586061bb7e119588bab7ba8fb3b1f67ef7bec..93516c141ad13888e47b57d72d247c897c2ebdbf 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -15,8 +15,11 @@ namespace corsika::process {
namespace sibyll {
class Decay : public corsika::process::DecayProcess<Decay> {
+ mutable int fCount = 0;
+
public:
Decay() {}
+ ~Decay() { std::cout << "Sibyll::Decay n=" << fCount << std::endl; }
void Init() {
setHadronsUnstable();
setTrackedParticlesStable();
@@ -142,6 +145,7 @@ namespace corsika::process {
void DoDecay(Particle& p, Stack& s) const {
using corsika::geometry::Point;
using namespace corsika::units::si;
+ fCount++;
SibStack ss;
ss.Clear();
// copy particle to sibyll stack
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fbaf7c5f519391af653da149d7c9aa3d173d6b0e..eb45a570dd45dd679ff99437948a484aa7663c75 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -15,9 +15,11 @@ namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
+ mutable int fCount = 0;
+
public:
Interaction() {}
- ~Interaction() {}
+ ~Interaction() { std::cout << "Sibyll::Interaction n=" << fCount << std::endl; }
void Init() {
@@ -73,8 +75,8 @@ namespace corsika::process::sibyll {
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
+ std::cout << "Interaction: LambdaInt: \n"
+ << " input energy: " << p.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kBeam << endl
<< " beam XS code:" << kBeam << endl
<< " beam pid:" << p.GetPID() << endl
@@ -102,7 +104,6 @@ namespace corsika::process::sibyll {
<< "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
GrammageType int_length = kTarget * nucleon_mass / sig;
- // pick random step lenth
std::cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length << std::endl;
@@ -110,17 +111,6 @@ namespace corsika::process::sibyll {
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-
- /*
- what are the units of the output? slant depth or 3space length?
-
- */
- //
- // int a = 0;
- // const double next_step = -int_length * log(s_rndm_(a));
- // std::cout << "Interaction: "
- // << "next step (g/cm2): " << next_step << std::endl;
- // return next_step;
}
template <typename Particle, typename Stack>
@@ -137,9 +127,7 @@ namespace corsika::process::sibyll {
<< "DoInteraction: " << 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 =
+ const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point pOrig = p.GetPosition();
@@ -157,14 +145,12 @@ namespace corsika::process::sibyll {
const int kTarget = corsika::particles::Oxygen::
GetNucleusA(); // env.GetTargetParticle().GetPID();
- cout << "defining target momentum.." << endl;
// FOR NOW: target is always at rest
- const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
+ const EnergyType Etarget = 0_GeV + corsika::particles::Proton::GetMass();
const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
- cout << "target momentum (GeV/c): "
- << pTarget.GetComponents() / 1_GeV * constants::c << endl;
- cout << "beam momentum (GeV/c): "
- << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
+ cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
+ cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
+ << endl;
cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "time: " << tOrig << endl;
@@ -182,8 +168,7 @@ namespace corsika::process::sibyll {
// total momentum
MomentumVector Ptot = p.GetMomentum();
// invariant mass, i.e. cm. energy
- EnergyType Ecm =
- sqrt(Etot * Etot - Ptot.squaredNorm()); // sqrt( 2. * E * 0.93827_GeV );
+ EnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
/*
get transformation between Stack-frame and SibStack-frame
for EAS Stack-frame is lab. frame, could be different for CRMC-mode
@@ -207,8 +192,9 @@ namespace corsika::process::sibyll {
<< " DoInteraction: should have dropped particle.." << std::endl;
// p.Delete(); delete later... different process
} else {
+ fCount++;
// Sibyll does not know about units..
- double sqs = Ecm / 1_GeV;
+ const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, kTarget, sqs);
// running decays
@@ -228,7 +214,8 @@ namespace corsika::process::sibyll {
// SibStack does not know about momentum yet so we need counter to access
// momentum array in Sibyll
int i = -1;
- MomentumVector Ptot_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ EnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
++i;
// skip particles that have decayed in Sibyll
@@ -263,10 +250,14 @@ namespace corsika::process::sibyll {
pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
- Ptot_final += pnew.GetMomentum();
+ Plab_final += pnew.GetMomentum();
+ E_final += en_lab;
+ Ecm_final += psib.GetEnergy();
}
- // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
- // * constants::c << endl;
+ std::cout << "conservation (all GeV): E_final=" << E_final / 1_GeV
+ << ", Ecm_final=" << Ecm_final / 1_GeV
+ << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
+ << std::endl;
}
}
return process::EProcessReturn::eOk;
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index 604988a1da2a1d1b72d6bd1b4beb819924a7afc3..9bd417da2398d6eaec29b24d85242c1ad0d867ab 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -35,7 +35,7 @@ namespace corsika::process::sibyll {
SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) {
return static_cast<SibyllCode>(
- corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
+ corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
}
corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) {
diff --git a/Processes/Sibyll/sibyll2.3c.cc b/Processes/Sibyll/sibyll2.3c.cc
index fefea6d0b1acbac71e804cc83234bfe5bb335248..de48630e207e9a1831bd4ec9f34c4603333f8cfe 100644
--- a/Processes/Sibyll/sibyll2.3c.cc
+++ b/Processes/Sibyll/sibyll2.3c.cc
@@ -6,7 +6,7 @@
double s_rndm_(int&) {
static corsika::random::RNG& rng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-
+
std::uniform_real_distribution<double> dist;
return dist(rng);
}
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index f92ee8eee67b1b242b5e248ba07b815e4ee0bbb3..6fdd4f6624edc1c522352e954ed70a9c52b4b8d4 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -12,6 +12,7 @@
#define _include_corsika_processes_TrackingLine_h_
#include <corsika/geometry/Point.h>
+#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
@@ -75,12 +76,19 @@ namespace corsika::process {
void Init() {}
auto GetTrack(Particle const& p) {
+ using std::cout;
+ using std::endl;
using namespace corsika::units::si;
+ using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
+ std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
+ << std::endl;
+ std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+
// to do: include effect of magnetic field
geometry::Line line(currentPosition, velocity);
@@ -125,6 +133,8 @@ namespace corsika::process {
min = *minIter;
}
+ std::cout << " t-intersect: " << min << std::endl;
+
return geometry::Trajectory<corsika::geometry::Line>(line, min);
}
};
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index f07bac485f36e6b1144430d249178759dda695a4..481f8d48dc1abab991086e63d1f3b0962313844f 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -120,7 +120,7 @@ namespace corsika::stack {
void SetPosition(const int i, const corsika::geometry::Point& v) {
fPosition[i] = v;
}
-
+
void SetTime(const int i, const corsika::units::si::TimeType& v) { fTime[i] = v; }
corsika::particles::Code GetPID(const int i) const { return fDataPID[i]; }
diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc
index 092baf5f9259ad658392a78b9a753c37ad74d3e4..dbb585a237d2ec201bd99449046046ea6afa7ecc 100644
--- a/Stack/SuperStupidStack/testSuperStupidStack.cc
+++ b/Stack/SuperStupidStack/testSuperStupidStack.cc
@@ -45,10 +45,9 @@ TEST_CASE("SuperStupidStack", "[stack]") {
auto pout = s.GetNextParticle();
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() == stack::super_stupid::MomentumVector(dummyCS, {1_GeV,
+ // 1_GeV, 1_GeV})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 *
+ // meter, 1 * meter}));
REQUIRE(pout.GetTime() == 100_s);
}