IAP GITLAB

Commit fd16e22a authored by Ralf Ulrich's avatar Ralf Ulrich

made many small details more robust, made a lot of testing

parent dbf6e995
Pipeline #2913 failed with stages
in 0 seconds
......@@ -101,7 +101,7 @@ int main() {
// create "world" as infinite sphere filled with protons
auto world = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
Point{rootCS, 0_m, 0_m, 0_m}, 100_km);
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
......
......@@ -57,6 +57,8 @@ int main() {
logging::SetLevel(logging::level::info);
C8LOG_INFO("vertical_EAS");
std::cout << "cascade_example" << std::endl;
const LengthType height_atmosphere = 112.8_km;
......@@ -72,7 +74,7 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto world = setup::Environment::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
......@@ -110,7 +112,7 @@ int main() {
rootCS, 0_m, 0_m,
height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -5000_km}, env};
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
......
......@@ -10,6 +10,7 @@
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
......@@ -17,6 +18,7 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/ShowerAxis.h>
#include <corsika/geometry/Sphere.h>
......@@ -59,7 +61,7 @@ using namespace corsika::units::si;
//
int main() {
logging::SetLevel(logging::level::trace);
logging::SetLevel(logging::level::info);
std::cout << "cascade_proton_example" << std::endl;
......@@ -73,20 +75,20 @@ int main() {
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto theMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
auto world = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
theMedium->SetModelProperties<MyHomogeneousModel>(
world->SetModelProperties<MyHomogeneousModel>(
environment::Medium::AirDry1Atm, geometry::Vector(rootCS, 0_T, 0_T, 1_T),
1_kg / (1_m * 1_m * 1_m),
NuclearComposition(std::vector<particles::Code>{particles::Code::Hydrogen},
std::vector<float>{(float)1.}));
universe.AddChild(std::move(theMedium));
universe.AddChild(std::move(world));
// setup particle stack, and add primary particle
setup::Stack stack;
......@@ -97,6 +99,7 @@ int main() {
double theta = 0.;
double phi = 0.;
Point injectionPos(rootCS, 0_m, 0_m, 0_m);
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
......@@ -112,11 +115,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, 0_m);
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});
}
// setup processes, decays and interactions
......@@ -130,20 +132,19 @@ int main() {
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay;
process::pythia::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV, true, true);
process::particle_cut::ParticleCut cut(60_GeV, true, true);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
process::track_writer::TrackWriter trackWriter("tracks.dat");
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
// assemble all processes into an ordered process list
// auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence = process::sequence(pythia, decay, cut, trackWriter, stackInspect);
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
auto sequence = process::sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
......
......@@ -71,13 +71,13 @@ int main(int argc, char** argv) {
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
auto builder = environment::make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 0_T, 1_T});
geometry::Vector{rootCS, 0_T, 50_uT, 0_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
......
......@@ -101,7 +101,7 @@ int main(int argc, char** argv) {
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 0_T, 1_T});
geometry::Vector{rootCS, 0_T, 50_uT, 0_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
......
......@@ -85,7 +85,7 @@ using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagnetic
int main(int argc, char** argv) {
logging::SetLevel(logging::level::trace);
logging::SetLevel(logging::level::info);
C8LOG_INFO("vertical_EAS");
......@@ -110,7 +110,7 @@ int main(int argc, char** argv) {
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 5000_mT, 0_T});
geometry::Vector{rootCS, 0_T, 50_uT, 0_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
......@@ -256,9 +256,9 @@ int main(int argc, char** argv) {
EnergySwitch(55_GeV));
auto decaySequence = process::sequence(decayPythia, decaySibyll);
stack_inspector::StackInspector<setup::Stack> stackInspect(1000, false, E0);
auto sequence =
process::sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
proposalCounted, em_continuous, cut, trackWriter, observationLevel, longprof);
auto sequence = process::sequence(stackInspect, hadronSequence, reset_particle_mass,
decaySequence, proposalCounted, em_continuous, cut,
trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
......
......@@ -39,6 +39,8 @@ namespace corsika::process {
* \function LeapFrogStep
*
* Performs one leap-frog step consistent of two halve-steps with steplength/2
* The step is caluculated analytically precisely to reach to the next volume
*boundary.
**/
template <typename TParticle>
auto LeapFrogStep(const TParticle& particle,
......@@ -50,7 +52,7 @@ namespace corsika::process {
} // charge of the particle
const int chargeNumber = particle.GetChargeNumber();
auto const* currentLogicalVolumeNode = particle.GetNode();
auto magneticfield =
MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->GetModelProperties().GetMagneticField(
particle.GetPosition());
geometry::Vector<SpeedType::dimension_type> velocity =
......@@ -85,6 +87,9 @@ namespace corsika::process {
class Tracking : public corsika::process::tracking::Intersect<Tracking> {
public:
Tracking()
: straightTracking_{tracking_line::Tracking()} {}
template <typename TParticle>
auto GetTrack(TParticle const& particle) {
using namespace corsika::units::si;
......@@ -111,23 +116,34 @@ namespace corsika::process {
// maximum step-length since we need to follow curved
// trajectories segment-wise -- at least if we don't employ concepts as "Helix
// Trajectories" or similar
const auto& magneticfield =
MagneticFieldVector const& magneticfield =
volumeNode.GetModelProperties().GetMagneticField(position);
const auto magnitudeB = magneticfield.norm();
const int chargeNumber = particle.GetChargeNumber();
auto const momentumVerticalMag =
particle.GetMomentum() -
particle.GetMomentum().parallelProjectionOnto(magneticfield);
corsika::units::si::MagneticFluxType const magnitudeB = magneticfield.norm();
int const chargeNumber = particle.GetChargeNumber();
bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
if (no_deflection) { return GetLinearTrajectory(particle); }
HEPMomentumType const pAlongB_delta =
(particle.GetMomentum() -
particle.GetMomentum().parallelProjectionOnto(magneticfield))
.norm();
if (pAlongB_delta == 0_GeV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
C8LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
return GetLinearTrajectory(particle);
}
LengthType const gyroradius =
(chargeNumber == 0 || magnitudeB == 0_T
? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
: momentumVerticalMag.norm() * 1_V /
(corsika::units::constants::c * abs(chargeNumber) * magnitudeB *
1_eV));
(pAlongB_delta * 1_V /
(corsika::units::constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
const double maxRadians = 0.01;
const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
const TimeType steplimit_time = steplimit / initialVelocity.norm();
C8LOG_DEBUG("gyroradius {}, steplimit: {} m = {} s", gyroradius, steplimit,
C8LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
steplimit_time);
// traverse the environment volume tree and find next
......@@ -148,39 +164,48 @@ namespace corsika::process {
const corsika::geometry::Sphere& sphere,
const TMedium& medium) {
using namespace corsika::units::si;
if (sphere.GetRadius() == 1_km * std::numeric_limits<double>::infinity()) {
return geometry::Intersections();
}
const int chargeNumber = particle.GetChargeNumber();
const auto& position = particle.GetPosition();
const auto& magneticfield = medium.GetMagneticField(position);
MagneticFieldVector const& magneticfield = medium.GetMagneticField(position);
const geometry::Vector<SpeedType::dimension_type> velocity =
particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
const geometry::Vector<dimensionless_d> directionBefore =
velocity.normalized(); // determine steplength to next volume
if (chargeNumber == 0 || magneticfield.norm() == 0_T) {
auto const projectedDirection = directionBefore.cross(magneticfield);
auto const projectedDirectionSqrNorm = projectedDirection.GetSquaredNorm();
bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
if (chargeNumber == 0 || magneticfield.norm() == 0_T || isParallel) {
return tracking_line::Tracking::Intersect(particle, sphere, medium);
}
bool const numericallyInside = sphere.Contains(particle.GetPosition());
const geometry::Vector<SpeedType::dimension_type> velocity =
particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
const auto absVelocity = velocity.norm();
auto energy = particle.GetEnergy();
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(absVelocity * energy * 1_V);
const geometry::Vector<dimensionless_d> directionBefore =
velocity.normalized(); // determine steplength to next volume
auto const denom =
(directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k;
const double a =
((directionBefore.cross(magneticfield)).dot(position - sphere.GetCenter()) *
k +
1) *
4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
4 / (1_m * 1_m * denom);
const double b = directionBefore.dot(position - sphere.GetCenter()) * 8 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k *
k * 1_m * 1_m * 1_m);
(denom * 1_m * 1_m * 1_m);
const double c = ((position - sphere.GetCenter()).GetSquaredNorm() -
(sphere.GetRadius() * sphere.GetRadius())) *
4 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k *
k * 1_m * 1_m * 1_m * 1_m);
4 / (denom * 1_m * 1_m * 1_m * 1_m);
C8LOG_TRACE("denom={}, a={}, b={}, c={}", denom, a, b, c);
std::complex<double>* solutions = solve_quartic(0, a, b, c);
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
......@@ -252,6 +277,38 @@ namespace corsika::process {
throw std::runtime_error(
"The Volume type provided is not supported in Intersect(particle, node)");
}
protected:
/**
* Use internally stored class tracking_line::Tracking to
* perform a straight line tracking, if no magnetic bendig was
* detected.
*
*/
template <typename TParticle>
auto GetLinearTrajectory(TParticle& particle) {
using namespace corsika::units::si;
// perform simple linear tracking
auto [straightTrajectory, minNode] = straightTracking_.GetTrack(particle);
// return as leap-frog trajectory
return std::make_tuple(
geometry::LeapFrogTrajectory(
straightTrajectory.GetLine().GetR0(),
straightTrajectory.GetLine().GetV0(),
MagneticFieldVector(particle.GetPosition().GetCoordinateSystem(), 0_T,
0_T, 0_T),
square(0_m) / (square(1_s) * 1_V),
straightTrajectory.GetDuration()), // trajectory
minNode); // next volume node
}
protected:
tracking_line::Tracking
straightTracking_; ///! we want this for neutral and B=0T tracks
}; // namespace tracking_leapfrog_curved
} // namespace tracking_leapfrog_curved
......
......@@ -52,6 +52,19 @@ namespace corsika::process {
class Tracking : public tracking_line::Tracking {
public:
/**
* \param firstFraction fraction of first leap-frog halve step
* relative to full linear step to next volume boundary. This
* should not be less than 0.5, otherwise you risk that
* particles will never travel from one volume to the next
* one. A cross should be possible (even likely). If
* firstFraction is too big (~1) the resulting calculated error
* will be largest.
*
*/
Tracking(double firstFraction = 0.55)
: firstFraction_(firstFraction) {}
template <typename Particle>
auto GetTrack(Particle& particle) {
using namespace corsika::units::si;
......@@ -72,39 +85,9 @@ namespace corsika::process {
typedef decltype(particle.GetNode()) node_type;
const node_type volumeNode = particle.GetNode();
auto magneticfield =
volumeNode->GetModelProperties().GetMagneticField(initialPosition);
// charge of the particle
const int chargeNumber = particle.GetChargeNumber();
const auto magnitudeB = magneticfield.GetNorm();
C8LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
magneticfield.GetComponents() / 1_uT, chargeNumber, magnitudeB / 1_T);
// we need to limit maximum step-length since we absolutely
// need to follow strongly curved trajectories segment-wise,
// at least if we don't employ concepts as "Helix
// Trajectories" or similar
auto const momentumVerticalMag =
particle.GetMomentum() -
particle.GetMomentum().parallelProjectionOnto(magneticfield);
bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
LengthType const gyroradius =
(no_deflection ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
: momentumVerticalMag.norm() * 1_V /
(corsika::units::constants::c * abs(chargeNumber) *
magnitudeB * 1_eV));
const double maxRadians = 0.01;
const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
C8LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
// calculate first halve step for "steplimit"
const auto initialMomentum = particle.GetMomentum();
const auto absMomentum = initialMomentum.norm();
const auto absVelocity = initialVelocity.norm();
const geometry::Vector<dimensionless_d> direction = initialVelocity.normalized();
// check if particle is moving at all
const auto absVelocity = initialVelocity.norm();
if (absVelocity * 1_s == 0_m) {
return std::make_tuple(
geometry::LineTrajectory(geometry::Line(initialPosition, initialVelocity),
......@@ -112,6 +95,15 @@ namespace corsika::process {
volumeNode);
}
// charge of the particle, and magnetic field
const int chargeNumber = particle.GetChargeNumber();
auto magneticfield =
volumeNode->GetModelProperties().GetMagneticField(initialPosition);
const auto magnitudeB = magneticfield.GetNorm();
C8LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
magneticfield.GetComponents() / 1_uT, chargeNumber, magnitudeB / 1_T);
bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
// check, where the first halve-step direction has geometric intersections
const auto [initialTrack, initialTrackNextVolume] =
tracking_line::Tracking::GetTrack(particle);
......@@ -128,9 +120,39 @@ namespace corsika::process {
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
HEPMomentumType const pAlongB_delta =
(particle.GetMomentum() -
particle.GetMomentum().parallelProjectionOnto(magneticfield))
.norm();
if (pAlongB_delta == 0_GeV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
C8LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
LengthType const gyroradius =
(pAlongB_delta * 1_V /
(corsika::units::constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
// we need to limit maximum step-length since we absolutely
// need to follow strongly curved trajectories segment-wise,
// at least if we don't employ concepts as "Helix
// Trajectories" or similar
const double maxRadians = 0.01;
const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
C8LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
// calculate first halve step for "steplimit"
const auto initialMomentum = particle.GetMomentum();
const auto absMomentum = initialMomentum.norm();
const geometry::Vector<dimensionless_d> direction = initialVelocity.normalized();
// avoid any intersections within first halve steplength
LengthType const firstHalveSteplength =
std::min(steplimit, initialTrackLength) / 2;
std::min(steplimit, initialTrackLength * firstFraction_);
C8LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
firstHalveSteplength, steplimit, initialTrackLength);
......@@ -203,6 +225,9 @@ namespace corsika::process {
new_direction_normalized * absVelocity), // trajectory
(switch_volume ? finalTrackNextVolume : volumeNode));
}
protected:
double firstFraction_;
};
} // namespace tracking_leapfrog_straight
......
......@@ -35,6 +35,7 @@ namespace corsika::process {
class Tracking : public corsika::process::tracking::Intersect<Tracking> {
public:
template <typename TParticle>
auto GetTrack(TParticle const& particle) {
using namespace corsika::units::si;
......@@ -115,6 +116,7 @@ namespace corsika::process {
1_s
: n.dot(delta) / c);
}
};
} // namespace tracking_line
......
......@@ -20,14 +20,39 @@
namespace corsika::setup {
// Note: Tracking and Trajectory must fit together !
/**
Note/Warning: Tracking and Trajectory must fit together !
tracking_leapfrog_curved::Tracking is the result of the Bachelor
thesis of Andre Schmidt, KIT. This is a leap-frog algorithm with
an analytical, precise calculation of volume intersections. This
algorithm needs a LeapFrogTrajectory.
tracking_leapfrog_straight::Tracking is a more simple and direct
leap-frog implementation. The two halve steps are coded explicitly
as two straight segments. Intersections with other volumes are
calculate only on the straight segments. This algorithm is based
on LineTrajectory.
tracking_line::Tracking is a pure straight tracker. It is based on
LineTrajectory.
*/
typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking;
// tracking_leapfrog_straight::Tracking tracking;
//typedef corsika::process::tracking_leapfrog_straight::Tracking Tracking;
//typedef corsika::process::tracking_line::Tracking Tracking;
/// definition of Trajectory base class, to be used in tracking and cascades
// typedef corsika::geometry::LineTrajectory Trajectory;
//typedef corsika::geometry::LineTrajectory Trajectory;
typedef corsika::geometry::LeapFrogTrajectory Trajectory;
/**
The following section is for unit testing only. Eventually it should
be moved to "tests".
*/
namespace testing {
template <typename TTrack>
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment