IAP GITLAB

Commit 90f2ff8d authored by Andre Schmidt's avatar Andre Schmidt

sth

parent 4f4bbe0d
Pipeline #2525 failed with stages
in 0 seconds
......@@ -246,7 +246,7 @@ 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
......
......@@ -64,13 +64,18 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
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-5) {
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);
if (direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
(plane_.GetNormal().dot(trajectory.GetR0() - plane_.GetCenter()) *
plane_.GetNormal().dot(direction.cross(magneticfield)) * 2 * k) < 0)
{
return std::numeric_limits<double>::infinity() * 1_m;
}
LengthType MaxStepLength1 =
( sqrt(direction.dot(plane_.GetNormal()) * direction.dot(plane_.GetNormal()) -
......
......@@ -430,9 +430,13 @@ namespace corsika::process {
<< std::endl;
//if the used Steplengths are too big, min gets false and we do another Step
std::cout << Steplength1 << Steplength2 << std::endl;
std::cout << Steplength1 << " " << Steplength2 << std::endl;
std::cout << intersection << " " << Steplimit << std::endl;
if ((Steplength1 > Steplimit || Steplength1 == 0_m) && Steplength2 > Steplimit) {
std::cout << "Test: (1) " << (Steplength1 > Steplimit || Steplength1 < 1e-6_m) << std::endl;
std::cout << "Test: (2) " << (Steplength2 > Steplimit) << std::endl;
if ((Steplength1 > Steplimit || Steplength1 < 1e-6_m) && Steplength2 > Steplimit) {
min = Steplimit / velocity.norm();
auto lineWithB = MagneticStep(p, line, Steplimit, true);
......@@ -443,6 +447,8 @@ namespace corsika::process {
histSlog(lineWithB.ArcLength(0_s, min) / 1_m);
//histLS(L-lineWithB.ArcLength(0_s,min) / 1_m);
std::cout << "is it here? TrackingLine.h (2)" << std::endl;
if (chargeNumber != 0) {
if (L * 1_m / lineWithB.ArcLength(0_s, min) - 1 > 0) {
histLSrelmag(L * 1_m / lineWithB.ArcLength(0_s, min) - 1);
......@@ -477,8 +483,6 @@ namespace corsika::process {
//histB(lineWithB.ArcLength(0_s,min) / 1_m);
histBlog(lineWithB.ArcLength(0_s, min) / 1_m);
}
std::cout << "Kontrolle: " << 1_m / 0 - 1_m << std::endl;
int pdg = static_cast<int>(particles::GetPDG(p.GetPID()));
if (abs(pdg) == 13) {
......@@ -505,6 +509,8 @@ namespace corsika::process {
Steplimit * 1.0001, Steplimit, minIter->second);
}
std::cout << "is it here? TrackingLine.h (3)" << std::endl;
auto lineWithB = MagneticStep(p, line, velocity.norm() * min, true);
//histL(L);
......@@ -530,18 +536,22 @@ namespace corsika::process {
(corsika::units::constants::cSquared * abs(chargeNumber) *
magneticfield.GetNorm() * 1_eV);
std::cout << "Schrittlaenge " << velocity.norm() * min << " gyroradius " << gyroradius << std::endl;
LengthType B = 2 * gyroradius * asin(lineWithB.ArcLength(0_s, min) / 2 / gyroradius);
std::cout << "Bogenlaenge" << B << std::endl;
//histB(B / 1_m);
histBlog(B / 1_m);
if (L - B / 1_m > 0) {
//histLB(L-B/1_m);
//histBS(B/1_m - lineWithB.ArcLength(0_s,min) / 1_m);
if (B != 0_m)
histLBrelgeo(L * 1_m / B - 1);
}
else {
Bbig.push_back(B / 1_m);
//irgendetwas geht schief: Schrittlnge viel grer als gyroradius
if (lineWithB.ArcLength(0_s, min) < gyroradius) {
LengthType B = 2 * gyroradius * asin(lineWithB.ArcLength(0_s, min) / 2 / gyroradius);
std::cout << "Bogenlaenge" << B << std::endl;
//histB(B / 1_m);
histBlog(B / 1_m);
if (L - B / 1_m > 0) {
//histLB(L-B/1_m);
//histBS(B/1_m - lineWithB.ArcLength(0_s,min) / 1_m);
if (B != 0_m)
histLBrelgeo(L * 1_m / B - 1);
}
else {
Bbig.push_back(B / 1_m);
}
}
}
else {
......
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