diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 36c0fc0d19d18fa2918c272f88403f43e79d1de5..f439cf2ab3fd895457b49610720c46da1ddd4bd8 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 cff0ff5aee8cf19ff79471ccfe750b6bda21733b..94c686ff2d4186582afa2120d73549965df5c340 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&) {