diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 1b41d0a68fe8f8226815e0ef223741ab57458e20..9896259cb0800160e7c01640c4ca0a68e1164581 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -408,7 +408,7 @@ namespace corsika::cascade { // determine steplength for the magnetic field // because Steplength should not be min_distance - /* + int chargeNumber; if (corsika::particles::IsNucleus(vParticle.GetPID())) { chargeNumber = vParticle.GetNuclearZ(); @@ -419,16 +419,21 @@ namespace corsika::cascade { auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition()); geometry::Vector<dimensionless_d> const directionBefore = vParticle.GetMomentum().normalized(); auto c = directionBefore.cross(magneticfield) * chargeNumber * corsika::units::constants::c * 1_eV / - (vParticle.GetMomentum().norm() * 1_V) + (vParticle.GetMomentum().norm() * 1_V); LengthType Steplength = min_distance; if (chargeNumber != 0) { Steplength = sqrt(2 / c.squaredNorm() * (sqrt(c.squaredNorm() * min_distance * min_distance + 1) -1)); std::cout << "Steplength " << Steplength << std::endl; - } - */ + } + if (Steplength == 0_m) { + Steplength = min_distance; + } + // This formula hasnt been tested - auto [position, direction, L2] = fTracking.MagneticStep(vParticle, min_distance); + //auto [position, direction, L2] = fTracking.MagneticStep(vParticle, min_distance); + auto [position, direction, L2] = fTracking.MagneticStep(vParticle, Steplength); + //histL2(L2); histLlog2(L2); int pdg = static_cast<int>(particles::GetPDG(vParticle.GetPID())); @@ -481,7 +486,7 @@ namespace corsika::cascade { - int chargeNumber = 0; + //int chargeNumber = 0; if (corsika::particles::IsNucleus(vParticle.GetPID())) { chargeNumber = vParticle.GetNuclearZ(); } else { diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc index 036bbcacfaf32379ba9633035ce1d5ea7d7cf19a..ff70be5fe18012723d1195b03842fbc16832432c 100644 --- a/Processes/ObservationPlane/ObservationPlane.cc +++ b/Processes/ObservationPlane/ObservationPlane.cc @@ -71,7 +71,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition()); auto direction = trajectory.GetV0().normalized(); - if (chargeNumber != 0 && 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-6) { auto const* currentLogicalVolumeNode = vParticle.GetNode(); auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition()); auto k = chargeNumber * corsika::units::constants::c * 1_eV / (vParticle.GetMomentum().norm() * 1_V);