Commit 0d20b954 authored by Alessandro Mirone's avatar Alessandro Mirone

modif per fittare spettro abscalc

parent 1c62a640
......@@ -16,8 +16,11 @@ res_param ={
6:[0.6845,2.477,3.617],
7:[0.6662,2.517,4.353],
8:[0.5879,2.575,3.878],
9:[0.6421,2.639,3.893]
9: [0.6421,2.639,3.893],
100:[0.01,0.005,0.02],
}
#Temperature (important : floating type is mandatory)
T = 297.0
DOIXS=0
\ No newline at end of file
......@@ -34,7 +34,7 @@ from PyMca import Gefit
from numpy import array
DOIXS=1
#--------------------------------------------------------------------------------------------------------
# Basic function Model
......@@ -180,7 +180,7 @@ class Model:
npars = contribution.nofMyParams()
parnames = contribution.parNames()
fact=self.fact4fft*mask[icontribution]
if icontribution==0:
if icontribution==0 or DOIXS==0:
value_and_deri = contribution.ft_and_derivatives( self.Sch, self.xmin , Stokes=0 )
result=result+value_and_deri [0]*fact # so far we exploit only function itself not derivatives..
el_center = contribution.get_Center()
......@@ -234,7 +234,7 @@ def print_logo():
def get_xy(event):
if event.button == 1:
if event.button == 2:
if event.inaxes:
global xy,cid
print 'Peak added -> ', event.xdata, event.ydata
......@@ -435,7 +435,7 @@ def read_configuration_file(cfgfn,allowed_keys={} ):
cfgfn is the filename of the configuration file.
the function return an object containing information from configuration file (cf inside cfg file).
"""
global DOIXS
try:
s=open(cfgfn,"r")
except:
......@@ -449,7 +449,9 @@ def read_configuration_file(cfgfn,allowed_keys={} ):
exec(s)
cfg = Config()
if hasattr(cfg, "DOIXS"):
DOIXS=cfg.DOIXS
for key in allowed_keys.keys() :
if key not in dir(cfg):
raise Exception , ("Key not found in config file : %s"%key)
......@@ -649,6 +651,7 @@ class Params_and_Functions:
File.write('--------------------------------------------------------\n')
def get_dotstripped_path_name( name ):
global HDFD
posslash=name.rfind("/")
posdot =name.rfind(".")
if posdot>posslash:
......@@ -656,24 +659,53 @@ def get_dotstripped_path_name( name ):
else:
pass
if(posslash>-1):
return name, name[posslash+1:]
if HDFD: return name, name[posslash+1:]
else: return name[:posslash] , name[posslash+1:]
else:
return name, name
if HDFD: return name, name
else : raise Exception
def extract_data_from(fn, itercount,lowlim, highlim ):
nome=fn+str(itercount)
f=open(nome,"r")
res=[]
for l in f:
nl=string.split( l )
if(len(nl)==2):
nl=map(string.atof,nl)
res.append(nl)
res=np.array(res)
x1 = res[:,0]
x2 = res[:,1]
sigmas=x2*0+0.0001
mask = np.less(lowlim,x1)*np.less(x1,highlim)
return x1[mask], x2[mask],sigmas[mask], nome
#--------------------------------------------------------------------------------------------------------
# Main
#--------------------------------------------------------------------------------------------------------
def main(argv, SHOW, BLOCK):
global RECORD, REPLAY
global RECORD, REPLAY, HDFD
print_logo()
fn = argv[1]
cfg_filename = argv[2]
hdf = h5py.File(fn,'r')
try:
hdf = h5py.File(fn,'r')
cfg_filename = argv[2]
HDFD=1
except:
HDFD=0
cfg_filename = argv[2]
orbstart = string.atoi( argv[3] )
orbend = string.atoi( argv[4] )
lowlim = string.atof( argv[5] )
highlim = string.atof( argv[6] )
# the allowed keys will be available as cfg members after reading parameter file
allowed_keys={"res_param":DictType,"T":FloatType}
allowed_keys={"res_param":DictType,"T":FloatType, "DOIXS":IntType}
# EXAMPLE OF CFG FILE :
""""
# Parameters for resolution function
......@@ -692,11 +724,19 @@ def main(argv, SHOW, BLOCK):
mod=None
const=None
CONVOLUTION_METHOD="PSEUDOVOIGT"
itercount=0
while( HDFD or ( itercount <=(orbend-orbstart) ) ):
while(1):
if interactive_Entry:
( scan_num , detect_num,
Ene_array ,Intens_array, Intens_Err) = interactive_extract_data_from_h5(hdf)
if not HDFD:
Ene_array ,Intens_array, Intens_Err, filename = extract_data_from(fn, itercount+orbstart, lowlim, highlim)
itercount+=1
scan_num , detect_num = itercount,100
else:
( scan_num , detect_num,
Ene_array ,Intens_array, Intens_Err) = interactive_extract_data_from_h5(hdf)
if CONVOLUTION_METHOD=="PSEUDOVOIGT":
# we build hera a pseudo_voigt for convolution, based on configuration parameters peculiar to the detector
mu,gaussian_w,lorentz_w = cfg.res_param[int(detect_num)]
......@@ -709,7 +749,7 @@ def main(argv, SHOW, BLOCK):
xy, noel = interactive_GUI_get_init_peak_params(Ene_array,Intens_array)
const=None
skip = (xy == []) # xy is a list : [ e0, height0, e1, height....]
if noel : # means : energy range was not containing zero , and elastic peak has not been set by
if noel and DOIXS: # means : energy range was not containing zero , and elastic peak has not been set by
# the above GUI routine. We are going to ask for it now and prepend Ec, Ael to xy
if not REPLAY==0:
exec(getinstruction(REPLAY))
......@@ -730,7 +770,7 @@ def main(argv, SHOW, BLOCK):
# setting up parameter list : ( position, height, width, position, height.... )
param_list = np.zeros([len(xy),3 ],"d")
param_list[:,:2]=xy
wel,wj =0.1,0.1 #widths of elastic and excitation peaks (initial guess)
wel,wj =0.01,0.01 #widths of elastic and excitation peaks (initial guess)
param_list[0,2]=wel
param_list[1:,2]=wj
......@@ -764,7 +804,7 @@ def main(argv, SHOW, BLOCK):
# Plot(mod,params_and_functions.par_array,Ene_array,Intens_array, Intens_Err, show_graph=1)
# plt.show()
const = default_build_constrains(params_and_functions ,position=3,intensity=2,irange=[0.,params_and_functions.maxheight()*2],width=3)
const = default_build_constrains(params_and_functions ,position=3,intensity=2,irange=[-params_and_functions.maxheight(),params_and_functions.maxheight()*4],width=3)
refined_param, chisq, sigmapar = Gefit.LeastSquaresFit(mod.Ft_I ,params_and_functions.par_array ,
constrains=const ,xdata=Ene_array ,
......@@ -791,7 +831,11 @@ def main(argv, SHOW, BLOCK):
print '--------------------------------------------------------------'
print 'Output parameters :'
params_and_functions.print_params(cfg.T,sigmapar, File=sys.stdout) # on the screen
output_dir, output_stripped_name = get_dotstripped_path_name(hdf.filename)
if HDFD:
output_dir, output_stripped_name = get_dotstripped_path_name(hdf.filename)
else:
output_dir, output_stripped_name = get_dotstripped_path_name(filename)
if not os.path.exists(output_dir):
os.mkdir(output_dir)
out = open('%s/%s_%s_%s.param'%(output_dir,output_stripped_name ,scan_num,detect_num),'w')
......@@ -839,7 +883,8 @@ def main(argv, SHOW, BLOCK):
plt.close()
# now we exit from the main
hdf.close()
if(HDFD) :
hdf.close()
#--------------------------------------------------------------------------------------------------------
def getinstruction(REPLAY):
......@@ -854,15 +899,15 @@ if __name__ == '__main__':
SHOW=1
BLOCK=0
argv = sys.argv
if len(argv) in [3,4]:
if len(argv)==4:
if sys.argv[3]=="RECORD":
if len(argv) in [3,4,7,8]:
if len(argv) in [4,8]:
if sys.argv[len(argv)-1 ]=="RECORD":
open("interactive_session.log","w")
RECORD=1
elif "REPLAY" in sys.argv[3]:
print sys.argv[3]
exec(sys.argv[3])
if "REPLAYSHOW" in sys.argv[3] :
elif "REPLAY" in sys.argv[len(argv)-1]:
print sys.argv[len(argv)-1]
exec(sys.argv[len(argv)-1])
if "REPLAYSHOW" in sys.argv[len(argv)-1] :
REPLAY=open( REPLAYSHOW,"r")
SHOW=1
BLOCK=1
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment