From 27367abf52bd4a1e9b053cb5ab81988874fc9e7d Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Wed, 26 Feb 2025 22:53:57 -0500 Subject: [PATCH] snapshot --- Raphael/dwba_zr.py | 20 ++++++++++---------- Raphael/reactionData.py | 5 +++-- dwuck4/extractData.py | 21 ++++++++++++++++----- 3 files changed, 29 insertions(+), 17 deletions(-) diff --git a/Raphael/dwba_zr.py b/Raphael/dwba_zr.py index 8650f51..8e85ed7 100755 --- a/Raphael/dwba_zr.py +++ b/Raphael/dwba_zr.py @@ -107,7 +107,7 @@ class DWBA_ZR: #---------------------------------------- other constants print("========================================") - D0 = 1.55e+4 # for (d,p) + self.D0 = 1.55e+4 # for (d,p) mass_I = self.dwI.mu k_I = self.dwI.k @@ -126,15 +126,18 @@ class DWBA_ZR: # print(f"spin a : {self.spin_a}") # print(f"spin B : {self.spin_B}") - # self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.s +1) - self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1) + self.fmsq2mb = 10 + self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.s +1) + # self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1) # print(f" spin factor : {self.spinFactor}") - self.xsecScalingfactor = D0 * mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * self.spinFactor + self.xsecScalingfactor = mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * self.spinFactor * self.D0 self.radialInt = None + print(f" (JA, JB, s) : {self.spin_A}, {self.spin_B}, {self.s}") + print(f" Spin factor : {self.spinFactor:.6f}") print(f"Xsec Scaling factor : {self.xsecScalingfactor:.6f}") self.PreCalNineJ() @@ -215,7 +218,8 @@ class DWBA_ZR: wf2 = self.wfu_O[int(L2)][index2] pf1 = np.exp(1j*self.dwI.CoulombPhaseShift(L1)) pf2 = np.exp(1j*self.dwO.CoulombPhaseShift(L2)) - integral = simpson (bs*wf1*wf2, dx=self.boundState.dr) + min_length = min(len(bs), len(wf1), len(wf2)) + integral = simpson(bs[:min_length] * wf1[:min_length] * wf2[:min_length], dx=self.boundState.dr) indexL2 = int(L2 - L1 + self.l) # product = integral * pf1 * pf2 product = integral * pf1 * pf2 * self.massBoverMassA @@ -398,10 +402,6 @@ class DWBA_ZR: else: fact1 = pow(-1, m) * np.power(1j, L1-L2-self.l) * (2*L2+1) * np.sqrt((2*self.l+1)*(2*self.s+1)*(2*L1+1)*(2*J2+1)) fact2 = np.sqrt( quantum_factorial(L2-abs(m)) / quantum_factorial(L2 + abs(m)) ) - # fact3 = clebsch_gordan(J2, mb-m, self.j, m-mb+ma, J1, ma) - # fact4 = clebsch_gordan(L1, 0, self.spin_a, ma, J1, ma) - # fact5 = clebsch_gordan(L2, -m, self.spin_b, mb, J2, mb-m) - # fact6 = clebsch_gordan(L2, 0, self.l, 0, L1, 0) fact3 = self.GetPreCalCG(J2, mb-m, self.j, m-mb+ma, J1, ma) fact4 = self.GetPreCalCG(L1, 0, self.spin_a, ma, J1, ma) fact5 = self.GetPreCalCG(L2, -m, self.spin_b, mb, J2, mb-m) @@ -461,7 +461,7 @@ class DWBA_ZR: xsec += np.abs(haha)**2 # print(f"{ma:4.1f}, {mb:4.1f}, {m:4.0f}, {haha:.6f}, {np.abs(haha)**2:.6e}, {xsec:.6e}") - return xsec * self.xsecScalingfactor * 10 # factor 10 for fm^2 = 10 mb + return xsec * self.xsecScalingfactor * self.fmsq2mb def CalAngDistribution(self, angMin:float = 0, angMax:float = 180, angStep:float = 1): self.angMin = angMin diff --git a/Raphael/reactionData.py b/Raphael/reactionData.py index 17e5f02..bbed58f 100755 --- a/Raphael/reactionData.py +++ b/Raphael/reactionData.py @@ -126,13 +126,14 @@ class ReactionData: self.mass_O = self.dwO.mu - Eout2 = self.ELab + self.Q_value #this is incorrec, but used in ptolmey infileCreator + # Eout2 = self.ELab + self.Q_value #this is incorrec, but used in ptolmey infileCreator print("==================================================") print(self.reactionStr) print(f"Transfer Orbtial : {orbital}") print(f"Q-value : {self.Q_value:10.6f} MeV") print(f"Binding : {self.BindingEnergy:10.6f} MeV") - print(f" Eout : {self.Eout} MeV | {Eout2}") + # print(f" Eout : {self.Eout} MeV | {Eout2}") + print(f" Eout : {self.Eout} MeV ") return True diff --git a/dwuck4/extractData.py b/dwuck4/extractData.py index 4fa56c1..55a01ba 100755 --- a/dwuck4/extractData.py +++ b/dwuck4/extractData.py @@ -183,7 +183,7 @@ def extract_ElasticXsec(file_path=filename): #------------------------------------------------------- -def extract_Xsec(file_path=filename): +def extract_Xsec(factor:float = 10, file_path=filename): x_data, y_data = [], [] start_line = None @@ -203,7 +203,7 @@ def extract_Xsec(file_path=filename): columns = line.split() if len(columns) >= 2: x_data.append(float(columns[0])) # Convert to float - y_data.append(float(columns[1])*10) ## factor convert fm^2 to mb + y_data.append(float(columns[1])*factor) ## factor convert fm^2 to mb return [x_data, y_data] @@ -398,12 +398,23 @@ plot_SMatrix(sAmpOut, sb) # elXsec_data = extract_ElasticXsec() # plot_Xsec(elXsec_data) -xsec_data = extract_Xsec() +JA=2.5 +JB=0 +j=2.5 +D0 = 1.55 +spinFactor=(2*JB+1)/(2*JA+1)/(2*j+1) +scalingFactor=spinFactor*D0*10 + +print(f"spin factor : {spinFactor}") +print(f" D0 : {D0}") +print(f" scaling : {scalingFactor}") + +xsec_data = extract_Xsec(scalingFactor) plot_Xsec(xsec_data) x_data, y_data = xsec_data for i, r in enumerate(x_data): - if i % 5 != 0: - continue + # if i % 5 != 0: + # continue print(f"{{{r:7.3f}, {y_data[i]:10.7f}}},") def plot_RadialMatrix2(ma:float, mb:float, isPlot:bool=True):