diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 97acd45670ef02ca00f234d9bac8419be033ed4d..81c6fc535538d6789ee6170fbe44e6b2fb7fedd5 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -90,11 +90,14 @@ target_link_libraries (vertical_EAS
CORSIKAlogging
CORSIKArandom
ProcessSibyll
+ ProcessPythia
+ ProcessUrQMD
+ ProcessSwitch
CORSIKAcascade
ProcessEnergyLoss
+ ProcessObservationPlane
ProcessTrackWriter
ProcessTrackingLine
- ProcessHadronicElasticModel
ProcessParticleCut
ProcessStackInspector
CORSIKAprocesses
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 7b990016bbcb4e2a2a0794fddf16aa2a92399d37..edc1faf9b67dfb155684ad392fa3f5ec6a4be512 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -177,8 +177,4 @@ int main() {
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
-
- // basic check for unit-tests
- assert(cut.GetNumberEmParticles() == 29785);
- assert(cut.GetNumberInvParticles() == 26697);
}
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index d5a7097b17d0e6ca6112b1143ed7d9f9a8f191e1..f5c09343babcf62f87135a2fadb124c74a144f7e 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -103,6 +103,10 @@ int main() {
double theta = 0.;
double phi = 0.;
+ Point const injectionPos(
+ rootCS, 0_m, 0_m,
+ height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
+
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
@@ -118,12 +122,10 @@ int main() {
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- Point pos(rootCS, 0_m, 0_m,
- height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
- beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
+ beamCode, E0, plab, injectionPos, 0_ns, nuclA, nuclZ});
}
// setup processes, decays and interactions
@@ -139,7 +141,8 @@ int main() {
process::particle_cut::ParticleCut cut(80_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
- process::energy_loss::EnergyLoss eLoss;
+ process::energy_loss::EnergyLoss eLoss(
+ injectionPos, geometry::Vector<dimensionless_d>(rootCS, {0, 0, -1}));
// assemble all processes into an ordered process list
auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
diff --git a/Documentation/Examples/stopping_power.cc b/Documentation/Examples/stopping_power.cc
index edbd2bc3219017a7f05e94a429f51a7b3225050b..5e2423b0c63c557ccf18369677fb1bff5be5b1e5 100644
--- a/Documentation/Examples/stopping_power.cc
+++ b/Documentation/Examples/stopping_power.cc
@@ -40,7 +40,13 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
- process::energy_loss::EnergyLoss eLoss;
+ Point const injectionPos(
+ rootCS, 0_m, 0_m,
+ 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
+
+ Vector<dimensionless_d> showerAxis(rootCS, {0, 0, -1});
+
+ process::energy_loss::EnergyLoss eLoss(injectionPos, showerAxis);
setup::Stack stack;
@@ -68,12 +74,11 @@ int main() {
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- Point pos(rootCS, 0_m, 0_m,
- 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
+
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- beamCode, E0, plab, pos, 0_ns});
+ beamCode, E0, plab, injectionPos, 0_ns});
auto const p = stack.GetNextParticle();
HEPEnergyType dE = eLoss.TotalEnergyLoss(p, 1_g / square(1_cm));
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 0b3d4641a8fa6db56ea407049ceef2332d185228..78972df91f591a56a23c0bf6aeb6e8022a293549 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -10,24 +10,31 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/StackProcess.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
-#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
+#include <corsika/process/observation_plane/ObservationPlane.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
-#include <corsika/process/sibyll/Decay.h>
+//~ #include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
+#include <corsika/process/pythia/Decay.h>
+
+#include <corsika/process/urqmd/UrQMD.h>
+
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/track_writer/TrackWriter.h>
@@ -37,9 +44,7 @@
#include <corsika/utl/CorsikaFenv.h>
-#include <boost/type_index.hpp>
-using boost::typeindex::type_id_with_cvr;
-
+#include <iomanip>
#include <iostream>
#include <limits>
#include <typeinfo>
@@ -56,88 +61,125 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-//
-// The example main program for a particle cascade
-//
+void registerRandomStreams() {
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+ random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
+
+ random::RNGManager::GetInstance().SeedAll();
+}
+
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
- random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ registerRandomStreams();
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto& universe = *(env.GetUniverse());
- auto theMedium =
- EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
-
- // fraction of oxygen
- const float fox = 0.20946;
- using MyHomogeneousModel = HomogeneousMedium<environment::IMediumModel>;
- theMedium->SetModelProperties<MyHomogeneousModel>(
- 1_kg / (1_m * 1_m * 1_m),
+ auto theMedium = EnvType::CreateNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ auto constexpr temperature = 295_K; // AIRES default temperature for isothermal model
+ auto constexpr lambda =
+ -constants::R * temperature / (constants::g_sub_n * 28.966_g / mole);
+ // 1036 g/cm² taken from AIRES code
+ auto constexpr rho0 = -1036_g / square(1_cm) / lambda;
+ using FlatExp = environment::FlatExponential<environment::IMediumModel>;
+ theMedium->SetModelProperties<FlatExp>(
+ Point{rootCS, 0_m, 0_m, 0_m}, Vector<dimensionless_d>{rootCS, {0., 0., 1.}}, rho0,
+ lambda,
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
- std::vector<float>{(float)1. - fox, fox}));
-
- universe.AddChild(std::move(theMedium));
-
- const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ std::vector<float>{
+ 0.7847f,
+ 1.f - 0.7847f})); // values taken from AIRES manual, Ar removed for now
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const Code beamCode = Code::Nucleus;
- const int nuclA = 4;
- const int nuclZ = int(nuclA / 2.15 + 0.7);
- const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
- const HEPEnergyType E0 = nuclA * 10_TeV;
+ const Code beamCode = Code::Proton;
+ auto const mass = particles::GetMass(beamCode);
+ const HEPEnergyType E0 = 0.1_PeV;
double theta = 0.;
double phi = 0.;
- {
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
- HEPMomentumType P0 = elab2plab(E0, mass);
- auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
- return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
- -ptot * cos(theta));
- };
- auto const [px, py, pz] =
- momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
- auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
- cout << "input particle: " << beamCode << endl;
- cout << "input angles: theta=" << theta << " phi=" << phi << endl;
- cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- Point pos(rootCS, 0_m, 0_m,
- 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
- stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType, unsigned short, unsigned short>{
- beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
- }
+ Point const injectionPos(
+ rootCS, 0_m, 0_m, 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
+
+ // {
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, injectionPos, 0_ns});
+ // }
+
+ Line const line(injectionPos, plab.normalized() * 1_m * 1_Hz);
+ auto const velocity = line.GetV0().norm();
+
+ auto const observationHeight = 1.425_km;
+
+ setup::Trajectory const showerAxis(line, (112.8_km - observationHeight) / velocity);
+
+ auto const grammage = theMedium->GetModelProperties().IntegratedGrammage(
+ showerAxis, (112.8_km - observationHeight));
+ std::cout << "Grammage to ground: " << grammage / (1_g / square(1_cm)) << " g/cm²"
+ << std::endl;
+
+ universe.AddChild(std::move(theMedium));
// setup processes, decays and interactions
- tracking_line::TrackingLine tracking;
- stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
- random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ const std::vector<particles::Code> trackedHadrons = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
- process::sibyll::Decay decay;
- process::particle_cut::ParticleCut cut(20_GeV);
+ //~ process::sibyll::Decay decay(trackedHadrons);
+
+ process::pythia::Decay decay(trackedHadrons);
+ process::particle_cut::ParticleCut cut(5_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
- process::energy_loss::EnergyLoss eLoss;
+ process::energy_loss::EnergyLoss eLoss(showerAxis);
+
+ Plane const obsPlane(Point(rootCS, 0_m, 0_m, observationHeight),
+ Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
+ process::observation_plane::ObservationPlane observationLevel(obsPlane,
+ "particles.dat");
// assemble all processes into an ordered process list
- auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << stackInspect;
+
+ process::UrQMD::UrQMD urqmd;
+
+ auto sibyllSequence = sibyll << sibyllNuc;
+ process::switch_process::SwitchProcess switchProcess(urqmd, sibyllSequence, 55_GeV);
+ auto sequence = switchProcess << decay << eLoss << cut << observationLevel
+ << trackWriter;
// define air shower object, run simulation
+ tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
@@ -151,4 +193,7 @@ int main() {
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+
+ std::ofstream finish("finished");
+ finish << "run completed without error" << std::endl;
}
diff --git a/Environment/BaseExponential.h b/Environment/BaseExponential.h
index 462c73c26a0713233e7b3c37b30d25cb3f57bcb4..331a09c33b79ce5f40285b3bb91609b6cf1b9a9d 100644
--- a/Environment/BaseExponential.h
+++ b/Environment/BaseExponential.h
@@ -53,6 +53,8 @@ namespace corsika::environment {
units::si::GrammageType IntegratedGrammage(
geometry::Trajectory<geometry::Line> const& vLine, units::si::LengthType vL,
geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
+ if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
+
auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index 25cc2ffaac9c92553a8503b76e4d642856bec45c..b8cda38c3e2688fdd0f8096ddbaeba7884266fd6 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -27,7 +27,7 @@ namespace corsika::environment {
corsika::geometry::Point const&) const = 0;
// todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
- // approach for now, only lines are supported
+ // approach; for now, only lines are supported
virtual corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::LengthType) const = 0;
diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
index d97ba380880ddf018e685d010e688a9f639bd0ac..4c2a4f6705c4d1cf74403934acd39c151cc38cdf 100644
--- a/Environment/NuclearComposition.h
+++ b/Environment/NuclearComposition.h
@@ -15,6 +15,7 @@
#include <corsika/units/PhysicalUnits.h>
#include <cassert>
+#include <functional>
#include <numeric>
#include <random>
#include <stdexcept>
@@ -85,7 +86,6 @@ namespace corsika::environment {
auto WeightedSum(TFunction func) const {
using ResultQuantity = decltype(func(*fComponents.cbegin()));
- auto const sum = [](auto x, auto y) { return x + y; };
auto const prod = [&](auto const compID, auto const fraction) {
return func(compID) * fraction;
};
@@ -94,12 +94,12 @@ namespace corsika::environment {
return std::inner_product(
fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
ResultQuantity::zero(), // .zero() is defined for quantity types only
- sum, prod);
+ std::plus<ResultQuantity>(), prod);
} else {
return std::inner_product(
fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
ResultQuantity(0), // in other cases we have to use a bare 0
- sum, prod);
+ std::plus<ResultQuantity>(), prod);
}
}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 99c28b6c97819a7d6c2a92bb21590aa9244fcae6..0caa88ae6fd495258c0c1c473f2e29cf685c0cc4 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -122,7 +122,9 @@ namespace corsika::cascade {
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
auto pNext = fStack.GetNextParticle();
+ std::cout << "========= next: " << pNext.GetPID() << std::endl;
Step(pNext);
+ std::cout << "========= stack ============" << std::endl;
fProcessSequence.DoStack(fStack);
}
// do cascade equations, which can put new particles on Stack,
diff --git a/Framework/Geometry/Plane.h b/Framework/Geometry/Plane.h
index f6c201177989866ff7fbfec15253f911e592ec79..f82e883a5ebcf5c56f3dc16a8dc6c47b8567a93a 100644
--- a/Framework/Geometry/Plane.h
+++ b/Framework/Geometry/Plane.h
@@ -26,11 +26,16 @@ namespace corsika::geometry {
public:
Plane(Point const& vCenter, DimLessVec const& vNormal)
: fCenter(vCenter)
- , fNormal(vNormal) {}
+ , fNormal(vNormal.normalized()) {}
+
bool IsAbove(Point const& vP) const {
return fNormal.dot(vP - fCenter) > corsika::units::si::LengthType::zero();
}
+ units::si::LengthType DistanceTo(geometry::Point const& vP) const {
+ return (fNormal * (vP - fCenter).dot(fNormal)).norm();
+ }
+
Point const& GetCenter() const { return fCenter; }
DimLessVec const& GetNormal() const { return fNormal; }
};
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index c379e493bfc483481dbad4783f7504be9667eab4..35fa6dae14301213e428a4ddf2585535e6f810cd 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -341,7 +341,6 @@ namespace corsika::process {
/// marker to identify objectas ProcessSequence
template <typename A, typename B>
struct is_process_sequence<corsika::process::ProcessSequence<A, B>> : std::true_type {};
-
} // namespace corsika::process
#endif
diff --git a/Framework/Utilities/sgn.h b/Framework/Utilities/sgn.h
index b0084ee04827345548384b00cd73ea4c76209086..84fca4e550d80f3e9749c513566f80e66b5f15f9 100644
--- a/Framework/Utilities/sgn.h
+++ b/Framework/Utilities/sgn.h
@@ -1,4 +1,3 @@
-
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index febabc7764039da77c2012a627ae31de49daf87d..c9bf6afd5af77f67a64ec9f219be7813c48b22c6 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -15,7 +15,10 @@
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/geometry/Line.h>
+
#include <cmath>
+#include <fstream>
#include <iostream>
#include <limits>
@@ -32,11 +35,6 @@ namespace corsika::process::energy_loss {
return sqrt((Elab - m) * (Elab + m));
};
- EnergyLoss::EnergyLoss()
- : fEnergyLossTot(0_GeV)
- , fdX(10_g / square(1_cm)) // profile binning
- {}
-
/**
* PDG2018, passage of particles through matter
*
@@ -163,6 +161,7 @@ namespace corsika::process::energy_loss {
process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p,
SetupTrack const& t) {
if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
+
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
@@ -183,7 +182,7 @@ namespace corsika::process::energy_loss {
p.SetEnergy(Enew);
MomentumUpdate(p, Enew);
fEnergyLossTot += dE;
- GetXbin(p, t, dE);
+ FillProfile(p, t, dE);
return status;
}
@@ -211,37 +210,63 @@ namespace corsika::process::energy_loss {
vP.SetMomentum(pnew * Pnew / pnew.GetNorm());
}
-#include <corsika/geometry/CoordinateSystem.h>
-
- int EnergyLoss::GetXbin(SetupParticle const& vP, SetupTrack const& vTrack,
- const HEPEnergyType dE) {
+ void EnergyLoss::FillProfile(SetupParticle const& vP, SetupTrack const& vTrack,
+ const HEPEnergyType dE) {
using namespace corsika::geometry;
- CoordinateSystem const& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- Point const pos1(rootCS, 0_m, 0_m, 0_m);
- Point const pos2(rootCS, 0_m, 0_m, vTrack.GetPosition(0).GetCoordinates()[2]);
- auto const delta = (pos2 - pos1) / 1_s;
- Trajectory const t(Line(pos1, delta), 1_s);
+ auto const toStart = vTrack.GetPosition(0) - fInjectionPoint;
+ auto const toEnd = vTrack.GetPosition(1) - fInjectionPoint;
+
+ auto const v1 = (toStart * 1_Hz).dot(fShowerAxisDirection);
+ auto const v2 = (toEnd * 1_Hz).dot(fShowerAxisDirection);
+ geometry::Line const lineToStartBin(fInjectionPoint, fShowerAxisDirection * v1);
+ geometry::Line const lineToEndBin(fInjectionPoint, fShowerAxisDirection * v2);
+
+ SetupTrack const trajToStartBin(lineToStartBin, 1_s);
+ SetupTrack const trajToEndBin(lineToEndBin, 1_s);
+
+ GrammageType const grammageStart =
+ vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToStartBin,
+ trajToStartBin.GetLength());
+ GrammageType const grammageEnd =
+ vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToEndBin,
+ trajToEndBin.GetLength());
- GrammageType const grammage =
- vP.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
+ const int binStart = grammageStart / fdX;
+ const int binEnd = grammageEnd / fdX;
- const int bin = grammage / fdX;
+ std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
+ << grammageEnd << std::endl;
+
+ auto energyCount = HEPEnergyType::zero();
+
+ auto fill = [&](int bin, GrammageType weight) {
+ auto const increment = -dE * weight / (grammageEnd - grammageStart);
+ fProfile[bin] += increment;
+ energyCount += increment;
+
+ std::cout << "filling bin " << bin << " with weight " << weight << ": " << increment
+ << std::endl;
+ };
// fill longitudinal profile
- if (!fProfile.count(bin)) { cout << "EnergyLoss new x bin " << bin << endl; }
- fProfile[bin] += -dE / 1_GeV;
- return bin;
+ fill(binStart, (1 + binStart) * fdX - grammageStart);
+ fill(binEnd, grammageEnd - binEnd * fdX);
+
+ if (binStart == binEnd) { fill(binStart, -fdX); }
+
+ for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, fdX); }
+
+ std::cout << "total energy added to histogram: " << energyCount << std::endl;
}
void EnergyLoss::PrintProfile() const {
-
- cout << "EnergyLoss PrintProfile X-bin [g/cm2] dE/dX [GeV/g/cm2] " << endl;
+ std::ofstream file("EnergyLossProfile.dat");
+ cout << "# EnergyLoss PrintProfile X-bin [g/cm2] dE/dX [GeV/g/cm2] " << endl;
double const deltaX = fdX / 1_g * square(1_cm);
for (auto v : fProfile) {
- cout << v.first * deltaX << " " << v.second / deltaX << endl;
+ file << v.first * deltaX << " " << v.second / (deltaX * 1_GeV) << endl;
}
}
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index 4380a8ddcda43d3925df0ab839d010f3c0e0cad4..5061c7c616e34f6137693d2f16b04019ee389ddb 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -11,6 +11,8 @@
#ifndef _Processes_EnergyLoss_h_
#define _Processes_EnergyLoss_h_
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/units/PhysicalUnits.h>
@@ -29,13 +31,20 @@ namespace corsika::process::energy_loss {
void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew);
public:
- EnergyLoss();
+ template <typename TDim>
+ EnergyLoss(geometry::Point const& injectionPoint,
+ geometry::Vector<TDim> const& direction)
+ : fInjectionPoint(injectionPoint)
+ , fShowerAxisDirection(direction.normalized()) {}
+
+ EnergyLoss(setup::Trajectory const& trajectory)
+ : EnergyLoss(trajectory.GetPosition(0), trajectory.GetV0()){};
+
void Init() {}
process::EProcessReturn DoContinuous(setup::Stack::ParticleType&,
setup::Trajectory const&);
units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const&) const;
-
units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
void PrintProfile() const;
static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const&,
@@ -46,12 +55,19 @@ namespace corsika::process::energy_loss {
const units::si::GrammageType);
private:
- int GetXbin(setup::Stack::ParticleType const&, setup::Trajectory const&,
- units::si::HEPEnergyType);
+ void FillProfile(setup::Stack::ParticleType const&, setup::Trajectory const&,
+ units::si::HEPEnergyType);
+ // void FillProfileAbsorbed(setup::Stack::ParticleType const&, setup::Trajectory
+ // const&);
- units::si::HEPEnergyType fEnergyLossTot;
- units::si::GrammageType fdX; // profile binning
- std::map<int, double> fProfile; // longitudinal profile
+ units::si::HEPEnergyType fEnergyLossTot = units::si::HEPEnergyType::zero();
+ units::si::GrammageType const fdX = std::invoke([]() {
+ using namespace units::si;
+ return 10_g / square(1_cm);
+ }); // profile binning
+ std::map<int, units::si::HEPEnergyType> fProfile; // longitudinal profile
+ geometry::Point const fInjectionPoint;
+ geometry::Vector<units::si::dimensionless_d> const fShowerAxisDirection;
};
} // namespace corsika::process::energy_loss
diff --git a/Processes/ParticleCut/CMakeLists.txt b/Processes/ParticleCut/CMakeLists.txt
index d5a92febd68757f2978c97eee38022b38c22b53a..eed2c18f29e715c8794310f420ed58cc0dff03c1 100644
--- a/Processes/ParticleCut/CMakeLists.txt
+++ b/Processes/ParticleCut/CMakeLists.txt
@@ -29,6 +29,7 @@ target_link_libraries (
ProcessParticleCut
CORSIKAunits
CORSIKAparticles
+ CORSIKAprocesssequence
CORSIKAsetup
)
diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc
index 1c9ad8c671ade7fdfbc056be1d7f4eb8df38708b..5783c05fd32e31b691bb930931d9c98e49509dae 100644
--- a/Processes/ParticleCut/ParticleCut.cc
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -48,40 +48,16 @@ namespace corsika::process {
}
bool ParticleCut::ParticleIsInvisible(Code vCode) const {
- bool is_inv = false;
- // FOR NOW: switch
switch (vCode) {
case Code::NuE:
- is_inv = true;
- break;
case Code::NuEBar:
- is_inv = true;
- break;
case Code::NuMu:
- is_inv = true;
- break;
case Code::NuMuBar:
- is_inv = true;
- break;
- case Code::MuPlus:
- is_inv = true;
- break;
- case Code::MuMinus:
- is_inv = true;
- break;
-
- case Code::Neutron:
- is_inv = true;
- break;
-
- case Code::AntiNeutron:
- is_inv = true;
- break;
+ return true;
default:
- break;
+ return false;
}
- return is_inv;
}
EProcessReturn ParticleCut::DoSecondaries(corsika::setup::StackView& vS) {
diff --git a/Processes/ParticleCut/testParticleCut.cc b/Processes/ParticleCut/testParticleCut.cc
index c6948f27a4e865744c7efec9e6d276aebbc45c83..2375600c8bacb5908950dbeab0efdd0b46c19b3f 100644
--- a/Processes/ParticleCut/testParticleCut.cc
+++ b/Processes/ParticleCut/testParticleCut.cc
@@ -46,7 +46,7 @@ TEST_CASE("ParticleCut", "[processes]") {
particles::Code::Electron, particles::Code::MuPlus, particles::Code::NuE,
particles::Code::Neutron};
- SECTION("cut invisible") {
+ SECTION("cut on particle type") {
ParticleCut cut(20_GeV);
@@ -73,7 +73,7 @@ TEST_CASE("ParticleCut", "[processes]") {
cut.DoSecondaries(view);
- REQUIRE(view.GetSize() == 6);
+ REQUIRE(view.GetSize() == 8);
}
SECTION("cut low energy") {
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 20d2aaa5044e484b2727e5cdc03f011fe36be83b..eed8d5e9e0e9d345bb572b945721d540f69fa112 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -74,7 +74,7 @@ namespace corsika::process::pythia {
using namespace units::si;
HEPEnergyType E = p.GetEnergy();
- HEPMassType m = particles::GetMass(p.GetPID());
+ HEPMassType m = p.GetMass();
const double gamma = E / m;
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index 043b1d40524876c3ecdf3ad2003edcada07e0934..70565f0f82d793f028b520dcd849f99de9215b39 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -212,26 +212,12 @@ namespace corsika::process::pythia {
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
- // weighted sum
- int i = -1;
- si::CrossSectionType weightedProdCrossSection = 0_mb;
- // get weights of components from environment/medium
- const auto& w = mediumComposition.GetFractions();
- // loop over components in medium
- for (auto const targetId : mediumComposition.GetComponents()) {
- i++;
- cout << "Interaction: get interaction length for target: " << targetId << endl;
-
- auto const [productionCrossSection, elaCrossSection] =
- GetCrossSection(corsikaBeamId, targetId, ECoM);
- [[maybe_unused]] const auto& dummy_elaCrossSection = elaCrossSection;
-
- cout << "Interaction: IntLength: pythia return (mb): "
- << productionCrossSection / 1_mb << endl
- << "Interaction: IntLength: weight : " << w[i] << endl;
-
- weightedProdCrossSection += w[i] * productionCrossSection;
- }
+
+ auto const weightedProdCrossSection =
+ mediumComposition.WeightedSum([=](auto vTargetID) {
+ return std::get<0>(this->GetCrossSection(corsikaBeamId, vTargetID, ECoM));
+ });
+
cout << "Interaction: IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mb << endl
<< "Interaction: IntLength: average mass number: "
diff --git a/Processes/UrQMD/urqmd.f b/Processes/UrQMD/urqmd.f
index 9319dae63135f3fbe34fa422cd0ec0fc88f8a04c..b9951571ebedbad6d6c30f9833f2e4fbe1f5e5c7 100644
--- a/Processes/UrQMD/urqmd.f
+++ b/Processes/UrQMD/urqmd.f
@@ -338,6 +338,7 @@ cdh write(*,*)'(W) No collision in event ',event
c~ stop
endif
if (noc.ge.50000) then
+ print *,'UrQMD terminating...'
call exit(2) ! think of a better way to hand over the error
! to C++
endif
diff --git a/Tools/plot_tracks.sh b/Tools/plot_tracks.sh
index 135b410dc8d80523c8a5f340bf3c26634e35bb7f..16ecf4a6e961115ca11542e3d328d81f15ed4a9c 100755
--- a/Tools/plot_tracks.sh
+++ b/Tools/plot_tracks.sh
@@ -11,28 +11,34 @@
# with this script you can plot an animation of output of TrackWriter
track_dat=$1
-if [ -z "$track_dat" ]; then
- echo "usage: $0 <track.dat> [output.gif]" >&2
+muon_dat=$2
+
+if [ -z "$track_dat" ] || [ -z "$muon_dat" ]; then
+ echo "usage: $0 <hadron.dat> <muon.dat> [output.gif]" >&2
exit 1
fi
-output=$2
+output=$3
if [ -z "$output" ]; then
output="$track_dat.gif"
fi
cat <<EOF | gnuplot
-set term gif animate size 600,600
-set output "$output"
+set term png size 900,900
+#set output "$output"
+#set zrange [0:40e3]
+#set xrange [-10:10]
+#set yrange [-10:10]
set xlabel "x / m"
set ylabel "y / m"
set zlabel "z / m"
set title "CORSIKA 8 preliminary"
-do for [t=0:360:1] {
- set view 60, t
- splot "$track_dat" u 3:4:5:6:7:8 w vectors nohead t ""
+do for [t=0:359:1] {
+ set output sprintf("%03d_$output", t)
+ set view 90, t
+ splot "$muon_dat" u 3:4:5:6:7:8 w vectors nohead lt rgb "red" t "", "$track_dat" u 3:4:5:6:7:8 w vectors nohead lc rgb "black" t ""
}
EOF