IAP GITLAB

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

use ShowerAxis in EnergyLoss

parent baa8946a
No related branches found
No related tags found
1 merge request!206Shower axis
Hide whitespace changes
Inline Side-by-side
Processes/EnergyLoss/EnergyLoss.cc
+ 211
− 222
View file @ 781f4d46
...@@ -21,6 +21,7 @@ ...@@ -21,6 +21,7 @@
#include <fstream> #include <fstream>
#include <iostream> #include <iostream>
#include <limits> #include <limits>
#include <numeric>
using namespace std; using namespace std;
...@@ -29,248 +30,236 @@ using namespace corsika::units::si; ...@@ -29,248 +30,236 @@ using namespace corsika::units::si;
using SetupParticle = corsika::setup::Stack::ParticleType; using SetupParticle = corsika::setup::Stack::ParticleType;
using SetupTrack = corsika::setup::Trajectory; using SetupTrack = corsika::setup::Trajectory;
namespace corsika::process::energy_loss { using namespace corsika::process::energy_loss;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { EnergyLoss::EnergyLoss(environment::ShowerAxis const& shower_axis)
return sqrt((Elab - m) * (Elab + m)); : shower_axis_(shower_axis)
}; , profile_(int(shower_axis.maximumX() / dX_) + 1) {}
/**
* PDG2018, passage of particles through matter
*
* Note, that \f$I_{\mathrm{eff}}\f$ of composite media a determined from \f$ \ln I =
* \sum_i a_i \ln(I_i) \f$ where \f$ a_i \f$ is the fraction of the electron population
* (\f$\sim Z_i\f$) of the \f$i\f$-th element. This can also be used for shell
* corrections or density effects.
*
* The \f$I_{\mathrm{eff}}\f$ of compounds is not better than a few percent, if not
* measured explicitly.
*
* For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
*
*/
HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const dX) {
// all these are material constants and have to come through Environment
// right now: values for nitrogen_D
// 7 nitrogen_gas 82.0 0.49976 D E 0.0011653 0.0 1.7378 4.1323 0.15349 3.2125 10.54
auto Ieff = 82.0_eV;
[[maybe_unused]] auto Zmat = 7;
auto ZoverA = 0.49976_mol / 1_g;
const double x0 = 1.7378;
const double x1 = 4.1323;
const double Cbar = 10.54;
const double delta0 = 0.0;
const double aa = 0.15349;
const double sk = 3.2125;
// end of material constants
// this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
auto constexpr K = 0.307075_MeV / 1_mol * square(1_cm);
HEPEnergyType const E = p.GetEnergy();
HEPMassType const m = p.GetMass();
double const gamma = E / m;
int const Z = p.GetChargeNumber();
int const Z2 = Z * Z;
HEPMassType constexpr me = particles::Electron::GetMass();
auto const m2 = m * m;
auto constexpr me2 = me * me;
double const gamma2 = gamma * gamma;
double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2 (1-1/gamma)*(1+1/gamma);
// (gamma_2-1)/gamma_2 = (1-1/gamma2);
double constexpr c2 = 1; // HEP convention here c=c2=1
cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
[[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
HEPMassType const Wmax =
2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
// approx, but <<1% HEPMassType const Wmax = 2*me*c2*beta2*gamma2; for HEAVY
// PARTICLES Wmax ~ 2me v2 for non-relativistic particles
cout << "BetheBloch Wmax=" << Wmax << endl;
// Sternheimer parameterization, density corrections towards 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) {
delta = 2 * (log(10)) * x - Cbar;
} else if (x < x1 && x >= x0) {
delta = 2 * (log(10)) * x - Cbar + aa * pow((x1 - x), sk);
} 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, <~100MeV
// need more clarity about formulas and units
const double Cadj = 0;
/*
// https://www.nap.edu/read/20066/chapter/8#104
HEPEnergyType Iadj = 12_eV * Z + 7_eV; // Iadj<163eV
if (Iadj>=163_eV)
Iadj = 9.76_eV * Z + 58.8_eV * pow(Z, -0.19); // Iadj>=163eV
double const Cadj = (0.422377/eta2 + 0.0304043/(eta2*eta2) -
0.00038106/(eta2*eta2*eta2)) * 1e-6 * Iadj*Iadj + (3.858019/eta2 -
0.1667989/(eta2*eta2) + 0.00157955/(eta2*eta2*eta2)) * 1e-9 * Iadj*Iadj*Iadj;
*/
// Barkas correction O(Z3) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
// double A = 1;
// if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
// double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
// double const Llow = 0.01 * Erel;
// double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
// double const barkas = Z * Llow*Lhigh/(Llow+Lhigh); // RU, I think the Z was
// missing...
double const barkas = 1; // does not work yet
// Bloch correction for O(Z4) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
const double alpha = 1. / 137.035999173;
double const y2 = Z * Z * alpha * alpha / beta2;
double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
// cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch <<
// endl;
double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
return -K * Z2 * ZoverA / beta2 *
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
// radiation losses according to PDG 2018, ch. 33 ref. [5]
HEPEnergyType EnergyLoss::RadiationLosses(SetupParticle const& vP,
GrammageType const vDX) {
// simple-minded hard-coded value for b(E) inspired by data from
// http://pdg.lbl.gov/2018/AtomicNuclearProperties/ for N and O.
auto constexpr b = 3.0 * 1e-6 * square(1_cm) / 1_g;
return -vP.GetEnergy() * b * vDX;
}
HEPEnergyType EnergyLoss::TotalEnergyLoss(SetupParticle const& vP, auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
GrammageType const vDX) { return sqrt((Elab - m) * (Elab + m));
return BetheBloch(vP, vDX) + RadiationLosses(vP, vDX); };
}
process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, /**
SetupTrack const& t) { * PDG2018, passage of particles through matter
if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk; *
* Note, that \f$I_{\mathrm{eff}}\f$ of composite media a determined from \f$ \ln I =
GrammageType const dX = * \sum_i a_i \ln(I_i) \f$ where \f$ a_i \f$ is the fraction of the electron population
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength()); * (\f$\sim Z_i\f$) of the \f$i\f$-th element. This can also be used for shell
cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber() * corrections or density effects.
<< ", dX=" << dX / 1_g * square(1_cm) << "g/cm2" << endl; *
HEPEnergyType dE = TotalEnergyLoss(p, dX); * The \f$I_{\mathrm{eff}}\f$ of compounds is not better than a few percent, if not
auto E = p.GetEnergy(); * measured explicitly.
const auto Ekin = E - p.GetMass(); *
auto Enew = E + dE; * For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
cout << "EnergyLoss dE=" << dE / 1_MeV << "MeV, " *
<< " E=" << E / 1_GeV << "GeV, Ekin=" << Ekin / 1_GeV */
<< ", Enew=" << Enew / 1_GeV << "GeV" << endl; HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const dX) {
auto status = process::EProcessReturn::eOk;
if (-dE > Ekin) { // all these are material constants and have to come through Environment
dE = -Ekin; // right now: values for nitrogen_D
Enew = p.GetMass(); // 7 nitrogen_gas 82.0 0.49976 D E 0.0011653 0.0 1.7378 4.1323 0.15349 3.2125 10.54
status = process::EProcessReturn::eParticleAbsorbed; auto Ieff = 82.0_eV;
} [[maybe_unused]] auto Zmat = 7;
p.SetEnergy(Enew); auto ZoverA = 0.49976_mol / 1_g;
MomentumUpdate(p, Enew); const double x0 = 1.7378;
EnergyLossTot_ += dE; const double x1 = 4.1323;
FillProfile(p, t, dE); const double Cbar = 10.54;
return status; const double delta0 = 0.0;
const double aa = 0.15349;
const double sk = 3.2125;
// end of material constants
// this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
auto constexpr K = 0.307075_MeV / 1_mol * square(1_cm);
HEPEnergyType const E = p.GetEnergy();
HEPMassType const m = p.GetMass();
double const gamma = E / m;
int const Z = p.GetChargeNumber();
int const Z2 = Z * Z;
HEPMassType constexpr me = particles::Electron::GetMass();
auto const m2 = m * m;
auto constexpr me2 = me * me;
double const gamma2 = gamma * gamma;
double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2 (1-1/gamma)*(1+1/gamma);
// (gamma_2-1)/gamma_2 = (1-1/gamma2);
double constexpr c2 = 1; // HEP convention here c=c2=1
cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
[[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
HEPMassType const Wmax =
2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
// approx, but <<1% HEPMassType const Wmax = 2*me*c2*beta2*gamma2; for HEAVY
// PARTICLES Wmax ~ 2me v2 for non-relativistic particles
cout << "BetheBloch Wmax=" << Wmax << endl;
// Sternheimer parameterization, density corrections towards 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) {
delta = 2 * (log(10)) * x - Cbar;
} else if (x < x1 && x >= x0) {
delta = 2 * (log(10)) * x - Cbar + aa * pow((x1 - x), sk);
} else if (x < x0) { // and IF conductor (otherwise, this is 0)
delta = delta0 * pow(100, 2 * (x - x0));
} }
cout << "BetheBloch delta=" << delta << endl;
LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
SetupTrack const& vTrack) const { // with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
if (vParticle.GetChargeNumber() == 0) {
return units::si::meter * std::numeric_limits<double>::infinity(); // shell correction, <~100MeV
} // need more clarity about formulas and units
const double Cadj = 0;
auto constexpr dX = 1_g / square(1_cm); /*
auto const dE = -TotalEnergyLoss(vParticle, dX); // dE > 0 // https://www.nap.edu/read/20066/chapter/8#104
//~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass(); HEPEnergyType Iadj = 12_eV * Z + 7_eV; // Iadj<163eV
auto const maxLoss = 0.01 * vParticle.GetEnergy(); if (Iadj>=163_eV)
auto const maxGrammage = maxLoss / dE * dX; Iadj = 9.76_eV * Z + 58.8_eV * pow(Z, -0.19); // Iadj>=163eV
double const Cadj = (0.422377/eta2 + 0.0304043/(eta2*eta2) -
return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack, 0.00038106/(eta2*eta2*eta2)) * 1e-6 * Iadj*Iadj + (3.858019/eta2 -
maxGrammage) * 0.1667989/(eta2*eta2) + 0.00157955/(eta2*eta2*eta2)) * 1e-9 * Iadj*Iadj*Iadj;
1.0001; // to make sure particle gets absorbed when DoContinuous() is called */
// Barkas correction O(Z3) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
// double A = 1;
// if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
// double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
// double const Llow = 0.01 * Erel;
// double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
// double const barkas = Z * Llow*Lhigh/(Llow+Lhigh); // RU, I think the Z was
// missing...
double const barkas = 1; // does not work yet
// Bloch correction for O(Z4) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
const double alpha = 1. / 137.035999173;
double const y2 = Z * Z * alpha * alpha / beta2;
double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
// cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch <<
// endl;
double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
return -K * Z2 * ZoverA / beta2 *
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
// radiation losses according to PDG 2018, ch. 33 ref. [5]
HEPEnergyType EnergyLoss::RadiationLosses(SetupParticle const& vP,
GrammageType const vDX) {
// simple-minded hard-coded value for b(E) inspired by data from
// http://pdg.lbl.gov/2018/AtomicNuclearProperties/ for N and O.
auto constexpr b = 3.0 * 1e-6 * square(1_cm) / 1_g;
return -vP.GetEnergy() * b * vDX;
}
HEPEnergyType EnergyLoss::TotalEnergyLoss(SetupParticle const& vP,
GrammageType const vDX) {
return BetheBloch(vP, vDX) + RadiationLosses(vP, vDX);
}
process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack const& t) {
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 = TotalEnergyLoss(p, dX);
auto E = p.GetEnergy();
const auto Ekin = E - p.GetMass();
auto Enew = E + dE;
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.GetMass();
status = process::EProcessReturn::eParticleAbsorbed;
} }
p.SetEnergy(Enew);
void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP, MomentumUpdate(p, Enew);
corsika::units::si::HEPEnergyType Enew) { FillProfile(t, dE);
HEPMomentumType Pnew = elab2plab(Enew, vP.GetMass()); return status;
auto pnew = vP.GetMomentum(); }
vP.SetMomentum(pnew * Pnew / pnew.GetNorm());
LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
SetupTrack const& vTrack) const {
if (vParticle.GetChargeNumber() == 0) {
return units::si::meter * std::numeric_limits<double>::infinity();
} }
void EnergyLoss::FillProfile(SetupParticle const& vP, SetupTrack const& vTrack, auto constexpr dX = 1_g / square(1_cm);
const HEPEnergyType dE) { auto const dE = -TotalEnergyLoss(vParticle, dX); // dE > 0
//~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass();
auto const maxLoss = 0.01 * vParticle.GetEnergy();
auto const maxGrammage = maxLoss / dE * dX;
using namespace corsika::geometry; return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack,
maxGrammage) *
1.0001; // to make sure particle gets absorbed when DoContinuous() is called
}
auto const toStart = vTrack.GetPosition(0) - InjectionPoint_; void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,
auto const toEnd = vTrack.GetPosition(1) - InjectionPoint_; corsika::units::si::HEPEnergyType Enew) {
HEPMomentumType Pnew = elab2plab(Enew, vP.GetMass());
auto pnew = vP.GetMomentum();
vP.SetMomentum(pnew * Pnew / pnew.GetNorm());
}
auto const v1 = (toStart * 1_Hz).dot(ShowerAxisDirection_); void EnergyLoss::FillProfile(SetupTrack const& vTrack, const HEPEnergyType dE) {
auto const v2 = (toEnd * 1_Hz).dot(ShowerAxisDirection_);
geometry::Line const lineToStartBin(InjectionPoint_, ShowerAxisDirection_ * v1);
geometry::Line const lineToEndBin(InjectionPoint_, ShowerAxisDirection_ * v2);
SetupTrack const trajToStartBin(lineToStartBin, 1_s); GrammageType const grammageStart = shower_axis_.projectedX(vTrack.GetPosition(0));
SetupTrack const trajToEndBin(lineToEndBin, 1_s); GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1));
const auto deltaX = grammageEnd - grammageStart;
GrammageType const grammageStart = const int binStart = grammageStart / dX_;
vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToStartBin, const int binEnd = grammageEnd / dX_;
trajToStartBin.GetLength());
GrammageType const grammageEnd =
vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToEndBin,
trajToEndBin.GetLength());
const int binStart = grammageStart / dX_; std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
const int binEnd = grammageEnd / dX_; << grammageEnd << std::endl;
std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and " auto energyCount = HEPEnergyType::zero();
<< grammageEnd << std::endl;
auto energyCount = HEPEnergyType::zero(); auto fill = [&](int bin, GrammageType weight) {
if (deltaX > dX_threshold_) {
auto const increment = -dE * weight / deltaX;
profile_[bin] += increment;
energyCount += increment;
auto fill = [&](int bin, GrammageType weight) { std::cout << "filling bin " << bin << " with weight " << weight << ": " << increment
const auto dX = grammageEnd - grammageStart; << std::endl;
if (dX > dX_threshold_) { }
auto const increment = -dE * weight / (grammageEnd - grammageStart); };
Profile_[bin] += increment;
energyCount += increment;
std::cout << "filling bin " << bin << " with weight " << weight << ": " // fill longitudinal profile
<< increment << std::endl; fill(binStart, (1 + binStart) * dX_ - grammageStart);
} fill(binEnd, grammageEnd - binEnd * dX_);
};
// fill longitudinal profile if (binStart == binEnd) { fill(binStart, -dX_); }
fill(binStart, (1 + binStart) * dX_ - grammageStart);
fill(binEnd, grammageEnd - binEnd * dX_);
if (binStart == binEnd) { fill(binStart, -dX_); } for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); } std::cout << "total energy added to histogram: " << energyCount << std::endl;
}
std::cout << "total energy added to histogram: " << energyCount << std::endl; void EnergyLoss::PrintProfile() const {
} std::ofstream file("EnergyLossProfile.dat");
file << "# EnergyLoss profile" << std::endl
void EnergyLoss::PrintProfile() const { << "# lower X bin edge [g/cm2] dE/dX [GeV/g/cm2]" << endl;
std::ofstream file("EnergyLossProfile.dat"); double const deltaX = dX_ / 1_g * square(1_cm);
cout << "# EnergyLoss PrintProfile X-bin [g/cm2] dE/dX [GeV/g/cm2] " << endl; for (size_t i = 0; i < profile_.size(); ++i) {
double const deltaX = dX_ / 1_g * square(1_cm); file << std::scientific << std::setw(15) << i * deltaX << std::setw(15)
for (auto v : Profile_) { << profile_.at(i) / (deltaX * 1_GeV) << endl;
file << v.first * deltaX << " " << v.second / (deltaX * 1_GeV) << endl;
}
} }
}
} // namespace corsika::process::energy_loss HEPEnergyType EnergyLoss::GetTotal() const {
return std::accumulate(profile_.cbegin(), profile_.cend(), HEPEnergyType::zero());
}
Processes/EnergyLoss/EnergyLoss.h
+ 8
− 19
View file @ 781f4d46
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
#ifndef _Processes_EnergyLoss_h_ #ifndef _Processes_EnergyLoss_h_
#define _Processes_EnergyLoss_h_ #define _Processes_EnergyLoss_h_
#include <corsika/environment/ShowerAxis.h>
#include <corsika/geometry/Point.h> #include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h> #include <corsika/geometry/Vector.h>
#include <corsika/process/ContinuousProcess.h> #include <corsika/process/ContinuousProcess.h>
...@@ -31,21 +32,14 @@ namespace corsika::process::energy_loss { ...@@ -31,21 +32,14 @@ namespace corsika::process::energy_loss {
void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew); void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew);
public: public:
template <typename TDim> EnergyLoss(environment::ShowerAxis const& showerAxis);
EnergyLoss(geometry::Point const& injectionPoint,
geometry::Vector<TDim> const& direction)
: InjectionPoint_(injectionPoint)
, ShowerAxisDirection_(direction.normalized()) {}
EnergyLoss(setup::Trajectory const& trajectory)
: EnergyLoss(trajectory.GetPosition(0), trajectory.GetV0()){};
void Init() {} void Init() {}
process::EProcessReturn DoContinuous(setup::Stack::ParticleType&, process::EProcessReturn DoContinuous(setup::Stack::ParticleType&,
setup::Trajectory const&); setup::Trajectory const&);
units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&, units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const&) const; setup::Trajectory const&) const;
units::si::HEPEnergyType GetTotal() const { return EnergyLossTot_; } units::si::HEPEnergyType GetTotal() const;
void PrintProfile() const; void PrintProfile() const;
static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const&, static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const&,
const units::si::GrammageType); const units::si::GrammageType);
...@@ -55,22 +49,17 @@ namespace corsika::process::energy_loss { ...@@ -55,22 +49,17 @@ namespace corsika::process::energy_loss {
const units::si::GrammageType); const units::si::GrammageType);
private: private:
void FillProfile(setup::Stack::ParticleType const&, setup::Trajectory const&, void FillProfile(setup::Trajectory const&, units::si::HEPEnergyType);
units::si::HEPEnergyType);
// void FillProfileAbsorbed(setup::Stack::ParticleType const&, setup::Trajectory
// const&);
units::si::HEPEnergyType EnergyLossTot_ = units::si::HEPEnergyType::zero();
units::si::GrammageType const dX_ = std::invoke([]() { units::si::GrammageType const dX_ = std::invoke([]() {
using namespace units::si; using namespace units::si;
return 10_g / square(1_cm); return 10_g / square(1_cm);
}); // profile binning }); // profile binning
std::map<int, units::si::HEPEnergyType> Profile_; // longitudinal profile environment::ShowerAxis const& shower_axis_;
geometry::Point const InjectionPoint_; std::vector<units::si::HEPEnergyType> profile_; // longitudinal profile
geometry::Vector<units::si::dimensionless_d> const ShowerAxisDirection_;
}; };
const units::si::GrammageType dX_threshold_ = std::invoke([]() { units::si::GrammageType const dX_threshold_ = std::invoke([]() {
using namespace units::si; using namespace units::si;
return 0.0001_g / square(1_cm); return 0.0001_g / square(1_cm);
}); });
......
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