IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 4064b5d5 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan: Committed by Maximilian Reininghaus
Browse files

fixed bug in EnergyLoss binning

parent 101d7274
No related branches found
No related tags found
1 merge request!116Some improvements here and there
......@@ -147,7 +147,8 @@ namespace corsika::process::energy_loss {
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack& t) {
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());
......@@ -169,7 +170,7 @@ namespace corsika::process::energy_loss {
p.SetEnergy(Enew);
MomentumUpdate(p, Enew);
fEnergyLossTot += dE;
GetXbin(p, dE);
GetXbin(p, t, dE);
return status;
}
......@@ -199,17 +200,17 @@ namespace corsika::process::energy_loss {
#include <corsika/geometry/CoordinateSystem.h>
int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& vP,
int EnergyLoss::GetXbin(SetupParticle const& vP, SetupTrack const& vTrack,
const HEPEnergyType dE) {
using namespace corsika::geometry;
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point pos1(rootCS, 0_m, 0_m, 0_m);
Point pos2(rootCS, 0_m, 0_m, vP.GetPosition().GetCoordinates()[2]);
Vector delta = (pos2 - pos1) / 1_s;
Trajectory t(Line(pos1, delta), 1_s);
Point const pos1(rootCS, 0_m, 0_m, 0_m);
Point const pos2(rootCS, 0_m, 0_m, vTrack.GetPosition(0).GetCoordinates()[2]);
auto const delta = (pos2 - pos1) / 1_s;
Trajectory const t(Line(pos1, delta), 1_s);
GrammageType const grammage =
vP.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
......
......@@ -25,35 +25,30 @@ namespace corsika::process::energy_loss {
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));
units::si::square(1e-2 * units::si::meter));
void MomentumUpdate(corsika::setup::Stack::ParticleType&,
corsika::units::si::HEPEnergyType Enew);
void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew);
public:
EnergyLoss();
void Init() {}
process::EProcessReturn DoContinuous(setup::Stack::ParticleType&,
setup::Trajectory const&);
units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const&) const;
corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&,
corsika::setup::Trajectory&);
corsika::units::si::LengthType MaxStepLength(
corsika::setup::Stack::ParticleType const&,
corsika::setup::Trajectory const&) const;
corsika::units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
void PrintProfile() const;
private:
static corsika::units::si::HEPEnergyType BetheBloch(
corsika::setup::Stack::ParticleType const& p,
const corsika::units::si::GrammageType dX);
static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const& p,
const units::si::GrammageType dX);
int GetXbin(corsika::setup::Stack::ParticleType& p,
const corsika::units::si::HEPEnergyType dE);
int GetXbin(setup::Stack::ParticleType const&, setup::Trajectory const&, units::si::HEPEnergyType);
corsika::units::si::HEPEnergyType fEnergyLossTot;
corsika::units::si::GrammageType fdX; // profile binning
std::map<int, double> fProfile; // longitudinal profile
units::si::HEPEnergyType fEnergyLossTot;
units::si::GrammageType fdX; // profile binning
std::map<int, double> fProfile; // longitudinal profile
};
} // namespace corsika::process::energy_loss
......
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