IAP GITLAB

Skip to content
Snippets Groups Projects
Commit cd5772fb authored by ralfulrich's avatar ralfulrich
Browse files

style, remove Stack from ProcessInterface

parent 4a7d7f14
No related branches found
No related tags found
No related merge requests found
Showing
with 64 additions and 68 deletions
......@@ -4,7 +4,7 @@ add_subdirectory (Sibyll)
if (PYTHIA8_FOUND)
add_subdirectory (Pythia)
endif (PYTHIA8_FOUND)
add_subdirectory (StackInspector)
# add_subdirectory (StackInspector)
add_subdirectory (TrackWriter)
add_subdirectory (TrackingLine)
add_subdirectory (HadronicElasticModel)
......
......@@ -125,7 +125,7 @@ namespace corsika::process::EnergyLoss {
// Barkas correction O(Z3) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
double A = 1;
[[maybe_unused]] double A = 1;
if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
// double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
// double const Llow = 0.01 * Erel;
......@@ -148,7 +148,7 @@ namespace corsika::process::EnergyLoss {
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) {
process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t) {
if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
......@@ -178,16 +178,16 @@ namespace corsika::process::EnergyLoss {
return units::si::meter * std::numeric_limits<double>::infinity();
}
void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& p,
void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,
corsika::units::si::HEPEnergyType Enew) {
HEPMomentumType Pnew = elab2plab(Enew, p.GetMass());
auto pnew = p.GetMomentum();
p.SetMomentum(pnew * Pnew / pnew.GetNorm());
HEPMomentumType Pnew = elab2plab(Enew, vP.GetMass());
auto pnew = vP.GetMomentum();
vP.SetMomentum(pnew * Pnew / pnew.GetNorm());
}
#include <corsika/geometry/CoordinateSystem.h>
int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& p,
int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& vP,
const HEPEnergyType dE) {
using namespace corsika::geometry;
......@@ -195,12 +195,12 @@ namespace corsika::process::EnergyLoss {
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point pos1(rootCS, 0_m, 0_m, 0_m);
Point pos2(rootCS, 0_m, 0_m, p.GetPosition().GetCoordinates()[2]);
Point pos2(rootCS, 0_m, 0_m, vP.GetPosition().GetCoordinates()[2]);
Vector delta = (pos2 - pos1) / 1_s;
Trajectory t(Line(pos1, delta), 1_s);
GrammageType const grammage =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
vP.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
const int bin = grammage / fdX;
......
......@@ -35,8 +35,7 @@ namespace corsika::process::EnergyLoss {
void Init() {}
corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&,
corsika::setup::Trajectory&,
corsika::setup::Stack&);
corsika::setup::Trajectory&);
corsika::units::si::LengthType MaxStepLength(corsika::setup::Stack::ParticleType&,
corsika::setup::Trajectory&);
......
......@@ -83,7 +83,7 @@ namespace corsika::process::HadronicElasticModel {
}
template <>
process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p, Stack&) {
process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p) {
if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; }
using namespace units::si;
......
......@@ -65,8 +65,8 @@ namespace corsika::process::HadronicElasticModel {
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&);
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&);
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
};
} // namespace corsika::process::HadronicElasticModel
......
......@@ -24,7 +24,8 @@ using std::vector;
using namespace corsika;
using namespace corsika::setup;
using Particle = Stack::StackIterator; // ParticleType;
using Projectile = corsika::setup::StackView::ParticleType;
using Particle = corsika::setup::Stack::ParticleType;
using Track = Trajectory;
namespace corsika::process::sibyll {
......@@ -103,27 +104,27 @@ namespace corsika::process::sibyll {
}
template <>
units::si::TimeType Decay::GetLifetime(Particle const& p) {
units::si::TimeType Decay::GetLifetime(Particle const& vP) {
using namespace units::si;
HEPEnergyType E = p.GetEnergy();
HEPMassType m = p.GetMass();
HEPEnergyType E = vP.GetEnergy();
HEPMassType m = vP.GetMass();
const double gamma = E / m;
const TimeType t0 = particles::GetLifetime(p.GetPID());
const TimeType t0 = particles::GetLifetime(vP.GetPID());
auto const lifetime = gamma * t0;
const auto mkin =
(E * E - p.GetMomentum().squaredNorm()); // delta_mass(p.GetMomentum(), E, m);
cout << "Decay: code: " << p.GetPID() << endl;
(E * E - vP.GetMomentum().squaredNorm()); // delta_mass(vP.GetMomentum(), E, m);
cout << "Decay: code: " << vP.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV"
cout << "Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV << " GeV"
<< endl;
cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV << " "
<< m / 1_GeV * m / 1_GeV << endl;
auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
cout << "Decay: MinStep: tau: " << lifetime << endl;
......@@ -132,23 +133,23 @@ namespace corsika::process::sibyll {
}
template <>
void Decay::DoDecay(Particle& p, Stack&) {
void Decay::DoDecay(Projectile& vP) {
using geometry::Point;
using namespace units::si;
fCount++;
SibStack ss;
ss.Clear();
const particles::Code pCode = p.GetPID();
const particles::Code pCode = vP.GetPID();
// copy particle to sibyll stack
ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
p.GetMomentum(),
ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), vP.GetEnergy(),
vP.GetMomentum(),
// setting particle mass with Corsika values, may be inconsistent
// with sibyll internal values
particles::GetMass(pCode));
// remember position
Point const decayPoint = p.GetPosition();
TimeType const t0 = p.GetTime();
Point const decayPoint = vP.GetPosition();
TimeType const t0 = vP.GetTime();
// set all particles/hadrons unstable
// setHadronsUnstable();
SetUnstable(pCode);
......@@ -166,7 +167,7 @@ namespace corsika::process::sibyll {
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
if (psib.HasDecayed()) continue;
// add to corsika stack
p.AddSecondary(
vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()), psib.GetEnergy(),
......@@ -174,8 +175,6 @@ namespace corsika::process::sibyll {
}
// empty sibyll stack
ss.Clear();
// remove original particle from corsika stack
p.Delete();
}
} // namespace corsika::process::sibyll
......@@ -35,10 +35,10 @@ namespace corsika::process {
void SetHadronsUnstable();
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle const& p);
corsika::units::si::TimeType GetLifetime(Particle const&);
template <typename Particle, typename Stack>
void DoDecay(Particle& p, Stack&);
template <typename Projectile>
void DoDecay(Projectile&);
};
} // namespace sibyll
} // namespace corsika::process
......
......@@ -29,7 +29,8 @@ using std::tuple;
using namespace corsika;
using namespace corsika::setup;
using Particle = Stack::StackIterator; // ParticleType;
using Particle = Stack::StackIterator; // ParticleType;
using Projectile = StackView::StackIterator; // ParticleType;
using Track = Trajectory;
namespace corsika::process::sibyll {
......@@ -211,7 +212,7 @@ namespace corsika::process::sibyll {
*/
template <>
process::EProcessReturn Interaction::DoInteraction(Particle& p, Stack&) {
process::EProcessReturn Interaction::DoInteraction(Projectile& p) {
using namespace units;
using namespace utl;
......@@ -336,7 +337,6 @@ namespace corsika::process::sibyll {
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
throw std::runtime_error("energy too low for SIBYLL");
// p.Delete(); delete later... different process
} else {
fCount++;
// Sibyll does not know about units..
......@@ -350,9 +350,6 @@ namespace corsika::process::sibyll {
sib_list_(print_unit);
fNucCount += get_nwounded() - 1;
// delete current particle
p.Delete();
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
......
......@@ -65,8 +65,8 @@ namespace corsika::process::sibyll {
event is copied (and boosted) into the shower lab frame.
*/
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
private:
corsika::environment::Environment const& fEnvironment;
......
......@@ -31,7 +31,8 @@ using std::vector;
using namespace corsika;
using namespace corsika::setup;
using Particle = Stack::StackIterator; // ParticleType;
using Particle = Stack::ParticleType; // StackIterator; // ParticleType;
using Projectile = StackView::StackIterator; // StackView::ParticleType;
using Track = Trajectory;
namespace corsika::process::sibyll {
......@@ -312,7 +313,7 @@ namespace corsika::process::sibyll {
}
template <>
process::EProcessReturn NuclearInteraction::DoInteraction(Particle& p, Stack& s) {
process::EProcessReturn NuclearInteraction::DoInteraction(Projectile& p) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
......@@ -599,12 +600,9 @@ namespace corsika::process::sibyll {
PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
// create inelastic interaction
cout << "calling HadronicInteraction..." << endl;
fHadronicInteraction.DoInteraction(inelasticNucleon, s);
fHadronicInteraction.DoInteraction(inelasticNucleon);
}
// delete parent particle
p.Delete();
cout << "NuclearInteraction: DoInteraction: done" << endl;
return process::EProcessReturn::eOk;
......
......@@ -60,8 +60,8 @@ namespace corsika::process::sibyll {
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&);
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s);
template <typename Projectle>
corsika::process::EProcessReturn DoInteraction(Projectle& p);
private:
corsika::environment::Environment const& fEnvironment;
......
......@@ -122,12 +122,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction model(env);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
......@@ -146,13 +147,14 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction hmodel(env);
NuclearInteraction model(env, hmodel);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
......@@ -169,6 +171,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
const std::vector<particles::Code> particleList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
......@@ -177,8 +181,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
Decay model(particleList);
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
stack);
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
}
......@@ -36,8 +36,8 @@ template <typename Stack>
StackInspector<Stack>::~StackInspector() {}
template <typename Stack>
process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
Stack& s) {
process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&,
setup::Trajectory&) {
if (!fReport) return process::EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
......
......@@ -31,7 +31,7 @@ namespace corsika::process {
~StackInspector();
void Init();
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s);
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&);
corsika::units::si::LengthType MaxStepLength(Particle&,
corsika::setup::Trajectory&);
......
......@@ -33,13 +33,13 @@ namespace corsika::process::TrackWriter {
}
template <>
process::EProcessReturn TrackWriter::DoContinuous(Particle& p, Track& t, Stack&) {
process::EProcessReturn TrackWriter::DoContinuous(Particle& vP, Track& vT) {
using namespace units::si;
auto const start = t.GetPosition(0).GetCoordinates();
auto const delta = t.GetPosition(1).GetCoordinates() - start;
auto const& name = particles::GetName(p.GetPID());
auto const start = vT.GetPosition(0).GetCoordinates();
auto const delta = vT.GetPosition(1).GetCoordinates() - start;
auto const& name = particles::GetName(vP.GetPID());
fFile << name << " " << p.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
fFile << name << " " << vP.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
<< start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
<< delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
......
......@@ -28,8 +28,8 @@ namespace corsika::process::TrackWriter {
void Init();
template <typename Particle, typename Track, typename Stack>
corsika::process::EProcessReturn DoContinuous(Particle& p, Track& t, Stack&);
template <typename Particle, typename Track>
corsika::process::EProcessReturn DoContinuous(Particle&, Track&);
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle&, Track&);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment