IAP GITLAB

Commit 51fd68a7 authored by André Schmidt's avatar André Schmidt

Merge branch 'andre-master' into 'master'

solved transition error

See merge request !11
parents c8b31f62 f2e86927
Pipeline #2031 failed with stages
......@@ -246,7 +246,8 @@ int main(int argc, char** argv) {
55_GeV);
auto decaySequence = decayPythia << decaySibyll;
auto sequence = switchProcess << decaySequence << longprof << eLoss << cut
<< observationLevel << trackWriter;
<< observationLevel;
// << trackWriter
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
......
......@@ -211,18 +211,17 @@ 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;
}
std::cout << "Decay after " << distance_decay << std::endl;
std::cout << "Interaction after " << distance_interact << std::endl;
std::cout << " move particle by : " << min_distance << std::endl;
// 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();
......@@ -232,19 +231,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:
......@@ -323,7 +322,7 @@ namespace corsika::cascade {
assert(assertion()); // numerical and logical nodes don't match
} else { // boundary crossing, step is limited by volume boundary
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
std::cout << "boundary crossing! current node = " << currentLogicalNode <<" next node = " << nextVol << std::endl;
vParticle.SetNode(nextVol);
// DoBoundary may delete the particle (or not)
fProcessSequence.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
......
......@@ -65,7 +65,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 /
......@@ -82,12 +82,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;
}
}
......
......@@ -114,12 +114,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)
......@@ -178,12 +179,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)
......
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