Source Code for Module reflectometry.model1d.adaptor.reflTheory

  1  import os 
  2  import sys 
  3  import numpy 
  4   
  5  from park.theory.parkTheory import Theory 
  6  from park import XmlParameter 
  7  from park import XmlDataArray 
  8  from park import setDefault 
  9   
 10  from calcReflTheory    import calcReflTheory 
 11  from calcMagReflTheory import calcMagReflTheory, calcUnpolarizedMagTheory 
 12  from xmlProfile import PROFILE_TAG 
 13   
 14  TIMMING = False 
 15  if TIMMING: 
 16      import time 
 17   
 18   
 19 -class ReflParameter(XmlParameter): 
 20   
 21 -    def __init__(self): 
 22          super(ReflParameter, self).__init__() 
 23          self.setDefault() 
 24   
 25   
 26 -    def _postParse(self):  # set the default value 
 27          self.setDefault() 
 28   
 29   
 30 -    def setDefault(self): 
 31          setDefault(self, 'depth',   100.0) 
 32          setDefault(self, 'rho',     1.0) 
 33          setDefault(self, 'mu',      0.0) 
 34          setDefault(self, 'rough',   0.0) 
 35   
 36   
 37   
 38 -class ReflParameterMag(XmlParameter): 
 39   
 40 -    def __init__(self): 
 41          super(ReflParameterMag, self).__init__() 
 42          self.setDefault() 
 43   
 44   
 45 -    def _postParse(self):  # set the default value 
 46          self.setDefault() 
 47   
 48   
 49 -    def setDefault(self): 
 50          setDefault(self, 'depth',   100.0) 
 51          setDefault(self, 'rho',     1.0) 
 52          setDefault(self, 'mu',      0.0) 
 53          setDefault(self, 'phi',     0.0) 
 54          setDefault(self, 'theta',   0.0) 
 55          setDefault(self, 'rough',   0.0) 
 56   
 57   
 58   
 59   
 60 -class ReflTheory(Theory): 
 61      """ 
 62      Calculate the theory reflectivity of model. 
 63      """ 
 64 -    def __init__(self): 
 65          Theory.__init__(self) 
 66   
 67   
 68 -    def _getMeta4Dataset(self): 
 69          """ 
 70          Get metadat from first data of dataset 
 71          """ 
 72          try: 
 73              _meta = self._dataset.getChildren('data')[0].getChildren('meta')[0] 
 74          except: 
 75              _meta = None 
 76          return _meta 
 77   
 78   
 79 -    def _getProfile4Dataset(self): 
 80          """ 
 81          Get profile from dataset 
 82          """ 
 83          try: 
 84              _profile =  self._dataset.getChildren( PROFILE_TAG )[0] 
 85          except: 
 86              _profile = None 
 87   
 88          if _profile != None and hasattr(_profile, 'getDataArrayList'): 
 89              return _profile 
 90          else: 
 91              return None 
 92   
 93   
 94   
 95 -    def _getQdQ(self, data): 
 96          try: 
 97              Q  = data[0].getData() 
 98              dQ = data[1].getData() 
 99              array = XmlDataArray() 
