From 73da63d0eb1e76bb413ea614a5e8e986b91f547f Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Fri, 22 Feb 2019 23:03:49 +0100
Subject: [PATCH] simplified a bit. have to complify later.
---
Documentation/Examples/cascade_example.cc | 4 ++--
Processes/EnergyLoss/EnergyLoss.cc | 21 ++++++++++++---------
2 files changed, 14 insertions(+), 11 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 36c0fc0d1..f439cf2ab 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -277,8 +277,8 @@ int main() {
// assemble all processes into an ordered process list
// auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut << trackWriter;
// auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut << trackWriter;
- auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
- // auto sequence = stackInspect << sibyll << sibyllNuc << decay << cut << trackWriter;
+ //auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
+ auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index cff0ff5ae..94c686ff2 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -58,7 +58,7 @@ namespace corsika::process::EnergyLoss {
// 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;
- auto Zmat = 7;
+ [[maybe_unused]] auto Zmat = 7;
auto ZoverA = 0.49976_mol/1_g;
const double x0 = 1.7378;
const double x1 = 4.1323;
@@ -84,7 +84,7 @@ namespace corsika::process::EnergyLoss {
double const beta2 = (gamma2-1)/gamma2; // 1-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
cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
- double const eta2 = beta2/(1 - beta2);
+ [[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
@@ -99,7 +99,7 @@ namespace corsika::process::EnergyLoss {
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)
+ } else if (x<x0) { // and IF conductor (otherwise, this is 0)
delta = delta0*pow(100,2*(x-x0));
}
cout << "BetheBloch delta=" << delta << endl;
@@ -123,19 +123,22 @@ namespace corsika::process::EnergyLoss {
double A = 1;
if (p.GetPID() == particles::Code::Nucleus)
A = p.GetNuclearA();
- double const Erel = p.GetEnergy()/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 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;
+ 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;
}
process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) {
--
GitLab