#-----------------------------------------------------------------------------
# Name:        unstable.py
# Purpose:     Unstable macros! Use with caution!
#
# Author:      Andras Wacha
#
# Created:     2010/02/22
# RCS-ID:      $Id: unstable.py$
# Licence:     GPLv2
#-----------------------------------------------------------------------------

import numpy as np
import pylab
import utils
import types
import scipy.stats.stats
import scipy.signal
import matplotlib.widgets
import B1io
import guitools
import asamacros
import fitting
from c_unstable import indirectdesmear
HC=12398.419 #Planck's constant times speed of light, in eV*Angstrom units

#EXPERIMENTAL (DANGER ZONE)
def stackdata(tup):
"""Stack two or more scattering data dictionaries above each other.

Inputs:
tup: a tuple containing the dictionaries

Output:
a scattering data dictionary of the output

NOTE: EXPERIMENTAL!!!!
"""
print "BIG FAT WARNING: stackdata() is an EXPERIMENTAL function. You may not get what you expect!"
data={}
data['q']=np.vstack(tuple([t['q'].reshape(t['q'].size,1) for t in tup]))
data['Intensity']=np.vstack(tuple([t['Intensity'].reshape(t['Intensity'].size,1) for t in tup]))
data['Error']=np.vstack(tuple([t['Error'].reshape(t['Error'].size,1) for t in tup]))
tmp=np.vstack((data['q'].reshape(1,data['q'].size),data['Intensity'].reshape(1,data['Intensity'].size),data['Error'].reshape(1,data['Error'].size)))
tmp=tmp.transpose()
#tmp.sort(0)
data['q']=tmp[:,0]
data['Intensity']=tmp[:,1]
data['Error']=tmp[:,2]
#print data['q'].min()
#print data['q'].max()
return data
def selectasaxsenergies(f1f2,energymin,energymax,Nenergies=3,kT=1000,NITERS=30000,energydistance=0,stepsize=0.5):
"""Select energies for ASAXS measurement by minimizing the condition number of the ASAXS matrix.

Inputs:
f1f2: a numpy array of 3 columns and N rows. Each row should be: energy, f1, f2
energymin: smallest energy
energymax: largest energy
Nenergies: the number of energies
kT: the temperature times Boltzmann's constant, for the Metropolis algorithm
NITERS: how much iterations should we do
energydistance: if two energies are nearer than this value, they are considered
the same
stepsize: the step size for the energies. This should be larger than
energydistance.

Returns:
energies: array of the chosen energies

also two graphs will be plotted
"""

def matrixcond(f1f2,energies,atomicnumber=0):
"""Calculate the condition number of the ASAXS matrix

Inputs:
f1f2: matrix for the anomalous dispersion coefficients
energies: the chosen energies
atomicnumber: the atomic number of the resonant atom. Set it zero
if you want the evaluation according to Stuhrmann.

Outputs:
the condition number. If f1 denotes the column vector of the f1
values and f2 for the f2 values, then the ASAXS matrix is
calculated as:

B=inv(A^T.A).A^T

where

A=[1, 2* (Z+f1), (Z+f1)^2+f2^2]

and Z is the atomic  number.

The 2nd order (=euclidean) condition number of B will be returned.
The np.linalg.linalg.cond() function is used to determine this.  If
the matrix is non-square (ie. rectangular), this type of condition
number can still be determined from the singular value
decomposition.

"""
f1=np.interp(energies,f1f2[:,0],f1f2[:,1])
f2=np.interp(energies,f1f2[:,0],f1f2[:,2])
A=np.ones((len(energies),3));
A[:,1]=2*(f1+atomicnumber);
A[:,2]=(f1+atomicnumber)**2+f2**2;
B=np.dot(np.inv(np.dot(A.T,A)),A.T);
return np.linalg.linalg.cond(B)

np.random.seed()
f1f2=f1f2[f1f2[:,0]<=(energymax+100),:]
f1f2=f1f2[f1f2[:,0]>=(energymin-100),:]
energies=np.random.rand(Nenergies)*(energymax-energymin)+energymin
c=matrixcond(f1f2,energies)
ok=False
oldenergies=energies.copy()
oldc=c
cs=np.zeros(NITERS)
drops=0
eidx=0
sign=0
condmin=c
energiesmin=energies
print 'Initial energies: ',energies
for i in range(NITERS):
oldenergies=energies.copy()
oldc=c
ok=False
while not ok:
#which energy to modify?
eidx=int(np.random.rand()*Nenergies)
#modify energy in which direction?
sign=2*np.floor(np.random.rand()*2)-1
#modify energy
energies[eidx]=energies[eidx]+sign*stepsize
# if the modified energy is inside the bounds and the current
# energies are different, go on.
if energies.min()>=energymin and energies.max()<=energymax and len(utils.unique(energies,lambda a,b:(abs(a-b)<energydistance)))==Nenergies:
ok=True
else: # if not, drop this and re-calculate new energy
energies=oldenergies.copy()
#        print energies
#        print oldenergies
#calculate the condition number of the ASAXS eval. matrix with these energies.
try:
c=matrixcond(f1f2,energies)
except np.linalg.LinAlgError:
energies=oldenergies
c=oldc
drops=drops+1
else:
if c>oldc: #if the condition number is larger than the old one,
if np.random.rand()>(np.exp(c-oldc)/kT): # drop with some probability
energies=oldenergies
c=oldc
drops=drops+1
cs[i]=c # save the current value for the condition number
if np.mod(i,1000)==0: # printing is slow, only print every 1000th step
#            print i
pass
if c<condmin:
condmin=c;
energiesmin=energies.copy()
energies=energiesmin
f1end=np.interp(energies,f1f2[:,0],f1f2[:,1])
f2end=np.interp(energies,f1f2[:,0],f1f2[:,2])
pylab.semilogx(cs)
pylab.xlabel('Step number')
pylab.ylabel('Condition number of the matrix')
a=pylab.gca().axis()
pylab.gca().axis((a[0],a[1],a[2],cs[0]))
pylab.figure()
pylab.plot(f1f2[:,0],f1f2[:,1])
pylab.plot(f1f2[:,0],f1f2[:,2])
pylab.plot(energies,f1end,markersize=10,marker='o',linestyle='')
pylab.plot(energies,f2end,markersize=10,marker='o',linestyle='')
ax=pylab.gca().axis()
pylab.plot([energymin,energymin],[ax[2],ax[3]],color='black',linestyle='--')
pylab.plot([energymax,energymax],[ax[2],ax[3]],color='black',linestyle='--')
pylab.xlabel('Energy (eV)')
pylab.ylabel('f1 and f2')
pylab.title('f1 and f2 values from Monte Carlo simulation.\nkT=%f, N=%d, cond_opt=%f' % (kT,Nenergies,condmin))
print 'Drops: ',drops
print 'Energies: ',energies
print 'f1 values: ',f1end
print 'f2 values: ',f2end
print 'Optimal condition number: ',condmin
print 'Step size: ',stepsize
print 'kT: ',kT
return energies
"""Do radial histogramming on 2D scattering images, according to the idea
of Teemu Ikonen

Inputs:
data: the intensity matrix
energy: the (calibrated) beam energy (eV)
distance: the distance from the sample to the detector (mm)
res: pixel size in mm-s. Both x and y (row and column) direction can
be given if wished, in a list with two elements. A scalar value
means that the pixel size is equal in both directions
bcx: the coordinate of the beam center in the x (row) direction,
starting from ZERO
bcy: the coordinate of the beam center in the y (column) direction,
starting from ZERO
q: the q bins at which the histogram is requested. It should be
defined in 1/Angstroems.
I: the intensity bins

Output:
the histogram matrix
"""
if type(res)!=types.ListType:
res=[res,res];
if len(res)==1:
res=[res[0], res[0]]
if len(res)>2:
raise ValueError('res should be a scalar or a nonempty vector of length<=2')
raise ValueError('data and mask should be of the same shape')
M=data.shape[0] # number of rows
N=data.shape[1] # number of columns

