IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 0d71c250 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

reverted stuff added for debugging

parent 2e730f78
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 48
− 82
View file @ 0d71c250
......@@ -216,9 +216,9 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
template <bool deleteParticle>
struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
//~ environment::BaseNodeType const& fA, fB;
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
......@@ -226,12 +226,7 @@ struct MyBoundaryCrossingProcess
EProcessReturn DoBoundaryCrossing(Particle& p,
typename Particle::BaseNodeType const& from,
typename Particle::BaseNodeType const& to) {
std::cout << "boundary crossing! from: " << from.GetModelProperties().GetName()
<< "; to: " << to.GetModelProperties().GetName() << std::endl;
//~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) {
p.Delete();
//~ }
std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
auto const& name = particles::GetName(p.GetPID());
auto const start = p.GetPosition().GetCoordinates();
......@@ -239,25 +234,15 @@ struct MyBoundaryCrossingProcess
fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
<< start[2] / 1_m << '\n';
if constexpr (deleteParticle) { p.Delete(); }
return EProcessReturn::eOk;
}
std::ofstream fFile;
void Init() {}
};
template <class T>
class TheNameModel : public T {
std::string const fName;
public:
template <typename... Args>
TheNameModel(std::string const& name, Args&&... args)
: T(std::forward<Args>(args)...)
, fName(name) {}
std::string const& GetName() const override { return fName; }
private:
std::ofstream fFile;
};
//
......@@ -273,41 +258,29 @@ int main() {
EnvType env;
auto& universe = *(env.GetUniverse());
auto outerMedium =
EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen
const float fox = 0.20946;
outerMedium->SetModelProperties<
TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
"outer", 1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{(float)1. - fox, fox}));
auto innerMedium = EnvType::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, -500_m}, 5000_m);
innerMedium->SetModelProperties<
TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
"inner", 1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{(float)1. - fox, fox}));
auto const props =
outerMedium
->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{1.f - fox, fox}));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
innerMedium->SetModelProperties(props);
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
universe.SetModelProperties<
TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
"Universe", 0_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
......@@ -319,17 +292,11 @@ int main() {
process::sibyll::Decay decay;
ProcessCut cut(20_GeV);
random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic;
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
MyBoundaryCrossingProcess boundaryCrossing("crossings.dat");
MyBoundaryCrossingProcess<false> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence =
/*p0 <<*/ sibyll << sibyllNuc << decay << cut /* << hadronicElastic */
<< boundaryCrossing << trackWriter;
auto sequence = sibyll << sibyllNuc << decay << cut << boundaryCrossing << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
......@@ -339,30 +306,29 @@ int main() {
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) * particles::Neutron::GetMass();
const HEPEnergyType E0 = nuclA * 10_TeV;
double phi = 0;
for (double theta = 0; theta < 180; theta += 20) {
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
for (int l = 0; l < 1; ++l) {
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType>{beamCode, E0, plab, pos, 0_ns});
}
}
const HEPEnergyType E0 = nuclA * 100_TeV;
double const phi = 0;
double const theta = 0;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns});
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
......
Framework/Cascade/Cascade.h
+ 9
− 15
View file @ 0d71c250
......@@ -20,6 +20,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <cassert>
#include <cmath>
#include <iostream>
#include <type_traits>
......@@ -91,9 +92,6 @@ namespace corsika::cascade {
auto const* numericalNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
p.SetNode(numericalNode);
std::cout << "initial node " << p.GetNode()->GetModelProperties().GetName()
<< std::endl;
});
}
......@@ -229,19 +227,15 @@ namespace corsika::cascade {
} else { // step-length limitation within volume
std::cout << "step-length limitation" << std::endl;
}
auto const* numericalNodeAfterStep =
fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
std::cout << "nodes: " << currentLogicalNode->GetModelProperties().GetName()
<< " " << numericalNodeAfterStep->GetModelProperties().GetName()
<< std::endl;
auto const assertion = [&] {
auto const* numericalNodeAfterStep =
fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
return numericalNodeAfterStep == currentLogicalNode;
};
if (numericalNodeAfterStep != currentLogicalNode) {
std::cout << "position " << particle.GetPosition().GetCoordinates()
<< std::endl;
throw std::runtime_error("numerical and logical nodes don't match");
}
} else { // boundary crossing
assert(assertion()); // numerical and logical nodes don't match
} else { // boundary crossing
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
particle.SetNode(nextVol);
fProcessSequence.DoBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
......@@ -255,7 +249,7 @@ namespace corsika::cascade {
Stack& fStack;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
};
}; // namespace corsika::cascade
} // namespace corsika::cascade
......
Processes/Sibyll/testSibyll.cc
+ 3
− 0
View file @ 0d71c250
......@@ -99,6 +99,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
......@@ -123,6 +124,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
Interaction model;
......@@ -147,6 +149,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
particle.SetNode(nodePtr);
Interaction hmodel;
NuclearInteraction model(hmodel);
......
Processes/TrackingLine/TrackingLine.h
+ 13
− 9
View file @ 0d71c250
......@@ -43,7 +43,7 @@ namespace corsika::process {
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
auto GetTrack(Particle const& p) {
auto GetTrack(Particle const& p) {
using namespace corsika::units::si;
using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
......@@ -52,8 +52,10 @@ namespace corsika::process {
auto const currentPosition = p.GetPosition();
std::cout << "TrackingLine pid: " << p.GetPID()
<< " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates() << " ["
<< p.GetNode()->GetModelProperties().GetName() << "]\n";
std::cout << "TrackingLine pos: "
<< currentPosition.GetCoordinates()
// << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
......@@ -84,8 +86,10 @@ namespace corsika::process {
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
std::cout << "intersection times: " << t1 / 1_s << "; " << t2 / 1_s << " "
<< vtn.GetModelProperties().GetName() << std::endl;
std::cout << "intersection times: " << t1 / 1_s << "; "
<< t2 / 1_s
// << " " << vtn.GetModelProperties().GetName()
<< std::endl;
if (t1.magnitude() > 0)
intersections.emplace_back(t1, &vtn);
else if (t2.magnitude() > 0)
......@@ -119,10 +123,10 @@ namespace corsika::process {
min = minIter->first;
}
std::cout << " t-intersect: " << min << " "
<< minIter->second->GetModelProperties().GetName() << std::endl;
//~ std::cout << "point of intersection: " <<
//line.GetPosition(min).GetCoordinates() << std::endl;
std::cout << " t-intersect: "
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
velocity.norm() * min, minIter->second);
......
Setup/SetupEnvironment.h
+ 3
− 3
View file @ 0d71c250
......@@ -12,13 +12,13 @@
#ifndef _include_corsika_setup_environment_h_
#define _include_corsika_setup_environment_h_
#include <corsika/environment/Environment.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NameModel.h>
#include <corsika/environment/Environment.h>
namespace corsika::setup {
using IEnvironmentModel = environment::NameModel<environment::IMediumModel>;
using IEnvironmentModel = environment::IMediumModel;
using SetupEnvironment = environment::Environment<IEnvironmentModel>;
}
} // namespace corsika::setup
#endif
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

data privacy - impressum