Source Code for Module reflectometry.reduction.ratiocor

  1  # This program is public domain 
  2  """ 
  3  RatioIntensity estimates an intensity scan based on the measured 
  4  intensity from a known scatterer. 
  5   
  6  You need to perform this type of scan on reflectometers whose beam 
  7  profile changes as a function of angle, thus making it impossible 
  8  to measure the beam intensity directly. 
  9   
 10  Even with the ratio scan the reduction may need a scaling 
 11  step to recover the absolute intensity. 
 12   
 13  Note that we will assume that slit opening as a function of 
 14  angle is identical for the measured intensity and the measured 
 15  reflection. 
 16   
 17  We will also assume that the intensity scan and the measured 
 18  reflection are performed on blanks with identical geometery 
 19  so the footprint is the same on both.  No footprint correction 
 20  will be required. 
 21   
 22  Example:: 
 23   
 24      # Compute the intensity scan assuming that data measures 
 25      # the specular reflection off 85% water, 15% deuterated water. 
 26      data.apply(WaterIntensity(D2O=15)) 
 27   
 28  """ 
 29  from __future__ import with_statement 
 30   
 31  from math import * 
 32  import reflectometry.model1d as reflfit 
 33  from reflectometry.reduction.correction import Correction 
 34  from reflectometry.reduction.fresnel import Fresnel 
 35  import numpy 
 36   
 37  __all__ = ['RatioIntensity', 'WaterIntensity'] 
 38   
 39 -class RatioIntensity(Correction): 
 40      """ 
 41      Compute the incident intensity assuming that data measures a 
 42      reflection off known sample. 
 43   
 44      model = reflectometry.model1d.model 
 45   
 46      Returns a correction object that can be applied to data. 
 47      """ 
 48      model = None 
 49   
 50 -    def __init__(self, model): 
 51          self.model = model 
 52   
 53 -    def apply(self, data): 
 54          intensity_ratio(data=data,model=self.model) 
 55   
 56 -    def __str__(self): 
 57          return "RatioIntensity(%s)"%str(model) 
 58   
 59 -class WaterIntensity(Correction): 
 60      """ 
 61      Compute the incident intensity assuming that data measures a 
 62      reflection off water. 
 63   
 64      D2O = percentage is the portion that is deuterated, from 0 to 100. 
 65      probe = 'neutron' or 'xray' 
 66   
 67      Returns a correction object that can be applied to data. 
 68      """ 
 69 -    def _set_fresnel_coefficients(self): 
 70          # SLDs calculated with NBCU.XLS from the NCNR website. 
 71          if self.probe == 'neutron': 
 72              rho_H2O,mu_H2O = -0.560, 0.001 
 73              rho_D2O,mu_D2O = 5.756, 0.00 
 74          elif self.probe == 'xray': 
 75              rho_H2O,mu_H2O = 9.466, 0.098 
 76              rho_D2O,mu_D2O = 8.516, 0.089 
 77          else: 
 78              raise ValueError, "water_ratio needs probe 'neutron' or 'xray'" 
 79          D2O = self._D2O 
 80          self.model.rho = 1e-6*((100-D2O)*(rho_H2O) + D2O*rho_H2O)/100. 
 81          self.model.mu = 1e-6*((100-D2O)*(mu_H2O) + D2O*mu_H2O)/100. 
 82   
 83 -    def _setD2O(self, D2O): 
 84          self._D2O = D2O 
 85          self._set_fresnel_coefficients 
 86 -    def _getD2O(self): 
 87          return self._D2O 
 88 -    def _setprobe(self,probe): 
 89          self._probe = probe 
 90          self._set_fresnel_coefficients 
 91 -    def _getprobe(self): 
 92          return self._probe 
 93      D2O = property(_getD2O,_setD2O) 
 94      probe = property(_getprobe,_setprobe) 
 95   
 96 -    def __init__(self,D2O=0,probe=None,background=0): 
 97          self.model = Fresnel(sigma=3) 
 98          self._probe = probe 
 99          self._D2O = D2O 
100          self._set_fresnel_coefficients() 
101   
102 -    def apply(self,data): 
103          intensity_ratio(data=data,model=self.model) 
104   
105 -    def __str__(self): 
106          if self._D2O == 0: 
107              return "WaterIntensity(probe='%s')"%self._probe 
108          else: 
109              return "WaterIntensity(D2O=%d,probe='%s')"%(self._D2O,self._probe) 
110   
111 -def intensity_ratio(data=None,model=None): 
112      """ 
113      Compute the incident intensity assuming that data measures a 
114      reflection off the given model. 
115   
116      Returns a copy of the dataset that can be subsequently written 
117      to a data file for use in reflpak or other reduction software. 
118      """ 
119   
120      # TODO for now we are not using resolution information when 
121      # computing the data/theory ratio 
122      Q,L = data.Qz,data.detector.wavelength 
123      R = model.reflectivity(Q,L) 
124      data.R, data.dR = data.R/R, data.dR/R 
125      data.intent = 'intensity' 
126      return 
127   
128      # Compute the convolved reflectivity. 
129      # Find the theory points necessary to compute an accurate estimate of 
130      # the convolved model without over- or under-sampling. 
131      Q,dQ,L = data.Qz,data.dQ,data.wavelength 
132      Qi = model.fitQ(Q,dQ) 
133      Ri = model.reflectivity(Qi,L) 
134      R = reflfit.convolve(Qi,Ri,Q,dQ) 
135      data.R, data.dR = data.R/R, data.dR/R 
136   
137 -def demo(): 
138      import pylab 
139      from reflectometry.reduction import normalize 
140      from reflectometry.reduction.examples import ng7 as dataset 
141      spec = dataset.spec()[0] 
142      water = WaterIntensity(D2O=20,probe=spec.probe) 
143      spec.apply(normalize()) 
144      theory = water.model(spec.Qz,spec.detector.wavelength) 
145   
146      pylab.subplot(211) 
147      pylab.title('Data normalized to water scattering (%g%% D2O)'%water.D2O) 
148      pylab.xlabel('Qz (inv Ang)') 
149      pylab.ylabel('Reflectivity') 
150      pylab.semilogy(spec.Qz,theory,'-',label='expected') 
151      scale = theory[0]/spec.R[0] 
152      pylab.errorbar(spec.Qz,scale*spec.R,scale*spec.dR,fmt='.',label='measured') 
153   
154      spec.apply(water) 
155      pylab.subplot(212) 
156      #pylab.title('Intensity correction factor') 
157      pylab.xlabel('Slit 1 opening (mm)') 
158      pylab.ylabel('Incident intensity') 
159      pylab.yscale('log') 
160      pylab.errorbar(spec.slit1.x,spec.R,spec.dR,fmt='.',label='correction') 
161   
162      pylab.show() 
163   
164  if __name__ == "__main__": demo() 
165