#Copper

sheets, Counts1, Counts2, Counts3 = np.loadtxt('AttenuationCopper.csv', delimiter = ',', unpack = True)

sheetthicknessC = 0.66 #cm
thickC = sheets*sheetthicknessC #total thickness
errorsheetC = 0.005 #cm
errorthickC = sheets*errorsheetC
countsError = np.sqrt(AvgCountsC)

time = 20
Distance = 12.5 #cm

AvgCountsC = (Counts1+Counts2+Counts3)/3
RateC = AvgCountsC/time
errorRateC = np.sqrt((countsError/AvgCountsC)**2 + (timeError/time)**2)


def fitLine(p, x):
    '''
    Straight line
    '''
    f = p[0] + p[1]*x
    return f
#
def fitError(p, x, y, xerr, yerr):
    '''
    Error function for straight line fit
    '''
    e = (y - fitLine(p, x))/(np.sqrt(yerr**2 + p[1]**2*xerr**2))
    return e
#
# Set initial values of fit parameters, run fit

xDataC = thickC
yDataC = RateC
xErrorC = errorthickC
yErrorC = errorRateC

nPoints = 7

pInit = [1.0, -1.0]
out = least_squares(fitError, pInit, args=(xDataC, yDataC, xErrorC, yErrorC))
#
fitOK = out.success
#
# Test if fit failed
if not fitOK:
    print(" ")
    print("Fit failed")
else:
    #
    # get output
    pFinal = out.x
    cVal = pFinal[0]
    mVal = pFinal[1]
    #
    #   Calculate chis**2 per point, summed chi**2 and chi**2/NDF
    chiarr = fitError(pFinal, xDataC, yDataC, xErrorC, yErrorC)**2
    chisq = np.sum(chiarr)
    NDF = nPoints - 2
    redchisq = chisq/NDF
#
    np.set_printoptions(precision = 3)
    print(" ")
    print("Fit quality:")
    print("chisq per point = \n",chiarr)
    print("chisq = {:5.2f}, chisq/NDF = {:5.2f}.".format(chisq, redchisq))
    #
    # Compute covariance
    jMat = out.jac
    jMat2 = np.dot(jMat.T, jMat)
    detJmat2 = np.linalg.det(jMat2)
    #
    if detJmat2 < 1E-32:
        print("Value of determinat detJmat2",detJmat2)
        print("Matrix singular, error calculation failed.")
        print(" ")
        print("Parameters returned by fit:")
        print("Intercept = {:5.2f}".format(cVal))
        print("Gradient = {:5.2f}".format(mVal))
        print(" ")
        cErr = 0.0
        mErr = 0.0
    else:
        covar = np.linalg.inv(jMat2)
        cErr = np.sqrt(covar[0, 0])
        mErr = np.sqrt(covar[1, 1])
        #
        print(" ")
        print("Parameters returned by fit:")
        print("Intercept (Sin theta i) = {:5.2f} +- {:5.2f}".format(cVal, cErr))
        print("Gradient () = {:5.2f} +- {:5.2f}".format(mVal, mErr))
        print(" ")
    #
    # Calculate fitted function values
    fitDataC = fitLine(pFinal, xDataC)
    #


    # Plot data
fig = plt.figure(figsize = (8, 6))
plt.title('Data with fit')
plt.xlabel('Thickness of sheets')
plt.ylabel('Rate')
plt.errorbar(xDataC, yDataC, xerr = xErrorC, yerr = yErrorC, fmt='r', \
             linestyle = '', label = "Data")
plt.plot(xDataC, fitDataC, color = 'b', linestyle = '-', label = "Fit")
#plt.xlim(1.0, 7.0)
#plt.ylim(0.0, 3.0)
plt.grid(color = 'g')

plt.savefig("Copper.png")
plt.show()