import numpy as np
import pandas as pd
from pydiana.tools import pyroot_fitting as prf
import scipy as sc
import plotly.graph_objects as go
import pydiana.tools.dataframe_manip as dfman
import plotly.subplots as psb
import copy
from pydiana.shell.plotly_theme import plotly_color_palette
from pydiana.tools import graph_tools as gt
from pydiana.tools import basic_tools as bts
import plotly.express as px
#####################
# #
# Linearity Fit #
# #
#####################
[docs]
def linearity_fit(x,y,err_y=None,err_x=None,
unit_x = 'cycle',unit_y = 'mV',
a_lims:list=None,b_lims:list=None,
xlabel='Cycles',ylabel='Mean',
points2skip:list=None,
add_results:bool=True,
showplot:bool=True,
fit_options:str='QS',
sort:bool=True,
):
"""
Function for fitting a second degree polynomial to x and y and plot residuals plot.
Parameters
----------
x : array-like containing the number of cycles used for the LED
y : array-like containing the means
err_y : array-like containing the errors on the y
err_x : array-like containing the errors on the x (if left to None 0 errors will be applied)
unit_x : string for indicating the x units
unit_y : string for indicating the y units
a_lims : array-like containing min and max range for the "a" parameter
b_lims : array-like containing min and max range for the "b" parameter
xlabel : string containing the xlabel of the plot
ylabel : string containing the ylabel of the plot
points2skip : list with the indixes of the points to be skipped in the fit in increasing order of "x" (start from 0). If None everything is used
add_results : bool for enabling/disabling the writing of the results on the plot
showplot : switch for enabling the showing of the plot
Returns
-------
fit_results : pandas dataframe containing fit results
fig : figure instance
"""
#initializing fitting function
parlims = None
if a_lims is not None:
parlims = {0:a_lims}
if b_lims is not None:
if parlims is None:
parlims = {1:b_lims}
else:
parlims[1]=b_lims
xrange = (min(x),max(x))
linearity =prf.root_func(name="linearity",
expr = "[0]*x*(1+[1]*[0]*x)",
xrange = xrange,
parnames = {0:f"a[{unit_y}/{unit_x}]",1:f'b[1/{unit_y}]'},
parlims = parlims
)
results,fig=gt.make_fit_scatter(x=x,y=y,
func=linearity,
err_x=err_x,err_y=err_y,
xlab=f'{xlabel} [{unit_x}]',
ylab=f"{ylabel} [{unit_y}]",
title=f"Linearity measurement and Residuals <br> Fit function {ylabel}({xlabel}) = a*{xlabel}*(1+b*a*{xlabel})",
logx=False,
logy=False,
showplot=False,
points2skip=points2skip,
residuals=True,
summary=False,
sort=sort,
fit_options= fit_options,
get_detail_matrices=False,
)
print(f"Linearization valid up to: {ylabel} = ({-1/(4*results['Values']['b']):.2e} +/- {-results['Errors']['b']/(4*results['Values']['b']**2):.2e}) {unit_y}")
if add_results: #block part to add to plot a square with fit results and remove legend
res = copy.deepcopy(results)
res["Values"]["b"] = res["Values"]["b"]*100
res["Errors"]["b"] = res["Errors"]["b"]*100
res["Units"]["b"] = f"%/{unit_y}"
fig = gt.add_summary_legend(fig,res)
if showplot:
fig.show()
return results,fig
[docs]
def apply_linearity(df,params,varcol:str,sep_field:str=None):
"""
Function for applying the linearization procedure to the data
Parameters
----------
df : pandas dataframe with the data
params : pandas dataframe with the results of one linearization fit or dictionary of dataframes with various results from linearization fits (must specify sep_field)
varcol : String with the name of the dataframe column where to calculate the linearization
sep_field : Used only if params is a dict, string containing the name of the column of df that contains the value of param's keys for each row of df.
Returns
-------
1D array with the linearized version of the specified variable
"""
linearized = None
linearization = lambda amp,a,b: (-1 + np.sqrt(1+4*b*amp))/(2*b)
if isinstance(params,dict):
if sep_field is None:
print('Please specify a separation field if params is a dict')
return None
sep_field = df[sep_field]
linearized = [linearization(amp,
params[s]['Values']['a'] ,
params[s]['Values']['b'] ) if s in list(params.keys()) else np.nan for amp,s in zip(df[varcol],sep_field)]
linearized = np.array(linearized)
else:
linearized = linearization(df[varcol],
params['Values']['a'],
params['Values']['b']
)
return linearized
[docs]
def plot_linearity_uncertainty(params,xrange,varname :str = 'Variable',npoints:int=1000,showplot:bool=False,title=None):
"""
Function for applying the linearization procedure to the data
Parameters
----------
df : pandas dataframe with the data
params : pandas dataframe with the results of one linearization fit or dictionary of dataframes with various results from linearization fits (must specify sep_field)
varcol : String with the name of the dataframe column where to calculate the linearization
sep_field : Used only if params is a dict, string containing the name of the column of df that contains the value of param's keys for each row of df.
Returns
-------
1D array with the linearized version of the specified variable
"""
error = None
errorfunc = lambda amp,b,err_b: np.abs((-1 + np.sqrt(1+4*b*amp) -2*amp*b)/(2*b*b*np.sqrt(1+4*amp*b)))*np.abs(err_b)
linearized = None
linearization = lambda amp,a,b: (-1 + np.sqrt(1+4*b*amp))/(2*b)
df2plot = pd.DataFrame()
if title is None:
title = f'Uncertainty on linearization of {varname}'
if isinstance(params,dict):
varrange = [np.linspace(np.min(xrange[s]),
np.max(xrange[s]),npoints) for s in params]
error = [errorfunc(varrange[i],
params[s]['Values']['b'] ,
params[s]['Errors']['b'] ) for i,s in enumerate(list(params.keys())) ]
linearized = [linearization(varrange[i],
params[s]['Values']['a'] ,
params[s]['Values']['b'] ) for i,s in enumerate(list(params.keys())) ]
df2plot[varname] = np.concatenate(varrange)
df2plot[f"Linearized {varname}"] = np.concatenate(linearized)
df2plot['Keys'] = np.concatenate([[s]*len(error[i]) for i,s in enumerate(list(params.keys())) ])
df2plot['Error'] = np.concatenate(error)
else:
varrange = np.linspace(np.min(xrange),np.max(xrange),npoints)
error = errorfunc(varrange,
params['Values']['b'],
params['Errors']['b']
)
linearized = linearization(varrange,
params['Values']['a'],
params['Values']['b']
)
df2plot[varname] = varrange
df2plot[f"Linearized {varname}"] = linearized
df2plot['Error'] = error
df2plot = df2plot.dropna()
fig = px.line(df2plot, y=varname, x=f"Linearized {varname}", color='Keys' if isinstance(params,dict) else None,
)
x_left = df2plot[f"Linearized {varname}"]-df2plot['Error']
x_right = df2plot[f"Linearized {varname}"]+df2plot['Error']
x=np.concatenate([x_right,x_left[::-1]])
y =np.concatenate([df2plot[varname].values,df2plot[varname].values[::-1]])
fig.add_trace(go.Scatter(
y=y, # x, then x reversed
x=x, # upper, then lower reversed
fill='toself',
fillcolor=plotly_color_palette[1],
line=dict(color='rgba(255,255,255,0)'),
hoverinfo="skip",
name='+/- 1 stddev',
showlegend=True,
))
fig.add_trace(go.Scatter(y=df2plot[varname],x=df2plot[varname],line_dash='dash',line_color='black',showlegend=True,name='Linear Approx'))
fig.update_layout(title=title)
if showplot : fig.show()
return df2plot ,fig
#################################
# #
# Absolute Calibration Fit #
# #
#################################
[docs]
def absolute_calibration_fit(means:list,std_devs:list,
err_means:list,err_std_devs:list,
unit:str = 'mV',photon_wavelength:float = 248,
points2skip :list=None,hasnoise:bool=True,fixs0:bool=False,
photon_wavelength_error:float = 8, showplot:bool =True,
sort:bool=True,get_detail_matrices:bool=False,summary:bool=False,
legendpos_x=0.3,legendpos_y=0.65,fit_options:str='QS',useKeV:bool=False,
respORCalConst ='both', enableSquarePar=False):
"""
Function for performing the absolute calibration fit
UPDATE: now there is the possibility to perform a quadratic fit (instead of a linear one);
this (in principle) will correct eventual gain fluctuation (i.e. responsivity parameter instability)
Parameters
----------
means : array-like with the means for the absolute calibration
std_devs : array-like with the standard deviations for the absolute calibration
err_means : array-like with the errors on the means
err_std_devs : array-like with the errors on the standard deviations
unit : unit of means/std_devs
points2skip : list with the indixes of the points to be skipped in the fit in increasing order of mean amplitude (start from 0). If None everything is used. If sort is True then the order is given from the x axis.
hasnoise : boolean for indicating if the list of points provided contains a noise point
photon_wavelength : wavelength used for calibration in nm
photon_wavelength_error : error on the photon wavelength in nm
showplot : bool switch for enabling plotting
sort : bool for enable point sorting based on x axis
get_detail_matrices : bool for enabling the returning of the covariance and correlation matrices
summary : bool for enabling plot summary legend
legendpos_x : position of the summary legend on the x axis
legendpos_x : position of the summary legend on the y axis
enableSquarePar : bool which enable or not a parameter for the square dependence
Returns
-------
fit_results : pandas dataframe with relevant information on the calibration
fig : figure instance of the calibration
"""
if points2skip is None:
points2fit = range(len(means))[:]
else:
if isinstance(points2skip,int):
points2fit= [i for i in range(len(means)) if i != points2skip]
else:
points2fit = [i for i in range(len(means)) if i not in points2skip]
means = np.array(means)
std_devs = np.array(std_devs)
err_means = np.array(err_means)
err_std_devs = np.array(err_std_devs)
var = std_devs **2
err_var = 2*std_devs * err_std_devs
#photon energy conversion
nm2eV = sc.constants.c*sc.constants.h*1e9/sc.constants.e #ev*nm
photon_energy = None
if photon_wavelength is not None:
photon_energy = nm2eV/ photon_wavelength
photon_energy_error = photon_energy/photon_wavelength * photon_wavelength_error
photon_energy_unit='eV'
if useKeV:
photon_energy_unit='keV'
#initial conditions
responsivity = np.mean(np.diff(var)/np.diff(means))
var_0 = var[0] if hasnoise else np.abs(np.mean(var-responsivity*means))
parfix=None
if fixs0:
parfix = [1]
#initializing fitting function
if enableSquarePar == False:
absol_cal = prf.root_func(expr = '[1]+[0]*x',
name = 'absolute_cal',
parnames={0:f'r[{unit}/photon]',1:f'var_0[{unit}^2]'},
parvals={0:responsivity,1:var_0},
parfix=parfix,
xrange=(min(means),max(means)),
)
fit_result,matrices,fig=gt.make_fit_scatter(x=means,
y=var,
func=absol_cal,
err_x=err_means,
err_y=err_var,
ylab=f"Variance [{unit}<sup>2</sup>]",
xlab=f"Mean [{unit}]",
title="Absolute Calibration<br>Fit function σ² = σ<sub>0</sub>² + r ·μ",
logx=False,
logy=False,
showplot=False,
points2skip=points2skip,
residuals=True,
summary=False,
sort=sort,
fit_options=fit_options,
get_detail_matrices=True,
)
else: ## so enableSquarePar == True
quadTerm=0. ##XXX: to check if it's a good initialization
absol_cal = prf.root_func(expr = '[1]+[0]*x+[2]*x*x',
name = 'absolute_cal_square_enabled',
parnames={0:f'r[{unit}/photon]',1:f'var_0[{unit}^2]',2:f'w'},
parvals={0:responsivity,1:var_0,2:quadTerm},
parfix=parfix,
xrange=(min(means),max(means)),
)
fit_result,matrices,fig=gt.make_fit_scatter(x=means,
y=var,
func=absol_cal,
err_x=err_means,
err_y=err_var,
ylab=f"Variance [{unit}<sup>2</sup>]",
xlab=f"Mean [{unit}]",
title="Absolute Calibration<br>Fit function σ² = σ<sub>0</sub>² + r ·μ + w ·μ²",
logx=False,
logy=False,
showplot=False,
points2skip=points2skip,
residuals=True,
summary=False,
sort=sort,
fit_options=fit_options,
get_detail_matrices=True,
)
#Extract interesting quantities
fit_result['Units']['MinX']=unit
fit_result['Units']['MaxX']=unit
r=fit_result['Values']["r"]
var_0=fit_result['Values']['var_0']
sigma_0 = np.sqrt(var_0)
err_r = fit_result['Errors']['r']
err_var_0 =fit_result['Errors']['var_0']
err_sigma_0 = 0.5*err_var_0/sigma_0
sig = pd.DataFrame()
sig['Variables']=['sigma_0']
sig['Values']=[sigma_0]
sig['Errors'] = [err_sigma_0]
sig['Units'] = [unit]
if hasnoise:
sig = pd.DataFrame()
idx = np.argmin(means)
stdnoise= std_devs[idx]
errstdnoise = err_std_devs[idx]
sig['Variables']=['sigma_0','Noise_s0']
sig['Values']=[sigma_0,stdnoise]
sig['Errors'] = [err_sigma_0,errstdnoise]
sig['Units'] = [unit,unit]
sig = sig.set_index('Variables')
fit_result = dfman.dataframe_insert(fit_result,sig,row='Chi2')
result=fit_result
if (photon_wavelength is not None) and (photon_wavelength_error is not None):
calconst = photon_energy/r
err_calconst = np.sqrt((err_r *calconst/r)**2 + (photon_energy_error/r)**2)
responsivity = r/photon_energy
err_responsivity =np.sqrt( (err_r/photon_energy)**2 + (responsivity*photon_energy_error/photon_energy)**2)
var_0_ev = var_0*calconst*calconst
err_var_0_ev = np.sqrt((err_var_0*calconst*calconst)**2 + (2*var_0*calconst*err_calconst)**2)
sigma_0_ev = sigma_0*calconst
err_sigma_0_ev = np.sqrt((err_sigma_0*calconst)**2 + (sigma_0*err_calconst)**2)
energy_fit_result=pd.DataFrame()
energy_fit_result['Variables'] = ('Calibration Constant',
'responsivity',
'var_0_eV',
'sigma_0_eV',
'photon_energy',
)
energy_fit_result['Values'] =((calconst/1000) if useKeV else calconst,
(responsivity*1000) if useKeV else responsivity,
var_0_ev,
sigma_0_ev,
photon_energy,
)
energy_fit_result['Errors'] = ((err_calconst/1000) if useKeV else err_calconst,
(err_responsivity*1000) if useKeV else err_responsivity,
err_var_0_ev,
err_sigma_0_ev,
photon_energy_error,
)
energy_fit_result['Units'] = (f"{photon_energy_unit}/{unit}",
f"{unit}/{photon_energy_unit}",
"eV^2",
"eV",
"eV",
)
if hasnoise:
idx = np.argmin(means)
stdnoise= std_devs[idx] * calconst
stdnoiseerr= np.sqrt((std_devs[idx] * err_calconst)**2 + (err_std_devs[idx] * calconst)**2)
energy_fit_result=pd.DataFrame()
energy_fit_result['Variables'] = ('Calibration Constant',
'responsivity',
'var_0_eV',
'sigma_0_eV',
'Noise_s0_eV',
'photon_energy',
)
energy_fit_result['Values'] =((calconst/1000) if useKeV else calconst,
(responsivity*1000) if useKeV else responsivity,
var_0_ev,
sigma_0_ev,
stdnoise,
photon_energy,
)
energy_fit_result['Errors'] = ((err_calconst/1000) if useKeV else err_calconst,
(err_responsivity*1000) if useKeV else err_responsivity,
err_var_0_ev,
err_sigma_0_ev,
stdnoiseerr,
photon_energy_error,
)
energy_fit_result['Units'] = (f"{photon_energy_unit}/{unit}",
f"{unit}/{photon_energy_unit}",
"eV^2",
"eV",
"eV",
"eV",
)
energy_fit_result = energy_fit_result.set_index('Variables')
results = dfman.dataframe_insert(fit_result,energy_fit_result,row='Chi2')
if summary:
vars2get = ['Calibration Constant', 'responsivity','Chi2','NDF','Probability','Reduced Chi2','MinX','MaxX']
if 'resp' in respORCalConst:
vars2get = ['responsivity','Chi2','NDF','Probability','Reduced Chi2','MinX','MaxX']
elif 'CalConst' in respORCalConst:
vars2get = ['Calibration Constant','Chi2','NDF','Probability','Reduced Chi2','MinX','MaxX']
if hasnoise:
vars2get = vars2get[:2]+['Noise_s0_eV']+vars2get[2:]
if enableSquarePar:
vars2get = ["w"]+vars2get
idx = [i for i,f in enumerate(results.index.to_list()) if f in vars2get]
print(results)
res2print = results.iloc[idx]#.reindex(vars2get)
print(res2print)
fig = gt.add_summary_legend(figin=fig,
df=res2print,
showplot=False,
y=legendpos_y,
x=legendpos_x,
)
if showplot:
fig.show()
results = [results,fig]
if get_detail_matrices:
results = [results,matrices,fig]
return results