IAP GITLAB

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

atmosphere set to 295 K isothermic flat

parent ff3c78bc
No related branches found
No related tags found
1 merge request!121Vertical EAS
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
#include <corsika/setup/SetupTrajectory.h> #include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h> #include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h> #include <corsika/environment/FlatExponential.h>
#include <corsika/environment/NuclearComposition.h> #include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h> #include <corsika/geometry/Sphere.h>
...@@ -84,10 +84,7 @@ public: ...@@ -84,10 +84,7 @@ public:
} else { } else {
// TODO: center-of-mass energy hard coded // TODO: center-of-mass energy hard coded
const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV); const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
if (p.GetEnergy() < fECut || Ecm < 10_GeV) return (p.GetEnergy() < fECut || Ecm < 10_GeV);
return true;
else
return false;
} }
} }
...@@ -124,18 +121,13 @@ public: ...@@ -124,18 +121,13 @@ public:
case Code::NuEBar: case Code::NuEBar:
is_inv = true; is_inv = true;
break; break;
case Code::NuMu: case Code::NuMu:
is_inv = true; is_inv = true;
break; break;
case Code::NuMuBar: case Code::NuMuBar:
is_inv = true; is_inv = true;
break; break;
case Code::MuPlus:
is_inv = true;
break;
case Code::MuMinus:
is_inv = true;
break;
case Code::Neutron: case Code::Neutron:
is_inv = true; is_inv = true;
...@@ -210,7 +202,6 @@ public: ...@@ -210,7 +202,6 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; } HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
}; };
// //
// The example main program for a particle cascade // The example main program for a particle cascade
// //
...@@ -225,19 +216,24 @@ int main() { ...@@ -225,19 +216,24 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem(); const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto& universe = *(env.GetUniverse()); auto& universe = *(env.GetUniverse());
auto theMedium = auto theMedium = EnvType::CreateNode<Sphere>(
EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen auto constexpr temperature = 295_K; // AIRES default temperature for isothermal model
auto constexpr lambda =
-constants::R * temperature / (constants::g_sub_n * 28.966_g / mole);
// 1036 g/cm² taken from AIRES code
auto constexpr rho0 = -1036_g / square(1_cm) / lambda;
using FlatExp = environment::FlatExponential<environment::IMediumModel>; using FlatExp = environment::FlatExponential<environment::IMediumModel>;
theMedium->SetModelProperties<FlatExp>( theMedium->SetModelProperties<FlatExp>(
Point{rootCS, 0_m, 0_m, }, Point{rootCS, 0_m, 0_m, 0_m}, Vector<dimensionless_d>{rootCS, {0., 0., 1.}}, rho0,
lambda,
environment::NuclearComposition( environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen, std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen, particles::Code::Oxygen},
particles::Code::Argon}, std::vector<float>{
std::vector<float>{0.7847,0.2105,0.005})); // values taken from AIRES manual 0.7847f,
1.f - 0.7847f})); // values taken from AIRES manual, Ar removed for now
universe.AddChild(std::move(theMedium)); universe.AddChild(std::move(theMedium));
...@@ -273,11 +269,9 @@ int main() { ...@@ -273,11 +269,9 @@ int main() {
// setup particle stack, and add primary particle // setup particle stack, and add primary particle
setup::Stack stack; setup::Stack stack;
stack.Clear(); stack.Clear();
const Code beamCode = Code::Nucleus; const Code beamCode = Code::Proton;
const int nuclA = 4; auto const mass = particles::GetMass(beamCode);
const int nuclZ = int(nuclA / 2.15 + 0.7); const HEPEnergyType E0 = 10_TeV;
const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
const HEPEnergyType E0 = nuclA * 10_TeV;
double theta = 0.; double theta = 0.;
double phi = 0.; double phi = 0.;
...@@ -298,10 +292,10 @@ int main() { ...@@ -298,10 +292,10 @@ int main() {
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl; cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, Point pos(rootCS, 0_m, 0_m,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, stack.AddParticle(
corsika::stack::MomentumVector, geometry::Point, std::tuple<particles::Code, units::si::HEPEnergyType,
units::si::TimeType, unsigned short, unsigned short>{ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ}); beamCode, E0, plab, pos, 0_ns});
} }
// define air shower object, run simulation // define air shower object, run simulation
......
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