change Steplength

Documentation/Examples/vertical_EAS.cc

@@ -34,6 +34,7 @@
 #include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/proposal/ContinuousProcess.h>
 #include <corsika/process/proposal/Interaction.h>
+#include <corsika/process/track_writer/TrackWriter.h>
 #include <corsika/process/pythia/Decay.h>
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
@@ -226,6 +227,7 @@ int main(int argc, char** argv) {
 
   process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
 
+  process::track_writer::TrackWriter trackWriter("tracks.dat");
   process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
 
   Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));

Framework/Cascade/Cascade.h

@@ -214,7 +214,7 @@ namespace corsika::cascade {
                                         vParticle.GetEnergy() * units::constants::c;
                                     
       // determine geometric tracking
-      auto [step, geomMaxLength, nextVol, magMaxLength, directionBefore, directionAfter] = fTracking.GetTrack(vParticle);
+      auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
       [[maybe_unused]] auto const& dummy_nextVol = nextVol;
       
       // convert next_step from grammage to length
@@ -226,18 +226,46 @@ namespace corsika::cascade {
       LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
       std::cout << "distance_max=" << distance_max << std::endl;
 
-      // take minimum of geometry, interaction, decay for next step
+      // take minimum of geometry, interaction, decay, magnetic field for next step
       auto const min_distance = std::min(
-          {distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength});
+          {distance_interact, distance_decay, distance_max, geomMaxLength});
 
       C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
 
+      //determine displacement by the magnetic field
+	    auto const* currentLogicalVolumeNode = vParticle.GetNode();
+      int chargeNumber;
+      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);
+    		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;
+    		// Change of direction by magnetic field
+    		geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
+                                                                Steplength * k; 
+    		// Second Movement
+    		position = position + directionAfter * Steplength / 2;
+    		geometry::Vector<dimensionless_d> const direction = (position - vParticle.GetPosition()) / 
+    									                                      (position - vParticle.GetPosition()).GetNorm();
+        vParticle.SetMomentum(direction * vParticle.GetMomentum().GetNorm());
+      }
+
       // here the particle is actually moved along the trajectory to new position:
       // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
       vParticle.SetPosition(step.PositionFromArclength(min_distance));
       // .... also update time, momentum, direction, ...  
-	  vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) + 
-	  	directionAfter * min_distance /magMaxLength) * vParticle.GetMomentum().GetNorm());
       vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
 
       step.LimitEndTo(min_distance);
@@ -263,7 +291,7 @@ namespace corsika::cascade {
 
         TStackView secondaries(vParticle);
 
-        if (min_distance != distance_max && min_distance != magMaxLength) {
+        if (min_distance != distance_max) {
           /*
             Create SecondaryView object on Stack. The data container
             remains untouched and identical, and 'projectil' is identical

Framework/Utilities/CMakeLists.txt

@@ -20,6 +20,7 @@ set (
   UTILITIES_SOURCES  
   COMBoost.cc
   CorsikaData.cc
+  quartic.cpp
   ${CORSIKA_FENV})
 
 set (
@@ -31,6 +32,7 @@ set (
   sgn.h
   CorsikaFenv.h
   MetaProgramming.h
+  quartic.h
   )
 
 set (

Processes/TrackingLine/TrackingLine.h

@@ -15,8 +15,7 @@
 #include <corsika/geometry/Vector.h>
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
-#include <corsika/logging/Logging.h>
-
+#include <corsika/utl/quartic.h>
 #include <optional>
 #include <type_traits>
 #include <utility>
@@ -60,7 +59,7 @@ namespace corsika::process {
                   << std::endl;
         std::cout << "TrackingLine   E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
         
-        // determine velocity after adding magnetic field
+        // determine direction of the particle after adding magnetic field
         auto const* currentLogicalVolumeNode = p.GetNode();
         int chargeNumber;
         if(corsika::particles::IsNucleus(p.GetPID())) {
@@ -68,39 +67,51 @@ namespace corsika::process {
         } else {
         	chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
         }
-        geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
-        auto magMaxLength = 1_m/0;
-        auto directionAfter = directionBefore;
         if(chargeNumber != 0) {
         	auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
         	std::cout << "TrackingLine   B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
-        	geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
-        		velocity.parallelProjectionOnto(magneticfield);
-		LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V / 
-					      (corsika::units::constants::cSquared * abs(chargeNumber) * 
-					      magneticfield.GetNorm() * 1_eV);
-		//steplength depending on how exact it should be
-		LengthType const Steplength = 0.01 * gyroradius;
-		// First Movement
-		auto position = currentPosition + directionBefore * Steplength / 2;
-		// Change of direction by magnetic field at position
-		magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(position);
-		directionAfter = directionBefore + directionBefore.cross(magneticfield) * chargeNumber * 
-				 Steplength * corsika::units::constants::cSquared * 1_eV / 
-				 (p.GetEnergy() * velocity.GetNorm() * 1_V); 
-		// Second Movement
-		position = position + directionAfter * Steplength / 2;
-		magMaxLength = (position - currentPosition).GetNorm();
-		geometry::Vector<dimensionless_d> const direction = (position - currentPosition) / 
-									magMaxLength;
-		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();
+        	//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<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) {
+              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<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
+                                                                  Steplength * k;
+		      // Second Movement
+		      position = position + directionAfter * Steplength / 2;
+		      geometry::Vector<dimensionless_d> 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;
+                    << " GeV " << std::endl;
         }
-        
         std::cout << "TrackingLine   v: " << velocity.GetComponents() << std::endl;
         
         geometry::Line line(currentPosition, velocity);
@@ -166,8 +177,7 @@ namespace corsika::process {
         C8LOG_DEBUG(" t-intersect: {} ", min);
 
         return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
-                               velocity.norm() * min, minIter->second, magMaxLength, 
-                               directionBefore, directionAfter);
+                               velocity.norm() * min, minIter->second);
       }
     };
 

Processes/TrackingLine/testTrackingLine.cc

@@ -86,7 +86,7 @@ TEST_CASE("TrackingLine") {
                                                             0_m / second, 1_m / second);
     Line line(origin, v);
 
-    auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
+    auto const [traj, geomMaxLength, nextVol, magMaxLength, beforeDirection, afterDirection] = tracking.GetTrack(p);
     [[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
     [[maybe_unused]] auto dummy_nextVol = nextVol;