create_alpha_decay_file.py
Go to the documentation of this file.
1 # Author: M. Reichert
2 # Date : 15.03.2023
3 
4 # Creates a file with alpha decay half-lives. The half-lifes
5 # are calculated by the Viola-Seaborg formula. The fit parameters
6 # are partly taken from Dong&Ren 2005 and have been partly fitted
7 # to experimental available data of the reaclib.
8 # To create this file, a winvn file is needed.
9 import numpy as np
10 from tqdm import tqdm
11 
12 
13 # Path to a valid winvn file
14 path_winvn = "winvne_v2.0.dat"
15 
16 
17 def is_stable_nucleus(Z,N):
18  """
19  Check if the nucleus is a stable one.
20  """
21  Nstable = [ 0, 1, 1, 2, 3, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 10, 10, 10,
22  11, 12, 12, 12, 13, 14, 14, 14, 15, 16, 16, 16, 17, 18, 20, 18, 20, 18,
23  20, 22, 20, 22, 20, 22, 23, 24, 26, 24, 24, 25, 26, 27, 28, 28, 26, 28,
24  29, 30, 30, 28, 30, 31, 32, 32, 30, 32, 33, 34, 36, 34, 36, 34, 36, 37,
25  38, 40, 38, 40, 38, 40, 41, 42, 42, 40, 42, 43, 44, 46, 44, 46, 42, 44,
26  46, 47, 48, 50, 48, 46, 48, 49, 50, 50, 50, 51, 52, 54, 52, 50, 52, 53,
27  54, 55, 56, 52, 54, 55, 56, 57, 58, 60, 58, 56, 58, 59, 60, 62, 64, 60,
28  62, 58, 60, 62, 63, 64, 66, 64, 62, 64, 65, 66, 67, 68, 69, 70, 72, 74,
29  70, 72, 68, 70, 71, 72, 73, 74, 74, 70, 72, 74, 75, 76, 77, 78, 80, 78,
30  74, 76, 78, 79, 80, 81, 82, 82, 78, 80, 82, 84, 82, 82, 83, 85, 86, 88,
31  82, 87, 88, 90, 92, 90, 90, 91, 92, 93, 94, 96, 94, 90, 92, 94, 95, 96,
32  97, 98, 98, 94, 96, 98, 99, 100, 102, 100, 98, 100, 101, 102, 103, 104, 106, 104,
33  104, 105, 106, 107, 108, 108, 108, 109, 110, 112, 110, 111, 112, 113, 114, 116, 114, 116,
34  114, 116, 117, 118, 120, 118, 116, 118, 119, 120, 121, 122, 124, 122, 124, 122, 124, 125,
35  126]
36 
37  Zstable = [ 1, 1, 2, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 8, 8, 8, 9, 10, 10, 10, 11, 12, 12, 12,
38  13, 14, 14, 14, 15, 16, 16, 16, 16, 17, 17, 18, 18, 18, 19, 19, 20, 20, 20, 20, 20, 21, 22, 22,
39  22, 22, 22, 23, 24, 24, 24, 24, 25, 26, 26, 26, 26, 27, 28, 28, 28, 28, 28, 29, 29, 30, 30, 30,
40  30, 30, 31, 31, 32, 32, 32, 32, 33, 34, 34, 34, 34, 34, 35, 35, 36, 36, 36, 36, 36, 36, 37, 38,
41  38, 38, 38, 39, 40, 40, 40, 40, 41, 42, 42, 42, 42, 42, 42, 44, 44, 44, 44, 44, 44, 44, 45, 46,
42  46, 46, 46, 46, 46, 47, 47, 48, 48, 48, 48, 48, 48, 49, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50,
43  51, 51, 52, 52, 52, 52, 52, 52, 53, 54, 54, 54, 54, 54, 54, 54, 54, 55, 56, 56, 56, 56, 56, 56,
44  56, 57, 58, 58, 58, 58, 59, 60, 60, 60, 60, 60, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64,
45  65, 66, 66, 66, 66, 66, 66, 66, 67, 68, 68, 68, 68, 68, 68, 69, 70, 70, 70, 70, 70, 70, 70, 71,
46  72, 72, 72, 72, 72, 73, 74, 74, 74, 74, 75, 76, 76, 76, 76, 76, 77, 77, 78, 78, 78, 78, 78, 79,
47  80, 80, 80, 80, 80, 80, 80, 81, 81, 82, 82, 82, 82 ]
48 
49  Nstable = np.array(Nstable)
50  Zstable = np.array(Zstable)
51 
52  is_stable = np.any((Zstable == Z) & (Nstable == N))
53  return is_stable
54 
55 
56 
57 
58 def get_half_life(Z,N,Qval):
59  """
60  Calculate the half life of an alpha decay. The fit parameters
61  hlog_1, and fitpars_1 are from Dong & Ren 2005
62  (https://ui.adsabs.harvard.edu/abs/2005EPJA...26...69D/abstract).
63 
64  The function returns the a0 reaclib parameter as well as the half-life
65  """
66  hlog_1 = [0e0, 0.8937e0, 0.5720e0, 0.9380e0]
67  fitpars_1 = [1.64062e0, -8.54399e0, -0.19430e0, -33.9054e0]
68 
69  fitpars_2 = [1.71182818, -7.50480827, -0.25315377, -30.70284443]
70  hlog_2 = [0e0,0.047594788143095035,0.12139516604487433,0.39333238949044014]
71 
72  fitpars_3 = [1.70874516, -7.52264755, -0.25153138, -30.82453316]
73  hlog_3 = [0e0, 0.2140471879925886,0.06002092026985373,0.49989717651371746]
74 
75  fitpars_4 = [1.71370933, -7.34225699, -0.2497752 , -30.6826417]
76  hlog_4 = [0e0, -0.12418256402907601,1.179881795937427,0.7165978784278064]
77 
78  # Get correct fit parameters for the individual regions
79  if Z>=82 and N>126:
80  fitpars = fitpars_1
81  hlog = hlog_1
82  elif 82<N<=126 and Z>=82:
83  fitpars = fitpars_2
84  hlog = hlog_2
85  elif 82<N<=126 and Z<82:
86  fitpars = fitpars_3
87  hlog = hlog_3
88  elif 50<N<=82 and 50<=Z<82:
89  fitpars = fitpars_4
90  hlog = hlog_4
91  else:
92  return np.inf,np.inf
93 
94  if Z%2 == 0:
95  if N%2 == 0:
96  helper_idx = 0
97  else:
98  helper_idx = 1
99  else:
100  if N%2 == 0:
101  helper_idx = 2
102  else:
103  helper_idx = 3
104 
105  log10Ta = (fitpars[0]*Z + fitpars[1])*Qval**(-0.5) + (fitpars[2]*Z + fitpars[3]) + hlog[helper_idx]
106  try:
107  Ta = 10**log10Ta
108  a0 = np.log(np.log(2)/Ta)
109  except OverflowError:
110  a0 = np.inf
111  Ta = np.inf
112 
113  return a0,Ta
114 
115 
116 def get_name(Z,N,Zdict,Ndict):
117  """
118  Get the name of a nucleus with Z and N from the
119  dictionaries Zdict and Ndict.
120  """
121  names = np.array(list(Zdict.keys()))
122  Zarray = np.array([x[1] for x in np.array(list(Zdict.items()))]).astype(float)
123  Narray = np.array([x[1] for x in np.array(list(Ndict.items()))]).astype(float)
124  try:
125  idx = np.where((Zarray==Z) & (Narray==N))[0]
126  name = names[idx][0]
127  except IndexError:
128  name = None
129 
130  return name
131 
132 
133 
134 # Get mass excesses and other data from winvn
135 mexc = {}
136 Z_number = {}
137 N_number = {}
138 with open(path_winvn,'r') as f:
139  lines = f.readlines()
140  for line in lines:
141  if len(line.split())==7:
142  mexc[line.split()[0].strip()] = float(line.split()[5].strip())
143  Z_number[line.split()[0].strip()] = float(line.split()[2].strip())
144  N_number[line.split()[0].strip()] = float(line.split()[3].strip())
145 
146 # Prepare output string
147 out = ""
148 
149 # Get mass excess of He4
150 mexc_he4 = mexc['he4']
151 
152 for nuc_name in tqdm(mexc.keys()):
153  if Z_number[nuc_name]<=50:
154  continue
155 
156  if is_stable_nucleus(Z_number[nuc_name],N_number[nuc_name]):
157  continue
158 
159  # Get mass excesses of parent
160  mexc_parent = mexc[nuc_name]
161  # Get mass excess of daughter
162  daughter_name = get_name(Z_number[nuc_name]-2,N_number[nuc_name]-2,Z_number,N_number)
163  if not (daughter_name is None):
164  mexc_daughter = mexc[daughter_name]
165  # Get the Q-value
166  Qval = mexc_parent - mexc_daughter - mexc_he4
167  # Ignore nuclei with negative Q-values
168  if Qval<0:
169  continue
170  # Get the half life and a0
171  a0, Thalf = get_half_life(Z_number[nuc_name],N_number[nuc_name],Qval)
172 
173  # Save the result to a string
174  # Ignore very long lived rates
175  if Thalf>1e20:
176  continue
177  out += nuc_name.rjust(5)+" "+"{:.5e}".format(Thalf)+"\n"
178 
179 
180 # Save the rates
181 with open("alpha_decays.dat","w") as f:
182  f.write(out)
bin.create_alpha_decay_file.get_half_life
def get_half_life(Z, N, Qval)
Definition: create_alpha_decay_file.py:58
bin.create_alpha_decay_file.get_name
def get_name(Z, N, Zdict, Ndict)
Definition: create_alpha_decay_file.py:116
bin.create_alpha_decay_file.is_stable_nucleus
def is_stable_nucleus(Z, N)
Definition: create_alpha_decay_file.py:17