# Creating D matrix which is the distance of the sub-pixels from the origin.
Y,X=np.meshgrid(np.arange(data.shape[1]),np.arange(data.shape[0]));
D=np.sqrt((res[0]*(X-bcx))**2+
(res[1]*(Y-bcy))**2)
# Q-matrix is calculated from the D matrix
q1=4*np.pi*np.sin(0.5*np.arctan(D/float(distance)))*energy/float(HC)
q=np.array(q)
q1=q1[np.isfinite(data)]
data=data[np.isfinite(data)]
# initialize the output matrix
hist=np.zeros((len(I),len(q)))
# set the bounds of the q-bins in qmin and qmax
qmin=map(lambda a,b:(a+b)/2.0,q[1:],q[:-1])
qmin.insert(0,q[0])
qmin=np.array(qmin)
qmax=map(lambda a,b:(a+b)/2.0,q[1:],q[:-1])
qmax.append(q[-1])
qmax=np.array(qmax)
# go through every pixel
for l in range(len(q)):
indices=((q1<=qmax[l])&(q1>qmin[l])) # the indices of the pixels which belong to this q-bin
hist[:,l]=scipy.stats.stats.histogram2(data[indices],I)/np.sum(indices.astype('float'))
return hist
"""2d coloured plot of a matrix with tweaking in the colorscale.

Inputs:
A: the matrix
maxval: maximal value, see plot2dmatrix()
qs: qs see plot2dmatrix()
showqscale: see plot2dmatrix()
pmin: lower scaling limit (proportion, default=0)
pmax: upper scaling limit (proportion, default=1)
"""
f=pylab.figure()
f.donetweakplot=False
a2=pylab.axes((0.1,0.05,0.65,0.02))
a1=pylab.axes((0.1,0.08,0.65,0.02))
ab=pylab.axes((0.85,0.05,0.1,0.1))
ax=pylab.axes((0.1,0.15,0.8,0.75))
button=matplotlib.widgets.Button(ab,'OK')
def finish(a=None,fig=f):
f.donetweakplot=True
button.on_clicked(finish)
sl1=matplotlib.widgets.Slider(a1,'vmin',0,1,pmin)
sl2=matplotlib.widgets.Slider(a2,'vmax',0,1,pmax)
def redraw(tmp=None,ax=ax,sl1=sl1,sl2=sl2):
ax.cla()
pylab.gcf().show()
pylab.draw()
sl1.on_changed(redraw)
sl2.on_changed(redraw)
redraw()
while not f.donetweakplot:
f.waitforbuttonpress()
pylab.close(f)
print sl1.val,sl2.val
return (sl1.val,sl2.val)
def uglyui():
"""Ugly but usable user interface for SAXS and WAXS data treatment
"""
uiparams={'SAXS_bc':None,'SAXS_hperL':None,'WAXS_a':None,'WAXS_b':None,'wavelength':1.54}
choice=-1
try:
choice=int(raw_input("Select a number:"))
except:
choice=-1
return choice
def input_float(prompt='',low=-np.inf,high=np.inf):
val=None
while val is None:
val=raw_input(prompt)
try:
val=float(val)
if (val<low) or (val>high):
val=None
except:
val=None
def input_caseinsensitiveword(prompt='',list=[]):
word=None
while word is None:
word=raw_input(prompt)
if len(list)==0:
return word
word=word.upper()
for w in list:
if word==w.upper():
return w
word=None
def do_desmear(asa):
s=input_float('Smoothing parameter for desmearing (negative to select by hand):')
if s <0:
print "Setting up GUI for smoothing."
smoothlow=input_float('Lowest smoothing value: ',0)
smoothhigh=input_float('Highest smoothing value: ',0)
smoothmode=input_caseinsensitiveword('Mode of the smoothing scale bar (lin or log): ',['lin','log'])
s={'low':smoothlow,'high':smoothhigh,'mode':smoothmode,'val':0.5*(smoothlow+smoothhigh)}
p={}
p['pixelmin']=input_float('Lowest pixel to take into account (starting from 0):',0,len(asa['position']))
p['pixelmax']=input_float('Highest pixel to take into account (starting from 0):',p['pixelmin'],len(asa['position']))
tmp=input_caseinsensitiveword('Do you have a desmearing matrix saved to a file (y or n):',['y','n'])
if tmp=='y':
try:
except:
print "Could not load file. Falling back to manual selection"
tmp='n'
mat=None
if tmp=='n':
p['beamcenter']=input_float('Pixel coordinate of the beam center:')
p['pixelsize']=input_float('Pixel size in micrometers:',0)
p['lengthbaseh']=input_float('Length of the base of the horizontal beam trapezoid',0)
p['lengthtoph']=input_float('Length of the top of the horizontal beam trapezoid',0)
p['lengthbasev']=0
p['lengthtopv']=0
print "Desmearing..."
pixels,desmeared,smoothed,mat,params,smoothing=asamacros.directdesmear(asa['position'],s,p)
x=np.arange(len(asa['position']))
outname=raw_input('Output file name:')
try:
f=open(outname,'wt')
f.write('# pixel\toriginal\tsmoothed\tdesmeared\n')
for i in range(len(pixels)):
f.write('%d\t%g\t%g\t%g\n' %(pixels[i],asa['position'][x==pixels[i]],smoothed[i],desmeared[i]))
f.close()
except:
print "Could not write file %s" % outname
tmp=input_caseinsensitiveword('Would you like to save the smearing matrix for later use (y or n):',['y','n'])
if tmp=='y':
outname=raw_input('File to save the matrix:')
try:
np.savetxt(outname,mat)
except:
print "Could not write file %s" % outname
return pixels,desmeared,smoothed,mat,params,smoothing
'Do Tripalmitine calibration','Desmear',
'Plot original dataset'
'q-calibration of SAXS data',
'q-calibration of WAXS data',
'Show parameters',
'Set parameters'])
if a==0:
return
elif a==1:
fname=raw_input('AgSt measurement file basename (without .P00 extension but may contain path):')
if asa is None:
print "Cannot find file %s.{p00,e00,inf} (If on Linux, check the case)" % fname
return
if raw_input('Desmear before picking peaks (y or n):').upper()=='Y':
pixels,desmeared,smoothed,mat,params,smoothing=do_desmear(asa)
x=pixels
y=desmeared
else:
x=np.arange(len(asa['position']))
y=asa['position']
npeaks=input_float('How many AgSt peaks do you have (at least 2):',2)
a,b,aerr,berr=asamacros.agstcalib(x,y,np.arange(npeaks),returnq=False,wavelength=uiparams['wavelength'])
uiparams['SAXS_bc']=b
uiparams['SAXS_hperL']=a
elif a==2:
fname=raw_input('Tripalmitine measurement file basename (without .P00 extension but may contain path):')
if asa is None:
print "Cannot find file %s.{p00,e00,inf} (If on Linux, check the case)" % fname
return
x=np.arange(len(asa['position']))
y=asa['position']
npeaks=input_float('How many AgSt peaks do you have (at least 2):',2)
a,b,aerr,berr=asamacros.tripcalib(x,y,returnq=False)
uiparams['WAXS_b']=b
uiparams['WAXS_a']=a
elif a==3:
fname=raw_input('Measurement file basename (without .P00 extension but may contain path):')
if asa is None:
print "Cannot find file %s.{p00,e00,inf} (If on Linux, check the case)" % fname
return
pixels,desmeared,smoothed,mat,params,smoothing=do_desmear(asa)
elif a==4:
fname=raw_input('Measurement file basename (without .P00 extension but may contain path):')
if asa is None:
print "Cannot find file %s.{p00,e00,inf} (If on Linux, check the case)" % fname
return
guitools.plotasa(asa)
elif a==5:
if ((uiparams['SAXS_bc'] is None) or (uiparams['SAXS_hperL'] is None) or
(uiparams['wavelength'] is None)):
print """Parameters for SAXS calibration are not yet set. Please set them
via the "Set parameters" or "AgSt calibration" menu items!"""
return
fname=raw_input('Measurement file basename (without .P00 extension but may contain path):')
if asa is None:
print "Cannot find file %s.{p00,e00,inf} (If on Linux, check the case)" % fname
return
x=4*np.pi*np.sin(0.5*np.arctan((np.arange(len(asa['position']))-uiparams['SAXS_bc'])*uiparams['SAXS_hperL']))/uiparams['wavelength']
outfile=raw_input('Output filename:')
pass
def tripcalib2(xdata,ydata,peakmode='Lorentz',wavelength=1.54,qvals=2*np.pi*np.array([0.21739,0.25641,0.27027]),returnq=True):
pcoord=[]
peaks=range(len(qvals))
for p in peaks:
tmp=guitools.findpeak(xdata,ydata,
('Zoom to peak %d (q = %f) and press ENTER' % (p,qvals[p])),
peakmode,scaling='lin')
pcoord.append(tmp)
pcoord=np.array(pcoord)
n=np.array(peaks)
a,b,aerr,berr=fitting.linfit(pcoord,qvals)
if returnq:
return a*xdata+b
else:
return a,b,aerr,berr