100          except: 
101              Q  = data[0] 
102              dQ = data[1] 
103              array = None 
104   
105          return (Q,dQ,array) 
106   
107   
108 -    def _getMagQdQ(self, data): 
109          Qs=[]; dQs=[]; arrays=[] 
110          try: 
111              for i in xrange(4): 
112                  Qs.append(  data[0+i*4].getData() ) 
113                  dQs.append( data[1+i*4].getData() ) 
114                  arrays.append( XmlDataArray() ) 
115          except: 
116              for i in xrange(4): 
117                  Qs.append(  data[0+i*4] ) 
118                  dQs.append( data[1+i*4] ) 
119                  arrays.append( None ) 
120   
121          return (Qs, dQs, arrays) 
122   
123   
124 -    def _getNonMagFx(self, data, params, _profile, _meta): 
125          """ 
126          Calculate the function for nonmagnetic case 
127          """ 
128          (Q, dQ, array) = self._getQdQ(data) 
129   
130          if TIMMING: t0 = time.clock() 
131          R = calcReflTheory( Q, params, _meta, _profile, dQ ) 
132          if TIMMING: print "Calc Theory time=", time.clock()-t0 
133   
134          if array is None: 
135              return (R,) 
136          else: 
137              array.setData(R) 
138              return (array, ) 
139   
140   
141 -    def _getMagFx(self, data, params, _profile, _meta): 
142          """ 
143          Calculate the function for magnetic case 
144          """ 
145          (Qs, dQs, arrays) = self._getMagQdQ(data) 
146   
147          if TIMMING: t0 = time.clock() 
148          R = calcMagReflTheory( Qs, params, _meta, _profile, dQs ) 
149          if TIMMING: print "Calc Theory time=", time.clock()-t0 
150   
151          if  arrays[0] == None: 
152              return R 
153          else: 
154              for i in xrange(4): 
155                  if R == None:  arrays[i].setData( None ) 
156                  else:          arrays[i].setData( R[i] ) 
157   
158              return arrays 
159   
160   
161 -    def _getUnpolarizedMagneticFx(self, data, params, _profile, _meta): 
162          """ 
163          Calculate the function for Unpolarized magnetic case 
164          """ 
165          (Q, dQ, array) = self._getQdQ(data) 
166   
167          if TIMMING: t0 = time.clock() 
168          R = calcUnpolarizedMagTheory( Q, params, _meta, _profile, dQ ) 
169          if TIMMING: print "Calc Theory time=", time.clock()-t0 
170   
171          if array is None: 
172              return (R,) 
173          else: 
174              array.setData(R) 
175              return (array, ) 
176   
177   
178 -    def isPolarized(self, data): 
179          """ 
180          Decide the data is polarized or not? 
181          """ 
182          if   len(data)==4:  return False 
183          elif len(data)==16: return True 
184          else:    pass 
185   
186   
187 -    def _getFx(self, data, params): 
188          """ 
189          Calculate the function. data is the data 
190          object return from the reduction data of dataset. 
191          The parameters are list of parameter. 
192          """ 
193          _meta    = self._getMeta4Dataset() 
194          _profile = self._getProfile4Dataset() 
195   
196          # When no model paremeters, just return 
197          if len(params) < 1: 
198              return 
199   
200          # The model is magnetic case or not? 
201          isMagnetic = False 
202          if hasattr(params[0], 'phi'): 
203              isMagnetic = True 
204   
205          # Double-check the model is magnetic or not? 
206          if _profile is not None: 
207              if _profile.isMagnetic() != isMagnetic: 
208                  raise TypeError, 'The model Type is not consistant' 
209   
210          # Get the theory curve 
211          if not isMagnetic: 
212              return self._getNonMagFx(data, params, _profile, _meta) 
213   
214          elif self.isPolarized(data): 
215              return self._getMagFx(data, params, _profile, _meta) 
216   
217          else: 
218              return self._getUnpolarizedMagneticFx(data,params, _profile, _meta) 
219   
220   
221 -    def _calcObjective(self, R, dR, calcR): 
222          R0 = numpy.ones(len(dR))*1.0e-10 
223          dR = numpy.choose( dR<1.0e-10, (dR,R0) ) 
224   
225          if dR == None or len(dR)==0:  dy =  R - calcR 
226          else:                         dy = (R - calcR)/dR 
227   
228          if len(R) >=1:  return numpy.dot(dy, dy)/len(R) 
229          else:           return numpy.dot(dy, dy) 
230   
231   
232 -    def _getObjectiveNonMagFx(self, data, params): 
233          """ 
234          Calculate the chisq for nonmagnetic case. 
235          """ 
236          (Q, dQ, array) = self._getQdQ(data) 
237   
238          _meta    = self._getMeta4Dataset() 
239          _profile = None 
240   
241          calcR = calcReflTheory( Q, params, _meta, _profile, dQ ) 
242   
243          try:    R = data[2].getData() 
244          except: R = data[2] 
245   
246          try:    dR = data[3].getData() 
247          except: dR = None 
248   
249          return self._calcObjective(R, dR, calcR) 
250   
251   
252 -    def _getObjectiveUnpolarizedMagFx(self, data, params): 
253          """ 
254          Calculate the chisq for Unpolarized magnetic case. 
255          """ 
256          (Q, dQ, array) = self._getQdQ(data) 
257   
258          _meta    = self._getMeta4Dataset() 
259          _profile = None 
260   
261          calcR = calcUnpolarizedMagTheory( Q, params, _meta, _profile, dQ ) 
262   
263          try:    R = data[2].getData() 
264          except: R = data[2] 
265   
266          try:    dR = data[3].getData() 
267          except: dR = None 
268   
269          return self._calcObjective(R, dR, calcR) 
270   
271   
272 -    def _getObjectiveMagFx(self, data, params): 
273          """ 
274          Calculate the chisq for magnetic case. 
275          """ 
276          (Qs, dQs, arrays) = self._getMagQdQ(data) 
277   
278          _meta    = self._getMeta4Dataset() 
279          _profile = None 
280   
281          calcR = calcMagReflTheory( Qs, params, _meta, _profile, dQs ) 
282   
283          R = [] 
284          try: 
285              for i in xrange(4): R.append ( data[2+i*4].getData() ) 
286          except: 
287              for i in xrange(4): R.append ( data[2+i*4] ) 
288   
289          dR = [] 
290          try: 
291              for i in xrange(4): dR.append ( data[3+i*4].getData() ) 
292          except: 
293              for i in xrange(4): dR.append ( data[3+i*4] ) 
294   
295   
296          _sum = 0.0 
297          for i in xrange(4): 
298              dy = (R[i] - calcR[i])/dR[i] 
299              _sum += numpy.dot(dy,dy) 
300   
301          return _sum/ ( len(R[0]) + len(R[1]) + len(R[2]) + len(R[3]) ) 
302   
303   
304   
305 -    def _getObjectiveFx(self, data, params): 
306          """ 
307          Calculate the chisq. 
308          """ 
309           # When no model paremeters, just return 
310          if len(params) < 1: 
311              return 
312   
313          # The model is magnetic case or not? 
314          isMagnetic = False 
315          if hasattr(params[0], 'phi'): 
316              isMagnetic = True 
317   
318          # Get the theory curve 
319          if not isMagnetic: 
320              return self._getObjectiveNonMagFx(data, params) 
321   
322          elif self.isPolarized(data): 
323              return self._getObjectiveMagFx(data, params) 
324   
325          else: 
326              return self._getObjectiveUnpolarizedMagFx(data, params) 
327   
328   
329   
330 -    def _getNonMagResidual(self, data, params): 
331   
332          (Q, dQ, array) = self._getQdQ(data) 
333   
334          _meta    = self._getMeta4Dataset() 
335          _profile = None 
336   
337          calcR = calcReflTheory( Q, params, _meta, dQ ) 
338   
339          try:    R = data[2].getData() 
340          except: R = data[2] 
341   
342          return  (R-calcR) 
343   
344   
345 -    def _getUnpolarizedMagResidual(self, data, params): 
346   
347          (Q, dQ, array) = self._getQdQ(data) 
348   
349          _meta    = self._getMeta4Dataset() 
350          _profile = None 
351   
352          calcR =  calcUnpolarizedMagTheory( Q, params, _meta, _profile, dQ ) 
353   
354          try:    R = data[2].getData() 
355          except: R = data[2] 
356   
357          return  (R-calcR) 
358   
359   
360 -    def _getMagResidual(self, data, params): 
361          (Qs, dQs, arrays) = self._getMagQdQ(data) 
362   
363          _meta    = self._getMeta4Dataset() 
364          _profile  = None 
365   
366          calcR = calcMagReflTheory( xdata, params, _meta ) 
367   
368          try:    R = data[2].getData() 
369          except: R = data[2] 
370   
371          dy = R - calcR[0] 
372   
373          return dy 
374   
375   
376 -    def _getResidual(self, data, params): 
377          """ 
378          Calculate the chisq. 
379          """ 
380          if   len(data) == 4:  return self._getNonMagResidual(data, params) 
381          elif len(data) == 16: return self._getMagResidual(data, params) 
382          else:                 return 0.0 
383