diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 621cab01d1facaa687c12f287c12ed2d54a6501f..f8f1152a6c195182936dd0c8b5745d1438b65dde 100644 --- a/Documentation/Examples/CMakeLists.txt +++ b/Documentation/Examples/CMakeLists.txt @@ -24,14 +24,16 @@ CORSIKA_ADD_TEST (stack_example) add_executable (cascade_example cascade_example.cc) target_compile_options(cascade_example PRIVATE -g) # do not skip asserts target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging - CORSIKArandom - ProcessSibyll - CORSIKAcascade + CORSIKArandom +# CORSIKAprocesses + ProcessSibyll + ProcessEnergyLoss ProcessStackInspector ProcessTrackWriter ProcessTrackingLine ProcessHadronicElasticModel CORSIKAprocesses + CORSIKAcascade CORSIKAparticles CORSIKAgeometry CORSIKAenvironment diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index f9d1e95a83aa38ed0475c06850698d2b867f6f99..a875053845c0c81df2e454a14655a18e052ce64c 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -11,6 +11,7 @@ #include <corsika/cascade/Cascade.h> #include <corsika/process/ProcessSequence.h> +#include <corsika/process/energy_loss/EnergyLoss.h> #include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h> #include <corsika/process/stack_inspector/StackInspector.h> #include <corsika/process/tracking_line/TrackingLine.h> @@ -262,10 +263,12 @@ int main() { // hadronicElastic(env); process::TrackWriter::TrackWriter trackWriter("tracks.dat"); + process::EnergyLoss::EnergyLoss eLoss(2_MeV / 1_g * square(1_cm)); // assemble all processes into an ordered process list // auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter; - auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter; + auto sequence = p0 << sibyll << sibyllNuc << decay << eLoss << cut << trackWriter; + // auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter; // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() // << "\n"; @@ -276,17 +279,16 @@ int main() { const Code beamCode = Code::Nucleus; const int nuclA = 56; const int nuclZ = int(nuclA / 2.15 + 0.7); - const HEPMassType mass = particles::Proton::GetMass() * nuclZ + - (nuclA - nuclZ) * particles::Neutron::GetMass(); + const HEPMassType mass = GetNucleusMass(nuclA, nuclZ); const HEPEnergyType E0 = nuclA * - 100_GeV; // 1_PeV crashes with bad COMboost in second interaction (crash later) + 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later) double theta = 0.; double phi = 0.; { auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt(Elab * Elab - m * m); + return sqrt((Elab - m) * (Elab + m)); }; HEPMomentumType P0 = elab2plab(E0, mass); auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) { @@ -315,6 +317,9 @@ int main() { cut.ShowResults(); const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy(); - cout << "total energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl; + cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl + << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; + cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl + << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl; + eLoss.SaveSave(); } diff --git a/Environment/DensityFunction.h b/Environment/DensityFunction.h index 32ee30adedfcbe9fccf1ebea2f09c8b9877f6015..92ba841be89e63d7edc49f6eec5cf5bafaa6fae0 100644 --- a/Environment/DensityFunction.h +++ b/Environment/DensityFunction.h @@ -1,4 +1,5 @@ -/** + +/* * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * * See file AUTHORS for a list of contributors. diff --git a/Environment/InhomogeneousMedium.h b/Environment/InhomogeneousMedium.h index feb6497cfcd82a28b8eafa183c068d106e187520..0b1a0c1b9866ce1a3c64547c06004e6797cb41c7 100644 --- a/Environment/InhomogeneousMedium.h +++ b/Environment/InhomogeneousMedium.h @@ -1,4 +1,5 @@ -/** + +/* * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * * See file AUTHORS for a list of contributors. diff --git a/Environment/LinearApproximationIntegrator.h b/Environment/LinearApproximationIntegrator.h index 910d2397189e7aa9b1456e00df94f72cfd49ff4f..90e14d99a472b9933d094087baefef802ab9d7a9 100644 --- a/Environment/LinearApproximationIntegrator.h +++ b/Environment/LinearApproximationIntegrator.h @@ -1,4 +1,5 @@ -/** + +/* * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * * See file AUTHORS for a list of contributors. diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h index 32063061058c446430cd1ba52e13924c20685eff..e8a1a58bf0d7cdaa8a6fab5836a592c674ad3a7d 100644 --- a/Framework/Geometry/Trajectory.h +++ b/Framework/Geometry/Trajectory.h @@ -37,20 +37,21 @@ namespace corsika::geometry { Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); } corsika::units::si::TimeType GetDuration() const { return fTimeLength; } + corsika::units::si::LengthType GetLength() const { return GetDistance(fTimeLength); } corsika::units::si::LengthType GetDistance(corsika::units::si::TimeType t) const { - assert(t > fTimeLength); + assert(t <= fTimeLength); assert(t >= 0 * corsika::units::si::second); - return T::ArcLength(0, t); + return T::ArcLength(0 * corsika::units::si::second, t); } void LimitEndTo(corsika::units::si::LengthType limit) { fTimeLength = T::TimeFromArclength(limit); } - + auto NormalizedDirection() const { - static_assert(std::is_same_v<T, corsika::geometry::Line>); - return T::GetV0().normalized(); + static_assert(std::is_same_v<T, corsika::geometry::Line>); + return T::GetV0().normalized(); } }; diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h index 5ac1f118b789b500ad26a23f840bda954b2c9800..737b7830821df05cdf0fa5fd93918741f0cf731e 100644 --- a/Framework/Geometry/Vector.h +++ b/Framework/Geometry/Vector.h @@ -75,12 +75,14 @@ namespace corsika::geometry { * think about whether squaredNorm() might be cheaper for your computation. */ auto norm() const { return BaseVector<dim>::qVector.norm(); } + auto GetNorm() const { return BaseVector<dim>::qVector.norm(); } /*! * returns the squared norm of the Vector. Before using this method, * think about whether norm() might be cheaper for your computation. */ auto squaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); } + auto GetSquaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); } /*! * returns a Vector \f$ \vec{v}_{\parallel} \f$ which is the parallel projection diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc index d73ae99b2f21590984aae1e90c683777549b7b32..dae272be1bc052b3174ad42ae0b2a3f50777322b 100644 --- a/Framework/Geometry/testGeometry.cc +++ b/Framework/Geometry/testGeometry.cc @@ -177,8 +177,10 @@ TEST_CASE("Trajectories") { CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates()); CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3)); - - CHECK((base.NormalizedDirection().GetComponents(rootCS) - QuantityVector<dimensionless_d>{1, 0, 0}).eVector.norm() == Approx(0).margin(absMargin)); + + CHECK((base.NormalizedDirection().GetComponents(rootCS) - + QuantityVector<dimensionless_d>{1, 0, 0}) + .eVector.norm() == Approx(0).margin(absMargin)); } SECTION("Helix") { diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h index dcbbe6398e84a0ad3f148eabf36ad01f380deb12..456711b4983766b21428766e70b1d25805212883 100644 --- a/Framework/Particles/ParticleProperties.h +++ b/Framework/Particles/ParticleProperties.h @@ -44,8 +44,8 @@ namespace corsika::particles { using PDGCodeType = std::underlying_type<PDGCode>::type; // forward declarations to be used in GeneratedParticleProperties - int16_t constexpr GetElectricChargeNumber(Code const); - corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const); + int16_t constexpr GetChargeNumber(Code const); + corsika::units::si::ElectricChargeType constexpr GetCharge(Code const); corsika::units::si::HEPMassType constexpr GetMass(Code const); PDGCode constexpr GetPDG(Code const); constexpr std::string const& GetName(Code const); @@ -72,17 +72,18 @@ namespace corsika::particles { } /*! - * returns electric charge of particle / (e/3), e.g. return 3 for a proton. + * returns electric charge number of particle return 1 for a proton. */ - int16_t constexpr GetElectricChargeNumber(Code const p) { - return detail::electric_charges[static_cast<CodeIntType>(p)]; + int16_t constexpr GetChargeNumber(Code const p) { + // electric_charges stores charges in units of (e/3), e.g. 3 for a proton + return detail::electric_charges[static_cast<CodeIntType>(p)] / 3; } /*! * returns electric charge of particle, e.g. return 1.602e-19_C for a proton. */ - corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) { - return GetElectricChargeNumber(p) * (corsika::units::constants::e * (1. / 3.)); + corsika::units::si::ElectricChargeType constexpr GetCharge(Code const p) { + return GetChargeNumber(p) * (corsika::units::constants::e); } constexpr std::string const& GetName(Code const p) { @@ -112,6 +113,14 @@ namespace corsika::particles { std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p); Code ConvertFromPDG(PDGCode); + + /** + * Get mass of nucleus + **/ + corsika::units::si::HEPMassType constexpr GetNucleusMass(const int vA, const int vZ) { + return Proton::GetMass() * vZ + (vA - vZ) * Neutron::GetMass(); + } + } // namespace corsika::particles #endif diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py index 946940f7dcf54551e2888e438c6eeeae0e1659e1..84ef0ee981a77d1ac68768ac548d02ec06ac19f3 100755 --- a/Framework/Particles/pdxml_reader.py +++ b/Framework/Particles/pdxml_reader.py @@ -418,8 +418,8 @@ def gen_classes(particle_db): string += " public:\n" string += " static constexpr Code GetCode() { return Type; }\n" string += " static constexpr corsika::units::si::HEPMassType GetMass() { return corsika::particles::GetMass(Type); }\n" - string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetElectricCharge(Type); }\n" - string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n" + string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetCharge(Type); }\n" + string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetChargeNumber(Type); }\n" string += " static std::string const& GetName() { return corsika::particles::GetName(Type); }\n" string += " static constexpr Code GetAntiParticle() { return AntiType; }\n" string += " static constexpr bool IsNucleus() { return corsika::particles::IsNucleus(Type); }\n" diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc index 8bb73e7ac252fafe334d12c93a22dc05e195c4a5..3fd99ea2c0e24f52ebd61cd978db97ce54fe3ac7 100644 --- a/Framework/Particles/testParticles.cc +++ b/Framework/Particles/testParticles.cc @@ -43,7 +43,7 @@ TEST_CASE("ParticleProperties", "[Particles]") { SECTION("Charges") { REQUIRE(Electron::GetCharge() / constants::e == Approx(-1)); REQUIRE(Positron::GetCharge() / constants::e == Approx(+1)); - REQUIRE(GetElectricCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1)); + REQUIRE(GetCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1)); } SECTION("Names") { diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h index 786044aef3055b1472e0150dd451fa5050628f4b..6bc71d708316fdfb22656b22b21b462104e50cd0 100644 --- a/Framework/Units/PhysicalUnits.h +++ b/Framework/Units/PhysicalUnits.h @@ -45,9 +45,11 @@ namespace corsika::units::si { using hepmomentum_d = phys::units::hepenergy_d; using hepmass_d = phys::units::hepenergy_d; - /// defining cross section + /// defining cross section as area using sigma_d = phys::units::area_d; + // constexpr quantity<area_d> barn{Rep(1e-28L) * square(meter)}; + /// add the unit-types using LengthType = phys::units::quantity<phys::units::length_d, double>; using TimeType = phys::units::quantity<phys::units::time_interval_d, double>; diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index 26a44fe6e1ea2ddf806a4f9fd1e2b0552153a64d..e23f660bc65bb9a86612d732f2ebb52429ba1b99 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -5,6 +5,7 @@ add_subdirectory (StackInspector) add_subdirectory (TrackWriter) add_subdirectory (TrackingLine) add_subdirectory (HadronicElasticModel) +add_subdirectory (EnergyLoss) #add_custom_target(CORSIKAprocesses) add_library (CORSIKAprocesses INTERFACE) @@ -12,3 +13,4 @@ add_dependencies(CORSIKAprocesses ProcessNullModel) add_dependencies(CORSIKAprocesses ProcessSibyll) add_dependencies(CORSIKAprocesses ProcessStackInspector) add_dependencies(CORSIKAprocesses ProcessTrackingLine) +add_dependencies(CORSIKAprocesses ProcessEnergyLoss) diff --git a/Processes/EnergyLoss/CMakeLists.txt b/Processes/EnergyLoss/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..f108ba859df8cf1addf5f9690879dc0cfddfd7ba --- /dev/null +++ b/Processes/EnergyLoss/CMakeLists.txt @@ -0,0 +1,64 @@ + +set ( + MODEL_SOURCES + EnergyLoss.cc + ) + +set ( + MODEL_HEADERS + EnergyLoss.h + ) + +set ( + MODEL_NAMESPACE + corsika/process/energy_loss + ) + +add_library (ProcessEnergyLoss STATIC ${MODEL_SOURCES}) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessEnergyLoss ${MODEL_NAMESPACE} ${MODEL_HEADERS}) + +set_target_properties ( + ProcessEnergyLoss + PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION 1 +# PUBLIC_HEADER "${MODEL_HEADERS}" + ) + +# target dependencies on other libraries (also the header onlys) +target_link_libraries ( + ProcessEnergyLoss + CORSIKAunits + CORSIKAparticles + CORSIKAgeometry + CORSIKAsetup + ) + +target_include_directories ( + ProcessEnergyLoss + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) + +install ( + TARGETS ProcessEnergyLoss + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib +# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE} + ) + + +# -------------------- +# code unit testing +#add_executable (testNullModel testNullModel.cc) + +#target_link_libraries ( +# testNullModel ProcessNullModel +# CORSIKAsetup +# CORSIKAgeometry +# CORSIKAunits +# CORSIKAthirdparty # for catch2 +# ) +# CORSIKA_ADD_TEST(testNullModel) + diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc new file mode 100644 index 0000000000000000000000000000000000000000..8646b4d56dbfa6a351d22c4618ab7adf16c0e1a9 --- /dev/null +++ b/Processes/EnergyLoss/EnergyLoss.cc @@ -0,0 +1,108 @@ + +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <corsika/process/energy_loss/EnergyLoss.h> + +#include <corsika/particles/ParticleProperties.h> + +#include <corsika/setup/SetupStack.h> +#include <corsika/setup/SetupTrajectory.h> + +#include <cmath> +#include <iostream> +#include <limits> + +using namespace std; + +using namespace corsika; +using namespace corsika::units::si; +using namespace corsika::setup; +using Particle = Stack::ParticleType; +using Track = Trajectory; + +namespace corsika::process::EnergyLoss { + + auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { + return sqrt((Elab - m) * (Elab + m)); + }; + + EnergyLoss::EnergyLoss(MeVgcm2 const vdEdX) + : fdEdX(vdEdX) + , fEnergyLossTot(0_GeV) {} + + process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) { + GrammageType const dX = + p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength()); + HEPEnergyType dE = -dX * fdEdX * pow(p.GetChargeNumber(), 2); + auto E = p.GetEnergy(); + const auto Ekin = E - p.GetParticleMass(); + auto Enew = E + dE; + cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber() + << ", dX=" << dX / 1_g * square(1_cm) << "g/cm2, dE=" << dE / 1_MeV << "MeV, " + << " E=" << E / 1_GeV << "GeV, Ekin=" << Ekin / 1_GeV + << ", Enew=" << Enew / 1_GeV << "GeV" << endl; + auto status = process::EProcessReturn::eOk; + if (-dE > Ekin) { + dE = -Ekin; + Enew = p.GetParticleMass(); + status = process::EProcessReturn::eParticleAbsorbed; + } + p.SetEnergy(Enew); + MomentumUpdate(p, Enew); + fEnergyLossTot += dE; + GetXbin(p, dE); + return status; + } + + units::si::LengthType EnergyLoss::MaxStepLength(Particle&, Track&) { + return units::si::meter * std::numeric_limits<double>::infinity(); + } + + void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& p, + corsika::units::si::HEPEnergyType Enew) { + HEPMomentumType Pnew = elab2plab(Enew, p.GetParticleMass()); + auto pnew = p.GetMomentum(); + p.SetMomentum(pnew * Pnew / pnew.GetNorm()); + } + +#include <corsika/geometry/CoordinateSystem.h> + + int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& p, + const HEPEnergyType dE) { + + using namespace corsika::geometry; + + const GrammageType deltaX = 10_g / square(1_cm); // binning + + CoordinateSystem const& rootCS = + RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); + Point pos1(rootCS, 0_m, 0_m, 0_m); + Point pos2(rootCS, 0_m, 0_m, p.GetPosition().GetCoordinates()[2]); + Vector delta = (pos2 - pos1) / 1_s; + Trajectory t(Line(pos1, delta), 1_s); + + GrammageType const grammage = + p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength()); + + const int bin = grammage / deltaX; + + if (!fSave.count(bin)) { cout << "EnergyLoss new x bin " << bin << endl; } + fSave[bin] += -dE / 1_GeV; + return bin; + } + + void EnergyLoss::SaveSave() { + + cout << "EnergyLoss Save " << endl; + for (auto v : fSave) { cout << v.first << " " << v.second << endl; } + } + +} // namespace corsika::process::EnergyLoss diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h new file mode 100644 index 0000000000000000000000000000000000000000..3af10c01748de1e66eba1251d3c859d3703fce1d --- /dev/null +++ b/Processes/EnergyLoss/EnergyLoss.h @@ -0,0 +1,57 @@ + +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#ifndef _Processes_EnergyLoss_h_ +#define _Processes_EnergyLoss_h_ + +#include <corsika/process/ContinuousProcess.h> +#include <corsika/units/PhysicalUnits.h> + +#include <corsika/setup/SetupStack.h> +#include <corsika/setup/SetupTrajectory.h> + +#include <map> + +namespace corsika::process::EnergyLoss { + + class EnergyLoss : public corsika::process::ContinuousProcess<EnergyLoss> { + + using MeVgcm2 = decltype(1e6 * units::si::electronvolt / units::si::gram * + corsika::units::si::square(1e-2 * units::si::meter)); + + void MomentumUpdate(corsika::setup::Stack::ParticleType&, + corsika::units::si::HEPEnergyType Enew); + + public: + EnergyLoss(MeVgcm2 const vdEdX); + void Init() {} + + corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&, + corsika::setup::Trajectory&, + corsika::setup::Stack&); + corsika::units::si::LengthType MaxStepLength(corsika::setup::Stack::ParticleType&, + corsika::setup::Trajectory&); + + corsika::units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; } + void SaveSave(); + + private: + int GetXbin(corsika::setup::Stack::ParticleType& p, + const corsika::units::si::HEPEnergyType dE); + + MeVgcm2 fdEdX; + corsika::units::si::HEPEnergyType fEnergyLossTot; + std::map<int, double> fSave; + }; + +} // namespace corsika::process::EnergyLoss + +#endif diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h index f03963872164cdd4d48b9c09ccc3053348f93ea8..7f37905de1672f740f1ea7ec44051c5a730d1037 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.h +++ b/Processes/HadronicElasticModel/HadronicElasticModel.h @@ -44,8 +44,8 @@ namespace corsika::process::HadronicElasticModel { using SquaredHEPEnergyType = decltype(corsika::units::si::HEPEnergyType() * corsika::units::si::HEPEnergyType()); - using eV2 = decltype(units::si::detail::static_pow<2>(units::si::electronvolt)); - using inveV2 = decltype(units::si::detail::static_pow<-2>(units::si::electronvolt)); + using eV2 = decltype(units::si::square(units::si::electronvolt)); + using inveV2 = decltype(1 / units::si::square(units::si::electronvolt)); corsika::random::RNG& fRNG = corsika::random::RNGManager::GetInstance().GetRandomStream( diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index 5ccc1b743ce35f1e72a1ebd1e5cccba95a8d25e9..3c27e4d9755180a45c70a42c1313a7540a1b5e6a 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/Processes/Sibyll/Interaction.cc @@ -277,8 +277,8 @@ namespace corsika::process::sibyll { sigEla; // to avoid not used warning in array binding } - const auto targetCode = mediumComposition.SampleTarget( - cross_section_of_components, fRNG); + const auto targetCode = + mediumComposition.SampleTarget(cross_section_of_components, fRNG); cout << "Interaction: target selected: " << targetCode << endl; /* FOR NOW: allow nuclei with A<18 or protons only. diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc index 6a37363235f66c455206e16c0134a9b623cf4ff4..575903e019dea219b4b94b39e5925f3243b76c58 100644 --- a/Processes/Sibyll/NuclearInteraction.cc +++ b/Processes/Sibyll/NuclearInteraction.cc @@ -356,7 +356,8 @@ namespace corsika::process::sibyll { [[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS } - const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, fRNG); + const auto targetCode = + mediumComposition.SampleTarget(cross_section_of_components, fRNG); cout << "Interaction: target selected: " << targetCode << endl; /* FOR NOW: allow nuclei with A<18 or protons only. diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h index 11923aba8497123b6ac9175b73ab6b1b27aa9725..d6d545238e9510d21ccc787f903b267e0c46c70d 100644 --- a/Stack/NuclearStackExtension/NuclearStackExtension.h +++ b/Stack/NuclearStackExtension/NuclearStackExtension.h @@ -159,7 +159,24 @@ namespace corsika::stack { int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); } /// @} - // int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); } + /** + * Overwrite normal GetParticleMass function with nuclear version + */ + corsika::units::si::HEPMassType GetParticleMass() const { + if (InnerParticleInterface<StackIteratorInterface>::GetPID() == + corsika::particles::Code::Nucleus) + return corsika::particles::GetNucleusMass(GetNuclearA(), GetNuclearZ()); + return InnerParticleInterface<StackIteratorInterface>::GetParticleMass(); + } + /** + * Overwirte normal GetChargeNumber function with nuclear version + **/ + int16_t GetChargeNumber() const { + if (InnerParticleInterface<StackIteratorInterface>::GetPID() == + corsika::particles::Code::Nucleus) + return GetNuclearZ(); + return InnerParticleInterface<StackIteratorInterface>::GetChargeNumber(); + } protected: void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); } diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h index 365a45aa963c235b34b5e762cd79192af40b56c2..f28b51322425d24de8e293a2c33c76d5fc87c995 100644 --- a/Stack/SuperStupidStack/SuperStupidStack.h +++ b/Stack/SuperStupidStack/SuperStupidStack.h @@ -117,10 +117,22 @@ namespace corsika::stack { corsika::units::si::TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); } + /** + * @name derived quantities + * + * @{ + */ corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection() const { return GetMomentum() / GetEnergy(); } + corsika::units::si::HEPMassType GetParticleMass() const { + return corsika::particles::GetMass(GetPID()); + } + int16_t GetChargeNumber() const { + return corsika::particles::GetChargeNumber(GetPID()); + } + ///@} }; /**