q=a*xdata+b
bc=(0-b)/float(a)
alpha=60*np.pi/180.0
h=52e-3
l=150
def xtoq(x,bc,alpha,h,l,wavelength=wavelength):
twotheta=np.arctan((x-bc)*h*np.sin(alpha)/(l-(x-bc)*h*np.cos(alpha)))
return 4*np.pi*np.sin(0.5*twotheta)/wavelength
def costfunc(p,x,y):
pass

"""Create matrix for detector sensitivity correction

Inputs:
fsnrange: FSN range of the sensitivity measurement
energypre: apparent (uncalibrated) energy of the pre-edge measurements
energypost: apparent (uncalibrated) energy of the after-edge measurements
title: title of the sensitivity foil measurements. If the measurements
were recorded with more titles, give a list of strings.
fsnDC: a single number or a list of FSN-s for dark current
measurements
energymeas: apparent energies for energy calibration
energycalib: calibrated energies for energy calibration
energyfluorescence: energy of the fluorescence
origx, origy: the centers of the beamstop (row, column), starting from 1.
a makemask() window will pop up for user interaction.
savefile: a .npz file to save results to. None to skip saving.

Outputs: sensdict
sens: a dictionary of the following fields:
sens: the sensitivity matrix of the 2D detector, by which all
measured data should be divided pointwise. The matrix is
normalized to 1 on the average.
errorsens: the calculated error of the sensitivity matrix
chia, chim, alpha, beta, S1S: these values (4 matrices and
a scalar) are needed for calculation of the correction
terms (taking the dependence of a_0... and S into
account)
"""
# Watch out: this part is VERY UGLY. If you want to understand what
# this does (and you better want it ;-)) please take a look at the
# PDF file attached to the source
global _B1config

pixelsize=_B1config['pixelsize']

# these are controls for non-physical tampering. Their theoretical
# value is in [brackets]
factor=1 # tune this if the scattering is not completely subtracted
# from the fluorescence [1]
t0scaling=1 # multiply the dark-current measurement time by this [1]
hackDCsub=True # set nonpositive elements of DC-subtracted PSD counts to their smallest positive value. [False]
transmerrors=True #if the error of the transmission should be accounted for [True]

# some housekeeping...
if energypost<energypre: # if the two energies were input as swapped
tmp=energypost
energypost=energypre
energypre=tmp
if type(title)==type(''): # make a list of "title"
title=[title]

EB1=[]; hEB1=[]
EB2=[]; hEB2=[]
F1=[]; hF1=[]
F2=[]; hF2=[]
print "makesensitivity: summarizing"
# now summarize...
print "list of sample titles:",title
print "pre-edge energy:",energypre
print "post-edge energy:",energypost

for i in range(len(data)):
EB1.append(data[i])
F1.append(data[i])
else:
print "Unknown title (neither empty beam, nor sample): ","*%s*" % header[i]['Title']
EB2.append(data[i])
F2.append(data[i])
else:
print "Unknown title (neither empty beam, nor sample): ","*%s*" % header[i]['Title']
else:
print "Unknown energy (neither preedge, nor post-edge): ",header[i]['Energy']
print "Neglecting this file."
# go on, this is not an error, just a warning.
#print a short summary
print "Empty beam measurements before the edge:",len(EB2)
print "Empty beam measurements after the edge:",len(EB1)
print "Foil measurements before the edge:",len(F2)
print "Foil measurements after the edge:",len(F1)
print "Energy precision:",_B1config['energyprecision']
#summarize the scattering matrices
D1=sum(F1) # the builtin sum function
D2=sum(F2)
E1=sum(EB1)
E2=sum(EB2)
# summarize the measurement times
t1=sum([h['MeasTime'] for h in hF1])
t2=sum([h['MeasTime'] for h in hF2])
te1=sum([h['MeasTime'] for h in hEB1])
te2=sum([h['MeasTime'] for h in hEB2])
# summarize the anode counts
a1=sum([h['Anode'] for h in hF1])
a2=sum([h['Anode'] for h in hF2])
ae1=sum([h['Anode'] for h in hEB1])
ae2=sum([h['Anode'] for h in hEB2])
# summarize the monitor counts
m1=sum([h['Monitor'] for h in hF1])
m2=sum([h['Monitor'] for h in hF2])
me1=sum([h['Monitor'] for h in hEB1])
me2=sum([h['Monitor'] for h in hEB2])
# calculate the transmissions
T1=np.array([h['Transm'] for h in hF1]).mean()
T2=np.array([h['Transm'] for h in hF2]).mean()
if transmerrors:
print "Taking error of transmission into account"
dT1=np.array([h['Transm'] for h in hF1]).std()
dT2=np.array([h['Transm'] for h in hF2]).std()
else:
print "NOT taking error of transmission into account"
dT1=0
dT2=0
print "Transmission before the edge:",T2,"+/-",dT2
print "Transmission after the edge:",T1,"+/-",dT1
print "monitor_sample_above: ",m1
print "monitor_sample_below: ",m2
print "monitor_empty_above: ",me1
print "monitor_empty_below: ",me2
print "monitor_dark: ",m0

