diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f5c09343babcf62f87135a2fadb124c74a144f7e..c219b3d49fe2afb893c2ef7a095e79c5dde92dc8 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -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});
}
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index c9bf6afd5af77f67a64ec9f219be7813c48b22c6..749e574a25e64d6626c45bca7649b5b355e1d746 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -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;
}
}
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index 5061c7c616e34f6137693d2f16b04019ee389ddb..ee4a707cf4890e529a763b5ebffd467dd6231de5 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -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
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index b982cdb7275d2526dbb4a8b91ff0dd952a3626f0..720ccc96fbafa71c8c8865650f550ca3c085d9df 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -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>
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 28d6f04a12ce271e8920f88100848ffdbb93f727..7b6b3aba6bc89475c9f25efb0c1f46df32594064 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -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