| 1 | #GSASII powder calculation module |
|---|
| 2 | ########### SVN repository information ################### |
|---|
| 3 | # $Date: 2011-03-31 21:41:13 +0000 (Thu, 31 Mar 2011) $ |
|---|
| 4 | # $Author: vondreele $ |
|---|
| 5 | # $Revision: 264 $ |
|---|
| 6 | # $URL: trunk/GSASIIpeak.py $ |
|---|
| 7 | # $Id: GSASIIpeak.py 264 2011-03-31 21:41:13Z vondreele $ |
|---|
| 8 | ########### SVN repository information ################### |
|---|
| 9 | import sys |
|---|
| 10 | import math |
|---|
| 11 | import wx |
|---|
| 12 | import time |
|---|
| 13 | import numpy as np |
|---|
| 14 | import numpy.linalg as nl |
|---|
| 15 | import GSASIIpath |
|---|
| 16 | import pypowder as pyp #assumes path has been amended to include correctr bin directory |
|---|
| 17 | import GSASIIplot as G2plt |
|---|
| 18 | |
|---|
| 19 | # trig functions in degrees |
|---|
| 20 | sind = lambda x: math.sin(x*math.pi/180.) |
|---|
| 21 | asind = lambda x: 180.*math.asin(x)/math.pi |
|---|
| 22 | tand = lambda x: math.tan(x*math.pi/180.) |
|---|
| 23 | atand = lambda x: 180.*math.atan(x)/math.pi |
|---|
| 24 | atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi |
|---|
| 25 | cosd = lambda x: math.cos(x*math.pi/180.) |
|---|
| 26 | acosd = lambda x: 180.*math.acos(x)/math.pi |
|---|
| 27 | rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p) |
|---|
| 28 | #numpy versions |
|---|
| 29 | npsind = lambda x: np.sin(x*np.pi/180.) |
|---|
| 30 | npasind = lambda x: 180.*np.arcsin(x)/math.pi |
|---|
| 31 | npcosd = lambda x: np.cos(x*math.pi/180.) |
|---|
| 32 | nptand = lambda x: np.tan(x*math.pi/180.) |
|---|
| 33 | npatand = lambda x: 180.*np.arctan(x)/np.pi |
|---|
| 34 | npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
|---|
| 35 | |
|---|
| 36 | def factorize(num): |
|---|
| 37 | ''' Provide prime number factors for integer num |
|---|
| 38 | Returns dictionary of prime factors (keys) & power for each (data) |
|---|
| 39 | ''' |
|---|
| 40 | factors = {} |
|---|
| 41 | orig = num |
|---|
| 42 | |
|---|
| 43 | # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1 |
|---|
| 44 | i, sqi = 2, 4 |
|---|
| 45 | while sqi <= num: |
|---|
| 46 | while not num%i: |
|---|
| 47 | num /= i |
|---|
| 48 | factors[i] = factors.get(i, 0) + 1 |
|---|
| 49 | |
|---|
| 50 | sqi += 2*i + 1 |
|---|
| 51 | i += 1 |
|---|
| 52 | |
|---|
| 53 | if num != 1 and num != orig: |
|---|
| 54 | factors[num] = factors.get(num, 0) + 1 |
|---|
| 55 | |
|---|
| 56 | if factors: |
|---|
| 57 | return factors |
|---|
| 58 | else: |
|---|
| 59 | return {num:1} #a prime number! |
|---|
| 60 | |
|---|
| 61 | def makeFFTsizeList(nmin=1,nmax=1023,thresh=15): |
|---|
| 62 | ''' Provide list of optimal data sizes for FFT calculations |
|---|
| 63 | Input: |
|---|
| 64 | nmin: minimum data size >= 1 |
|---|
| 65 | nmax: maximum data size > nmin |
|---|
| 66 | thresh: maximum prime factor allowed |
|---|
| 67 | Returns: |
|---|
| 68 | list of data sizes where the maximum prime factor is < thresh |
|---|
| 69 | ''' |
|---|
| 70 | plist = [] |
|---|
| 71 | nmin = max(1,nmin) |
|---|
| 72 | nmax = max(nmin+1,nmax) |
|---|
| 73 | for p in range(nmin,nmax): |
|---|
| 74 | if max(factorize(p).keys()) < thresh: |
|---|
| 75 | plist.append(p) |
|---|
| 76 | return plist |
|---|
| 77 | |
|---|
| 78 | def Transmission(Geometry,Abs,Diam): |
|---|
| 79 | #Calculate sample transmission |
|---|
| 80 | # Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate' |
|---|
| 81 | # Abs: absorption coeff in cm-1 |
|---|
| 82 | # Diam: sample thickness/diameter in mm |
|---|
| 83 | if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample |
|---|
| 84 | MuR = Abs*Diam/5.0 |
|---|
| 85 | if MuR <= 3.0: |
|---|
| 86 | T0 = 16/(3.*math.pi) |
|---|
| 87 | T1 = -0.045780 |
|---|
| 88 | T2 = -0.02489 |
|---|
| 89 | T3 = 0.003045 |
|---|
| 90 | T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4 |
|---|
| 91 | if T < -20.: |
|---|
| 92 | return 2.06e-9 |
|---|
| 93 | else: |
|---|
| 94 | return math.exp(T) |
|---|
| 95 | else: |
|---|
| 96 | T1 = 1.433902 |
|---|
| 97 | T2 = 0.013869+0.337894 |
|---|
| 98 | T3 = 1.933433+1.163198 |
|---|
| 99 | T4 = 0.044365-0.04259 |
|---|
| 100 | T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4 |
|---|
| 101 | return T/100. |
|---|
| 102 | |
|---|
| 103 | def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0): |
|---|
| 104 | #Calculate sample absorption |
|---|
| 105 | # Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate' |
|---|
| 106 | # Abs: absorption coeff in cm-1 |
|---|
| 107 | # Diam: sample thickness/diameter in mm |
|---|
| 108 | # Tth: 2-theta scattering angle - can be numpy array |
|---|
| 109 | # Phi: flat plate tilt angle - future |
|---|
| 110 | # Psi: flat plate tilt axis - future |
|---|
| 111 | MuR = Abs*Diam/5.0 |
|---|
| 112 | Sth2 = npsind(Tth/2.0)**2 |
|---|
| 113 | Cth2 = 1.-Sth2 |
|---|
| 114 | if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample |
|---|
| 115 | if MuR < 3.0: |
|---|
| 116 | T0 = 16.0/(3*np.pi) |
|---|
| 117 | T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \ |
|---|
| 118 | 0.109561*np.sqrt(Sth2)-26.04556 |
|---|
| 119 | T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \ |
|---|
| 120 | 1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3 |
|---|
| 121 | T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2 |
|---|
| 122 | Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4 |
|---|
| 123 | return np.exp(Trns) |
|---|
| 124 | else: |
|---|
| 125 | T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \ |
|---|
| 126 | 10.02088*Sth2**3-3.36778*Sth2**4 |
|---|
| 127 | T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \ |
|---|
| 128 | (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039 |
|---|
| 129 | T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \ |
|---|
| 130 | 0.13576*np.sqrt(Sth2)+1.163198 |
|---|
| 131 | T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202 |
|---|
| 132 | Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4 |
|---|
| 133 | return Trns/100. |
|---|
| 134 | elif 'Bragg' in Geometry: |
|---|
| 135 | return 1.0 |
|---|
| 136 | elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam |
|---|
| 137 | # and only defined for 2theta < 90 |
|---|
| 138 | T1 = np.exp(-MuR) |
|---|
| 139 | T2 = np.exp(-MuR/(1.-2.*Sth2)) |
|---|
| 140 | Tb = -2.*Abs*Sth2 |
|---|
| 141 | return (T1-T2)/Tb |
|---|
| 142 | elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector |
|---|
| 143 | cth = npcosd(Tth/2.0) |
|---|
| 144 | return (Diam/cth)*np.exp(-MuR/cth) |
|---|
| 145 | |
|---|
| 146 | def Polarization(Pola,Azm,Tth): |
|---|
| 147 | # Calculate x-ray polarization correction |
|---|
| 148 | # Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized |
|---|
| 149 | # Azm: azimuthal angle e.g. 0.0 in plane of polarization(?) |
|---|
| 150 | # Tth: 2-theta scattering angle - can be numpy array |
|---|
| 151 | pass |
|---|
| 152 | |
|---|
| 153 | |
|---|
| 154 | |
|---|
| 155 | def ValEsd(value,esd=0,nTZ=False): #NOT complete - don't use |
|---|
| 156 | # returns value(esd) string; nTZ=True for no trailing zeros |
|---|
| 157 | # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal |
|---|
| 158 | #get the 2 significant digits in the esd |
|---|
| 159 | edig = lambda esd: int(round(10**(math.log10(esd) % 1+1))) |
|---|
| 160 | #get the number of digits to represent them |
|---|
| 161 | epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd))) |
|---|
| 162 | |
|---|
| 163 | mdec = lambda esd: -int(math.log10(abs(esd))) |
|---|
| 164 | ndec = lambda esd: int(1.545-math.log10(abs(esd))) |
|---|
| 165 | if esd > 0: |
|---|
| 166 | fmt = '"%.'+str(ndec(esd))+'f(%d)"' |
|---|
| 167 | print fmt,ndec(esd),esd*10**(mdec(esd)+1) |
|---|
| 168 | return fmt%(value,int(esd*10**(mdec(esd)+1))) |
|---|
| 169 | elif esd < 0: |
|---|
| 170 | return str(round(value,mdec(esd))) |
|---|
| 171 | else: |
|---|
| 172 | text = "%F"%(value) |
|---|
| 173 | if nTZ: |
|---|
| 174 | return text.rstrip('0') |
|---|
| 175 | else: |
|---|
| 176 | return text |
|---|
| 177 | |
|---|
| 178 | |
|---|
| 179 | #GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model |
|---|
| 180 | |
|---|
| 181 | def DoPeakFit(peaks,background,limits,inst,data): |
|---|
| 182 | |
|---|
| 183 | def backgroundPrint(background,sigback): |
|---|
| 184 | if background[1]: |
|---|
| 185 | print 'Background coefficients for',background[0],'function' |
|---|
| 186 | ptfmt = "%12.5f" |
|---|
| 187 | ptstr = 'values:' |
|---|
| 188 | sigstr = 'esds :' |
|---|
| 189 | for i,back in enumerate(background[3:]): |
|---|
| 190 | ptstr += ptfmt % (back) |
|---|
| 191 | sigstr += ptfmt % (sigback[i+3]) |
|---|
| 192 | print ptstr |
|---|
| 193 | print sigstr |
|---|
| 194 | else: |
|---|
| 195 | print 'Background not refined' |
|---|
| 196 | |
|---|
| 197 | def instPrint(instVal,siginst,insLabels): |
|---|
| 198 | print 'Instrument Parameters:' |
|---|
| 199 | ptfmt = "%12.6f" |
|---|
| 200 | ptlbls = 'names :' |
|---|
| 201 | ptstr = 'values:' |
|---|
| 202 | sigstr = 'esds :' |
|---|
| 203 | for i,value in enumerate(instVal): |
|---|
| 204 | ptlbls += "%s" % (insLabels[i].center(12)) |
|---|
| 205 | ptstr += ptfmt % (value) |
|---|
| 206 | if siginst[i]: |
|---|
| 207 | sigstr += ptfmt % (siginst[i]) |
|---|
| 208 | else: |
|---|
| 209 | sigstr += 12*' ' |
|---|
| 210 | print ptlbls |
|---|
| 211 | print ptstr |
|---|
| 212 | print sigstr |
|---|
| 213 | |
|---|
| 214 | def peaksPrint(peaks,sigpeaks): |
|---|
| 215 | print 'Peak list=' |
|---|
| 216 | |
|---|
| 217 | begin = time.time() |
|---|
| 218 | Np = len(peaks[0]) |
|---|
| 219 | DataType = inst[1][0] |
|---|
| 220 | instVal = inst[1][1:] |
|---|
| 221 | Insref = inst[2][1:] |
|---|
| 222 | insLabels = inst[3][1:] |
|---|
| 223 | Ka2 = False |
|---|
| 224 | Ioff = 3 |
|---|
| 225 | if len(instVal) == 12: |
|---|
| 226 | lamratio = instVal[1]/instVal[0] |
|---|
| 227 | Ka2 = True |
|---|
| 228 | Ioff = 5 |
|---|
| 229 | insref = Insref[len(Insref)-7:-1] #just U,V,W,X,Y,SH/L |
|---|
| 230 | for peak in peaks: |
|---|
| 231 | dip = [] |
|---|
| 232 | dip.append(tand(peak[0]/2.0)**2) |
|---|
| 233 | dip.append(tand(peak[0]/2.0)) |
|---|
| 234 | dip.append(1.0/cosd(peak[0]/2.0)) |
|---|
| 235 | dip.append(tand(peak[0]/2.0)) |
|---|
| 236 | peak.append(dip) |
|---|
| 237 | B = background[2] |
|---|
| 238 | bcof = background[3:3+B] |
|---|
| 239 | Bv = 0 |
|---|
| 240 | if background[1]: |
|---|
| 241 | Bv = B |
|---|
| 242 | x,y,w,yc,yb,yd = data #these are numpy arrays! |
|---|
| 243 | V = [] |
|---|
| 244 | A = [] |
|---|
| 245 | swobs = 0.0 |
|---|
| 246 | smin = 0.0 |
|---|
| 247 | first = True |
|---|
| 248 | GoOn = True |
|---|
| 249 | Go = True |
|---|
| 250 | dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \ |
|---|
| 251 | style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT) |
|---|
| 252 | screenSize = wx.DisplaySize() |
|---|
| 253 | Size = dlg.GetSize() |
|---|
| 254 | dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0)) |
|---|
| 255 | try: |
|---|
| 256 | i = 0 |
|---|
| 257 | for xi in x : |
|---|
| 258 | Go = dlg.Update(i)[0] |
|---|
| 259 | if GoOn: |
|---|
| 260 | GoOn = Go |
|---|
| 261 | if limits[0] <= xi <= limits[1]: |
|---|
| 262 | yb[i] = 0.0 |
|---|
| 263 | dp = [] |
|---|
| 264 | for j in range(B): |
|---|
| 265 | t = (xi-limits[0])**j |
|---|
| 266 | yb[i] += t*bcof[j] |
|---|
| 267 | if background[1]: |
|---|
| 268 | dp.append(t) |
|---|
| 269 | yc[i] = yb[i] |
|---|
| 270 | Iv = 0 |
|---|
| 271 | for j in range(6): |
|---|
| 272 | if insref[j]: |
|---|
| 273 | dp.append(0.0) |
|---|
| 274 | Iv += 1 |
|---|
| 275 | for peak in peaks: |
|---|
| 276 | dip = peak[-1] |
|---|
| 277 | f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],instVal[-2],0.0) |
|---|
| 278 | yc[i] += f[0]*peak[2] |
|---|
| 279 | if f[0] > 0.0: |
|---|
| 280 | j = 0 |
|---|
| 281 | if insref[0]: #U |
|---|
| 282 | dp[Bv+j] += f[3]*dip[0] |
|---|
| 283 | j += 1 |
|---|
| 284 | if insref[1]: #V |
|---|
| 285 | dp[Bv+j] += f[3]*dip[1] |
|---|
| 286 | j += 1 |
|---|
| 287 | if insref[2]: #W |
|---|
| 288 | dp[Bv+j] += f[3] |
|---|
| 289 | j += 1 |
|---|
| 290 | if insref[3]: #X |
|---|
| 291 | dp[Bv+j] += f[4]*dip[2] |
|---|
| 292 | j += 1 |
|---|
| 293 | if insref[4]: #Y |
|---|
| 294 | dp[Bv+j] += f[4]*dip[3] |
|---|
| 295 | j += 1 |
|---|
| 296 | if insref[5]: #SH/L |
|---|
| 297 | dp[Bv+j] += f[5] |
|---|
| 298 | if Ka2: |
|---|
| 299 | pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0)) |
|---|
| 300 | f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],instVal[-2],0.0) |
|---|
| 301 | yc[i] += f2[0]*peak[2]*instVal[3] |
|---|
| 302 | if f[0] > 0.0: |
|---|
| 303 | j = 0 |
|---|
| 304 | if insref[0]: #U |
|---|
| 305 | dp[Bv+j] += f2[3]*dip[0]*instVal[3] |
|---|
| 306 | j += 1 |
|---|
| 307 | if insref[1]: #V |
|---|
| 308 | dp[Bv+j] += f2[3]*dip[1]*instVal[3] |
|---|
| 309 | j += 1 |
|---|
| 310 | if insref[2]: #W |
|---|
| 311 | dp[Bv+j] += f2[3]*instVal[3] |
|---|
| 312 | j += 1 |
|---|
| 313 | if insref[3]: #X |
|---|
| 314 | dp[Bv+j] += f2[4]*dip[2]*instVal[3] |
|---|
| 315 | j += 1 |
|---|
| 316 | if insref[4]: #Y |
|---|
| 317 | dp[Bv+j] += f2[4]*dip[3]*instVal[3] |
|---|
| 318 | j += 1 |
|---|
| 319 | if insref[5]: #SH/L |
|---|
| 320 | dp[Bv+j] += f2[5]*instVal[3] |
|---|
| 321 | for j in range(0,Np,2): |
|---|
| 322 | if peak[j+1]: dp.append(f[j/2+1]) |
|---|
| 323 | yd[i] = y[i]-yc[i] |
|---|
| 324 | swobs += w[i]*y[i]**2 |
|---|
| 325 | t2 = w[i]*yd[i] |
|---|
| 326 | smin += t2*yd[i] |
|---|
| 327 | if first: |
|---|
| 328 | first = False |
|---|
| 329 | M = len(dp) |
|---|
| 330 | A = np.zeros(shape=(M,M)) |
|---|
| 331 | V = np.zeros(shape=(M)) |
|---|
| 332 | A,V = pyp.buildmv(t2,w[i],M,dp,A,V) |
|---|
| 333 | i += 1 |
|---|
| 334 | finally: |
|---|
| 335 | dlg.Destroy() |
|---|
| 336 | Rwp = smin/swobs |
|---|
| 337 | Rwp = math.sqrt(Rwp)*100.0 |
|---|
| 338 | norm = np.diag(A) |
|---|
| 339 | for i,elm in enumerate(norm): |
|---|
| 340 | if elm <= 0.0: |
|---|
| 341 | print norm |
|---|
| 342 | return False,0,0,0,False |
|---|
| 343 | for i in xrange(len(V)): |
|---|
| 344 | norm[i] = 1.0/math.sqrt(norm[i]) |
|---|
| 345 | V[i] *= norm[i] |
|---|
| 346 | a = A[i] |
|---|
| 347 | for j in xrange(len(V)): |
|---|
| 348 | a[j] *= norm[i] |
|---|
| 349 | A = np.transpose(A) |
|---|
| 350 | for i in xrange(len(V)): |
|---|
| 351 | a = A[i] |
|---|
| 352 | for j in xrange(len(V)): |
|---|
| 353 | a[j] *= norm[i] |
|---|
| 354 | b = nl.solve(A,V) |
|---|
| 355 | A = nl.inv(A) |
|---|
| 356 | sig = np.diag(A) |
|---|
| 357 | for i in xrange(len(V)): |
|---|
| 358 | b[i] *= norm[i] |
|---|
| 359 | sig[i] *= norm[i]*norm[i] |
|---|
| 360 | sig[i] = math.sqrt(abs(sig[i])) |
|---|
| 361 | sigback = [0,0,0] |
|---|
| 362 | for j in range(Bv): |
|---|
| 363 | background[j+3] += b[j] |
|---|
| 364 | sigback.append(sig[j]) |
|---|
| 365 | backgroundPrint(background,sigback) |
|---|
| 366 | k = 0 |
|---|
| 367 | delt = [] |
|---|
| 368 | if Ka2: |
|---|
| 369 | siginst = [0,0,0,0,0] |
|---|
| 370 | else: |
|---|
| 371 | siginst = [0,0,0] |
|---|
| 372 | for j in range(6): |
|---|
| 373 | if insref[j]: |
|---|
| 374 | instVal[j+Ioff] += b[Bv+k] |
|---|
| 375 | siginst.append(sig[Bv+k]) |
|---|
| 376 | delt.append(b[Bv+k]) |
|---|
| 377 | k += 1 |
|---|
| 378 | else: |
|---|
| 379 | delt.append(0.0) |
|---|
| 380 | siginst.append(0.0) |
|---|
| 381 | delt.append(0.0) #dummies for azm |
|---|
| 382 | siginst.append(0.0) |
|---|
| 383 | instPrint(instVal,siginst,insLabels) |
|---|
| 384 | inst[1] = [DataType,] |
|---|
| 385 | for val in instVal: |
|---|
| 386 | inst[1].append(val) |
|---|
| 387 | B = Bv+Iv |
|---|
| 388 | for peak in peaks: |
|---|
| 389 | del peak[-1] # remove dip from end |
|---|
| 390 | delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2] |
|---|
| 391 | delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0) |
|---|
| 392 | for j in range(0,len(peak[:-1]),2): |
|---|
| 393 | if peak[j+1]: |
|---|
| 394 | peak[j] += b[B] |
|---|
| 395 | B += 1 |
|---|
| 396 | peak[4] += delsig |
|---|
| 397 | if peak[4] < 0.0: |
|---|
| 398 | print 'ERROR - negative sigma' |
|---|
| 399 | return False,0,0,0,False |
|---|
| 400 | peak[6] += delgam |
|---|
| 401 | if peak[6] < 0.0: |
|---|
| 402 | print 'ERROR - negative gamma' |
|---|
| 403 | return False,0,0,0,False |
|---|
| 404 | runtime = time.time()-begin |
|---|
| 405 | data = [x,y,w,yc,yb,yd] |
|---|
| 406 | return True,smin,Rwp,runtime,GoOn |
|---|
| 407 | |
|---|
| 408 | def ComputePDF(data,xydata): |
|---|
| 409 | for key in data: |
|---|
| 410 | print key,data[key] |
|---|
| 411 | #subtract backgrounds - if any |
|---|
| 412 | xydata['Sample corrected'] = xydata['Sample'] |
|---|
| 413 | if 'Sample Bkg.' in xydata: |
|---|
| 414 | xydata['Sample corrected'][1][1] -= (xydata['Sample Bkg.'][1][1]+ |
|---|
| 415 | data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult'] |
|---|
| 416 | if 'Container' in xydata: |
|---|
| 417 | xydata['Sample corrected'][1][1] -= (xydata['Container'][1][1]+ |
|---|
| 418 | data['Container']['Add'])*data['Container']['Mult'] |
|---|
| 419 | if 'Container Bkg.' in xydata: |
|---|
| 420 | xydata['Sample corrected'][1][1] += (xydata['Container Bkg.'][1][1]+ |
|---|
| 421 | data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult'] |
|---|
| 422 | |
|---|
| 423 | |
|---|
| 424 | |
|---|
| 425 | return xydata['Sample corrected'],[] |
|---|
| 426 | |
|---|
