From 8753b89c759b040649748e96c3d4ea367fb6c5e7 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 19 Feb 2019 11:24:55 +0100
Subject: [PATCH] fixed zero mass bug for nuclei and photons
---
Framework/Particles/ParticleProperties.h | 6 ++++
Processes/EnergyLoss/EnergyLoss.cc | 28 ++++++++++++-------
Processes/Sibyll/Decay.cc | 2 +-
.../NuclearStackExtension.h | 4 +--
Stack/SuperStupidStack/SuperStupidStack.h | 2 +-
5 files changed, 28 insertions(+), 14 deletions(-)
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index 456711b49..891e11403 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -61,6 +61,8 @@ namespace corsika::particles {
* returns mass of particle in natural units
*/
corsika::units::si::HEPMassType constexpr GetMass(Code const p) {
+ if (p == Code::Nucleus)
+ throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
return detail::masses[static_cast<CodeIntType>(p)];
}
@@ -75,6 +77,8 @@ namespace corsika::particles {
* returns electric charge number of particle return 1 for a proton.
*/
int16_t constexpr GetChargeNumber(Code const p) {
+ if (p == Code::Nucleus)
+ throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
// electric_charges stores charges in units of (e/3), e.g. 3 for a proton
return detail::electric_charges[static_cast<CodeIntType>(p)] / 3;
}
@@ -83,6 +87,8 @@ namespace corsika::particles {
* returns electric charge of particle, e.g. return 1.602e-19_C for a proton.
*/
corsika::units::si::ElectricChargeType constexpr GetCharge(Code const p) {
+ if (p == Code::Nucleus)
+ throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
return GetChargeNumber(p) * (corsika::units::constants::e);
}
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 6164f8833..148f37feb 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -71,7 +71,7 @@ namespace corsika::process::EnergyLoss {
// this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
auto const K = 0.307075_MeV / 1_mol * square(1_cm);
HEPEnergyType const E = p.GetEnergy();
- HEPMassType const m = particles::GetMass(p.GetPID());
+ HEPMassType const m = p.GetMass();
double const gamma = E / m;
double const Z = p.GetChargeNumber();
double const Z2 = pow(Z, 2);
@@ -79,14 +79,18 @@ namespace corsika::process::EnergyLoss {
auto const m2 = m * m;
auto const me2 = me * me;
double const gamma2 = pow(gamma, 2);
+ cout << "gamma2=" << gamma2 << endl;
+
double const beta2 = (gamma2-1)/gamma2; // (1-1/gamma)*(1+1/gamma); (gamma_2-1)/gamma_2 = (1-1/gamma2);
double const c2 = 1; // HEP convention here c=c2=1
- double const eta2 = beta2/(1-beta2);
-
+ cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
+ double const eta2 = beta2/(1 - beta2);
HEPMassType const Wmax = 2*me*c2*beta2*gamma2 / (1 + 2*gamma*me/m + me2/m2);
-
+ cout << "BetheBloch Wmax=" << Wmax << endl;
+
// Sternheimer parameterization, high energies
// NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization -> MISSING
+ cout << "BetheBloch p.GetMomentum().GetNorm()/m=" << p.GetMomentum().GetNorm()/m << endl;
double const x = log10(p.GetMomentum().GetNorm()/m);
double delta = 0;
if (x>=x1) {
@@ -96,7 +100,8 @@ namespace corsika::process::EnergyLoss {
} else if (x<x0) { // AND IF CONDUCTOR (otherwise, this is 0)
delta = delta0*pow(100,2*(x-x0));
}
-
+ cout << "BetheBloch delta=" << delta << endl;
+
// with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
// shell correction
@@ -123,20 +128,23 @@ namespace corsika::process::EnergyLoss {
}
process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) {
+ if (p.GetChargeNumber()==0)
+ return process::EProcessReturn::eOk;
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
+ cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
+ << ", dX=" << dX / 1_g * square(1_cm) << "g/cm2" << endl;
HEPEnergyType dE = BetheBloch(p, dX);
auto E = p.GetEnergy();
- const auto Ekin = E - p.GetParticleMass();
+ const auto Ekin = E - p.GetMass();
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, "
+ cout << "EnergyLoss 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();
+ Enew = p.GetMass();
status = process::EProcessReturn::eParticleAbsorbed;
}
p.SetEnergy(Enew);
@@ -152,7 +160,7 @@ namespace corsika::process::EnergyLoss {
void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& p,
corsika::units::si::HEPEnergyType Enew) {
- HEPMomentumType Pnew = elab2plab(Enew, p.GetParticleMass());
+ HEPMomentumType Pnew = elab2plab(Enew, p.GetMass());
auto pnew = p.GetMomentum();
p.SetMomentum(pnew * Pnew / pnew.GetNorm());
}
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 470ea1819..6816a350b 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -108,7 +108,7 @@ namespace corsika::process::sibyll {
using namespace units::si;
HEPEnergyType E = p.GetEnergy();
- HEPMassType m = particles::GetMass(p.GetPID());
+ HEPMassType m = p.GetMass();
const double gamma = E / m;
diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h
index d6d545238..6c1fd2b80 100644
--- a/Stack/NuclearStackExtension/NuclearStackExtension.h
+++ b/Stack/NuclearStackExtension/NuclearStackExtension.h
@@ -162,11 +162,11 @@ namespace corsika::stack {
/**
* Overwrite normal GetParticleMass function with nuclear version
*/
- corsika::units::si::HEPMassType GetParticleMass() const {
+ corsika::units::si::HEPMassType GetMass() const {
if (InnerParticleInterface<StackIteratorInterface>::GetPID() ==
corsika::particles::Code::Nucleus)
return corsika::particles::GetNucleusMass(GetNuclearA(), GetNuclearZ());
- return InnerParticleInterface<StackIteratorInterface>::GetParticleMass();
+ return InnerParticleInterface<StackIteratorInterface>::GetMass();
}
/**
* Overwirte normal GetChargeNumber function with nuclear version
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index f28b51322..1610aede7 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -126,7 +126,7 @@ namespace corsika::stack {
const {
return GetMomentum() / GetEnergy();
}
- corsika::units::si::HEPMassType GetParticleMass() const {
+ corsika::units::si::HEPMassType GetMass() const {
return corsika::particles::GetMass(GetPID());
}
int16_t GetChargeNumber() const {
--
GitLab