print "anode_sample_above: ",a1
print "anode_sample_below: ",a2
print "anode_empty_above: ",ae1
print "anode_empty_below: ",ae2
print "anode_dark: ",a0

print "time_sample_above: ",t1
print "time_sample_below: ",t2
print "time_empty_above: ",te1
print "time_empty_below: ",te2
print "time_dark: ",t0

print "makesensitivity: calculating auxiliary values"
# error values of anode counts
da1=np.sqrt(a1);    da2=np.sqrt(a2);    dae1=np.sqrt(ae1);
dae2=np.sqrt(ae2);    da0=np.sqrt(a0)

# errors of monitor counts
dm1=np.sqrt(m1);    dm2=np.sqrt(m2);    dme1=np.sqrt(me1);
dme2=np.sqrt(me2);    dm0=np.sqrt(m0)

# errors of 2d detector counts
dD1=np.sqrt(D1);    dD2=np.sqrt(D2);    dE1=np.sqrt(E1);
dE2=np.sqrt(E2);    dD0=np.sqrt(D0)

# Dark current correction: abc -> abcx
D1x=D1-t1/t0*D0 # scattering matrices...
D2x=D2-t2/t0*D0
E1x=E1-te1/t0*D0
E2x=E2-te2/t0*D0
a1x=a1-t1/t0*a0 # anode counts...
a2x=a2-t2/t0*a0
ae1x=ae1-te1/t0*a0
ae2x=ae2-te2/t0*a0
m1x=m1-t1/t0*m0 # monitor...
m2x=m2-t2/t0*m0
me1x=me1-te1/t0*m0
me2x=me2-te2/t0*m0

if hackDCsub:
#In principle, correction for dark current should not make
# scattering images negative. This hack corrects that problem.
print "Tampering with dark current subtraction: setting negative values to positive!"
D1x[D1x<=0]=np.nanmin(D1x[D1x>0])
D2x[D2x<=0]=np.nanmin(D2x[D2x>0])
E1x[E1x<=0]=np.nanmin(E1x[E1x>0])
E2x[E2x<=0]=np.nanmin(E2x[E2x>0])
print "Tampering done."

print "Scattering images, corrected by dark current:"
print utils.matrixsummary(D1x,"D1x")
print utils.matrixsummary(E1x,"E1x")
print utils.matrixsummary(D2x,"D2x")
print utils.matrixsummary(E2x,"E2x")

#two-theta for the pixels
tth=np.arctan(utils2d.calculateDmatrix(D1,pixelsize,origx,origy)/hF1[0]['Dist'])

# angle-dependent correction matrices
C0=gasabsorptioncorrectiontheta(energyfluorescence,tth)
C1,dC1=absorptionangledependenttth(tth,T1,diffaswell=True)
C2,dC2=absorptionangledependenttth(tth,T2,diffaswell=True)
print "Angle-dependent correction matrices:"
print utils.matrixsummary(C1,"C1")
print utils.matrixsummary(C2,"C2")
print utils.matrixsummary(dC1,"dC1")
print utils.matrixsummary(dC2,"dC2")

# some auxiliary variables:
P1=D1x*a1x/(T1*m1x*D1x.sum())*C0*C1
Q1=-E1x*ae1x/(me1x*E1x.sum())*C0
P2=-factor*D2x*a2x/(T2*m2x*D2x.sum())*C0*C2
Q2=factor*E2x*ae2x/(me2x*E2x.sum())*C0

#    pylab.figure();    pylab.imshow(C0);    pylab.title('C0');    pylab.gcf().show(); pylab.colorbar()
#    pylab.figure();    pylab.imshow(C1);    pylab.title('C1');    pylab.gcf().show(); pylab.colorbar()
#    pylab.figure();    pylab.imshow(C2);    pylab.title('C2');    pylab.gcf().show(); pylab.colorbar()
#    pylab.figure();    pylab.imshow(P1);    pylab.title('P1');    pylab.gcf().show(); pylab.colorbar()
#    pylab.figure();    pylab.imshow(P2);    pylab.title('P2');    pylab.gcf().show(); pylab.colorbar()
#    pylab.figure();    pylab.imshow(Q1);    pylab.title('Q1');    pylab.gcf().show(); pylab.colorbar()
#    pylab.figure();    pylab.imshow(Q2);    pylab.title('Q2');    pylab.gcf().show(); pylab.colorbar()

# the unnormalized, unmasked sensitivity matrix.
S1=P1+Q1+P2+Q2
#    pylab.figure();    pylab.imshow(S1);    pylab.title('S1');    pylab.gcf().show(); pylab.colorbar()
pylab.figure()

# multiply the matrix by the mask: masked areas will be zeros.
S1S=S1.sum()
S=S1/S1S # normalize.

# now the S matrix is ready.
print "makesensitivity: calculating error terms"

#we put together the error of S. The terms are denoted ET_<variable>,
# each corresponding to (\frac{d S_{ij}}{d <variable>})^2*Delta^2 <variable>,
# with exception of the matrices, where ET_<mx> corresponds to
# \sum_{mn}(\frac{d S_{ij}}{d <mx>_{mn}})^2Delta^2<mx>_{mn}.

diff=P1/a1x
ET_a1=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*da1**2
diff=P2/a2x
ET_a2=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*da2**2
diff=Q1/ae1x
ET_ae1=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dae1**2
diff=Q2/ae2x
ET_ae2=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dae2**2

diff=P1/m1x
ET_m1=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dm1**2
diff=P2/m2x
ET_m2=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dm2**2
diff=Q1/me1x
ET_me1=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dme1**2
diff=Q2/me2x
ET_me2=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dme2**2

diff=-(t1/t0*P1/a1x+te1/t0*Q1/ae1x+t2/t0*P2/a2x+te2/t0*Q2/ae2x)
ET_a0=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*da0**2
chia=-diff # save it for returning

diff=(t1/t0*P1/m1x+te1/t0*Q1/me1x+t2/t0*P2/m2x+te2/t0*Q2/me2x)
ET_m0=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dm0**2
chim=diff # save it for returning

# the error of the transmissions also have this form, only diff differs.
diff=-P1/T1+P1/C1*dC1
ET_T1=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dT1**2

diff=-P2/T2+P2/C2*dC2
ET_T2=1/S1S**2*(diff**2+S**2*diff.sum()**2-2*S*diff*diff.sum())*dT2**2

