IAP GITLAB

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

Merge branch 'andre-master' into 'master'

Andre master

See merge request !8
parents 0b3ca6f3 a3fc176b
Pipeline #2020 failed with stages
......@@ -97,7 +97,7 @@ int main(int argc, char** argv) {
node1->SetModelProperties<UniformMagneticField<HomogeneousMedium<IEnvModel>>>(
B, 1_g / (1_cm * 1_cm * 1e9_cm), composition);
/*auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
std::make_unique<geometry::Sphere>(center, seaLevel + 100.0_km));
node2->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
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
......
......@@ -196,62 +196,67 @@ 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
std::cout << "Interaction: " << distance_interact << std::endl;
std::cout << "Decay: " << distance_decay << std::endl;
std::cout << "ObsPlane: " << distance_max << std::endl;
std::cout << "Transition: " << geomMaxLength << std::endl;
auto const min_distance = std::min(
{distance_interact, distance_decay, distance_max, geomMaxLength});
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();
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();
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());
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
auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
min_distance * k;
// Second Movement
position = position + directionAfter * min_distance / 2;
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
vParticle.SetPosition(position);
} else {
vParticle.SetPosition(step.PositionFromArclength(min_distance));
}
// .... also update time, momentum, direction, ...
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
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
auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
min_distance * k;
// Second Movement
position = position + directionAfter * min_distance / 2;
auto distance = position - vParticle.GetPosition();
//distance.norm() != min_distance for distance_interact, distance_decay if q != 0
//small error can be neglected
velocity = distance.normalized() * velocity.norm();
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
geometry::Line line(vParticle.GetPosition(), velocity);
geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / velocity.GetNorm());
vParticle.SetPosition(position);
vParticle.SetTime(vParticle.GetTime() + distance.norm() / units::constants::c);
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
step.LimitEndTo(min_distance);
// apply all continuous processes on particle + track
process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, stepNew);
if (status == process::EProcessReturn::eParticleAbsorbed) {
std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
......
......@@ -53,11 +53,46 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
}
}
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vParticle,
setup::Trajectory const& trajectory) {
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);
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return std::numeric_limits<double>::infinity() * 1_m;
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
return MaxStepLength2 * 1.0001;
} else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
return MaxStepLength1 * 1.0001;
}
}
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
......
......@@ -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