IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 0df2e42d authored by Felix Riehn's avatar Felix Riehn
Browse files

Merge branch... 

Merge branch '243-fpe-due-to-zero-track-length-in-energy-loss-for-very-short-lived-particles' into 'master'

Resolve "FPE due to zero track length in energy loss for very short lived particles"

Closes #243

See merge request AirShowerPhysics/corsika!176
parents 91cb1961 62a2d2ba
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 1
− 1
View file @ 0df2e42d
......@@ -124,7 +124,7 @@ int main() {
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
units::si::TimeType, unsigned short, unsigned short>{
beamCode, E0, plab, injectionPos, 0_ns, nuclA, nuclZ});
}
......
Processes/EnergyLoss/EnergyLoss.cc
+ 24
− 21
View file @ 0df2e42d
......@@ -181,7 +181,7 @@ namespace corsika::process::energy_loss {
}
p.SetEnergy(Enew);
MomentumUpdate(p, Enew);
fEnergyLossTot += dE;
EnergyLossTot_ += dE;
FillProfile(p, t, dE);
return status;
}
......@@ -215,13 +215,13 @@ namespace corsika::process::energy_loss {
using namespace corsika::geometry;
auto const toStart = vTrack.GetPosition(0) - fInjectionPoint;
auto const toEnd = vTrack.GetPosition(1) - fInjectionPoint;
auto const toStart = vTrack.GetPosition(0) - InjectionPoint_;
auto const toEnd = vTrack.GetPosition(1) - InjectionPoint_;
auto const v1 = (toStart * 1_Hz).dot(fShowerAxisDirection);
auto const v2 = (toEnd * 1_Hz).dot(fShowerAxisDirection);
geometry::Line const lineToStartBin(fInjectionPoint, fShowerAxisDirection * v1);
geometry::Line const lineToEndBin(fInjectionPoint, fShowerAxisDirection * v2);
auto const v1 = (toStart * 1_Hz).dot(ShowerAxisDirection_);
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);
SetupTrack const trajToEndBin(lineToEndBin, 1_s);
......@@ -233,8 +233,8 @@ namespace corsika::process::energy_loss {
vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToEndBin,
trajToEndBin.GetLength());
const int binStart = grammageStart / fdX;
const int binEnd = grammageEnd / fdX;
const int binStart = grammageStart / dX_;
const int binEnd = grammageEnd / dX_;
std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
<< grammageEnd << std::endl;
......@@ -242,21 +242,24 @@ namespace corsika::process::energy_loss {
auto energyCount = HEPEnergyType::zero();
auto fill = [&](int bin, GrammageType weight) {
auto const increment = -dE * weight / (grammageEnd - grammageStart);
fProfile[bin] += increment;
energyCount += increment;
std::cout << "filling bin " << bin << " with weight " << weight << ": " << increment
<< std::endl;
const auto dX = grammageEnd - grammageStart;
if (dX > dX_threshold_) {
auto const increment = -dE * weight / (grammageEnd - grammageStart);
Profile_[bin] += increment;
energyCount += increment;
std::cout << "filling bin " << bin << " with weight " << weight << ": "
<< increment << std::endl;
}
};
// fill longitudinal profile
fill(binStart, (1 + binStart) * fdX - grammageStart);
fill(binEnd, grammageEnd - binEnd * fdX);
fill(binStart, (1 + binStart) * dX_ - grammageStart);
fill(binEnd, grammageEnd - binEnd * dX_);
if (binStart == binEnd) { fill(binStart, -fdX); }
if (binStart == binEnd) { fill(binStart, -dX_); }
for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, fdX); }
for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
std::cout << "total energy added to histogram: " << energyCount << std::endl;
}
......@@ -264,8 +267,8 @@ namespace corsika::process::energy_loss {
void EnergyLoss::PrintProfile() const {
std::ofstream file("EnergyLossProfile.dat");
cout << "# EnergyLoss PrintProfile X-bin [g/cm2] dE/dX [GeV/g/cm2] " << endl;
double const deltaX = fdX / 1_g * square(1_cm);
for (auto v : fProfile) {
double const deltaX = dX_ / 1_g * square(1_cm);
for (auto v : Profile_) {
file << v.first * deltaX << " " << v.second / (deltaX * 1_GeV) << endl;
}
}
......
Processes/EnergyLoss/EnergyLoss.h
+ 12
− 8
View file @ 0df2e42d
......@@ -34,8 +34,8 @@ namespace corsika::process::energy_loss {
template <typename TDim>
EnergyLoss(geometry::Point const& injectionPoint,
geometry::Vector<TDim> const& direction)
: fInjectionPoint(injectionPoint)
, fShowerAxisDirection(direction.normalized()) {}
: InjectionPoint_(injectionPoint)
, ShowerAxisDirection_(direction.normalized()) {}
EnergyLoss(setup::Trajectory const& trajectory)
: EnergyLoss(trajectory.GetPosition(0), trajectory.GetV0()){};
......@@ -45,7 +45,7 @@ namespace corsika::process::energy_loss {
setup::Trajectory const&);
units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const&) const;
units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
units::si::HEPEnergyType GetTotal() const { return EnergyLossTot_; }
void PrintProfile() const;
static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const&,
const units::si::GrammageType);
......@@ -60,16 +60,20 @@ namespace corsika::process::energy_loss {
// void FillProfileAbsorbed(setup::Stack::ParticleType const&, setup::Trajectory
// const&);
units::si::HEPEnergyType fEnergyLossTot = units::si::HEPEnergyType::zero();
units::si::GrammageType const fdX = std::invoke([]() {
units::si::HEPEnergyType EnergyLossTot_ = units::si::HEPEnergyType::zero();
units::si::GrammageType const dX_ = std::invoke([]() {
using namespace units::si;
return 10_g / square(1_cm);
}); // profile binning
std::map<int, units::si::HEPEnergyType> fProfile; // longitudinal profile
geometry::Point const fInjectionPoint;
geometry::Vector<units::si::dimensionless_d> const fShowerAxisDirection;
std::map<int, units::si::HEPEnergyType> Profile_; // longitudinal profile
geometry::Point const InjectionPoint_;
geometry::Vector<units::si::dimensionless_d> const ShowerAxisDirection_;
};
const units::si::GrammageType dX_threshold_ = std::invoke([]() {
using namespace units::si;
return 0.0001_g / square(1_cm);
});
} // namespace corsika::process::energy_loss
#endif
Processes/StackInspector/StackInspector.cc
+ 12
− 9
View file @ 0df2e42d
......@@ -32,9 +32,9 @@ template <typename TStack>
StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack,
const HEPEnergyType vE0)
: StackProcess<StackInspector<TStack>>(vNStep)
, fReportStack(vReportStack)
, fE0(vE0)
, fStartTime(std::chrono::system_clock::now()) {}
, ReportStack_(vReportStack)
, E0_(vE0)
, StartTime_(std::chrono::system_clock::now()) {}
template <typename TStack>
StackInspector<TStack>::~StackInspector() {}
......@@ -47,7 +47,7 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
for (const auto& iterP : vS) {
HEPEnergyType E = iterP.GetEnergy();
Etot += E;
if (fReportStack) {
if (ReportStack_) {
geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(); // for printout
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
......@@ -60,12 +60,15 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
}
auto const now = std::chrono::system_clock::now();
const std::chrono::duration<double> elapsed_seconds = now - fStartTime;
const std::chrono::duration<double> elapsed_seconds = now - StartTime_;
std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
double const progress = (fE0 - Etot) / fE0;
auto const dE = E0_ - Etot;
if (dE < dE_threshold_) return process::EProcessReturn::eOk;
double const progress = dE / E0_;
double const eta_seconds = elapsed_seconds.count() / progress;
std::time_t const eta_time = std::chrono::system_clock::to_time_t(
fStartTime + std::chrono::seconds((int)eta_seconds));
StartTime_ + std::chrono::seconds((int)eta_seconds));
cout << "StackInspector: "
<< " time=" << std::put_time(std::localtime(&now_time), "%T")
......@@ -79,8 +82,8 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
template <typename TStack>
void StackInspector<TStack>::Init() {
fReportStack = false;
fStartTime = std::chrono::system_clock::now();
ReportStack_ = false;
StartTime_ = std::chrono::system_clock::now();
}
#include <corsika/cascade/testCascade.h>
......
Processes/StackInspector/StackInspector.h
+ 8
− 4
View file @ 0df2e42d
......@@ -39,12 +39,16 @@ namespace corsika::process {
/**
* To set a new E0, for example when a new shower event is started
*/
void SetE0(const corsika::units::si::HEPEnergyType vE0) { fE0 = vE0; }
void SetE0(const corsika::units::si::HEPEnergyType vE0) { E0_ = vE0; }
private:
bool fReportStack;
corsika::units::si::HEPEnergyType fE0;
decltype(std::chrono::system_clock::now()) fStartTime;
bool ReportStack_;
corsika::units::si::HEPEnergyType E0_;
const corsika::units::si::HEPEnergyType dE_threshold_ = std::invoke([]() {
using namespace units::si;
return 1_eV;
});
decltype(std::chrono::system_clock::now()) StartTime_;
};
} // namespace stack_inspector
......
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