#we should take extra care here. Some elements of D1x can be zero.
# This shouldn't be a problem, as zero elements are invalid, so they
# should be under the mask. Numpy signals float divisions by zero
# by setting the result to NaN (not-a-number). Summarizing a matrix
# however, if it contains NaN elements, renders the sum to be NaN
# as well. Therefore we use np.nansum() instead of sum().
alpha=P1/D1x
beta=P1/np.nansum(D1x)**2
dD=dD1
ET_D1=1/S1S**2*(alpha**2-2*alpha*beta-2*S*alpha*(alpha-np.nansum(beta)))*dD**2+ \
1/S1S**2*(beta**2+S**2*np.nansum(beta)**2-2*S*beta*np.nansum(beta))*np.nansum(dD**2) + \
1/S1S**2*S*np.nansum(alpha**2*dD**2)+ \
1/S1S**2*(2*S*beta-2*S**2*np.nansum(beta))*np.nansum(alpha*dD**2)
alpha=Q1/E1x
beta=Q1/np.nansum(E1x)**2
dD=dE1
ET_E1=1/S1S**2*(alpha**2-2*alpha*beta-2*S*alpha*(alpha-np.nansum(beta)))*dD**2+ \
1/S1S**2*(beta**2+S**2*np.nansum(beta)**2-2*S*beta*np.nansum(beta))*np.nansum(dD**2) + \
1/S1S**2*S*np.nansum(alpha**2*dD**2)+ \
1/S1S**2*(2*S*beta-2*S**2*np.nansum(beta))*np.nansum(alpha*dD**2)
alpha=P2/D2x
beta=P2/np.nansum(D2x)**2
dD=dD2
ET_D2=1/S1S**2*(alpha**2-2*alpha*beta-2*S*alpha*(alpha-np.nansum(beta)))*dD**2+ \
1/S1S**2*(beta**2+S**2*np.nansum(beta)**2-2*S*beta*np.nansum(beta))*np.nansum(dD**2) + \
1/S1S**2*S*np.nansum(alpha**2*dD**2)+ \
1/S1S**2*(2*S*beta-2*S**2*np.nansum(beta))*np.nansum(alpha*dD**2)
alpha=Q2/E2x
beta=Q2/np.nansum(E2x)**2
dD=dE2
ET_E2=1/S1S**2*(alpha**2-2*alpha*beta-2*S*alpha*(alpha-np.nansum(beta)))*dD**2+ \
1/S1S**2*(beta**2+S**2*np.nansum(beta)**2-2*S*beta*np.nansum(beta))*np.nansum(dD**2) + \
1/S1S**2*S*np.nansum(alpha**2*dD**2)+ \
1/S1S**2*(2*S*beta-2*S**2*np.nansum(beta))*np.nansum(alpha*dD**2)

alpha=P1/D1x+P2/D2x+Q1/E1x+Q2/E2x
beta=P1/np.nansum(D1x)**2+P2/np.nansum(D2x)**2+Q1/np.nansum(E1x)**2+Q2/np.nansum(E2x)**2
dD=dD0
ET_D0=1/S1S**2*(alpha**2-2*alpha*beta-2*S*alpha*(alpha-np.nansum(beta)))*dD**2+ \
1/S1S**2*(beta**2+S**2*np.nansum(beta)**2-2*S*beta*np.nansum(beta))*np.nansum(dD**2) + \
1/S1S**2*S*np.nansum(alpha**2*dD**2)+ \
1/S1S**2*(2*S*beta-2*S**2*np.nansum(beta))*np.nansum(alpha*dD**2)
#the last alpha and beta are returned!!!

# the error matrix
dS=np.sqrt(ET_a1+ET_a2+ET_ae1+ET_ae2+ET_a0+ET_m1+ET_m2+ET_me1+ET_me2+ET_m0+ \
ET_D1+ET_D2+ET_E1+ET_E2+ET_D0+ET_T1+ET_T2)
print "The error terms:"
print utils.matrixsummary(ET_a1,"ET_a1")
pylab.hist(ET_a1.flatten(),bins=100,log=True)
print utils.matrixsummary(ET_a2,"ET_a2")
print utils.matrixsummary(ET_ae1,"ET_ae1")
print utils.matrixsummary(ET_ae2,"ET_ae2")
print utils.matrixsummary(ET_a0,"ET_a0")
print utils.matrixsummary(ET_m1,"ET_m1")
print utils.matrixsummary(ET_m2,"ET_m2")
print utils.matrixsummary(ET_me1,"ET_me1")
print utils.matrixsummary(ET_me2,"ET_me2")
print utils.matrixsummary(ET_m0,"ET_m0")
print utils.matrixsummary(ET_T1,"ET_T1")
print utils.matrixsummary(ET_T2,"ET_T2")
print utils.matrixsummary(ET_D1,"ET_D1")
print utils.matrixsummary(ET_D2,"ET_D2")
print utils.matrixsummary(ET_E1,"ET_E1")
print utils.matrixsummary(ET_E2,"ET_E2")
print utils.matrixsummary(ET_D0,"ET_D0")
print "----------------"
print utils.matrixsummary(dS**2,"dS^2")
print utils.matrixsummary(dS,"dS")
print utils.matrixsummary(S,"S")
print utils.matrixsummary(dS/S,"dS/S")
# set nans to zero
pylab.figure()
dS[np.isnan(dS)]=0
pylab.colorbar()
pylab.title('Sensitivity')
pylab.gcf().show()
pylab.figure()
pylab.colorbar()
pylab.title('Error of sensitivity')
pylab.gcf().show()
pylab.figure()
pylab.colorbar()
pylab.title('Relative error of sensitivity')
pylab.gcf().show()

#take care of zeros (mask them and set them to 1 afterwards to avoid division by zero)
S[S<=0]=1

if savefile is not None:
if savefile[-4:].upper()=='.MAT':
scipy.io.savemat(savefile,result,appendmat=False)
else:
np.savez(savefile,**result)
return result
def makesensitivity(fsn1,fsn2,fsnend,fsnDC,energymeas,energycalib,energyfluorescence,origx,origy):
"""Create matrix for detector sensitivity correction

Inputs:
fsn1: file sequence number for first measurement of foil at energy E1
fsn2: file sequence number for first measurement of foil at energy E2
fsnend: last FSN in the sensitivity measurement sequence
fsnDC: a single number or a list of FSN-s for dark current
measurements
energymeas: apparent energies for energy calibration
energycalib: calibrated energies for energy calibration
energyfluorescence: energy of the fluorescence
origx, origy: the centers of the beamstop (row,column), starting from 1

sens: the sensitivity matrix of the 2D detector, by which all
measured data should be divided pointwise. The matrix is
normalized to 1 on the average.
errorsens: the error of the sensitivity matrix.
the sensitivity is invalid. This should be used as a base
"""
global _B1config

pixelsize=_B1config['pixelsize']

fsns=range(min(fsn1,fsn2),fsnend+1) # the fsn range of the sensitivity measurement

E1header=[h for h in header if h['FSN']==fsn1 ][0] # the header of the first measurement at E1
# select the fsns for the measurements at E1

E2header=[h for h in header if h['FSN']==fsn2 ][0] # the header of the first measurement at E2
# select the fsns for the measurements at E2
# read data and header structures for E1 and E2 measurements (not empty beams)

# find all the empty beam FSNs for E1
# find all the empty beam FSNs for E2
# read data and header structures for empty beam measurements at E1 and E2

#subtract background, and correct for transmission (and sensitivity :-): to override this correction, ones() and zeros() are given)

#theta for pixels
#transmissions below and above the edge

