diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index 7684de0b4f36d8b52287954cac20c489eed7580d..4f55633f6554b5686e3a1c17b52bdd37f2db4ab9 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -97,7 +97,7 @@ int main(int argc, char** argv) { node1->SetModelProperties>>( B, 1_g / (1_cm * 1_cm * 1e9_cm), composition); - auto node2 = std::make_unique>( + /*auto node2 = std::make_unique>( std::make_unique(center, seaLevel + 100.0_km)); node2->SetModelProperties>>( B, center, 540.1778_g / (1_cm * 1_cm) / 772170.16_cm, -772170.16_cm, composition, @@ -124,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 diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 9e0ae52b5e76496f003c9b67ffb30dadb51bd7e7..28492691feda56a3dfb25962ca5a4e7b0bc40c13 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -229,26 +229,24 @@ namespace corsika::cascade { geometry::Vector 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; + auto position = vParticle.GetPosition() + directionBefore * min_distance / 2; // Change of direction by magnetic field geometry::Vector const directionAfter = directionBefore + directionBefore.cross(magneticfield) * - Steplength * k; + min_distance * k; // Second Movement - position = position + directionAfter * Steplength / 2; - geometry::Vector const direction = (position - vParticle.GetPosition()) / - (position - vParticle.GetPosition()).GetNorm(); - vParticle.SetMomentum(direction * vParticle.GetMomentum().GetNorm()); + position = position + directionAfter * min_distance / 2; + // here the particle is actually moved along the trajectory to new position: + // std::visit(setup::ParticleUpdate{vParticle}, step); + vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm()); + vParticle.SetPosition(position); + } else { + vParticle.SetPosition(step.PositionFromArclength(min_distance)); } - - // here the particle is actually moved along the trajectory to new position: - // std::visit(setup::ParticleUpdate{vParticle}, step); - vParticle.SetPosition(step.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); diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h index 3e4aad42d0ae8ef94d7ce6cd5e8bd76016966c00..b90c5d3116ad3a3ad75130bc0b70b4073b7b3d62 100644 --- a/Processes/TrackingLine/TrackingLine.h +++ b/Processes/TrackingLine/TrackingLine.h @@ -61,76 +61,36 @@ namespace corsika::process { // << " [" << p.GetNode()->GetModelProperties().GetName() << "]" << std::endl; std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl; - - // determine direction of the particle after adding magnetic field - auto const* currentLogicalVolumeNode = p.GetNode(); - int chargeNumber; - if(corsika::particles::IsNucleus(p.GetPID())) { - chargeNumber = p.GetNuclearZ(); - } else { - chargeNumber = corsika::particles::GetChargeNumber(p.GetPID()); - } - if(chargeNumber != 0) { - auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition); - std::cout << "TrackingLine B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl; - auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V); - geometry::Vector 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* solutions = solve_quartic(0, a, b, c); - std::vector 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 const directionAfter = directionBefore + directionBefore.cross(magneticfield) * - Steplength * k; - // Second Movement - position = position + directionAfter * Steplength / 2; - geometry::Vector 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; - } + std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV + << " GeV " << std::endl; std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl; - geometry::Line line(currentPosition, velocity); - - //auto const* currentLogicalVolumeNode = p.GetNode(); + auto const* currentLogicalVolumeNode = p.GetNode(); //~ auto const* currentNumericalVolumeNode = //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition); auto const numericallyInside = currentLogicalVolumeNode->GetVolume().Contains(currentPosition); - std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false"); + std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false") << std::endl; auto const& children = currentLogicalVolumeNode->GetChildNodes(); auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes(); std::vector> intersections; + + //charge of the particle + int chargeNumber; + if(corsika::particles::IsNucleus(p.GetPID())) { + chargeNumber = p.GetNuclearZ(); + } else { + chargeNumber = corsika::particles::GetChargeNumber(p.GetPID()); + } + auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition); + std::cout << " Magnetic Field: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl; + auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V); + geometry::Vector const directionBefore = velocity.normalized(); + geometry::Vector velocity1 = velocity; + geometry::Vector velocity2 = velocity; // for entering from outside auto addIfIntersects = [&](auto const& vtn) { @@ -138,6 +98,47 @@ namespace corsika::process { auto const& sphere = dynamic_cast( 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 + 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 / + ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m); + double c = ((currentPosition - sphere.GetCenter()).GetSquaredNorm() - + (sphere.GetRadius() * sphere.GetRadius())) * 4 / + ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m); + std::complex* solutions = solve_quartic(0, a, b, c); + std::vector tmp; + 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) { + 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? + } + + // First Movement + //assuming magnetic field does not change during movement + auto position = currentPosition + directionBefore * Steplength / 2; + // Change of direction by magnetic field + geometry::Vector const directionAfter = directionBefore + directionBefore.cross(magneticfield) * + Steplength * k; + // Second Movement + position = position + directionAfter * Steplength / 2; + geometry::Vector 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 + geometry::Line line(currentPosition, velocity1); if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) { auto const [t1, t2] = *opt; @@ -160,6 +161,47 @@ namespace corsika::process { currentLogicalVolumeNode->GetVolume()); // 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 + 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 / + ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m); + double c = ((currentPosition - sphere.GetCenter()).GetSquaredNorm() - + (sphere.GetRadius() * sphere.GetRadius())) * 4 / + ((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m); + std::complex* solutions = solve_quartic(0, a, b, c); + std::vector tmp; + 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) { + 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? + } + + // First Movement + //assuming magnetic field does not change during movement + auto position = currentPosition + directionBefore * Steplength / 2; + // Change of direction by magnetic field + geometry::Vector const directionAfter = directionBefore + directionBefore.cross(magneticfield) * + Steplength * k; + // Second Movement + position = position + directionAfter * Steplength / 2; + geometry::Vector const direction = (position - currentPosition) / + (position - currentPosition).GetNorm(); + velocity2 = direction * velocity.GetNorm(); + }//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); [[maybe_unused]] auto dummy_t1 = t1; intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent()); @@ -183,6 +225,20 @@ namespace corsika::process { << min // << " " << minIter->second->GetModelProperties().GetName() << std::endl; + + // determine direction of the particle after adding magnetic field + //assuming magnetic field does not change during movement + // First Movement + auto position = currentPosition + velocity * min / 2; + // Change of direction by magnetic field + geometry::Vector const directionAfter = directionBefore + velocity.cross(magneticfield) * + min * k; + // Second Movement + position = position + directionAfter * velocity.norm() * min / 2; + geometry::Vector const direction = (position - currentPosition) / + (position - currentPosition).GetNorm(); + velocity = direction * velocity.GetNorm(); + geometry::Line line(currentPosition, velocity); return std::make_tuple(geometry::Trajectory(line, min), velocity.norm() * min, minIter->second); diff --git a/Processes/TrackingLine/TrackingLine_interpolation.h b/Processes/TrackingLine/TrackingLine_interpolation.h index 228855cfe2c0ba93e69f0f1aebc1d890fa5f9442..ea0c39d21b027fe6978e8a4283e1d3ba17be1ba9 100644 --- a/Processes/TrackingLine/TrackingLine_interpolation.h +++ b/Processes/TrackingLine/TrackingLine_interpolation.h @@ -115,6 +115,11 @@ namespace corsika::process { 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 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 << " GeV " << std::endl; @@ -143,6 +148,8 @@ namespace corsika::process { auto const& sphere = dynamic_cast( volume); // for the moment we are a bit bold here and assume // everything is a sphere, crashes with exception if not + + std::cout << "Mittelpunkt: " << sphere.GetCenter().GetCoordinates() << std::endl; if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) { auto const [t1, t2] = *opt;