data reduction / profile fitting on area detector
Kay Diederichs
dikay at sun1.ruf.uni-freiburg.de
Fri Jul 10 06:56:09 EST 1992
Dear colleagues,
I wrote a little program to test the usefulness of profile fitting in
area detector data reduction. Profile fitting as opposed to straight
summation is supposed to reduce the variance in the data and thus to
decrease r-factors. The effect should be strongest for reflections
with low signal/noise ratio, and negligible for strong reflections.
To my knowledge, profile fitting is used in MOSFLM, DENZO and PROCESS
(that is the program that reduces R-AXIS II frames) but not in XDS.
( I don't know about OSC, MADNES, BUDDHA ...)
One question that I wanted to answer to myself is whether profile fitting
could compensate for the higher background of image plate detectors (with
oscillation ranges around 1 degree) when compared to gas-filled multiwire
detectors (oscillation range around 1/5 degree).
To summarize, here are a few numbers I get (using thresh=0.02,width=1.5
and always under the assumption that the average background and reflection
profiles are known exactly):
No offset of profile box -
R-int (%) / av. Intensity
countrate theoretical Straight profile-fitted
background peak Intensity summation non-weighted proper weights
50 10 71. 45.9 / 68. 34.6 / 69. 34.6 / 69.
50 50 353. 9.6 / 352. 7.9 / 351. 7.8 / 352.
50 250 1767. 2.6 / 1769. 2.3 / 1770. 2.3 / 1769.
50 1000 7068. 1.1 / 7068. 1.1 / 7060. 1.1 / 7063.
250 10 71. 91.2 / 67. 74.7 / 66. 74.8 / 66.
250 50 353. 19.3 / 349. 14.0 / 349. 14.0 / 349.
250 250 1767. 4.2 / 1763. 3.5 / 1764. 3.5 / 1764.
250 1000 7068. 1.3 / 7066. 1.3 / 7066. 1.2 / 7066.
If we compare the numbers, we see
1) properly weighted profile fitting is always better than not weighted
2) the internal R factors look a good deal better for the profile fitted
intensities, but profile fitting cannot make up for reduced background.
Now, if I offset the assumed-known profile by 1 pixel (which probably
occurs for many reflections, even if on the average the displacement is less):
50 10 71. 50.0 / 64. 40.3 / 54. 40.4 / 55.
50 50 353. 10.7 / 332. 9.5 / 280. 9.5 / 289.
50 250 1767. 2.6 / 1666. 2.9 / 1415. 2.6 / 1535.
50 1000 7068. 1.1 / 6654. 1.3 / 5656. 1.1 / 6471.
250 10 71. 94.3 / 65. 84.4 / 57. 84.4 / 57.
250 50 353. 19.7 / 331. 18.1 / 281. 18.1 / 283.
250 250 1767. 4.4 / 1662. 4.2 / 1411. 4.1 / 1452.
250 1000 7068. 1.3 / 6656. 1.5 / 5656. 1.4 / 6080.
3) the summed intensities are much more robust against displacement of
the profile box. The larger error in the profile-fitted intensities does
not show up in the R-factors. The intensities obtained by straight
summation are about 6 - 10% lower, those of weighted profile fitting
14 - 23 % lower than the correct values. If the decrease in intensities
of reflections with displaced profile boxes are translated into R-scale
values, the profile fitted intensities perform much worse than the summed
ones (my program prints out the values, but they are not in the table).
So, is profile fitting really worth the extra effort? Could there be a
weighting function that is as robust against displacement as straight
summation and yet reduces the variance?
Here is the program; try yourself. I would like to get your comments and
suggestions.
Kay
C cut here, it should compile and run on any computer
parameter (ng=4,ng2=2*ng+1,nprof=1000)
integer ir(-ng:ng,-ng:ng)
real prof(-ng:ng,-ng:ng),vali(nprof),valp(nprof),valp2(nprof)
real profn(-ng:ng,-ng:ng)
logical first
common iseed
1001 format(11i7)
1002 format(35x,a)
1 print*,'bkg per pix, peak per pix, threshold, width,',
& ' offset (e.g. 10,100,.02,1.5,0.) ?'
read(5,*,end=99) ibkg,hei,thresh,del,disali
iseed=1234567
first=.true.
sump=0.
sumq=0
sumq2=0
do 100 i=-ng,ng
do 100 j=-ng,ng
g=exp(-(i**2+j**2)/del**2)
prof(j,i)=g
sump=sump+g
g=exp(-((i+disali)**2+j**2)/del**2)
sumq=sumq+g
profn(j,i)=g
100 continue
C normalize profile
do 110 i=-ng,ng
do 110 j=-ng,ng
profn(j,i)=profn(j,i)/sumq
sumq2=sumq2+profn(j,i)**2
110 continue
C generate profiles
do 150 k=1,nprof
do 130 i=-ng,ng
do 130 j=-ng,ng
call rnpoi(1,ibkg+hei*prof(j,i),ir(j,i))
130 continue
if (first) then
first=.false.
do 200 i=-ng,ng
200 print 1001,(ir(j,i),j=-ng,ng)
print*,'weighting function used in profile fitting:'
rmax=0.
do 250 i=-ng,ng
do 250 j=-ng,ng
250 rmax=max(rmax,profn(j,i)*(hei*profn(j,i))/
& (profn(j,i)*sump*hei+ibkg))
do 260 i=-ng,ng
260 print 1001,(nint(100.*profn(j,i)*(hei*profn(j,i))/
& (profn(j,i)*sump*hei+ibkg)/rmax),j=-ng,ng)
endif
call integ(ir,profn,ng2,val,ibkg,thresh/sumq)
vali(k)=val
call profi(ir,profn,ng2,val,ibkg)
valp(k)=val/sumq2
call profi2(ir,profn,ng2,val,ibkg)
valp2(k)=val
150 continue
sump=sump*hei
print*,'theoretical intensity,sigma:',sump,sqrt(sump+ng**2*ibkg)
print 1002,'mean I s.dev r-scale(%) r-int(%)'
call eval(vali,nprof,sump, ' integrated intensities: ')
call eval(valp,nprof,sump, ' profile-fitted intensities: ')
call eval(valp2,nprof,sump,' profile-fitted intensities 2: ')
goto 1
99 stop
end
C
subroutine integ(ir,prof,ng2,val,ibkg,thresh)
C Summation of counts is done for profile points which exceed a given threshold
C The threshold should probably be a few percent of the maximum of the profile.
C The XDS program uses thresh=0.02
integer ir(ng2,ng2)
real prof(ng2,ng2)
val=0.
sum1=0.
sum2=0.
do 100 i=1,ng2
do 100 j=1,ng2
if (prof(j,i).gt.thresh) then
val=val+ir(j,i)-ibkg
sum1=sum1+prof(j,i)
else
sum2=sum2+prof(j,i)
endif
100 continue
C scale up intensity slightly to account for missing profile part
val=val*(sum1+sum2)/sum1
return
end
C
subroutine profi(ir,prof,ng2,val,ibkg)
C profile fitting without using the variance
integer ir(ng2,ng2)
real prof(ng2,ng2)
val=0.
do 100 i=1,ng2
do 100 j=1,ng2
100 val=val+(ir(j,i)-ibkg)*prof(j,i)
return
end
C
subroutine profi2(ir,prof,ng2,val,ibkg)
C correct profile fitting
integer ir(ng2,ng2)
real prof(ng2,ng2)
sumz=0.
sumn=0.
C first get estimate of intensity (needed for variance)
val=0.
do 50 i=1,ng2
do 50 j=1,ng2
if (prof(j,i).gt.0.)val=val+ir(j,i)-ibkg
50 continue
do 100 i=1,ng2
do 100 j=1,ng2
if (prof(j,i).le.0.) goto 100
sumz=sumz+prof(j,i)*(ir(j,i)-ibkg)/(prof(j,i)*val+ibkg)
sumn=sumn+prof(j,i)**2/(prof(j,i)*val+ibkg)
100 continue
val=sumz/sumn
C now we have a good estimate of intensity
sumz=0.
sumn=0.
do 200 i=1,ng2
do 200 j=1,ng2
if (prof(j,i).le.0.) goto 200
sumz=sumz+prof(j,i)*(ir(j,i)-ibkg)/(prof(j,i)*val+ibkg)
sumn=sumn+prof(j,i)**2/(prof(j,i)*val+ibkg)
200 continue
val=sumz/sumn
return
end
C
subroutine eval(val,np,theo,string)
C compute r-factors. NB r-scale is actually computed without scaling!
real val(np)
character*(*) string
1001 format(a,4f10.2)
sum=0.
sumq=0.
sumz=0.
sumn=0.
do 100 i=1,np
sumz=sumz+abs(val(i)-theo)
sumn=sumn+theo
sum=sum+val(i)
sumq=sumq+val(i)**2
100 continue
sumzi=0.
sumni=0.
do 200 i=1,np
sumzi=sumzi+abs(val(i)-sum/np)
sumni=sumni+val(i)
200 continue
print 1001,string,sum/np,sqrt(max(0.,sumq/np-sum**2/np**2)),
& 100*sumz/sumn,100*sumzi/sumni
return
end
C-----------------------------------------------------------------------
C IMSL Name: RNPOI (Single precision version)
C
C Computer: vax/SINGLE
C
C Revised: August 20, 1985
C
C Purpose: Generate pseudorandom numbers from a Poisson
C distribution.
C
C
More information about the Xtal-log
mailing list