#angle-dependent absorption
transmcorr1=absorptionangledependenttth(tth,transm1.mean())/transm1.mean()
transmcorr2=absorptionangledependenttth(tth,transm2.mean())/transm2.mean()

#subtract empty beam
B1=(A1/transm1.mean()-eb1)*transmcorr1
B2=(A2/transm2.mean()-eb2)*transmcorr2
errB1=np.sqrt(errA1**2/(transm1.mean())**2+erreb1**2)*transmcorr1
errB2=np.sqrt(errA2**2/(transm2.mean())**2+erreb2**2)*transmcorr2

factor=1 #empirical compensation factor to rule out small-angle scattering completely
C=B1-factor*B2
Cerr=np.sqrt(errB1**2+factor**2*errB2**2)
else:
C=B2-factor*B1
Cerr=np.sqrt(errB2**2+factor**2*errB1**2)
C=C*gasabsorptioncorrectiontheta(energyfluorescence,tth)
cc=C.mean()
sens=C/cc
errorsens=Cerr/cc
#taking care of zeros
sens[sens==0]=1
pylab.imshow(sens)
pylab.colorbar()
pylab.axis('equal')
pylab.gcf().show()
"""Integrate, normalize, do absolute calibration... on a sequence

Inputs:
fsn1: a single number or a list of FSNs.
thicknesses: either one thickness in cm or a dictionary,
containing thickness values for each sample. For example, if
you have measured My_Samplename_1, Al2O3, DPPC_CdS with
thicknesses 1, 1.2 and 0.5 mm, respectively, then you should
give: thicknesses={'My_Samplename_1':0.1,'Al2O3':0.12,
'DPPC_CdS':0.05}
orifsn: which element of the sequence should be used for
determining the origin. 0 means empty beam...
Or you can give a tuple or a list of two: in this case these
will be the coordinates of the origin and no auto-searching
will be performed.
fsndc: one FSN or a list of FSNS corresponding to the empty beam
measurements. If sens is a sensitivity dict, this value will
be disregarded.
sens: sensitivity matrix, or a sensitivity dict.
errorsens: error of the sensitivity matrix. If sens is a sensitivity
dict, this value will be disregarded.
energymeas: list of apparent energies, for energy calibration
energycalib: list of theoretical energies, corresponding to the
apparent energies.
distminus: this will be subtracted from the sample-to-detector
distance read from the measurement files, but only for
samples, not for references
detshift: this will be subtracted from the sample-to-detector
distance read from all measurement files, including
references!
orig: first guess for the origin. A list of two.

Outputs:
qs: q-scales of intensities
ints: intensities
errs: errors
"""
try:
ni=float(thicknesses)
print "Using thickness %f cm for all samples except references" % (ni)
except:
pass

if (len(energycalib)!=len(energymeas)) or len(energycalib)<2:
print "Stopping. Variables energycalib and energymeas should contain equal amount of\npoints and at least two points to be able to make the energy calibration."
return [],[],[],[]
#find the reference measurement
referencenumber=None
for k in range(len(ints)):
referencenumber=k
if k is None:
print "No reference measurements found! At least one file with title\n 'Reference_on_GC_holder_before_sample_sequence' should be part of the sequence!"
return [],[],[],[]
GCdata=None
for r in _B1config['refdata']:
refthick=r['thick']
if GCdata is None:
print "No calibration data exists with ref. position %.2f +/- %.2f." % (header[referencenumber]['PosRef'],_B1config['refposprecision'])
return [],[],[],[]
print "FSN %d: Using GLASSY CARBON REFERENCE with nominal thickness %.f micrometers." % (header[referencenumber]['FSN'],refthick*1e4)

#re-integrate GC measurement to the same q-bins
print "Re-integrating GC data to the same bins at which the reference is defined"
Aerrs[referencenumber],
GCdata[:,0])
print "Re-integration done."
GCdata=GCdata[(AGC>=_B1config['GCareathreshold']) & (intGC>=_B1config['GCintsthreshold']),:]
errGC=errGC[(AGC>=_B1config['GCareathreshold']) & (intGC>=_B1config['GCintsthreshold'])]
qGC=qGC[(AGC>=_B1config['GCareathreshold']) & (intGC>=_B1config['GCintsthreshold'])]
intGC=intGC[(AGC>=_B1config['GCareathreshold']) & (intGC>=_B1config['GCintsthreshold'])]

if len(intGC)<2:
print "ERROR: re-integrated reference does not have enough points! Saving not normalized!"
mult=1
errmult=0
else:
intGC=intGC/refthick
errGC=errGC/refthick

mult,errmult=utils.multfactor(qGC,GCdata[:,1],GCdata[:,2],intGC,errGC)

print "Factor for GC normalization: %.2g +/- %.2f %%" % (mult,errmult/mult*100)
pylab.clf()
pylab.plot(GCdata[:,0],GCdata[:,1],'.',label='Calibrated reference')
pylab.xlabel('q (1/Angstrom)')
pylab.ylabel('Scattering cross-section (1/cm)')
pylab.title('Reference FSN %d multiplied by %.2e, error percentage %.2f' %(header[referencenumber]['FSN'],mult,(errmult/mult*100)))
#pylab.xscale('log')
pylab.yscale('log')
#pause('on')
utils.pause()

#normalize all data to 1/cm
for k in range(len(ints)):
if k!=referencenumber:
thick=None
try:
except:
pass
if thick==None:
try:
thick=float(thicknesses)
except:
print "Did not find thickness for sample %s. Stopping." % header[k]['Title']
print "Using thickness %f cm for sample %s" % (thick,header[k]['Title'])
errs[k]=np.sqrt((mult*errs[k])**2+(errmult*ints[k])**2)/thick
ints[k]=mult*ints[k]/thick
Aerrs[k]=np.sqrt((mult*Aerrs[k])**2+(errmult*As[k])**2)/thick
As[k]=mult*As[k]/thick
else:
if len(intGC)<2:
norm=False
else:
norm=True
if norm:
else:

"""Integrate a sequence

Inputs:
fsn1: range of fsns. The first should be the empty beam
measurement.
fsndc: one FSN or a list of FSNS corresponding to the empty beam
measurements. If sens is a sensitivity dict, this value will
be disregarded.
sens: sensitivity matrix, or a sensitivity dict.
errorsens: error of the sensitivity matrix. If sens is a
sensitivity dict, this value will be disregarded.
orifsn: which element of the sequence should be used for
determining the origin. 0 means empty beam...
Or you can give a tuple or a list of two: in this case these
will be the coordinates of the origin and no auto-searching
will be performed.
energymeas: list of apparent energies, for energy calibration
energycalib: list of theoretical energies, corresponding to the
apparent energies.
distminus: this will be subtracted from the sample-to-detector
distance read from the measurement files, but only for
samples, not for references
detshift: this will be subtracted from the sample-to-detector
distance read from all measurement files, including
references!
orig: helper data for the beam finding procedures. You have
several possibilities:
A) a vector/list/tuple of two: findbeam_sector() will be
tried. In this case this is the initial value of the
beam center
B) a vector/list/tuple of four: xmin,xmax,ymin,ymax:
the borders of the rectangle, around the beam, if a
semitransparent beam-stop was used. In this case
findbeam_semitransparent() will be tried, and the beam
center will be determined for each measurement,
independently (disregarding orifsn).
C) a vector/list/tuple of five: Ntheta,dmin,dmax,bcx,bcy:
findbeam_azimuthal will be used. Ntheta, dmin and dmax
are the respective parameters for azimintpix(), while
bcx and bcy are the x and y coordinates for the origin
at the first guess.
same size as that of the measurement data. In this
case findbeam_gravity() will be used.
transm: you can give this if you know the transmission of the
sample from another measurement. Leave it None to use the
measured transmission.

qs: q-range. Automatically determined from the mask, the energy
and the sample-to-detector distance
ints: intensities corresponding to the q values
errs: calculated absolute errors
Areas: effective areas during integration
As: list of corrected 2D data
Aerrs: list of corrected 2D errors
"""
global _B1config
distancetoreference=_B1config['distancetoreference'] #mm. The references are nearer to the detector than the samples
pixelsize=_B1config['pixelsize']

