IAP GITLAB

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

fixing problems

parent cc13c82c
No related branches found
No related tags found
1 merge request!278Magnetic Tracking
...@@ -70,8 +70,10 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa ...@@ -70,8 +70,10 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
auto const* currentLogicalVolumeNode = vParticle.GetNode(); auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition()); auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
auto direction = trajectory.GetV0().normalized(); auto direction = trajectory.GetV0().normalized();
std::cout << " " << plane_.GetNormal().dot(trajectory.GetV0().cross(magneticfield)) << std::endl;
if (chargeNumber != 0 && abs(plane_.GetNormal().dot(trajectory.GetV0().cross(magneticfield))) * 1_s / 1_m / 1_T > 1e-6) { if (chargeNumber != 0 && abs(plane_.GetNormal().dot(trajectory.GetV0().cross(magneticfield))) * 1_s / 1_m / 1_T > 1e-5) {
auto const* currentLogicalVolumeNode = vParticle.GetNode(); auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition()); auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
auto k = chargeNumber * corsika::units::constants::c * 1_eV / (vParticle.GetMomentum().norm() * 1_V); auto k = chargeNumber * corsika::units::constants::c * 1_eV / (vParticle.GetMomentum().norm() * 1_V);
...@@ -93,11 +95,11 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa ...@@ -93,11 +95,11 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) { if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return std::numeric_limits<double>::infinity() * 1_m; return std::numeric_limits<double>::infinity() * 1_m;
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) { } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
std::cout << " distance to obs plane : " << MaxStepLength2 << std::endl; std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
return MaxStepLength2 * 1.0001; return MaxStepLength2 * (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2).norm() * 1.001;
} else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) { } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
std::cout << " distance to obs plane : " << MaxStepLength1 << std::endl; std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
return MaxStepLength1 * 1.0001; return MaxStepLength1 * (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2).norm() * 1.001;
} }
} }
TimeType const timeOfIntersection = TimeType const timeOfIntersection =
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment