1
2 """
3 Correct from counts per pixel to counts per unit area on the detector.
4
5 Usage
6 =====
7
8 Example::
9
10 # construct the correction based on a detector measurement
11 from reflectometry.reduction import load, measured_area_correction
12 floodfill = load('detector measurement')
13 areacor = measured_area_correction(floodfill)
14
15 # apply the correction
16 data = load(d).apply(areacor)
17
18 If you want to the pixels to have equal area then you need to include
19 a rebin=True flag::
20
21 areacor = measured_area_correction(floodfill,rebin=True)
22
23 You can also set the areacor.rebin attribute directly.
24
25 If you know the pixel widths in x and y then you don't have to extract
26 them from a measurement, but can instead use::
27
28 areacor = area_correction(wx,wy,source="provenance")
29
30 Certain instruments have predefined area corrections, which are
31 returned by the area_correction() method of the refldata object.
32
33 Limitations
34 ===========
35
36 The correction applies independently in x and y using the sum of the pixels
37 across that detector dimension. Small detector rotations and misaligned
38 wires will lead to artificial broadening.
39
40 The correction does not apply if there are other detector artifacts, such as
41 pixel efficiency variation. Sliding the detector across a narrow slit and
42 summing the total counts across the detector will show any efficiency
43 variation. We do not attempt to correct for this.
44
45 This theory also assumes the floodfill is uniform, which will depend on
46 how the measurement is performed. For example, incoherent scattering
47 off water will yield higher intensities closer to the beam, and measurements
48 of a point source with a flat detector will yield higher counts where the
49 detector is nearer the point. Data not taken in the ideal condition of
50 sliding the detector behind a narrow slit while maintaining uniform intensity
51 from the beam will need to scaled appropriately before estimating pixel widths.
52
53 Even after correcting for pixel widths in x and y, there may be a residual
54 signal on the detector. This can be viewed by applying the detector
55 correction to the floodfill itself. As a second order pixel area effect due
56 to variations in wire position across the rows and columns of the detector,
57 the pixels could be further scaled by this variation. The scale factor
58 should be S = nx*ny * F/sum(F), where F is the corrected detector image.
59 This is not implemented here.
60
61 Theory
62 ======
63
64 Assuming the pixels have none uniform area, we can estimate the area
65 of each pixel by looking at the counts across the detector in a floodfill
66 environment.
67
68 First lets define the detector. We will do so in one dimension, but
69 it is a separable problem which can be applied independently in the
70 x and y directions.
71
72 Detector has varying width pixels:
73
74 |......|....|.....|....|
75 w_1 w_2 w_n
76
77 L = sum(w) = length of the detector
78 C_i = measured counts in cell i
79 D = sum(C)/L = average flux on the detector
80
81 Given a uniform flux on the detector, what is the probability of seeing
82 the particular measured set of counts P(C|w). We will assume a Poisson
83 distribution for the counts in the individual pixels.
84
85 D = flux in counts/unit area
86 lambda_i = D w_i = expected counts in cell i
87 P_i(C_i|D w_i) = exp(-D w_i) (D w_i)**C_i / C_i!
88 P(C|w) = prod P_i(C_i|D w_i)
89
90 Minimize the log probability by setting the derivative to zero:
91
92 log P = sum log P_i
93 log P_i = -D w_i + C_i log D w_i - log C_i!
94
95 d/dw_i log P = d/dw_i log P_i
96 = d/dw_i (-D w_i) + C_i d/dw_i log D w_i
97 = -D + C_i D / (D w_i)
98 d/dw_i log P = 0
99 => D = C_i/w_i
100 => w_i = C_i/D = L * C_i / C
101 """
102
103 import numpy
104
105 from reflectometry.reduction.rebin import rebin2d
106
108 """
109 Given a detector measurement (either a flood fill or a main
110 beam scan) return an area correction based on the pixel widths
111 estimated from the data.
112
113 The resulting correction will scale the pixels to counts per
114 unit area, so must be applied all datasets for the absolute
115 reflection amplitude to be calculated.
116
117 If rebin is true, then instead of scaling pixels to counts
118 per unit area, the pixel boundaries will be adjusted so that
119 they have equal area using standard rebinning techniques.
120
121 For main beam scans, the data will need to be normalized to
122 the same number of incident neutrons on all frames, and all
123 frames summed together. Time of flight data will need to be
124 combined into a single file, with the time channels collapsed
125 to a single channel containing all the wavelengths of interest
126 for each measurement, and the resulting frames normalized and
127 summed. In any case this preprocessing should happen once
128 for each beam scan, resulting in a small file that can be applied
129 to multiple datasets.
130 """
131
132 Cxy = numpy.sum(data.detector.counts,axis=0)
133
134 C = numpy.sum(Cxy)
135
136 Cx = numpy.sum(Cxy,axis=0)
137 Cy = numpy.sum(Cxy,axis=1)
138
139 Lx,Ly = self.detector.size
140 wx = Cx / C * Lx
141 wy = Cy / C * Ly
142
143 return AreaCorrection(wx,wy,rebin=rebin,source=data.name)
144
146 """
147 Convert detector counts from counts per pixel to counts per unit area.
148 """
149 properties = ["wx","wy","source","rebin"]
150
151 - def __init__(self, wx, wy, rebin=False, source="unknown"):
152 """
153 Create a pixel area correction function.
154
155 wx,wy is the solid angle of the pixels as measured on
156 the detector. This function will normalize the counts
157 on the detector by pixel area. This will change the
158 pixel widths in the data file.
159
160 If rebin is True then adjust pixel boundaries so the
161 pixels have equal area.
162
163 Source is a string to report in the log as the origin
164 of the correction data.
165 """
166 self.wx = numpy.array([wx])
167 self.wy = numpy.array([wy]).T
168 self.source = source
169 self.rebin = rebin
170
172 """Apply the area correction to the data"""
173 nx,ny = self.wx.shape[0],self.wy.shape[1]
174 assert data.detector.dims == [nx,ny], \
175 "area correction size does not match detector size"
176 if self.rebin:
177
178 x = numpy.concatenate([(0.),numpy.cumsum(wx)])
179 y = numpy.concatenate([(0.),numpy.cumsum(wy)])
180 Lx,Ly = numpy.sum(wx),numpy.sum(wy)
181 xo = numpy.linspace(0,Lx,nx+1)
182 yo = numpy.linspace(0,Ly,ny+1)
183 Io = numpy.empty((nx,ny),'d')
184
185
186 for i in xrange(data.detector.counts.shape[0]):
187 frame = data.detector.counts[i]
188 frame[:] = rebin2d(x,y,frame,xo,yo,Io)
189
190
191
192
193 fail.data.detector.widths_x = numpy.zeros(nx,'d')+Lx/nx
194 frame.data.detector.widths_y = numpy.zeros(ny,'d')+Ly/ny
195 else:
196
197 data.detector.widths_x = wx[0,:]
198 data.detector.widths_y = wy[:,0]
199
200
201 fail.data.detector.counts /= self.wx
202 data.detector.counts /= self.wy
203
205 return "AreaCorrection('%s')"%self.source
206