solved transition error

Framework/Cascade/Cascade.h

@@ -229,18 +229,15 @@ namespace corsika::cascade {
       // take minimum of geometry, interaction, decay for next step
       auto min_distance = std::min(
           {distance_interact, distance_decay, distance_max, geomMaxLength});
-      if (min_distance == geomMaxLength) {
-        min_distance = 1.001 * geomMaxLength;
-      } 
 
       C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
 
       // determine displacement by the magnetic field
-	    auto const* currentLogicalVolumeNode = vParticle.GetNode();
+	  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();
+  	  geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
       int chargeNumber;
       if (corsika::particles::IsNucleus(vParticle.GetPID())) {
         chargeNumber = vParticle.GetNuclearZ();
@@ -250,19 +247,19 @@ namespace corsika::cascade {
       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) *
+	  // 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;
+  	  // Second Movement
+  	  position = position + directionAfter * min_distance / 2;
       auto distance = position - vParticle.GetPosition();
       // distance.norm() != min_distance if q != 0
       // small error can be neglected
       if (distance.norm() != 0_m) {
-      velocity = distance.normalized() * velocity.norm();
+		velocity = distance.normalized() * velocity.norm();
       } // no velocity update for very small steps
       
       // here the particle is actually moved along the trajectory to new position:

Processes/ObservationPlane/ObservationPlane.cc

@@ -71,7 +71,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
   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) {
+  if (chargeNumber != 0 && plane_.GetNormal().dot(velocity.cross(magneticfield)) * 1_s / 1_m / 1_T > 0.0001) {
     auto const* currentLogicalVolumeNode = vParticle.GetNode();
     auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
     auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / 
@@ -88,12 +88,13 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
       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) {
+	  std::cout << " distance to obs plane : " << MaxStepLength2 << std::endl;
       return MaxStepLength2 * 1.0001;
     } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
+	  std::cout << " distance to obs plane : " << MaxStepLength1 << std::endl;
       return MaxStepLength1 * 1.0001;
     }
   } 

Processes/TrackingLine/TrackingLine.h

@@ -111,12 +111,13 @@ namespace corsika::process {
             for (int i = 0; i < 4; i++) {
               if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
                 tmp.push_back(solutions[i].real());
+                std::cout << "Solutions for next Volume: " << solutions[i].real() << std::endl;
               }
             }
             LengthType Steplength;
             if (tmp.size() > 0) {
               Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
-              std::cout << "s = " << Steplength << std::endl;
+              std::cout << "Steplength to next volume = " << Steplength << std::endl;
             } else {
               std::cout << "no intersection (1)!" << std::endl;
               // what to do when this happens? (very unlikely)
@@ -172,12 +173,18 @@ namespace corsika::process {
             for (int i = 0; i < 4; i++) {
               if (solutions[i].imag() == 0 && solutions[i].real() > 0) {
                 tmp.push_back(solutions[i].real());
+                std::cout << "Solutions for current Volume: " << solutions[i].real() << std::endl;
               }
             }
             LengthType Steplength;
             if (tmp.size() > 0) {
-              Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
-              std::cout << "s = " << Steplength << std::endl;
+				Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
+				if (numericallyInside == false) {
+					int p = std::min_element(tmp.begin(),tmp.end()) - tmp.begin();
+					tmp.erase(tmp.begin() + p);
+					Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
+				}
+				std::cout << "steplength out of current volume = " << Steplength << std::endl;
             } else {
               std::cout << "no intersection (2)!" << std::endl;
               // what to do when this happens? (very unlikely)