#finding beamcenter

print "B1integrate: doing energy calibration and correction for reference distance"
for k in range(len(data)):
print "Corrected sample-detector distance for fsn %d (ref. before)." % header[k]['FSN']
print "Corrected sample-detector distance for fsn %d (ref. after)." % header[k]['FSN']
else:
# subtract the background and the dark current, normalise by sensitivity and transmission

qs=[]
ints=[]
errs=[]
Areas=[]
As=[]
Aerrs=[]
print "Integrating data. Press Return after inspecting the images."
for k in range(len(Asub)):
continue
As.append(Asub[k])
Aerrs.append(errAsub[k])
pylab.clf()
#pylab.gcf().show()
pylab.draw()
#now do the integration
print "Now integrating..."
spam=time.time()
qs.append(q)
ints.append(I)
errs.append(e)
Areas.append(A)
print "...done. Integration took %f seconds" % (time.time()-spam)
utils.pause() # we put pause here, so while the user checks the 2d data, the integration is carried out.
pylab.clf()
pylab.subplot(121)
pylab.cla()
#print qs[-1]
#print ints[-1]
valid=np.isfinite(ints[-1])&np.isfinite(errs[-1])
pylab.errorbar(qs[-1][valid],ints[-1][valid],errs[-1][valid])
pylab.axis('tight')
pylab.xlabel('q (1/Angstrom)')
pylab.ylabel('Intensity (arb. units)')
pylab.xscale('log')
pylab.yscale('log')
pylab.subplot(122)
pylab.cla()
pylab.plot(qs[-1],Areas[-1],'.')
pylab.xlabel('q (1/Angstrom)')
pylab.ylabel('Effective area (pixels)')
pylab.gcf().show()
utils.pause()
"""Subtract dark current and empty beam from the measurements and
carry out corrections for detector sensitivity, dead time and beam
flux (monitor counter).

Inputs:
data: list of scattering image matrices
fsndc: FSN for the dark current measurements. Can be a single
integer number or a list of numbers. In the latter case the
DC data are summed up. However, if sens is a sensitivity
dict, this value will be disregarded.
sens: sensitivity matrix, or a sensitivity dict.
senserr: error of the sensitivity matrix. If sens is a
sensitivity dict, this value will be disregarded.
transm: if given, disregard the measured transmission of the
sample.
oldversion: do an old-style error propagation (an approximation
only, not a true and mathematically exact error propagation)
Asub: the corrected matrix
errAsub: the error of the corrected matrix
injectionEB: 'y' if an injection between the empty beam and
sample measurement occured. 'n' otherwise
"""
global _B1config
hackDCsub=1

if type(sens)==type({}): # if sens is a new-type sensitivity dict
md=sens['m0']
td=sens['t0']
D=sens['D0']
S=sens['sens']
senserr=sens['errorsens']
else: # the original way
#sum up darkcurrent measurements, if more.
# summarize transmission, anode, monitor and measurement time data
# for dark current files
S=sens
dmd=np.sqrt(md)
dD=np.sqrt(D)
# initialize the result lists to emptys
Asub=[]
errAsub=[]
injectionEB=[]

for k in range(len(data)): # for each measurement
if header[k]['Title']==_B1config['ebtitle']: # if this is an empty beam measurement, skip.
continue
ebindex=None
for i in range(len(header)): # find empty beam in sequence
ebindex=i
if ebindex is None:
print 'Cannot find all empty beam measurements.\nWhere is the empty FSN %d belonging to FSN %d? Ignoring.'% (header[k]['FSNempty'],header[k]['FSN'])
continue

if oldversion:
# subtract dark current and normalize by sensitivity and transmission (1 in case of empty beam)
# subtract dark current from scattering patterns and normalize by sensitivity and transmission

#two-theta
#angle-dependent corrections:
Ca=absorptionangledependenttth(tth,Ta)
# subtract background, but first check if an injection occurred
K=A2*C0*Ca-Abg*C0
dK=np.sqrt(A2err**2*C0**2*Ca**2+Abgerr**2*C0**2)
else: # new version
A=data[k]
dA=np.sqrt(A)
B=data[ebindex]
dB=np.sqrt(B)
daa=np.sqrt(aa)
dab=np.sqrt(ab)
dma=np.sqrt(ma)
dmb=np.sqrt(mb)
dTa=0 # if the error of the transmission is not zero, please set it here.
if (Ta<=0): # older measurements on the glassy carbon did not save the transmission.
print ""
print ""
print "---------------------------------------------------------------------"
print "VERY BIG, VERY FAT WARNING!!!"
print "The transmission of this sample (",header[k]['Title'],") is nonpositive!"
print "ASSUMING IT TO BE %f." % (_B1config['GCtransmission'])
print "Note that this may foul the calibration into absolute intensity units!"
print "----------------------------------------------------------------------"
print ""
utils.pause()
Ta=_B1config['GCtransmission']
# <anything>x will be the DC corrected version of <anything>
Ax=A-D*ta/td
max=ma-md*ta/td
Bx=B-D*tb/td
mbx=mb-md*tb/td
if hackDCsub:
print "Tampering with DC subtraction!"
Ax[Ax<=0]=np.nanmin(Ax[Ax>0])
Bx[Bx<=0]=np.nanmin(Bx[Bx>0])

#            print "DC corrected counts:"
#            print "ma:",max
#            print "mb:",mbx
#            print "aa:",aax
#            print "ab:",abx
#            print utils.matrixsummary(Ax,"A")
#            print utils.matrixsummary(Bx,"B")
#two-theta for the pixels

#angle-dependent corrections:
Ca,dCa=absorptionangledependenttth(tth,Ta,diffaswell=True)

#auxiliary variables:
P=Ax/(Ta*max*S)*aax/np.nansum(Ax)*C0*Ca
Q=-Bx/(mbx*S)*abx/np.nansum(Bx)*C0
#            print "Auxiliary matrices:"
#            print utils.matrixsummary(P,"P")
#            print utils.matrixsummary(Q,"Q")
# K is the resulting matrix (corrected for dark current,
# lost anode counts (dead time), sensitivity, monitor,
# transmission, background, and various angle-dependent errors)
K=P+Q
#now calculate the different error terms: ET_x is the contribution
# to the error of K from x.
ET_Ta=(-P/Ta+dCa/Ca)**2*dTa**2
ET_ma=(P/max)**2*dma**2
ET_mb=(Q/mbx)**2*dmb**2
ET_md=(ta/td*P/max+tb/td*Q/mbx)**2*dmd**2
ET_aa=(P/aax)**2*daa**2
ET_ab=(Q/abx)**2*dab**2
ET_S=(K/S)**2*senserr**2

