1 """
2 Polarized data corrections
3
4 Polarized neutron measurements select a particular neutron spin state
5 for the incident and reflected neutrons. Experiments which probe the
6 interaction between spin state and sample must first determine the
7 efficiency with which each spin state can be selected and apply that
8 effect to the data before presenting it to the user.[1]
9
10 The way that spin states are selected will vary between instruments.
11 Some instruments will have a polarizing supermirror+flipper
12 arrangement, with the mirror either transmitting or reflecting
13 a particular polarization state, and a magnetic field tuned to
14 flip the polarization state when a current is applied or otherwise
15 leave it in the selected state. TOF sources will require a
16 time-varying field to adjust for the different transit times
17 for different wavelengths through the flipper field. He3 polarizers
18 can select spin up or spin down for transmission and so do not
19 require a flipper. In any case, the efficiency will vary with
20 Q, either due to increased divergence or wavelength dependence.
21
22
23 1. Prepare the intensity data
24 ===============================
25
26 To perform a polarization data reduction you must first have measured
27 beam data under different polarization conditions::
28
29 pp spin up incident, spin up reflected
30 pm spin up incident, spin down measured
31 mp spin down incident, spin up measured
32 mm spin down incident, spin down measured
33
34 These data are passed into the correction using a PolarizedData container::
35
36 I = PolarizedData()
37 I.xlabel, P.xunits = 'Slit 1', 'mm' # On TOF, maybe 'wavelength','nm'
38 I.ylabel, P.yunits = 'Intensity', 'counts'
39 I.pp.set(x=Spp, y=Ipp, variance=Ipp)
40 I.pm.set(x=Spm, y=Ipm, variance=Ipm)
41 I.mp.set(x=Smp, y=Imp, variance=Imp)
42 I.mm.set(x=Smm, y=Imm, variance=Imm)
43
44 The intensities given should be 1-dimensional. For 1D and 2D detectors, the
45 integrated beam image should be used. Time-of-flight data will be
46 organized by time channel while monochromatic data is organized by slit
47 opening. Any ordering will do so long as the intensity is well behaved
48 between points.
49
50 We will assume the correction is uniform across the detector since it
51 is impractical to measure the efficiency as a function of position.
52 In fact the slightly different path lengths for the various positions
53 on the detector will lead to a decrease in efficiency as you move away
54 from the center of the beam, and so it may have a subtle influence on
55 off-specular polarized reflectometry.
56
57 The different cross sections must have corresponding ordinate axes. On TOF
58 machines this will usually not be a problem since all time channels are
59 measured simultaneously. On scanning instruments, there are a variety
60 of factors which may result in points that don't correspond well, such
61 as user error, or undersampling in the spin-flip cross sections.
62
63 Even if the data are already aligned across the four cross sections,
64 it is still desirable to smooth the polarized data prior to estimating
65 because the efficiency estimate is unstable.
66
67 Use alignment or smoothing corrections to make your data regular::
68
69 I.apply(align) # use linear interpolation to regrid the data
70 I.apply(smooth) # a weighted irregular Savitsky-Golay filter
71
72
73 2. Compute polarization efficiencies
74 ====================================
75
76 The polarization efficiency correction is controlled by the
77 PolarizationEfficiency class::
78
79 eff = PolarizationEfficiency(I)
80
81 Once the class is constructed and the resulting efficiency attributes
82 are available::
83
84 eff.fp front polarizer efficiency
85 eff.ff front flipper efficiency
86 eff.rp rear polarizer efficiency
87 eff.rf rear flipper efficiency
88 eff.Ic overall beam intensity
89
90 The formalism from Majkrzak, et al. (see attached PDF) defines the
91 following, which are attributes to the efficiency object::
92
93 eff.F = fp
94 eff.R = rp
95 eff.x = 1 - 2*ff
96 eff.y = 1 - 2*rf
97 eff.beta = Ic/2
98
99 There are a few attributes of the class that will change how the
100 efficiencies are calculated::
101
102 eff.min_intensity = 1e-2
103 eff.min_efficiency = 0.7
104 eff.FRbalance = 0.5
105
106 FRbalance determines the relative front-back weighting of the
107 polarization efficiency. From the data we can only estimate the
108 product FR of the front and rear polarization efficiencies, and
109 the user must decide how to distribute the remainder.
110 The FRbalance should vary between 0 (front polarizer is 100% efficient)
111 through 0.5 (distribute inefficiency equally) to 1 (rear polarizer
112 is 100% efficient). The particular formula used is::
113
114 F = (F*R)^FRbalance
115 R = (F*R)/F
116
117 The min_intensity and min_efficiency numbers are used to bound
118 the range of the correction. Normally you won't need to set
119 them. These are class attributes, so you can set them globally
120 in PolarizationEfficiency instead of eff, but still override
121 them for particular instances by setting them in eff.
122
123
124 3. Apply the correction to the data
125 ===================================
126
127 You must first prepare your data for the polarization correction by
128 placing it in a PolarizedData container, and possibly subtracting
129 the background measurements (the polarization correction is linear so it
130 doesn't matter mathematically if background subtraction happens before
131 or after the polarization correction, but the polarization correction
132 is slow enough that it should probably only be done once).
133
134 Once again the measurements need to be aligned, but in this case
135 smoothing should not be used. The data will need to be normalized
136 to the same monitor as the intensity measurement before adding them
137 to the container.
138
139 For scanning instruments, the data will need to retain the slit
140 settings used for the measurement so that the correct intensity
141 measurement can be used for the correction.
142
143
144 TODO: create a complete data reduction doc including footprint, etc.
145 TODO: use array masks to identify interpolated data.
146 TODO: incorporate time for He3 polarizer
147 TODO: extend to 1D and 2D detectors
148 TODO: cross check results against Asterix algorithm
149 TODO: implement alignment and smoothing
150 TODO: consider applying the efficiency correction to the theory rather than
151 the data --- if the inversion is unstable, this may be more reliable.
152 TODO: with both spin up and spin down neutrons in sample simultaneously
153 in some proportion, do the cross sections need to be coherently to account
154 for the theory?
155
156 [1] C.F. Majkrzak (1996). Physica B 221, 342-356.
157 """
158
159 import numpy as N
160 from reflectometry.reduction.correction import Correction
161 from reflectometry.reduction.data import PolarizedData
162 from reflectometry.reduction.wsolve import wsolve
166 """
167 Polarization efficiency correction object. Create a correction object
168 from a polarized direct beam measurement and apply it to measured data.
169
170 E.g.,
171 eff = PolarizationEfficiency(beam)
172 data.apply(eff)
173
174 """
175
176 properties = ['FRbalance','min_efficiency','min_intensity',
177 'clip','spinflip']
178 min_efficiency = 0.7
179 """Minimum efficiency cutoff"""
180 min_intensity = 1e-2
181 """Minimum intensity cutoff"""
182 clip = True
183 """Perform clipping to keep efficiencies physical"""
184 spinflip = True
185 """Correct both spinflip and non-spinflip data if available"""
186 _beam = PolarizedData()
187 _FRbalance = 0.5
188
189
192 if not beam.isaligned():
193 raise ValueError, "polarization efficiency correction needs aligned intensity measurements"
194 self._beam = beam
195 self.__update()
196 beam = property(_getbeam,_setbeam,
197 doc="measured beam intensity")
198
201 if val > 1: val = 1
202 elif val < 0: val = 0
203 self._FRbalance = val
204 self.__update()
205 FRbalance = property(_getFRbalance,_setFRbalance,
206 doc="relative balance of front to back efficiency")
207
208 @property
209 - def x(self): return (1-2*self.ff)
210 @property
211 - def y(self): return (1-2*self.rf)
212 @property
213 - def F(self): return self.fp
214 @property
215 - def R(self): return self.rp
216 @property
217 - def beta(self): return 0.5*self.Ic
218
220 """
221 Define the polarization efficiency correction for the beam
222
223 Keywords:
224 beam: measured beam intensity for all four cross sections
225 FRbalance: portion of efficiency to assign to front versus rear
226 clip: [True|False] force efficiencies below 100%
227 spinflip: [True|False] compute spinflip correction
228 """
229 self.__dirty = False
230 self.set(**kw)
231
232 - def set(self, **kw):
233 """
234 Update a set of correction parameters.
235 """
236 self.__initializing = True
237 for k,v in kw.iteritems():
238 assert hasattr(self,k), "No %s in %s"%(k,self.__class__.__name__)
239 setattr(self,k,v)
240 self.__initializing = False
241 if self.__dirty: self.__update()
242
243
250
252 return "PolarizationEfficiency('%s')"%self.beam.name
253
255 """
256 Compute polarizer and flipper efficiencies from the intensity data.
257
258 If clip is true, reject points above or below particular efficiencies.
259 The minimum intensity is 1e-10. The minimum efficiency is 0.9.
260
261 The computed values are systematically related to the efficiencies:
262 Ic: intensity is 2*beta
263 fp: front polarizer efficiency is F
264 rp: rear polarizer efficiency is R
265 ff: front flipper efficiency is (1-x)/2
266 rf: rear flipper efficiency is (1-y)/2
267 reject is the indices of points which are clipped because they
268 are below the minimum efficiency or intensity.
269
270 See PolarizationEfficiency.pdf for details on the calculation.
271 """
272 if self.__initializing:
273 self.__dirty = True
274 return
275 else:
276 self.__dirty = False
277
278
279 beam = self._beam
280 assert beam.isaligned(), "need aligned data"
281 pp,pm,mp,mm = beam.pp.v,beam.pm.v,beam.mp.v,beam.mm.v
282 Ic = (pp*mm-pm*mp) / (pp+mm-pm-mp)
283 reject = (Ic!=Ic)
284 if self.clip: reject |= clip(Ic, self.min_intensity, N.inf)
285
286
287
288
289
290 FR = pp/Ic - 1
291 if self.clip: reject |= clip(FR, self.min_efficiency**2, 1)
292 fp = FR ** self.FRbalance
293 rp = FR / fp
294
295
296
297
298
299
300
301 x = (pm/Ic - 1)/FR
302 y = (mp/Ic - 1)/FR
303 ff = (1-x)/2
304 rf = (1-y)/2
305
306 if self.clip: reject |= clip(ff, self.min_efficiency, 1)
307 if self.clip: reject |= clip(rf, self.min_efficiency, 1)
308
309 self.Ic, self.fp, self.ff, self.rp, self.rf = Ic,fp,ff,rp,rf
310 self.reject = reject
311
312 -def clip(field,lo,hi,nanval=0.):
313 """clip the values to the range, returning the indices of the values
314 which were clipped. Note that this modifies field in place. NaN
315 values are clipped to the nanval default.
316 """
317 idx = N.isnan(field); field[idx] = nanval; reject = idx
318 idx = field<lo; field[idx] = lo; reject |= idx
319 idx = field>hi; field[idx] = hi; reject |= idx
320 return reject
321
323 """Apply the efficiency correction in eff to the data."""
324
325 F = eff.F
326 R = eff.R
327 Fx = eff.F*eff.x
328 Ry = eff.R*eff.y
329 beta = eff.beta
330
331
332 if spinflip and data.pm and data.mp:
333 spinflip = False
334
335 if spinflip:
336 Y = N.vstack([ data.pp.v, data.pm.v, data.mp.v, data.mm.v ]) / beta
337 dY = N.vstack([ data.pp.dv, data.pm.dv, data.mp.dv, data.mm.dv ]) / beta
338
339
340
341
342
343 H = N.array([
344 [(1+F)*(1+R), (1+Fx)*(1+R), (1+F)*(1+Ry), (1+Fx)*(1+Ry)],
345 [(1-F)*(1+R), (1-Fx)*(1+R), (1-F)*(1+Ry), (1-Fx)*(1+Ry)],
346 [(1+F)*(1-R), (1+Fx)*(1-R), (1+F)*(1-Ry), (1+Fx)*(1-Ry)],
347 [(1-F)*(1-R), (1-Fx)*(1-R), (1-F)*(1-Ry), (1-Fx)*(1-Ry)],
348 ])
349 else:
350 Y = N.array([ data.pp.v, data.mm.v ]) / beta
351 dY = N.array([ data.pp.dv, data.mm.dv ]) / beta
352 H = N.array([
353 [(1+F)*(1+R), (1+Fx)*(1+Ry)],
354 [(1-F)*(1-R), (1-Fx)*(1-Ry)],
355 ])
356
357
358
359
360
361
362
363 X = N.zeros(Y.shape)
364 dX = N.zeros(dY.shape)
365 reject = eff.reject&True
366 for i in xrange(H.shape[1]):
367 A = H[:,:,i]
368 y = Y[:,i]
369 dy = dY[:,i]
370 try:
371 s = wsolve(A,y,dy)
372 except ValueError:
373 reject[i] = True
374 x,dx = y,dy
375 else:
376 x,dx = s.x,s.std
377 X[:,i] = x
378 dX[:,i] = dx
379
380 if spinflip:
381 data.pp.v, data.pp.dv = X[0,:], dX[0,:]
382 data.pm.v, data.pm.dv = X[1,:], dX[1,:]
383 data.mp.v, data.mp.dv = X[2,:], dX[2,:]
384 data.mm.v, data.mm.dv = X[3,:], dX[3,:]
385 else:
386 data.pp.v, data.pp.dv = X[0,:], dX[0,:]
387 data.mm.v, data.mm.dv = X[1,:], dX[1,:]
388 if hasattr(data.reject):
389 data.reject |= reject
390 else:
391 data.reject = reject
392
400
401 if __name__ == "__main__": demo()
402