IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 56c63903 authored by ralfulrich's avatar ralfulrich
Browse files

first physical energy loss model

parent 3df59dda
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Processes/EnergyLoss/EnergyLoss.cc
+ 86
− 1
View file @ 56c63903
......@@ -38,10 +38,95 @@ namespace corsika::process::EnergyLoss {
: fdEdX(vdEdX)
, fEnergyLossTot(0_GeV) {}
/**
* PDG2018, passage of particles through matter
*
* Note, that Ieff of composite media a determined from lnI = sum
* a_i lnI_i where a_i is the fraction of the electron population
* (~Z_i) of the ith element. This can also be used for shell
* corrections or density effects.
*
* The Ieff of compounds is not bettre than a few percent, if not
* measured explicitly.
*
* For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
*
*/
HEPEnergyType EnergyLoss::BetheBloch(Particle& p, GrammageType const dX) {
// all these are material constants and have to come through Environment
// 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 ZoverA = 0.49976_mol/1_g;
const double x0 = 1.7378;
const double x1 = 4.1323;
const double Cbar = 10.54;
const double delta0 = 0.0;
const double aa = 0.15349;
const double sk = 3.2125;
// end of material constants
// this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
auto const K = 0.307075_MeV / 1_mol * square(1_cm);
HEPEnergyType const E = p.GetEnergy();
HEPMassType const m = particles::GetMass(p.GetPID());
double const gamma = E / m;
double const Z = p.GetChargeNumber();
double const Z2 = pow(Z, 2);
HEPMassType const me = particles::Electron::GetMass();
auto const m2 = m * m;
auto const me2 = me * me;
double const gamma2 = pow(gamma, 2);
double const beta2 = (gamma2-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
double const eta2 = beta2/(1-beta2);
HEPMassType const Wmax = 2*me*c2*beta2*gamma2 / (1 + 2*gamma*me/m + me2/m2);
// Sternheimer parameterization, high energies
// NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization -> MISSING
double const x = log10(p.GetMomentum().GetNorm()/m);
double delta = 0;
if (x>=x1) {
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)
delta = delta0*pow(100,2*(x-x0));
}
// with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
// shell correction
// need more clarity about formulas and units
const double Cadj = 0;
/*
// https://www.nap.edu/read/20066/chapter/8#104
HEPEnergyType Iadj = 12_eV * Z + 7_eV; // Iadj<163eV
if (Iadj>=163_eV)
Iadj = 9.76_eV * Z + 58.8_eV * pow(Z, -0.19); // Iadj>=163eV
double const Cadj = (0.422377/eta2 + 0.0304043/(eta2*eta2) - 0.00038106/(eta2*eta2*eta2)) * 1e-6 * Iadj*Iadj +
(3.858019/eta2 - 0.1667989/(eta2*eta2) + 0.00157955/(eta2*eta2*eta2)) * 1e-9 * Iadj*Iadj*Iadj;
*/
// Bloch correction for higher-order Born approximation
//z2 L2
// Barkas correction
//double barkass = 0;
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) * dX;
}
process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) {
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
HEPEnergyType dE = -dX * fdEdX * pow(p.GetChargeNumber(), 2);
HEPEnergyType dE = BetheBloch(p, dX);
auto E = p.GetEnergy();
const auto Ekin = E - p.GetParticleMass();
auto Enew = E + dE;
......
Processes/EnergyLoss/EnergyLoss.h
+ 2
− 0
View file @ 56c63903
......@@ -44,6 +44,8 @@ namespace corsika::process::EnergyLoss {
void SaveSave();
private:
corsika::units::si::HEPEnergyType BetheBloch(corsika::setup::Stack::ParticleType& p, const corsika::units::si::GrammageType dX);
int GetXbin(corsika::setup::Stack::ParticleType& p,
const corsika::units::si::HEPEnergyType dE);
......
Processes/EnergyLoss/Properties8.dat 0 → 100644
+ 2907
− 0
View file @ 56c63903
This diff is collapsed.
Processes/EnergyLoss/ReadData.py 0 → 100644
+ 73
− 0
View file @ 56c63903
#!/usr/bin/python
import os
import sys
import string
data = open("Properties8.dat", "r")
nEl = 0
nOpt = 0
count = 0 # the data comes with flexible row-blocks, need to count
for line in data.readlines():
line = line.rstrip()
if line == "" or line[0] == "#":
continue
# print str(count) + " " + line
if count == 0:
index = int(line[0:4].strip()) # sternheimers index
name = line[5:11].strip() # tag
sig = int(line[11:14].strip()) # significant figures of atomic mass if element
weight = float(line[16:27].strip()) # atomic weight, if element
weight_error = int((line[31:32].strip() or "0"))# error in last place
ZoverA = float(line[34:42].strip()) # Z/A
rho = float(line[42:52].strip()) # Sternheimers density
rho_corr = float(line[52:63].strip())# Corrected density
state = line[64].strip() # Solid Liquid Gas Diatomicgas
nEl = int(line[66:69].strip()) # number of elements
nAtom = int(line[70:72].strip()) # atoms of el. 1
nOpt = int(line[73:74].strip()) # number of optional lines
type = line[75] # type: Element, Radiactive element, Inorganic compound, Organic compound, Polymer, Mixture, Biological
elif count == 1:
short = line[0:4].strip() # if element
long = line[4:100].strip()
elif count == 2:
desc = line # description and formula
elif count == 3:
Ieff = float(line[0:10].strip()) # ieff
Cbar = float(line[11:19].strip()) # cbar
x0 = float(line[21:28].strip()) # x0
x1 = float(line[30:37].strip()) # x1
aa = float(line[39:46].strip()) # a / aa
sk = float(line[48:55].strip()) # k / sk
delta0 = float(line[59:].strip()) # delta0 / dlt0
elif count == 4 and count < 4+nEl:
elZ = int(line[0:11].strip()) # Z
# num frac.
# weigh frac.
pass
elif count>=4+nEl and count<=4+nEl+nOpt:
# print "opt : " + str(count) + " " + str(4+nEl+nOpt) + " " + line
# property A25
# value F20
# if 1:5 is "Note:" following lines are extra comments
if (count == 4+nEl+nOpt):
count = -1
print str(elZ) + " " + long + " " + str(Ieff) + " " + str(ZoverA) + " " + state + " " + type + " " + str(rho_corr) + " " + str(delta0) + " " + str(x0) + " " + str(x1) + " " + str(aa) + " " + str(sk) + " " + str(Cbar)
count += 1
if (count != 0):
print ("Error " + str(count))
Processes/EnergyLoss/SummaryPropTable.dat 0 → 100644
+ 138
− 0
View file @ 56c63903
# RU Mo 18. Feb 14:07:40 CET 2019
# Data table obtained from http://pdg.lbl.gov/2018/AtomicNuclearProperties/summary_prop_table.dat
# see also http://pdg.lbl.gov/2018/AtomicNuclearProperties/expert.html
# Z A density dEmin Xo lam_I lam_pi_I E_c (elec) E_c (pstrn) state
1 1.00794 0.00 4.103 63.04 52.01 80.26 344.80 338.33 D
1 1.00794 0.07 4.034 63.04 52.01 80.26 278.02 271.50 L
1 2.01410 0.00 2.053 125.98 71.78 110.12 345.50 339.03 D
1 2.01410 0.17 2.019 125.98 71.78 110.12 278.02 271.50 L
2 4.00260 0.00 1.937 94.32 71.00 103.61 257.13 252.23 G
2 4.00260 0.12 1.936 94.32 71.00 103.61 208.25 203.34 L
3 6.94100 0.53 1.639 82.78 71.34 103.26 149.06 145.33 S
4 9.01218 1.85 1.595 65.19 77.80 109.90 113.70 110.68 S
5 10.81100 2.37 1.623 52.69 83.32 115.54 93.95 91.41 S
6 12.01070 2.00 1.749 42.70 85.81 117.79 82.08 79.85 S
6 12.01070 2.21 1.742 42.70 85.81 117.79 81.74 79.51 S
6 12.01070 2.27 1.745 42.70 85.81 117.79 81.67 79.44 S
6 12.01070 3.52 1.725 42.70 85.81 117.79 80.17 77.94 S
7 14.00670 0.00 1.825 37.99 89.68 121.69 91.73 89.71 D
7 14.00670 0.81 1.813 37.99 89.68 121.69 75.47 73.49 L
8 15.99940 0.00 1.801 34.24 90.17 121.91 81.45 79.62 D
8 15.99940 1.14 1.788 34.24 90.17 121.91 66.82 65.04 L
9 18.99840 0.00 1.676 32.93 97.42 127.24 73.15 71.48 D
9 18.99840 1.51 1.634 32.93 97.42 127.24 59.81 58.17 L
10 20.17970 0.00 1.724 28.93 98.99 128.66 67.02 65.47 G
10 20.17970 1.20 1.695 28.93 98.99 128.66 55.10 53.59 L
11 22.98977 0.97 1.639 27.74 102.55 132.22 51.38 49.99 S
12 24.30500 1.74 1.674 25.03 104.12 133.65 46.55 45.25 S
13 26.98154 2.70 1.615 24.01 107.16 136.67 42.70 41.48 S
14 28.08550 2.33 1.664 21.82 108.35 137.74 40.19 39.05 S
15 30.97376 2.20 1.613 21.21 111.36 140.72 37.92 36.84 S
16 32.06500 2.00 1.652 19.50 112.44 141.70 35.85 34.83 S
17 35.45300 0.00 1.630 19.28 115.70 144.94 40.05 39.04 D
17 35.45300 1.57 1.608 19.28 115.70 144.94 34.32 33.35 L
18 39.94800 0.00 1.519 19.55 119.73 148.99 38.03 37.06 G
18 39.94800 1.40 1.508 19.55 119.73 148.99 32.84 31.91 L
19 39.09830 0.86 1.623 17.32 118.98 148.08 31.62 30.72 S
20 40.07800 1.55 1.655 16.14 119.83 148.83 29.56 28.71 S
21 44.95591 2.99 1.522 16.55 123.90 152.91 27.61 26.79 S
22 47.86700 4.54 1.477 16.16 126.20 155.17 26.01 25.23 S
23 50.94150 6.11 1.436 15.84 128.54 157.46 24.67 23.91 S
24 51.99610 7.18 1.456 14.94 129.31 158.16 23.52 22.79 S
25 54.93804 7.44 1.428 14.64 131.44 160.24 22.59 21.89 S
26 55.84500 7.87 1.451 13.84 132.08 160.80 21.68 21.00 S
27 58.93319 8.90 1.419 13.62 134.23 162.90 20.82 20.16 S
28 58.69340 8.90 1.468 12.68 134.06 162.64 20.05 19.41 S
29 63.54600 8.96 1.403 12.86 137.30 165.87 19.42 18.79 S
30 65.38000 7.13 1.411 12.43 138.49 166.99 18.93 18.33 S
31 69.72300 5.90 1.379 12.47 141.22 169.69 18.57 17.98 S
32 72.64000 5.32 1.370 12.25 143.00 171.42 18.16 17.58 S
33 74.92159 5.73 1.370 11.94 144.36 172.73 17.65 17.09 S
34 78.97100 4.50 1.343 11.91 146.71 175.04 17.34 16.79 S
35 79.90400 0.01 1.388 11.42 147.24 175.51 19.16 18.60 D
35 79.90400 3.10 1.380 11.42 147.24 175.51 17.28 16.75 L
36 83.79800 0.00 1.357 11.37 149.42 177.65 18.61 18.06 G
36 83.79800 2.42 1.357 11.37 149.42 177.65 17.03 16.51 L
37 85.46780 1.53 1.356 11.03 150.34 178.51 16.61 16.10 S
38 87.62000 2.54 1.353 10.76 151.50 179.62 15.97 15.47 S
39 88.90584 4.47 1.359 10.41 152.18 180.26 15.30 14.81 S
40 91.22400 6.51 1.349 10.20 153.41 181.44 14.74 14.27 S
41 92.90637 8.57 1.343 9.92 154.28 182.26 14.23 13.77 S
42 95.95000 10.22 1.329 9.80 155.84 183.78 13.85 13.39 S
43 97.90722 11.50 1.325 9.58 156.83 184.72 13.46 13.01 S
44 101.07000 12.41 1.307 9.48 158.40 186.25 13.12 12.68 S
45 102.90550 12.41 1.310 9.27 159.30 187.11 12.84 12.41 S
46 106.42000 12.02 1.289 9.20 161.00 188.76 12.57 12.15 S
47 107.86820 10.50 1.299 8.97 161.68 189.41 12.36 11.94 S
48 112.41400 8.65 1.277 9.00 163.81 191.49 12.21 11.80 S
49 114.81800 7.31 1.278 8.85 164.91 192.55 12.06 11.65 S
50 118.71000 7.31 1.263 8.82 166.67 194.27 11.86 11.46 S
51 121.76000 6.69 1.259 8.73 168.02 195.58 11.70 11.31 S
52 127.60000 6.24 1.227 8.83 170.56 198.08 11.54 11.16 S
53 126.90447 4.93 1.263 8.48 170.25 197.74 11.42 11.04 S
54 131.29300 0.01 1.255 8.48 172.12 199.56 12.30 11.91 G
54 131.29300 2.95 1.255 8.48 172.12 199.56 11.66 11.28 L
55 132.90545 1.87 1.254 8.31 172.79 200.20 11.34 10.97 S
56 137.32700 3.50 1.231 8.31 174.62 201.98 10.98 10.62 S
57 138.90547 6.14 1.231 8.14 175.26 202.59 10.63 10.27 S
58 140.11600 6.77 1.234 7.96 175.75 203.04 10.40 10.04 S
59 140.90766 6.77 1.243 7.76 176.07 203.32 10.21 9.86 S
60 144.24200 7.01 1.231 7.71 177.40 204.62 10.02 9.68 S
61 144.91275 7.26 1.240 7.51 177.66 204.85 9.83 9.49 S
62 150.36000 7.52 1.210 7.57 179.79 206.95 9.66 9.32 S
63 151.96400 5.24 1.219 7.44 180.41 207.53 9.59 9.26 S
64 157.25000 7.90 1.188 7.48 182.42 209.50 9.34 9.01 S
65 158.92535 8.23 1.188 7.36 183.05 210.10 9.17 8.85 S
66 162.50000 8.55 1.175 7.32 184.38 211.39 9.02 8.70 S
67 164.93033 8.80 1.170 7.23 185.27 212.25 8.86 8.55 S
68 167.25900 9.03 1.168 7.14 186.12 213.07 8.73 8.42 S
69 168.93422 9.32 1.170 7.03 186.72 213.64 8.59 8.28 S
70 173.05400 6.90 1.159 7.02 188.19 215.08 8.52 8.22 S
71 174.96680 9.84 1.160 6.92 188.87 215.72 8.36 8.06 S
72 178.49000 13.31 1.152 6.89 190.10 216.93 8.22 7.93 S
73 180.94788 16.65 1.149 6.82 190.96 217.75 8.09 7.79 S
74 183.84000 19.30 1.145 6.76 191.95 218.71 7.97 7.68 S
75 186.20700 21.02 1.143 6.69 192.75 219.49 7.85 7.57 S
76 190.23000 22.57 1.132 6.68 194.11 220.81 7.75 7.47 S
77 192.21700 22.42 1.134 6.59 194.77 221.45 7.67 7.39 S
78 195.08400 21.45 1.128 6.54 195.72 222.37 7.59 7.31 S
79 196.96657 19.32 1.134 6.46 196.34 222.96 7.53 7.26 S
80 200.59000 13.55 1.130 6.44 197.52 224.11 7.53 7.26 L
81 204.38330 11.72 1.125 6.42 198.75 225.30 7.50 7.23 S
82 207.20000 11.35 1.122 6.37 199.64 226.17 7.43 7.16 S
83 208.98040 9.75 1.128 6.29 200.21 226.71 7.39 7.13 S
84 208.98243 9.32 1.141 6.16 200.21 226.69 7.33 7.06 S
85 209.98715 5.00 1.160 6.07 200.52 226.98 7.60 7.33 S
86 222.01758 0.01 1.116 6.28 204.25 230.66 7.71 7.43 G
87 223.01974 1.87 1.118 6.19 204.56 230.94 7.59 7.32 S
88 226.02541 5.00 1.111 6.15 205.47 231.83 7.19 6.94 S
89 227.02775 10.07 1.112 6.06 205.77 232.11 7.00 6.75 S
90 232.03770 11.72 1.098 6.07 207.26 233.57 6.91 6.66 S
91 231.03588 15.37 1.107 5.93 206.96 233.25 6.77 6.52 S
92 238.02891 18.95 1.081 6.00 209.02 235.28 6.65 6.41 S
93 237.04817 20.25 1.095 5.87 208.73 234.97 6.57 6.33 S
94 244.06420 19.84 1.071 5.93 210.77 236.98 6.49 6.25 S
95 243.06138 13.67 1.089 5.80 210.48 236.67 6.50 6.25 S
96 247.07035 13.51 1.082 5.79 211.63 237.79 6.44 6.20 S
97 247.07031 14.00 1.106 5.69 211.63 237.77 6.59 6.35 S
98 251.07959 14.00 1.097 5.68 212.77 238.88 6.53 6.29 S
99 252.08300 14.00 1.102 5.61 213.05 239.15 6.47 6.24 S
100 257.09510 14.00 1.090 5.62 214.45 240.52 6.42 6.19 S
101 258.09843 14.00 1.095 5.55 214.73 240.78 6.37 6.13 S
102 259.10100 14.00 1.099 5.48 215.00 241.04 6.31 6.08 S
103 262.11000 14.00 1.096 5.45 215.83 241.84 6.26 6.03 S
104 267.12200 14.00 1.084 5.47 217.20 243.19 6.21 5.98 S
105 268.12500 14.00 1.088 5.40 217.48 243.44 6.16 5.93 S
106 271.13300 14.00 1.085 5.38 218.29 244.23 6.11 5.89 S
107 270.13400 14.00 1.097 5.27 218.02 243.95 6.06 5.84 S
108 269.13410 14.00 1.110 5.17 217.75 243.66 6.01 5.79 S
109 276.15100 14.00 1.090 5.23 219.63 245.51 5.97 5.75 S
110 281.16200 14.00 1.079 5.24 220.96 246.82 5.92 5.71 S
111 280.16400 14.00 1.091 5.14 220.69 246.54 5.88 5.66 S
112 285.00000 14.00 1.080 5.16 221.96 247.78 5.83 5.61 S
113 284.17800 14.00 1.091 5.07 221.75 247.55 5.79 5.57 S
114 289.18700 14.00 1.080 5.08 223.05 248.83 5.75 5.53 S
115 288.19200 14.00 1.092 4.99 222.79 248.56 5.71 5.49 S
116 293.20000 14.00 1.081 5.00 224.08 249.83 5.67 5.45 S
117 294.00000 14.00 1.086 4.95 224.29 250.01 5.62 5.41 S
118 294.21500 0.01 1.092 4.88 224.34 250.05 5.67 5.45 G
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