ET_A=(P**2/Ax**2-2*P**2/(Ax*np.nansum(Ax)))*dA**2+P**2/np.nansum(Ax)**2*np.nansum(dA**2)
ET_B=(Q**2/Bx**2-2*Q**2/(Bx*np.nansum(Bx)))*dB**2+Q**2/np.nansum(Bx)**2*np.nansum(dB**2)
alpha=(ta/td*P/Ax+tb/td*Q/Bx)
beta=(ta/td*P/np.nansum(Ax)+tb/td*Q/np.nansum(Bx))

#            print utils.matrixsummary(alpha,"alpha_dark")
#            print utils.matrixsummary(beta,"beta_dark")
ET_D=alpha**2*dD**2+beta**2*np.nansum(dD**2)-2*alpha*beta*dD**2

print "error analysis for sample %s" %header[k]['Title']
print utils.matrixsummary(ET_Ta,"Transmission             ")
print utils.matrixsummary(ET_ma,"Monitor (sample)         ")
print utils.matrixsummary(ET_mb,"Monitor (empty beam)     ")
print utils.matrixsummary(ET_md,"Monitor (dark current)   ")
print utils.matrixsummary(ET_aa,"Anode (sample)           ")
print utils.matrixsummary(ET_ab,"Anode (empty beam)       ")
print utils.matrixsummary(ET_S,"Sensitivity              ")
print utils.matrixsummary(ET_A,"Scattering (sample)      ")
print utils.matrixsummary(ET_B,"Scattering (empty beam)  ")
print utils.matrixsummary(ET_D,"Scattering (dark current)")
print utils.matrixsummary(K,"Corrected matrix         ")
print utils.matrixsummary(dK,"Total error              ")
print utils.matrixsummary(dK**2,"Squared total error      ")
print utils.matrixsummary(dK/K,"Relative error           ")
if type(sens)==type({}): # if sens is a new-type sensitivity dict
#correction terms, accounting for the dependence of S and md,ad,D
CT_md=-2*K/S*(ta/t0*P/max+tb/t0*Q/mbx)*1/sens['S1S']* \
(sens['chim']-S*np.nansum(sens['chim']))*dmd**2
CT_D=2*K/(S*sens['S1S'])*(ta/t0*P/Ax+tb/t0*Q/Bx)*(S*np.nansum(sens['beta'])-sens['beta']+(1-S+ta/t0*P/np.nansum(Ax)+tb/t0*Q/np.nansum(Bx))*sens['alpha'])*dD**2-\
2*K/(S*sens['S1S'])*(ta/t0*P/np.nansum(Ax)+tb/t0*Q/np.nansum(Bx))*(S*np.nansum(sens['beta'])-sens['beta'])*np.nansum(dD**2) + \
2*K/(S*sens['S1S'])*(ta/t0*P/np.nansum(Ax)+tb/t0*Q/np.nansum(Bx))*S*np.nansum(sens['alpha']*dD**2)
print "Correction terms:"
print utils.matrixsummary(CT_md,"Monitor (dark current)   ")
print utils.matrixsummary(CT_D,"Scattering (dark current)")
print utils.matrixsummary(dK,"Total error")
print utils.matrixsummmary(dK**2,"Squared total error")

#normalize by beam size
Bx=header[k]['XPixel'] # the empty beam should be measured with the same settings...
K=K/(Bx*By)
dK=dK/(Bx*By)
# now K and dK are ready, either using the old or the new version.

# check if injection occurred between the empty beam and the dark current measurement
print "Possibly an injection between sample and its background:"
print "Current in DORIS at the end of empty beam measurement %.2f." % header[ebindex]['Current2']
print "Current in DORIS at the beginning of sample measurement %.2f." % header[k]['Current1']
injectionEB.append('y')
else:
injectionEB.append('n')
Asub.append(K)
errAsub.append(dK)
"""Carry out dead-time corrections (anode counts/sum) on 2D
scattering data, subtract dark current, normalize by monitor
counts and sensitivity, and finally by pixel size.

Inputs:
data: list of 2d scattering images, to be added
datadc: list of 2d scattering images for dark current
sens: sensitivity matrix
senserr: error matrix of the sensitivity data

Outputs:
the normalized data and its error matrix
"""
# summarize transmission, anode, monitor and measurement time data
# for measurement files
if type(data)!=type([]):
raise TypeError("Parameter <data> should be a list!")
raise TypeError("Parameter <header> should be a list!")
raise TypeError("Parameter <datadc> should be a list!")
raise TypeError("Parameter <headerdc> should be a list!")

if transm is None:
transmave=transm.mean() # average transmission
transmerr=transm.std() # standard deviation of the transmission
else:
transmave=transm
transmerr=0
if transmave<=0:
print ""
print ""
print "---------------------------------------------------------------------"
print "VERY BIG, VERY FAT WARNING!!!"
print "The transmission of this sample (",header[0]['Title'],") is nonpositive!"
print "ASSUMING IT TO BE 0.5."
print "Note that this may foul the calibration into absolute intensity units!"
print ""
print ""
utils.pause()
transmave=0.5
h['Transm']=transmave

print "Average transmission %.4f +/- %.4f" % (transmave,transmerr)

# summarize transmission, anode, monitor and measurement time data
# for dark current files

A1=sum(data) # do not use np.sum()

# correct monitor counts with its dark current
mo2=mo1-modc*meastime1/meastimedc
an2=an1-andc*meastime1/meastimedc

mo2err=np.sqrt(mo1+modc*(meastime1/meastimedc)**2)
an2err=np.sqrt(an1+andc*(meastime1/meastimedc)**2)

sumA2=A2.sum()
# error of sumA2, not sum of error of A2.
sumA2err=np.sqrt((A2err**2).sum())

# summarized scattering pattern, subtracted the dark current,
# normalized by the monitor counter and the sensitivity
A3=A2/(sens*mo2)*an2/sumA2/transmave

print "Sum/Total of dark current: %.2f %. Counts/s %.1f." % (100*Adc.sum()/andc,andc/meastimedc)
print "Sum/Total before dark current correction: %.2f %. Counts on anode %.1f cps. Monitor %.1f cps." %(100*A.sum()/an1,an1/meastime1,mo2/meastime1)
print "Sum/Total after dark current correction: %.2f %." % (100*sumA2/anA2)

A3err=np.sqrt((A2err/(sens*mo2*transmave)*an2/sumA2)**2+ \
(senserr*A2/(sens**2*mo2*transmave)*an2/sumA2)**2+ \
(mo2err*A2/(sens*transmave*mo2**2)*an2/sumA2)**2+ \
(an2err*A2/(sens*mo2*transmave)/sumA2)**2+ \
(sumA2err*A2/(sens*mo2*transmave)*an2/sumA2**2)**2+ \
(transmerr*A2/(sens*mo2*transmave**2)*an2/sumA2)**2)

#normalize by beam size