From 939b63a79f354795f6202c1eed0619f48504d579 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Fri, 14 Mar 2025 15:19:20 -0400 Subject: [PATCH] add (a,a) reaction --- Raphael/reactionData.py | 101 ++++++++++++++++-------- dwuck4/DWtest3.DAT | 34 ++++---- dwuck4/extractData.py | 16 ++-- dwuck4/inFileCreatorDW.py | 158 ++++++++++++++++++++++++++------------ dwuck4/test.dat | 26 ------- 5 files changed, 202 insertions(+), 133 deletions(-) delete mode 100644 dwuck4/test.dat diff --git a/Raphael/reactionData.py b/Raphael/reactionData.py index a14d894..579f309 100755 --- a/Raphael/reactionData.py +++ b/Raphael/reactionData.py @@ -8,6 +8,12 @@ from clebschGordan import obeys_triangle_rule from distortedWave import DistortedWave import opticalPotentials as op +from enum import Enum + +class ReactionType(Enum): + inelastic = 1 #also elastic + transfer = 2 + chargeExchange = 3 def approximate_to_half_integer(value): return round(value * 2) / 2 @@ -50,41 +56,76 @@ class ReactionData: # print("-------- spin_B",self.spin_B) - if self.A_a == 2 and self.Z_a == 1: - self.spin_a = 1.0 - self.spin_b = 0.5 - else: - self.spin_a = 0.5 - self.spin_b = 1.0 - #====== transfering nucleon - self.s = 1/2 # spin of x, neutron or proton - self.A_x = abs(self.A_a - self.A_b) - self.Z_x = abs(self.Z_a - self.Z_b) + #------------ check reaction type + self.reactionType = 0 + if self.A_a == self.A_b and self.Z_a == self.Z_b : + self.reactionType = ReactionType.inelastic - mass_x = iso.GetMassFromAZ(self.A_x, self.Z_x) + if self.A_a == self.A_b and self.Z_a != self.Z_b : + self.reactionType = ReactionType.chargeExchange + print("Charge exchaneg reaction does not support.") + return - #======== core - if self.A_A < self.A_B : # (d,p) - self.A_c = self.A_A - self.Z_c = self.Z_A - self.BindingEnergy = mass_B - mass_A - mass_x + ExB - else: #(p,d) - self.A_c = self.A_B - self.Z_c = self.Z_B - self.BindingEnergy = mass_A - mass_B - ExB - mass_x + if self.A_a != self.A_b : + self.reactionType = ReactionType.transfer + + # if self.A_a == 2 and self.Z_a == 1: + # self.spin_a = 1.0 + # self.spin_b = 0.5 + # else: + # self.spin_a = 0.5 + # self.spin_b = 1.0 - #=================== digest orbital - match = re.search(r'[a-zA-Z]', orbital) # Find first letter - if match: - index = match.start() # Get position of the first letter - - self.node = int(orbital[:index]) - self.l_sym = orbital[index:index+1] - j_sym = orbital[index+1:] - self.j = eval(j_sym) - self.l = op.ConvertLSym(self.l_sym) + self.spin_a = float(eval(re.sub(r'[+-]', '', iso.GetJpi(self.A_a, self.Z_a)))) + self.spin_b = float(eval(re.sub(r'[+-]', '', iso.GetJpi(self.A_b, self.Z_b)))) + if self.reactionType == ReactionType.transfer: + #====== transfering nucleon + self.s = 1/2 # spin of x, neutron or proton + self.A_x = abs(self.A_a - self.A_b) + self.Z_x = abs(self.Z_a - self.Z_b) + + mass_x = iso.GetMassFromAZ(self.A_x, self.Z_x) + + #======== core + if self.A_A < self.A_B : # (d,p) + self.A_c = self.A_A + self.Z_c = self.Z_A + self.BindingEnergy = mass_B - mass_A - mass_x + ExB + else: #(p,d) + self.A_c = self.A_B + self.Z_c = self.Z_B + self.BindingEnergy = mass_A - mass_B - ExB - mass_x + + #=================== digest orbital + match = re.search(r'[a-zA-Z]', orbital) # Find first letter + if match: + index = match.start() # Get position of the first letter + + self.node = int(orbital[:index]) + self.l_sym = orbital[index:index+1] + j_sym = orbital[index+1:] + self.j = eval(j_sym) + self.l = op.ConvertLSym(self.l_sym) + + if self.reactionType == ReactionType.inelastic : + self.s = 0 + self.A_x = 0 + self.Z_x = 0 + + mass_x = 0 + + self.A_c = self.A_A + self.Z_c = self.Z_A + self.BindingEnergy = 0 + + self.node = 0 + self.l_sym = "0" + self.l = 2 + self.j = 2 + + #======== regulate spin to half integer self.j = approximate_to_half_integer(self.j) self.s = approximate_to_half_integer(self.s) self.spin_a = approximate_to_half_integer(self.spin_a) diff --git a/dwuck4/DWtest3.DAT b/dwuck4/DWtest3.DAT index 9971539..27a48e0 100644 --- a/dwuck4/DWtest3.DAT +++ b/dwuck4/DWtest3.DAT @@ -1,24 +1,16 @@ 10001310500100000 16O(a,a)16O 6.971 2+, 40MeV +181. +00. +01.0 -+15+01+02+02 ++15+01+02+04 +00.10 +00.000 +15.000 -+40.00 +04. +02. +16. +08. +01.30 +00. -+01. -88.955 +01.149 +00.675 -02.348 +01.345 +00.603 -+02. +01.394 +00.687 +40.872 +01.394 +00.687 --04. -14.228 +00.972 +01.011 +01.562 +00.477 --6.971 +04. +02. +16. +08. +01.30 +00. -+01. -49.544 +01.146 +00.675 -02.061 +01.146 +00.675 -+02. +30.680 +01.302 +00.528 --04. -21.184 +00.934 +00.590 +00.424 +00.934 +00.590 - --00.00 +01. +00. +16. +08. +01.30 +00. -+02.000 -46.380 +01.250 +00.735 +00. +00.000 +01.250 +00.735 +00.000 --03.000 +00.000 +01.250 +00.437 +00.000 +61.400 +01.250 +00.437 +00.000 - -+00.000 +01.000 +00.000 +16.000 +08. +01.250 +00.000 +00.000 +00.000 -+11.000 -46.380 +01.250 +00.735 +00. +00.000 +01.250 +00.735 +00.000 -+00.10 +03. --12.000 +00.000 +01.250 +00.437 +00. +61.400 +01.250 +00.437 +00.000 -+00.10 +03. - -9 END OF DATA DWUCK4 test cases \ No newline at end of file ++40.00 +04. +02. +16. +08. +01.35 +00. ++01. -150.36 +01.342 +00.658 -01.619 +01.426 +00.658 +-02. +01.394 +00.687 +97.576 +01.293 +00.636 +-6.971 +04. +02. +16. +08. +01.35 +00. ++01. -154.37 +01.342 +00.658 -00.644 +01.426 +00.658 +-02. +99.780 +01.293 +00.636 ++00.000 +01.000 +00.000 +16.000 +08. +01.35 +00. ++11.000 -154.37 +01.342 +00.658 -00.644 +01.426 +00.658 ++00.10 +02. +-12. +99.780 +01.293 +00.636 ++00.10 +02. +9 END OF DATA DWUCK4 test cases \ No newline at end of file diff --git a/dwuck4/extractData.py b/dwuck4/extractData.py index 8094b05..52775e4 100755 --- a/dwuck4/extractData.py +++ b/dwuck4/extractData.py @@ -390,20 +390,22 @@ extract_LmaxSaSb() ## must be run first bs_data = extract_BoundState() # plot_BoundState(bs_data) -sAmpIn, sAmpOut = extract_ScatAmp() -plot_SMatrix(sAmpIn, sa) -plot_SMatrix(sAmpOut, sb) +# sAmpIn, sAmpOut = extract_ScatAmp() +# plot_SMatrix(sAmpIn, sa) +# plot_SMatrix(sAmpOut, sb) # elXsec_data = extract_ElasticXsec() # plot_Xsec(elXsec_data) -JA=2.5 -JB=0 -j=2.5 -D0 = 1.55 +JA=0 +JB=2 +j=2 +D0 = 1.0 spinFactor=(2*JB+1)/(2*JA+1)/(2*j+1) scalingFactor=spinFactor*D0*10 +scalingFactor = 5.74 + print(f"spin factor : {spinFactor}") print(f" D0 : {D0}") print(f" scaling : {scalingFactor}") diff --git a/dwuck4/inFileCreatorDW.py b/dwuck4/inFileCreatorDW.py index 7cd917c..c40014b 100755 --- a/dwuck4/inFileCreatorDW.py +++ b/dwuck4/inFileCreatorDW.py @@ -18,7 +18,7 @@ sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) sys.path.append(os.path.join(os.path.dirname(__file__), '../Raphael')) from IAEANuclearData import IsotopeClass import opticalPotentials as op -from reactionData import ReactionData +from reactionData import ReactionData, ReactionType ##################################################### @@ -27,6 +27,12 @@ from reactionData import ReactionData ##################################################### import re +def format_number(x): + formatted = f"{x:+08.4f}" # Default to 4 decimal places + if len(formatted) > 8: # If it exceeds 8 characters, reduce decimal places + formatted = f"{x:+08.{max(0, 8 - len(str(int(x))) - 1)}f}" + return formatted[:8] # Ensure final length is exactly 8 + #================== digest reaction nuclei = re.split(r'[(),]', reaction) @@ -38,6 +44,12 @@ nu_B = nuclei[3] reactionData = ReactionData(nu_A, nu_a, nu_b, JB_pi, orbital, Ex, ELab) +reactionType = reactionData.reactionType + +if reactionType == ReactionType.chargeExchange: + print("charger exchange reaction not supported") + exit() + sym_A = reactionData.sym_A A_A = reactionData.A_A Z_A = reactionData.Z_A @@ -57,6 +69,8 @@ l_sym = reactionData.l_sym spin_a = reactionData.spin_a spin_b = reactionData.spin_b +spin_B = reactionData.spin_B + l = reactionData.l j = reactionData.j @@ -64,9 +78,11 @@ Q_value = reactionData.Q_value BindingEnergy = reactionData.BindingEnergy #=================== outfile name -fileOutName = str(sym_A) + str(A_A) + "_" + str(nu_a) + str(nu_b) + "_" \ +if reactionType == ReactionType.transfer: + fileOutName = str(sym_A) + str(A_A) + "_" + str(nu_a) + str(nu_b) + "_" \ + str(node) + str(l_sym) + str(int(2*j)) + "_" + str(Ex) + "_" + str(ELab) + "_" + JB_pi +".in" - +else: + fileOutName = str(sym_A) + str(A_A) + "_" + str(nu_a) + str(nu_b) + "_" + str(Ex) + "_" + str(ELab) + "_" + JB_pi +".in" #=================== find the maximum L for partial wave mass_I = reactionData.mass_I # reduced mass of incoming channel @@ -89,7 +105,7 @@ with open(fileOutName, "w") as file: if A_a == 1 : pot = op.Koning(A_A, Z_A, A_a*ELab, Z_a) if A_a == 4 : - pot == op.SuAndHan(A_A, Z_A, A_a*ELab) + pot = op.SuAndHan(A_A, Z_A, A_a*ELab) file.write(f"{A_a*ELab:+08.4f}") file.write(f"{A_a:+08.4f}") @@ -102,31 +118,32 @@ with open(fileOutName, "w") as file: file.write(f"{2*spin_a:+08.4f}\n") # Woods-Saxon file.write(f"{1:+08.4f}") - file.write(f"{-pot.v:+08.4f}") # real + file.write(format_number(-pot.v)) # real file.write(f"{pot.r0:+08.4f}") # file.write(f"{pot.a:+08.4f}") # file.write(f"{'':8s}") # spin-orbit skipped file.write(f"{-pot.vi:+08.4f}") # imag file.write(f"{pot.ri0:+08.4f}") # file.write(f"{pot.ai:+08.4f}\n") # + if A_a != 4 and Z_a != 2 : + # Spin-Orbit + file.write(f"{4:+08.4f}") + file.write(f"{-4*pot.vso:+08.4f}") # real + file.write(f"{pot.rso0:+08.4f}") # + file.write(f"{pot.aso:+08.4f}") # + file.write(f"{'':8s}") # spin-orbit skipped + file.write(f"{-4*pot.vsoi:+08.4f}") # imag + file.write(f"{pot.rsoi0:+08.4f}") # + file.write(f"{pot.asoi:+08.4f}\n") # # Woods-Saxon surface - file.write(f"{2:+08.4f}") + file.write(f"{-2:+08.4f}") file.write(f"{'':8s}") # real file.write(f"{'':8s}") # file.write(f"{'':8s}") # file.write(f"{'':8s}") # spin-orbit skipped - file.write(f"{4*pot.vsi:+08.4f}") # imag + file.write(format_number(4*pot.vsi)) # real file.write(f"{pot.rsi0:+08.4f}") # file.write(f"{pot.asi:+08.4f}\n") # - # Spin-Orbit - file.write(f"{-4:+08.4f}") - file.write(f"{-4*pot.vso:+08.4f}") # real - file.write(f"{pot.rso0:+08.4f}") # - file.write(f"{pot.aso:+08.4f}") # - file.write(f"{'':8s}") # spin-orbit skipped - file.write(f"{-4*pot.vsoi:+08.4f}") # imag - file.write(f"{pot.rsoi0:+08.4f}") # - file.write(f"{pot.asoi:+08.4f}\n") # #===== Block 6 Eout = A_a*ELab + Q_value - Ex if A_b == 1 : @@ -147,52 +164,95 @@ with open(fileOutName, "w") as file: file.write(f"{2*spin_b:+08.4f}\n") # Woods-Saxon file.write(f"{1:+08.4f}") - file.write(f"{-pot.v:+08.4f}") # real + file.write(format_number(-pot.v)) # real file.write(f"{pot.r0:+08.4f}") # file.write(f"{pot.a:+08.4f}") # file.write(f"{'':8s}") # spin-orbit skipped file.write(f"{-pot.vi:+08.4f}") # imag file.write(f"{pot.ri0:+08.4f}") # file.write(f"{pot.ai:+08.4f}\n") # + if A_a != 4 and Z_a != 2 : + # Spin-Orbit + file.write(f"{4:+08.4f}") + file.write(f"{-4*pot.vso:+08.4f}") # real + file.write(f"{pot.rso0:+08.4f}") # + file.write(f"{pot.aso:+08.4f}") # + file.write(f"{'':8s}") # spin-orbit skipped + file.write(f"{-4*pot.vsoi:+08.4f}") # imag + file.write(f"{pot.rsoi0:+08.4f}") # + file.write(f"{pot.asoi:+08.4f}\n") # # Woods-Saxon surface - file.write(f"{2:+08.4f}") + file.write(f"{-2:+08.4f}") file.write(f"{'':8s}") # real file.write(f"{'':8s}") # file.write(f"{'':8s}") # file.write(f"{'':8s}") # spin-orbit skipped - file.write(f"{4*pot.vsi:+08.4f}") # imag + file.write(format_number(4*pot.vsi)) # imag file.write(f"{pot.rsi0:+08.4f}") # file.write(f"{pot.asi:+08.4f}\n") # - # Spin-Orbit - file.write(f"{-4:+08.4f}") - file.write(f"{-4*pot.vso:+08.4f}") # real - file.write(f"{pot.rso0:+08.4f}") # - file.write(f"{pot.aso:+08.4f}") # - file.write(f"{'':8s}") # spin-orbit skipped - file.write(f"{-4*pot.vsoi:+08.4f}") # imag - file.write(f"{pot.rsoi0:+08.4f}") # - file.write(f"{pot.asoi:+08.4f}\n") # #====== bound state - file.write(f"{BindingEnergy:+08.4f}") - file.write(f"{A_x:+08.4f}") - file.write(f"{Z_x:+08.4f}") - file.write(f"{A_c:+08.4f}") - file.write(f"{Z_c:+08.4f}") - file.write(f"{1.30:+08.4f}") # Coulomb radius - file.write(f"{'':8s}") # - file.write(f"{'':8s}") # - file.write(f"{1:+08.4f}\n") # neutron spin x 2 - # Woods-Saxon - file.write(f"{-1:+08.4f}") - file.write(f"{-1:+08.4f}") # real - file.write(f"{1.1:+08.4f}") # - file.write(f"{0.65:+08.4f}") # - file.write(f"{24:+8.4f}\n") # spin-orbit - # orbital - file.write(f"{node:+08.4f}") - file.write(f"{l:+08.4f}") - file.write(f"{2*j:+08.4f}") - file.write(f"{1:+08.4f}") # 2 x nuetron spin - file.write(f"{58:+08.4f}\n") + if A_a - A_b == 1: #transfer + file.write(f"{BindingEnergy:+08.4f}") + file.write(f"{A_x:+08.4f}") + file.write(f"{Z_x:+08.4f}") + file.write(f"{A_c:+08.4f}") + file.write(f"{Z_c:+08.4f}") + file.write(f"{1.30:+08.4f}") # Coulomb radius + file.write(f"{'':8s}") # + file.write(f"{'':8s}") # + file.write(f"{1:+08.4f}\n") # neutron spin x 2 + # Woods-Saxon + file.write(f"{-1:+08.4f}") + file.write(f"{-1:+08.4f}") # real + file.write(f"{1.1:+08.4f}") # + file.write(f"{0.65:+08.4f}") # + file.write(f"{24:+8.4f}\n") # spin-orbit + # orbital + file.write(f"{node:+08.4f}") + file.write(f"{l:+08.4f}") + file.write(f"{2*j:+08.4f}") + file.write(f"{1:+08.4f}") # 2 x nuetron spin + file.write(f"{58:+08.4f}\n") + if A_a - A_b == 0 : # inelastic + if A_a == 2 : + pot = op.AnCai(A_A, Z_A, A_a*ELab) + if A_a == 1 : + pot = op.Koning(A_A, Z_A, A_a*ELab, Z_a) + if A_a == 4 : + pot = op.SuAndHan(A_A, Z_A, A_a*ELab) + file.write(f"{BindingEnergy:+08.4f}") + file.write(f"{A_x:+08.4f}") + file.write(f"{Z_x:+08.4f}") + file.write(f"{A_c:+08.4f}") + file.write(f"{Z_c:+08.4f}") + file.write(f"{pot.rc0:+08.4f}") # Coulomb radius + file.write(f"{'':8s}") # + file.write(f"{'':8s}") # + file.write(f"{0:+08.4f}\n") + + file.write(f"{11:+08.4f}") + # file.write(f"{-pot.v:+08.4f}") # real + file.write(format_number(-pot.v)) # real + file.write(f"{pot.r0:+08.4f}") # + file.write(f"{pot.a:+08.4f}") # + file.write(f"{'':8s}") # spin-orbit skipped + file.write(f"{-pot.vi:+08.4f}") # imag + file.write(f"{pot.ri0:+08.4f}") # + file.write(f"{pot.ai:+08.4f}\n") # + + file.write(f"{0.01:+08.4f}") # + file.write(f"{spin_B:+08.4f}\n") # + + file.write(f"{-12:+08.4f}") + file.write(f"{'':8s}") # real + file.write(f"{'':8s}") # + file.write(f"{'':8s}") # + file.write(f"{'':8s}") # spin-orbit skipped + file.write(format_number(4*pot.vsi)) # imag + file.write(f"{pot.rsi0:+08.4f}") # + file.write(f"{pot.asi:+08.4f}\n") # + + file.write(f"{0.01:+08.4f}") # + file.write(f"{spin_B:+08.4f}\n") # #======== end of input file.write("9 end of input card") \ No newline at end of file diff --git a/dwuck4/test.dat b/dwuck4/test.dat deleted file mode 100644 index 28a4735..0000000 --- a/dwuck4/test.dat +++ /dev/null @@ -1,26 +0,0 @@ -1011000030000000 FE56(P,P)FE56* L=3- SPIN ORBIT = OPTION 4 -+37.0000+00.0000+05.0000 -+15+02+03+03 -+00.1000+00.0000+15.0000 - -+22.5000+01.0078+01.0000+56.0000+26.0000+01.2500+00.0000+00.0000+01.0000 -+04. -28.2 +01.25 +00.735 +00. +00. +01.25 +00.735 +00. -+01.0000-46.3800+01.2500+00.7350+00. +00.0000+01.2500+00.7350+00.0000 --02.0000+00.0000+01.2500+00.4370+00.0000+61.4000+01.2500+00.4370+00.0000 - --04.4999+01.0078+01.0000+56.0000+26.0000+01.2500+00.0000+00.0000+01.0000 -+04. -28.2 +01.25 +00.735 +00. +00. +01.25 +00.735 +00. -+01.0000-46.3800+01.2500+00.7350+00. +00.0000+01.2500+00.7350+00.0000 --02.0000+00.0000+01.2500+00.4370+00.0000+61.4000+01.2500+00.4370+00.0000 - -+00.0000+01.0000+00.0000+56.0000+26.0000+01.2500+00.0000+00.0000+00.0000 -+02.0000-46.3800+01.2500+00.7350+00. +00.0000+01.2500+00.7350+00.0000 --03.0000+00.0000+01.2500+00.4370+00.0000+61.4000+01.2500+00.4370+00.0000 - -+00.0000+01.0000+00.0000+56.0000+26.0000+01.2500+00.0000+00.0000+00.0000 -+11.0000-46.3800+01.2500+00.7350+00. +00.0000+01.2500+00.7350+00.0000 -+00.10 +03. --12.0000+00.0000+01.2500+00.4370+00.0000+61.4000+01.2500+00.4370+00.0000 -+00.10 +03. - -9 END OF DATA DWUCK4 test cases \ No newline at end of file