IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 15ae8ca8 authored by ralfulrich's avatar ralfulrich
Browse files

improve boundary_example.

parent b6a74928
No related branches found
No related tags found
1 merge request!269Resolve "Distinguish between alternative interaction processes and competing processes"
...@@ -29,6 +29,7 @@ target_link_libraries (boundary_example ...@@ -29,6 +29,7 @@ target_link_libraries (boundary_example
CORSIKAgeometry CORSIKAgeometry
CORSIKAenvironment CORSIKAenvironment
CORSIKAprocesssequence CORSIKAprocesssequence
C8::ext::boost # boost::histogram
) )
CORSIKA_ADD_EXAMPLE (cascade_example) CORSIKA_ADD_EXAMPLE (cascade_example)
...@@ -120,6 +121,7 @@ CORSIKA_ADD_EXAMPLE (stopping_power stopping_power) ...@@ -120,6 +121,7 @@ CORSIKA_ADD_EXAMPLE (stopping_power stopping_power)
target_link_libraries (stopping_power target_link_libraries (stopping_power
CORSIKAsetup CORSIKAsetup
CORSIKAunits CORSIKAunits
CORSIKAlogging
ProcessEnergyLoss ProcessEnergyLoss
CORSIKAparticles CORSIKAparticles
CORSIKAgeometry CORSIKAgeometry
......
...@@ -101,24 +101,23 @@ int main() { ...@@ -101,24 +101,23 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem(); const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>( // create "world" as infinite sphere filled with protons
auto world = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
auto const props = auto const props =
outerMedium world->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>( 1_kg / (1_m * 1_m * 1_m),
1_kg / (1_m * 1_m * 1_m), environment::NuclearComposition(
environment::NuclearComposition( std::vector<particles::Code>{particles::Code::Proton},
std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.f}));
std::vector<float>{1.f}));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km); // add a "target" sphere with 5km readius at 0,0,0
auto target = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
target->SetModelProperties(props);
innerMedium->SetModelProperties(props); world->AddChild(std::move(target));
universe.AddChild(std::move(world));
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
// setup processes, decays and interactions // setup processes, decays and interactions
tracking_line::TrackingLine tracking; tracking_line::TrackingLine tracking;
...@@ -129,7 +128,7 @@ int main() { ...@@ -129,7 +128,7 @@ int main() {
process::particle_cut::ParticleCut cut(50_GeV, true, true); process::particle_cut::ParticleCut cut(50_GeV, true, true);
process::track_writer::TrackWriter trackWriter("tracks.dat"); process::track_writer::TrackWriter trackWriter("boundary_tracks.dat");
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat"); MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list // assemble all processes into an ordered process list
...@@ -147,8 +146,8 @@ int main() { ...@@ -147,8 +146,8 @@ int main() {
std::mt19937 rng; std::mt19937 rng;
for (int i = 0; i < 100; ++i) { for (int i = 0; i < 100; ++i) {
auto const theta = distTheta(rng); double const theta = distTheta(rng);
auto const phi = distPhi(rng); double const phi = distPhi(rng);
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m)); return sqrt((Elab - m) * (Elab + m));
...@@ -166,6 +165,7 @@ int main() { ...@@ -166,6 +165,7 @@ int main() {
"input angles: theta={} phi={}" "input angles: theta={} phi={}"
"input momentum: {} GeV", "input momentum: {} GeV",
beamCode, theta, phi, plab.GetComponents() / 1_GeV); beamCode, theta, phi, plab.GetComponents() / 1_GeV);
// shoot particles from inside target out
Point pos(rootCS, 0_m, 0_m, 0_m); Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle( stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType, std::tuple<particles::Code, units::si::HEPEnergyType,
...@@ -180,8 +180,8 @@ int main() { ...@@ -180,8 +180,8 @@ int main() {
C8LOG_INFO("Result: E0={}GeV", E0 / 1_GeV); C8LOG_INFO("Result: E0={}GeV", E0 / 1_GeV);
cut.ShowResults(); cut.ShowResults();
const HEPEnergyType Efinal = [[maybe_unused]] const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy(); (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy());
C8LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV, C8LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV,
(Efinal / E0 - 1.) * 100); (Efinal / E0 - 1.) * 100);
} }
...@@ -64,10 +64,9 @@ namespace corsika::process { ...@@ -64,10 +64,9 @@ namespace corsika::process {
geometry::Line line(currentPosition, velocity); geometry::Line line(currentPosition, velocity);
auto const* currentLogicalVolumeNode = p.GetNode(); auto const* currentLogicalVolumeNode = p.GetNode();
auto const numericallyInside =
currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
C8LOG_DEBUG("numericallyInside = {} ", numericallyInside); C8LOG_DEBUG("numericallyInside = {} ",
currentLogicalVolumeNode->GetVolume().Contains(currentPosition));
auto const& children = currentLogicalVolumeNode->GetChildNodes(); auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes(); auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment