IAP GITLAB

Skip to content

GitLab

C
corsika
Commit 4f3abf30 authored by Andre Schmidt's avatar Andre Schmidt
Browse files
Options

update ObservationPlane

parent 974b0364
Pipeline #2018 failed with stages
              in 0 seconds
              Hide whitespace changes
              Inline Side-by-side
              Showing with 84 additions and 167 deletions
              Framework/Cascade/Cascade.h
              View file @ 4f3abf30
              ......@@ -196,16 +196,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
              ......@@ -214,23 +214,23 @@ namespace corsika::cascade {
              std::cout << " move particle by : " << min_distance << std::endl;
              //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) *
              ......@@ -242,16 +242,18 @@ 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
              process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
              process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, stepWithoutB);
              if (status == process::EProcessReturn::eParticleAbsorbed) {
              std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
              ......
              Framework/Geometry/oCoordinateSystem.h deleted 100644 → 0
              View file @ 974b0364
              /*
              * (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
              View file @ 4f3abf30
              ......@@ -53,16 +53,54 @@ 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);
              return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
              auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
              return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
              //why is it *1.0001? should i do that too?
              }
              }
              Processes/TrackingLine/TrackingLine.h
              View file @ 4f3abf30
              ......@@ -80,7 +80,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());
              ......@@ -99,9 +99,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 /
              ......@@ -111,22 +111,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) *
              ......@@ -136,8 +136,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()) {
              ......@@ -162,9 +162,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 /
              ......@@ -174,22 +174,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) *
              ......@@ -199,7 +199,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);
              ......@@ -225,9 +225,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
              ......@@ -238,9 +239,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);
              }
              };
              ......
              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

              mailto ikp-admins [at] lists [dot] kit [dot] edu