IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 7666c11b authored by Andre Schmidt's avatar Andre Schmidt Committed by ralfulrich
Browse files

update ObservationPlane

parent 0c5ad7a6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/Cascade/Cascade.h
+ 15
− 12
View file @ 7666c11b
......@@ -214,16 +214,16 @@ namespace corsika::cascade {
vParticle.GetEnergy() * units::constants::c;
// determine geometric tracking
auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
auto [stepWithoutB, stepWithB, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length
LengthType const distance_interact =
currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
currentLogicalNode->GetModelProperties().ArclengthFromGrammage(stepWithB,
next_interact);
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, stepWithoutB);
std::cout << "distance_max=" << distance_max << std::endl;
// take minimum of geometry, interaction, decay for next step
......@@ -232,23 +232,23 @@ namespace corsika::cascade {
C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
//determine displacement by the magnetic field
// determine displacement by the magnetic field
auto const* currentLogicalVolumeNode = vParticle.GetNode();
int chargeNumber;
if(corsika::particles::IsNucleus(vParticle.GetPID())) {
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());
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);
// First Movement
//assuming magnetic field does not change during movement
// assuming magnetic field does not change during movement
auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
......@@ -260,17 +260,20 @@ namespace corsika::cascade {
vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
vParticle.SetPosition(position);
} else {
vParticle.SetPosition(step.PositionFromArclength(min_distance));
vParticle.SetPosition(stepWithoutB.PositionFromArclength(min_distance));
}
// .... also update time, momentum, direction, ...
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
step.LimitEndTo(min_distance);
// is this necessary?
stepWithoutB.LimitEndTo(min_distance);
stepWithB.LimitEndTo(min_distance);
// apply all continuous processes on particle + track
if (process_sequence_.DoContinuous(vParticle, step) ==
process::EProcessReturn::eParticleAbsorbed) {
process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, stepWithoutB);
if (status == process::EProcessReturn::eParticleAbsorbed) {
C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV);
if (!vParticle.isDeleted()) vParticle.Delete();
......
Framework/Geometry/oCoordinateSystem.h deleted 100644 → 0
+ 0
− 125
View file @ 0c5ad7a6
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_COORDINATESYSTEM_H_
#define _include_COORDINATESYSTEM_H_
#include <corsika/geometry/QuantityVector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/sgn.h>
#include <Eigen/Dense>
#include <stdexcept>
typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
typedef Eigen::Translation<double, 3> EigenTranslation;
namespace corsika::geometry {
class RootCoordinateSystem;
template <typename T>
class Vector;
using corsika::units::si::length_d;
class CoordinateSystem {
CoordinateSystem const* reference = nullptr;
EigenTransform transf;
CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
: reference(&reference)
, transf(transf) {}
CoordinateSystem()
: // for creating the root CS
transf(EigenTransform::Identity()) {}
protected:
static auto CreateCS() { return CoordinateSystem(); }
friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
/// create ONE unique root CS
public:
static EigenTransform GetTransformation(CoordinateSystem const& c1,
CoordinateSystem const& c2);
auto& operator=(const CoordinateSystem& pCS) {
reference = pCS.reference;
transf = pCS.transf;
return *this;
}
auto translate(QuantityVector<length_d> vector) const {
EigenTransform const translation{EigenTranslation(vector.eVector)};
return CoordinateSystem(*this, translation);
}
template <typename TDim>
auto RotateToZ(Vector<TDim> vVec) const {
auto const a = vVec.normalized().GetComponents(*this).eVector;
auto const a1 = a(0), a2 = a(1);
auto const s = utl::sgn(a(2));
auto const c = 1 / (1 + s * a(2));
Eigen::Matrix3d A, B;
if (s > 0) {
A << 1, 0, a1, // comment to prevent clang-format
0, 1, a2, // .
-a1, -a2, 1; // .
B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-a1 * a2 * c, -a2 * a2 * c, 0, // .
0, 0, -(a1 * a1 + a2 * a2) * c; // .
} else {
A << 1, 0, a1, // .
0, -1, a2, // .
a1, -a2, -1; // .
B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
-a1 * a2 * c, +a2 * a2 * c, 0, // .
0, 0, (a1 * a1 + a2 * a2) * c; // .
}
return CoordinateSystem(*this, EigenTransform(A + B));
}
template <typename TDim>
auto rotate(QuantityVector<TDim> axis, double angle) const {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
return CoordinateSystem(*this, rotation);
}
template <typename TDim>
auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
QuantityVector<TDim> axis, double angle) {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
EigenTranslation(translation.eVector)};
return CoordinateSystem(*this, transf);
}
auto const* GetReference() const { return reference; }
auto const& GetTransform() const { return transf; }
};
} // namespace corsika::geometry
#endif
Processes/ObservationPlane/ObservationPlane.cc
+ 45
− 9
View file @ 7666c11b
......@@ -59,20 +59,56 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
}
}
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vParticle,
setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
int chargeNumber;
if (corsika::particles::IsNucleus(vParticle.GetPID())) {
chargeNumber = vParticle.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
}
auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
geometry::Vector<SpeedType::dimension_type> const velocity = trajectory.GetV0();
if (chargeNumber != 0 && plane_.GetNormal().dot(velocity.cross(magneticfield)) * 1_s / 1_m / 1_T != 0) {
auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(velocity.GetSquaredNorm() * vParticle.GetEnergy() * 1_V);
LengthType MaxStepLength1 =
( sqrt(velocity.dot(plane_.GetNormal()) * velocity.dot(plane_.GetNormal()) / velocity.GetSquaredNorm() -
(plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
(plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
LengthType MaxStepLength2 =
( - sqrt(velocity.dot(plane_.GetNormal()) * velocity.dot(plane_.GetNormal()) / velocity.GetSquaredNorm() -
(plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
(plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
std::cout << "Test: " << MaxStepLength1 << " " << MaxStepLength2 << std::endl;
if (MaxStepLength1 < 0_m && MaxStepLength2 < 0_m) {
return std::numeric_limits<double>::infinity() * 1_m;
} else if (MaxStepLength1 < 0_m || MaxStepLength2 < MaxStepLength1) {
return MaxStepLength2;
} else if (MaxStepLength2 < 0_m || MaxStepLength1 < MaxStepLength2) {
return MaxStepLength1;
}
} else {
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
auto dist = (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
C8LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m);
return dist;
auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
//why is it *1.0001? should i do that too?
}
}
void ObservationPlane::ShowResults() const {
......
Processes/TrackingLine/TrackingLine.h
+ 26
− 24
View file @ 7666c11b
......@@ -77,7 +77,7 @@ namespace corsika::process {
//charge of the particle
int chargeNumber;
if(corsika::particles::IsNucleus(p.GetPID())) {
if (corsika::particles::IsNucleus(p.GetPID())) {
chargeNumber = p.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
......@@ -96,9 +96,9 @@ namespace corsika::process {
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
//creating Line with magnetic field
if(chargeNumber != 0) {
//determine steplength to next volume
// creating Line with magnetic field
if (chargeNumber != 0) {
// determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
......@@ -108,22 +108,22 @@ namespace corsika::process {
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
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) {
for (int i = 0; i < 4; i++) {
if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType Steplength;
if(tmp.size() > 0) {
if (tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection (1)!" << std::endl;
//what to do when this happens?
// what to do when this happens? (very unlikely)
}
// First Movement
//assuming magnetic field does not change during 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) *
......@@ -133,8 +133,8 @@ namespace corsika::process {
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity1 = direction * velocity.GetNorm();
} //instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
//line has some errors
} // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
// line has some errors
geometry::Line line(currentPosition, velocity1);
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
......@@ -156,9 +156,9 @@ namespace corsika::process {
// for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
//creating Line with magnetic field
if(chargeNumber != 0) {
//determine steplength to next volume
// creating Line with magnetic field
if (chargeNumber != 0) {
// determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
......@@ -168,22 +168,22 @@ namespace corsika::process {
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
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) {
for (int i = 0; i < 4; i++) {
if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType Steplength;
if(tmp.size() > 0) {
if (tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection (2)!" << std::endl;
//what to do when this happens?
// what to do when this happens? (very unlikely)
}
// First Movement
//assuming magnetic field does not change during 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) *
......@@ -193,7 +193,7 @@ namespace corsika::process {
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity2 = direction * velocity.GetNorm();
}//instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
} // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
geometry::Line line(currentPosition, velocity2);
[[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
......@@ -219,9 +219,10 @@ namespace corsika::process {
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
geometry::Line lineWithoutB(currentPosition, velocity);
// determine direction of the particle after adding magnetic field
//assuming magnetic field does not change during movement
// assuming magnetic field does not change during movement
// First Movement
auto position = currentPosition + velocity * min / 2;
// Change of direction by magnetic field
......@@ -232,9 +233,10 @@ namespace corsika::process {
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity = direction * velocity.GetNorm();
geometry::Line line(currentPosition, velocity);
geometry::Line lineWithB(currentPosition, velocity);
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
return std::make_tuple(geometry::Trajectory<geometry::Line>(lineWithoutB, min),
geometry::Trajectory<geometry::Line>(lineWithB, min),
velocity.norm() * min, minIter->second);
}
};
......
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

data privacy - impressum