IAP GITLAB

Commit b4b94403 authored by Andre Schmidt's avatar Andre Schmidt

change Steplength

parent 74d04301
Pipeline #2004 canceled with stages
......@@ -26,6 +26,7 @@
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
......@@ -123,7 +124,7 @@ int main(int argc, char** argv) {
node4->AddChild(std::move(node5));
node3->AddChild(std::move(node4));
node2->AddChild(std::move(node3));
node1->AddChild(std::move(node2));
node1->AddChild(std::move(node2));
env.GetUniverse()->AddChild(std::move(node1));
// setup particle stack, and add primary particle
......@@ -227,7 +228,7 @@ int main(int argc, char** argv) {
decaySibyll.PrintDecayConfig();
process::particle_cut::ParticleCut cut{50_GeV};
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss(showerAxis);
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
......@@ -245,7 +246,7 @@ int main(int argc, char** argv) {
55_GeV);
auto decaySequence = decayPythia << decaySibyll;
auto sequence = switchProcess << decaySequence << longprof << eLoss << cut
<< observationLevel;
<< observationLevel << trackWriter;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
......
......@@ -196,7 +196,7 @@ namespace corsika::cascade {
vParticle.GetEnergy() * units::constants::c;
// determine geometric tracking
auto [step, geomMaxLength, nextVol, magMaxLength, directionBefore, directionAfter] = fTracking.GetTrack(vParticle);
auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length
......@@ -208,18 +208,46 @@ namespace corsika::cascade {
LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
std::cout << "distance_max=" << distance_max << std::endl;
// take minimum of geometry, interaction, decay for next step
// take minimum of geometry, interaction, decay, magnetic field for next step
auto const min_distance = std::min(
{distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength});
{distance_interact, distance_decay, distance_max, geomMaxLength});
std::cout << " move particle by : " << min_distance << std::endl;
//determine displacement by the magnetic field
auto const* currentLogicalVolumeNode = vParticle.GetNode();
int chargeNumber;
if(corsika::particles::IsNucleus(vParticle.GetPID())) {
chargeNumber = vParticle.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
}
if(chargeNumber != 0) {
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
geometry::Vector<SpeedType::dimension_type> velocity = vParticle.GetMomentum() / vParticle.GetEnergy() *
corsika::units::constants::c;
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
LengthType const steplength = min_distance /
(directionBefore + directionBefore.cross(magneticfield) * k / 2).GetNorm();
// First Movement
//assuming magnetic field does not change during movement
auto position = vParticle.GetPosition() + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - vParticle.GetPosition()) /
(position - vParticle.GetPosition()).GetNorm();
vParticle.SetMomentum(direction * vParticle.GetMomentum().GetNorm());
}
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ...
vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) +
directionAfter * min_distance /magMaxLength) * vParticle.GetMomentum().GetNorm());
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
step.LimitEndTo(min_distance);
......@@ -246,7 +274,7 @@ namespace corsika::cascade {
TStackView secondaries(vParticle);
if (min_distance != distance_max && min_distance != magMaxLength) {
if (min_distance != distance_max) {
/*
Create SecondaryView object on Stack. The data container
remains untouched and identical, and 'projectil' is identical
......
......@@ -20,6 +20,7 @@ endif ()
set (
UTILITIES_SOURCES
COMBoost.cc
quartic.cpp
${CORSIKA_FENV})
set (
......@@ -30,6 +31,7 @@ set (
sgn.h
CorsikaFenv.h
MetaProgramming.h
quartic.h
)
set (
......
......@@ -18,6 +18,7 @@
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/quartic.h>
#include <optional>
#include <type_traits>
#include <utility>
......@@ -61,7 +62,7 @@ namespace corsika::process {
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
// determine velocity after adding magnetic field
// determine direction of the particle after adding magnetic field
auto const* currentLogicalVolumeNode = p.GetNode();
int chargeNumber;
if(corsika::particles::IsNucleus(p.GetPID())) {
......@@ -69,39 +70,51 @@ namespace corsika::process {
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
}
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
auto magMaxLength = 1_m/0;
auto directionAfter = directionBefore;
if(chargeNumber != 0) {
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
std::cout << "TrackingLine B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
velocity.parallelProjectionOnto(magneticfield);
LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
(corsika::units::constants::cSquared * abs(chargeNumber) *
magneticfield.GetNorm() * 1_eV);
//steplength depending on how exact it should be
LengthType const Steplength = 0.01 * gyroradius;
// First Movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field at position
magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(position);
directionAfter = directionBefore + directionBefore.cross(magneticfield) * chargeNumber *
Steplength * corsika::units::constants::cSquared * 1_eV /
(p.GetEnergy() * velocity.GetNorm() * 1_V);
// Second Movement
position = position + directionAfter * Steplength / 2;
magMaxLength = (position - currentPosition).GetNorm();
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
magMaxLength;
velocity = direction * velocity.GetNorm();
std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
<< " GeV " << std::endl;
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
//determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - .GetCenter()) * k + 1) /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 * 1_m * 1_m / 4);
double b = directionBefore.dot(currentPosition - .GetCenter()) * 2 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 * 1_m / 4);
double c = ((currentPosition - .GetCenter()).GetSquaredNorm() - .GetRadius^2) /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 / 4);
std::complex<double>* solutions = solve_quartic(0, a, b, c);
std::vector<double> tmp;
for(int i = 0; i < 4; i++) {
if(solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType const Steplength;
if(tmp.size() > 0) {
Steplength = *std::min_element(tmp.begin(),tmp.end()) * 1_m;
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection with anything!" << std::endl;
//what to do when this happens?
}
// First Movement
//assuming magnetic field does not change during movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity = direction * velocity.GetNorm();
std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
<< " GeV " << std::endl;
} else {
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
<< " GeV " << std::endl;
}
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
geometry::Line line(currentPosition, velocity);
......@@ -172,8 +185,7 @@ namespace corsika::process {
<< std::endl;
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
velocity.norm() * min, minIter->second, magMaxLength,
directionBefore, directionAfter);
velocity.norm() * min, minIter->second);
}
};
......
......@@ -88,7 +88,7 @@ TEST_CASE("TrackingLine") {
0_m / second, 1_m / second);
Line line(origin, v);
auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
auto const [traj, geomMaxLength, nextVol, magMaxLength, beforeDirection, afterDirection] = tracking.GetTrack(p);
[[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
[[maybe_unused]] auto dummy_nextVol = nextVol;
......
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