IAP GITLAB

Commit a5f10ee3 authored by André Schmidt's avatar André Schmidt

Merge branch 'andre-master' into 'master'

save before change

See merge request !12
parents 51fd68a7 2f4918e1
Pipeline #2036 canceled with stages
......@@ -82,6 +82,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
(plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return std::numeric_limits<double>::infinity() * 1_m;
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
......
/*
* (c) Copyright 2019 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.
*/
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <fstream>
using namespace corsika::process::observation_plane;
using namespace corsika::units::si;
ObservationPlane::ObservationPlane(
geometry::Plane const& obsPlane,
geometry::Vector<units::si::dimensionless_d> const& x_axis,
std::string const& filename, bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
, deleteOnHit_(deleteOnHit)
, xAxis_(x_axis.normalized())
, yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;
}
corsika::process::EProcessReturn ObservationPlane::DoContinuous(
setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
return process::EProcessReturn::eOk;
}
auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
<< particle.GetEnergy() * (1 / 1_eV) << ' '
<< displacement.dot(xAxis_) / 1_m << ' ' << displacement.dot(yAxis_) / 1_m
<< ' ' << std::endl;
if (deleteOnHit_) {
return process::EProcessReturn::eParticleAbsorbed;
} else {
return process::EProcessReturn::eOk;
}
}
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const& trajectory) {
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;
}
auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
}
......@@ -54,20 +54,6 @@ namespace corsika::process {
geometry::Vector<SpeedType::dimension_type> velocity =
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
std::complex<double>* solutions = solve_quartic(1, 0, 1, -20);
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());
}
}
double s = *std::min_element(tmp.begin(),tmp.end());
std::cout << "s = " << s << std::endl;
std::cout << "x1 = " << (solutions[0].real()>=0. ? " " : "") << solutions[0].real(); if(solutions[0].imag()!=0.0) std::cout << " + i * " << solutions[0].imag(); std::cout << std::endl;
std::cout << "x2 = " << (solutions[1].real()>=0. ? " " : "") << solutions[1].real(); if(solutions[1].imag()!=0.0) std::cout << " - i * " << -solutions[1].imag(); std::cout << std::endl;
std::cout << "x3 = " << (solutions[2].real()>=0. ? " " : "") << solutions[2].real(); if(solutions[2].imag()!=0.0) std::cout << " + i * " << solutions[2].imag(); std::cout << std::endl;
std::cout << "x4 = " << (solutions[3].real()>=0. ? " " : "") << solutions[3].real(); if(solutions[3].imag()!=0.0) std::cout << " - i * " << -solutions[3].imag(); std::cout << std::endl;
auto const currentPosition = p.GetPosition();
std::cout << "TrackingLine pid: " << p.GetPID()
<< " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
......@@ -97,8 +83,11 @@ namespace corsika::process {
(corsika::units::constants::cSquared * abs(chargeNumber) *
magneticfield.GetNorm() * 1_eV);
//steplength should consider more things than just gyroradius
double maxAngle = 0.1;
double maxAngle = 1e-5;
LengthType const Steplength = 2 * sin(maxAngle * M_PI / 180) * gyroradius;
std::cout << "Test: " << Steplength << " " << gyroradius << std::endl;
// First Movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field at position
......@@ -114,11 +103,6 @@ namespace corsika::process {
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();
double test =((directionBefore.cross(magneticfield)).dot(position-currentPosition) * k + 1) / (1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
std::cout << "Test: " << test << k << std::endl;
} else {
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
......
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