simplified a bit. have to complify later.

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";

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&) {