<gh_stars>0
"""

Corrections to spectrum shape for beta minus decay

Based on Hayen et al. paper 10.1103/RevModPhys.90.015008

@author J.Kostensalo
@version Apr-2-2021
"""

# Input parameters
#A = 137
#Q = 4.1622  # Q in MeV
#Z = 55   # Final Z
#gA = 1.0    # Axial-vector coupling for radiative corr

#infile = "137Xe_spectrum.csv"
#outfile = 'corrected_spectrum.csv'

A = input("")
Q = input("")  # Q in MeV
Z = input("")   # Final Z
gA = input("")    # Axial-vector coupling for radiative corr
infile = input("")
outfile = input("")



import math
import numpy as np
import scipy.special as scsp
import pandas as pd

alpha = 0.0072973525693 # fine structure constant
me = 0.510998950
W0 = Q/me + 1.0  # W0 in m_ec^2alpha = 0.0072973525693 # fine structure constant
R = 1.2*A**(1/3) / (197.3) * 0.510998950   # R 

excoef = pd.read_csv('exchangetable.dat')



'''
Finite nucelear mass correction
'''
def RN(We):
    M = A*1822.8885
    r0 = W0**2/(2*M**2) - 11/(6*M**2)
    r1 = W0/(3*M**2) 
    r2 = 2/M-4*W0/(3*M**2)
    r3 = 16/(3*M**2)
    return 1+r0+r1/We+r2*We+r3/(We**2)

'''
Atomic mismatch effect
'''
def atomic_mismatch(We):
    return max(0,1-1/(W0-We+0.00001)*(44.2*Z**(0.41)+2.3196*10**(-7)*Z**4.45)/(me*1000000))

def screening(We):
    return math.exp(1.5*math.pi*7.297353e-3**3*Z**(7.0/3.0)/(We*We-1)**1.5*(1/math.tanh(-math.pi*Z*7.297353e-3/math.sqrt(We*We-1))+1))

def radiative(We):
    eps = 0     #1E-6 # For zero devision error problems
    beta = math.sqrt(We*We-1)/We
    g = 3*math.log(940/me)-0.75+4/beta*(-scsp.spence(1-2*beta/(1+beta))) \
        +4*(math.atanh(beta)/beta-1)*((W0-We)/(3*We)-1.5+math.log(2*(W0-We+eps))) \
        +math.atanh(beta+eps)/beta*(2*(1+beta**2)+(W0-We+eps)**2/(6*We**2) \
        -4*math.atanh(beta+eps))
    LAM = math.sqrt(6)/(1.27*A**(1/3))        
    L = 1.026725*(1-2*alpha/(3*math.pi)*math.log(2*W0))**(4.5)
    d2 = Z*alpha**2*( math.log(940/me)-2/3*math.log(2*We)+35/9+math.pi**2/6 \
                     -6*math.log(2) +math.log(LAM/(940/me))-0.395 \
                     - 3/(math.sqrt(10)*math.pi)*LAM/(940/me)* (0.5+(0.395+4/3\
                    ) +math.log((940/me)/LAM)) + 3/(2*math.sqrt(10)*math.pi) \
                    *LAM/(940/me)*(1-math.pi/(2*math.sqrt(10))*LAM/(940/me)) \
                    +3/(math.sqrt(10)*math.pi)*gA*4.706*LAM/(940/me) \
                    *(1/3+0.395+math.log((940/me)/LAM))+math.pi/(4*math.sqrt(10))\
                    *LAM/(940/me))
    a = math.pi/3-3/(2*math.pi)
    b = 4/(3*math.pi)*(11/4-(0.395+4/3-math.log(math.sqrt(10)))-math.pi**2/6)
    d3 = Z**2*alpha**3*(a*math.log(LAM/We)+b*(math.log(2*We)-5/6) \
        +4*math.pi/3*((0.5*(math.log(R)**2-math.log(2*We)**2))+\
        5/3*math.log(2*R*We))-0.649*math.log(2*W0))
    return (1+alpha/(2*math.pi)*(g-3*math.log(940/me/(2.0*W0))))*(L+d2+d3)


'''
Atomic exchange effects
'''
def exchange(We):
    eps=0 #1E-6 # For zero devision error problems
    Wp = We-1
    _,a,b,c,d,e,f,g,h,i = excoef.loc[Z-2,:]
    return max(0,1+a/(Wp+eps)+b/(Wp**2+eps*eps)+c*math.exp(-d*Wp)+e*math.sin((max(0,We-g))**f+h)/(We**i))


#raw_spectrum = pd.read_csv("dev_example.dat", sep=',', lineterminator=';',header=None)
raw_spectrum = pd.read_csv(infile, sep=',',header=None)

energies = np.array(raw_spectrum.loc[:,0])
spectrum = np.array(raw_spectrum.loc[:,1])

energies_we = 1.0 + energies*0.001/me


corrected_spectrum = spectrum*0 # Initialize the corrected spectrum vector

for i in range(0,len(spectrum)):
	corrected_spectrum[i] = spectrum[i]*RN(energies_we[i])*atomic_mismatch(energies_we[i])*screening(energies_we[i])*radiative(energies_we[i])*exchange(energies_we[i])



d = {'EnergykeV': energies, 'Spectrum': corrected_spectrum}
df = pd.DataFrame(data=d)
df.to_csv(outfile, header = False, index=False)