Changeset f60a8c2 in sasview for sansdataloader/src/sans/dataloader/manipulations.py
- Timestamp:
- Apr 27, 2012 10:42:24 AM (12 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- 7d6351e
- Parents:
- 10bfeb3
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansdataloader/src/sans/dataloader/manipulations.py
r47045e6 rf60a8c2 1 1 """ 2 Data manipulations for 2D data sets. 3 Using the meta data information, various types of averaging 4 are performed in Q-space 5 """ 2 6 ##################################################################### 3 7 #This software was developed by the University of Tennessee as part of the 4 8 #Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 5 #project funded by the US National Science Foundation. 9 #project funded by the US National Science Foundation. 6 10 #See the license text in license.txt 7 11 #copyright 2008, University of Tennessee 8 12 ###################################################################### 9 13 10 """11 Data manipulations for 2D data sets.12 Using the meta data information, various types of averaging13 are performed in Q-space14 """15 14 #TODO: copy the meta data from the 2D object to the resulting 1D object 16 15 import math … … 33 32 theta = 0.5 * math.atan(plane_dist/det_dist) 34 33 return (4.0 * math.pi/wavelength) * math.sin(theta) 34 35 35 36 36 def get_q_compo(dx, dy, det_dist, wavelength, compo=None): … … 55 55 return out 56 56 57 57 58 def flip_phi(phi): 58 59 """ … … 63 64 Pi = math.pi 64 65 if phi < 0: 65 phi_out = phi 66 phi_out = phi + (2 * Pi) 66 67 elif phi > (2 * Pi): 67 phi_out = phi 68 phi_out = phi - (2 * Pi) 68 69 else: 69 phi_out = phi 70 phi_out = phi 70 71 return phi_out 72 71 73 72 74 def reader2D_converter(data2d=None): … … 93 95 qy_data = new_y.flatten() 94 96 q_data = numpy.sqrt(qx_data*qx_data + qy_data*qy_data) 95 if data2d.err_data == None or numpy.any(data2d.err_data <= 0): 97 if data2d.err_data == None or numpy.any(data2d.err_data <= 0): 96 98 new_err_data = numpy.sqrt(numpy.abs(new_data)) 97 99 else: 98 100 new_err_data = data2d.err_data.flatten() 99 mask = numpy.ones(len(new_data), dtype=bool) 100 101 mask = numpy.ones(len(new_data), dtype=bool) 102 103 #TODO: make sense of the following two lines... 101 104 output = Data2D() 102 105 output = data2d … … 109 112 110 113 return output 114 111 115 112 116 class _Slab(object): … … 149 153 raise RuntimeError, msg 150 154 151 # Get data 155 # Get data 152 156 data = data2D.data[numpy.isfinite(data2D.data)] 153 157 q_data = data2D.q_data[numpy.isfinite(data2D.data)] 154 158 err_data = data2D.err_data[numpy.isfinite(data2D.data)] 155 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 159 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 156 160 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 157 161 … … 159 163 if maj == 'x': 160 164 if self.fold: 161 x_min = 0 162 else: x_min = self.x_min 163 nbins = int(math.ceil((self.x_max - x_min)/self.bin_width)) 165 x_min = 0 166 else: 167 x_min = self.x_min 168 nbins = int(math.ceil((self.x_max - x_min) / self.bin_width)) 164 169 qbins = self.bin_width * numpy.arange(nbins) + x_min 165 170 elif maj == 'y': 166 if self.fold: y_min = 0 167 else: y_min = self.y_min 171 if self.fold: 172 y_min = 0 173 else: 174 y_min = self.y_min 168 175 nbins = int(math.ceil((self.y_max - y_min)/self.bin_width)) 169 qbins = self.bin_width * numpy.arange(nbins) + y_min 176 qbins = self.bin_width * numpy.arange(nbins) + y_min 170 177 else: 171 178 raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj) 172 179 173 x 174 y 180 x = numpy.zeros(nbins) 181 y = numpy.zeros(nbins) 175 182 err_y = numpy.zeros(nbins) 176 183 y_counts = numpy.zeros(nbins) 177 184 178 # Average pixelsize in q space 179 for npts in range(len(data)): 180 # default frac 185 # Average pixelsize in q space 186 for npts in range(len(data)): 187 # default frac 181 188 frac_x = 0 182 189 frac_y = 0 … … 191 198 continue 192 199 # binning: find axis of q 193 if maj == 'x': 200 if maj == 'x': 194 201 q_value = qx_data[npts] 195 min = x_min 196 if maj == 'y': 197 q_value = qy_data[npts] 202 min = x_min 203 if maj == 'y': 204 q_value = qy_data[npts] 198 205 min = y_min 199 206 if self.fold and q_value < 0: 200 q_value = -q_value 207 q_value = -q_value 201 208 # bin 202 i_q = int(math.ceil((q_value - min) /self.bin_width)) - 1209 i_q = int(math.ceil((q_value - min) / self.bin_width)) - 1 203 210 204 211 # skip outside of max bins … … 206 213 continue 207 214 208 # give it full weight209 #frac = 1210 211 215 #TODO: find better definition of x[i_q] based on q_data 212 x[i_q] += frac * q_value #min + (i_q + 1) * self.bin_width / 2.0216 x[i_q] += frac * q_value # min + (i_q + 1) * self.bin_width / 2.0 213 217 y[i_q] += frac * data[npts] 214 218 215 219 if err_data == None or err_data[npts] == 0.0: 216 if data[npts] < 0: data[npts] = -data[npts] 220 if data[npts] < 0: 221 data[npts] = -data[npts] 217 222 err_y[i_q] += frac * frac * data[npts] 218 223 else: 219 224 err_y[i_q] += frac * frac * err_data[npts] * err_data[npts] 220 y_counts[i_q] 225 y_counts[i_q] += frac 221 226 222 # Average the sums 227 # Average the sums 223 228 for n in range(nbins): 224 229 err_y[n] = math.sqrt(err_y[n]) 225 230 226 231 err_y = err_y / y_counts 227 y 228 x 229 idx = (numpy.isfinite(y) & numpy.isfinite(x)) 232 y = y / y_counts 233 x = x / y_counts 234 idx = (numpy.isfinite(y) & numpy.isfinite(x)) 230 235 231 236 if not idx.any(): 232 msg = "Average Error: No points inside ROI to average..." 237 msg = "Average Error: No points inside ROI to average..." 233 238 raise ValueError, msg 234 239 #elif len(y[idx])!= nbins: … … 237 242 return Data1D(x=x[idx], y=y[idx], dy=err_y[idx]) 238 243 244 239 245 class SlabY(_Slab): 240 246 """ … … 251 257 return self._avg(data2D, 'y') 252 258 259 253 260 class SlabX(_Slab): 254 261 """ … … 264 271 265 272 """ 266 return self._avg(data2D, 'x') 267 273 return self._avg(data2D, 'x') 274 275 268 276 class Boxsum(object): 269 277 """ … … 282 290 def __call__(self, data2D): 283 291 """ 284 Perform the sum in the region of interest 292 Perform the sum in the region of interest 285 293 286 294 :param data2D: Data2D object … … 293 301 # Average the sums 294 302 counts = 0 if y_counts == 0 else y 295 error 303 error = 0 if y_counts == 0 else math.sqrt(err_y) 296 304 297 305 return counts, error … … 299 307 def _sum(self, data2D): 300 308 """ 301 Perform the sum in the region of interest 302 303 :param data2D: Data2D object 304 305 :return: number of counts, 309 Perform the sum in the region of interest 310 311 :param data2D: Data2D object 312 313 :return: number of counts, 306 314 error on number of counts, number of entries summed 307 315 … … 311 319 msg += "of detectors: %g" % len(data2D.detector) 312 320 raise RuntimeError, msg 313 # Get data 321 # Get data 314 322 data = data2D.data[numpy.isfinite(data2D.data)] 315 323 q_data = data2D.q_data[numpy.isfinite(data2D.data)] 316 324 err_data = data2D.err_data[numpy.isfinite(data2D.data)] 317 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 325 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 318 326 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 319 327 320 y 328 y = 0.0 321 329 err_y = 0.0 322 330 y_counts = 0.0 323 331 324 # Average pixelsize in q space 325 for npts in range(len(data)): 326 # default frac 327 frac_x = 0 328 frac_y = 0 332 # Average pixelsize in q space 333 for npts in range(len(data)): 334 # default frac 335 frac_x = 0 336 frac_y = 0 329 337 330 338 # get min and max at each points … … 337 345 if self.y_min <= qy and self.y_max > qy: 338 346 frac_y = 1 339 #Find the fraction along each directions 347 #Find the fraction along each directions 340 348 frac = frac_x * frac_y 341 349 if frac == 0: … … 348 356 else: 349 357 err_y += frac * frac * err_data[npts] * err_data[npts] 350 y_counts += frac 358 y_counts += frac 351 359 return y, err_y, y_counts 352 360 353 361 354 355 362 class Boxavg(Boxsum): 356 363 """ … … 363 370 def __call__(self, data2D): 364 371 """ 365 Perform the sum in the region of interest 372 Perform the sum in the region of interest 366 373 367 374 :param data2D: Data2D object … … 373 380 374 381 # Average the sums 375 counts = 0 if y_counts == 0 else y /y_counts376 error = 0 if y_counts == 0 else math.sqrt(err_y)/y_counts382 counts = 0 if y_counts == 0 else y / y_counts 383 error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts 377 384 378 385 return counts, error 379 386 387 380 388 def get_pixel_fraction_square(x, xmin, xmax): 381 389 """ 382 Return the fraction of the length 390 Return the fraction of the length 383 391 from xmin to x.:: 384 392 … … 437 445 # Get the dq for resolution averaging 438 446 if data2D.dqx_data != None and data2D.dqy_data != None: 439 # The pinholes and det. pix contribution present 447 # The pinholes and det. pix contribution present 440 448 # in both direction of the 2D which must be subtracted when 441 449 # converting to 1D: dq_overlap should calculated ideally at 442 # q = 0. Note This method works on only pinhole geometry. 450 # q = 0. Note This method works on only pinhole geometry. 443 451 # Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. 444 452 z_max = max(data2D.q_data) 445 453 z_min = min(data2D.q_data) 446 454 x_max = data2D.dqx_data[data2D.q_data[z_max]] 447 x_min = data2D.dqx_data[data2D.q_data[z_min]] 455 x_min = data2D.dqx_data[data2D.q_data[z_min]] 448 456 y_max = data2D.dqy_data[data2D.q_data[z_max]] 449 y_min = data2D.dqy_data[data2D.q_data[z_min]] 457 y_min = data2D.dqy_data[data2D.q_data[z_min]] 450 458 # Find qdx at q = 0 451 459 dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) … … 453 461 if dq_overlap_x > min(data2D.dqx_data): 454 462 dq_overlap_x = min(data2D.dqx_data) 455 dq_overlap_x *= dq_overlap_x 463 dq_overlap_x *= dq_overlap_x 456 464 # Find qdx at q = 0 457 465 dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) … … 462 470 dq_overlap_y *= dq_overlap_y 463 471 464 dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) /2.0)472 dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 465 473 # Final protection of dq 466 474 if dq_overlap < 0: … … 470 478 # def; dqx_data = dq_r dqy_data = dq_phi 471 479 # Convert dq 2D to 1D here 472 dqx = dqx_data * dqx_data 480 dqx = dqx_data * dqx_data 473 481 dqy = dqy_data * dqy_data 474 482 dq_data = numpy.add(dqx, dqy) … … 484 492 qbins = self.bin_width * numpy.arange(nbins) + self.r_min 485 493 486 x 487 y 494 x = numpy.zeros(nbins) 495 y = numpy.zeros(nbins) 488 496 err_y = numpy.zeros(nbins) 489 497 err_x = numpy.zeros(nbins) 490 498 y_counts = numpy.zeros(nbins) 491 499 492 for npt in range(len(data)): 500 for npt in range(len(data)): 493 501 494 502 if ismask and not mask_data[npt]: 495 continue 503 continue 496 504 497 505 frac = 0 498 506 499 507 # q-value at the pixel (j,i) 500 q_value = q_data[npt] 501 data_n = data[npt] 508 q_value = q_data[npt] 509 data_n = data[npt] 502 510 503 511 ## No need to calculate the frac when all data are within range 504 512 if self.r_min >= self.r_max: 505 raise ValueError, "Limit Error: min > max" 513 raise ValueError, "Limit Error: min > max" 506 514 507 515 if self.r_min <= q_value and q_value <= self.r_max: 508 frac = 1 516 frac = 1 509 517 if frac == 0: 510 518 continue 511 i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) 512 513 # Take care of the edge case at phi = 2pi. 514 if i_q == nbins: 515 i_q = nbins - 1519 i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) 520 521 # Take care of the edge case at phi = 2pi. 522 if i_q == nbins: 523 i_q = nbins - 1 516 524 y[i_q] += frac * data_n 517 525 # Take dqs from data to get the q_average … … 524 532 err_y[i_q] += frac * frac * err_data[npt] * err_data[npt] 525 533 if dq_data != None: 526 # To be consistent with dq calculation in 1d reduction, 527 # we need just the averages (not quadratures) because 528 # it should not depend on the number of the q points 534 # To be consistent with dq calculation in 1d reduction, 535 # we need just the averages (not quadratures) because 536 # it should not depend on the number of the q points 529 537 # in the qr bins. 530 538 err_x[i_q] += frac * dq_data[npt] 531 539 else: 532 540 err_x = None 533 y_counts[i_q] 534 535 # Average the sums 541 y_counts[i_q] += frac 542 543 # Average the sums 536 544 for n in range(nbins): 537 if err_y[n] < 0: err_y[n] = -err_y[n] 545 if err_y[n] < 0: 546 err_y[n] = -err_y[n] 538 547 err_y[n] = math.sqrt(err_y[n]) 539 548 #if err_x != None: … … 541 550 542 551 err_y = err_y / y_counts 543 err_y[err_y ==0] = numpy.average(err_y)544 y 545 x 546 idx = (numpy.isfinite(y)) & (numpy.isfinite(x)) 552 err_y[err_y == 0] = numpy.average(err_y) 553 y = y / y_counts 554 x = x / y_counts 555 idx = (numpy.isfinite(y)) & (numpy.isfinite(x)) 547 556 548 557 if err_x != None: … … 551 560 d_x = None 552 561 553 if not idx.any(): 554 msg = "Average Error: No points inside ROI to average..." 562 if not idx.any(): 563 msg = "Average Error: No points inside ROI to average..." 555 564 raise ValueError, msg 556 565 … … 567 576 around the ring as a function of phi. 568 577 569 Phi_min and phi_max should be defined between 0 and 2*pi 578 Phi_min and phi_max should be defined between 0 and 2*pi 570 579 in anti-clockwise starting from the x- axis on the left-hand side 571 580 """ … … 601 610 q_data = data2D.q_data[numpy.isfinite(data2D.data)] 602 611 err_data = data2D.err_data[numpy.isfinite(data2D.data)] 603 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 612 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 604 613 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 605 614 … … 612 621 phi_err = numpy.zeros(self.nbins_phi) 613 622 614 for npt in range(len(data)): 623 for npt in range(len(data)): 615 624 frac = 0 616 625 # q-value at the point (npt) 617 626 q_value = q_data[npt] 618 data_n = data[npt] 627 data_n = data[npt] 619 628 620 629 # phi-value at the point (npt) … … 622 631 623 632 if self.r_min <= q_value and q_value <= self.r_max: 624 frac = 1 633 frac = 1 625 634 if frac == 0: 626 635 continue … … 628 637 i_phi = int(math.floor((self.nbins_phi) * phi_value / (2 * Pi))) 629 638 630 # Take care of the edge case at phi = 2pi. 631 if i_phi == self.nbins_phi: 632 i_phi = self.nbins_phi - 1 639 # Take care of the edge case at phi = 2pi. 640 if i_phi == self.nbins_phi: 641 i_phi = self.nbins_phi - 1 633 642 phi_bins[i_phi] += frac * data[npt] 634 643 … … 646 655 phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i + 0.5) 647 656 648 idx = (numpy.isfinite(phi_bins)) 657 idx = (numpy.isfinite(phi_bins)) 649 658 650 659 if not idx.any(): 651 msg = "Average Error: No points inside ROI to average..." 660 msg = "Average Error: No points inside ROI to average..." 652 661 raise ValueError, msg 653 662 #elif len(phi_bins[idx])!= self.nbins_phi: … … 656 665 return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) 657 666 667 658 668 def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11): 659 669 """ 660 670 Returns the fraction of the pixel defined by 661 the four corners (q_00, q_01, q_10, q_11) that 671 the four corners (q_00, q_01, q_10, q_11) that 662 672 has q < qmax.:: 663 673 … … 714 724 # this pixel and the constant-q ring. We only need 715 725 # to know if we have to include it or exclude it. 716 elif (q_00 + q_01 + q_10 + q_11) /4.0 <qmax:726 elif (q_00 + q_01 + q_10 + q_11) / 4.0 < qmax: 717 727 frac_max = 1.0 718 728 719 729 return frac_max 730 720 731 721 732 def get_intercept(q, q_0, q_1): … … 727 738 728 739 729 A B 740 A B 730 741 +-----------+--------+ <--- pixel size 731 0 1 742 0 1 732 743 Q_0 -------- Q ----- Q_1 <--- equivalent Q range 733 744 if Q_1 > Q_0, A is returned … … 738 749 if q_1 > q_0: 739 750 if (q > q_0 and q <= q_1): 740 return (q - q_0) /(q_1 - q_0)751 return (q - q_0) / (q_1 - q_0) 741 752 else: 742 753 if (q > q_1 and q <= q_0): 743 return (q - q_1) /(q_0 - q_1)754 return (q - q_1) / (q_0 - q_1) 744 755 return None 745 756 757 746 758 class _Sector: 747 759 """ 748 760 Defines a sector region on a 2D data set. 749 761 The sector is defined by r_min, r_max, phi_min, phi_max, 750 and the position of the center of the ring 762 and the position of the center of the ring 751 763 where phi_min and phi_max are defined by the right 752 764 and left lines wrt central line 753 and phi_max could be less than phi_min. 765 and phi_max could be less than phi_min. 754 766 755 Phi is defined between 0 and 2*pi in anti-clockwise 767 Phi is defined between 0 and 2*pi in anti-clockwise 756 768 starting from the x- axis on the left-hand side 757 769 """ … … 763 775 self.nbins = nbins 764 776 765 766 777 def _agv(self, data2D, run='phi'): 767 778 """ … … 781 792 q_data = data2D.q_data[numpy.isfinite(data2D.data)] 782 793 err_data = data2D.err_data[numpy.isfinite(data2D.data)] 783 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 794 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 784 795 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 785 796 dq_data = None … … 787 798 # Get the dq for resolution averaging 788 799 if data2D.dqx_data != None and data2D.dqy_data != None: 789 # The pinholes and det. pix contribution present 800 # The pinholes and det. pix contribution present 790 801 # in both direction of the 2D which must be subtracted when 791 802 # converting to 1D: dq_overlap should calculated ideally at 792 # q = 0. 803 # q = 0. 793 804 # Extrapolate dqy(perp) at q = 0 794 805 z_max = max(data2D.q_data) 795 806 z_min = min(data2D.q_data) 796 807 x_max = data2D.dqx_data[data2D.q_data[z_max]] 797 x_min = data2D.dqx_data[data2D.q_data[z_min]] 808 x_min = data2D.dqx_data[data2D.q_data[z_min]] 798 809 y_max = data2D.dqy_data[data2D.q_data[z_max]] 799 y_min = data2D.dqy_data[data2D.q_data[z_min]] 810 y_min = data2D.dqy_data[data2D.q_data[z_min]] 800 811 # Find qdx at q = 0 801 812 dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) … … 803 814 if dq_overlap_x > min(data2D.dqx_data): 804 815 dq_overlap_x = min(data2D.dqx_data) 805 dq_overlap_x *= dq_overlap_x 816 dq_overlap_x *= dq_overlap_x 806 817 # Find qdx at q = 0 807 818 dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) … … 812 823 dq_overlap_y *= dq_overlap_y 813 824 814 dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 825 dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 815 826 if dq_overlap < 0: 816 827 dq_overlap = y_min … … 819 830 # def; dqx_data = dq_r dqy_data = dq_phi 820 831 # Convert dq 2D to 1D here 821 dqx = dqx_data * dqx_data 832 dqx = dqx_data * dqx_data 822 833 dqy = dqy_data * dqy_data 823 834 dq_data = numpy.add(dqx, dqy) … … 827 838 x = numpy.zeros(self.nbins) 828 839 y = numpy.zeros(self.nbins) 829 y_err = numpy.zeros(self.nbins) 830 x_err = numpy.zeros(self.nbins) 840 y_err = numpy.zeros(self.nbins) 841 x_err = numpy.zeros(self.nbins) 831 842 y_counts = numpy.zeros(self.nbins) 832 843 … … 837 848 q_data_max = numpy.max(q_data) 838 849 839 for n in range(len(data)): 850 for n in range(len(data)): 840 851 frac = 0 841 852 … … 848 859 849 860 # phi-value of the pixel (j,i) 850 phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi 861 phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi 851 862 852 863 ## No need to calculate the frac when all data are within range 853 864 if self.r_min <= q_value and q_value <= self.r_max: 854 frac = 1 865 frac = 1 855 866 if frac == 0: 856 867 continue 857 868 #In case of two ROIs (symmetric major and minor regions)(for 'q2') 858 if run.lower() =='q2':859 ## For minor sector wing 869 if run.lower() == 'q2': 870 ## For minor sector wing 860 871 # Calculate the minor wing phis 861 872 phi_min_minor = flip_phi(phi_min - Pi) … … 869 880 phi_value < phi_max_minor) 870 881 871 #For all cases(i.e.,for 'q', 'q2', and 'phi') 872 #Find pixels within ROI 873 if phi_min > phi_max: 882 #For all cases(i.e.,for 'q', 'q2', and 'phi') 883 #Find pixels within ROI 884 if phi_min > phi_max: 874 885 is_in = is_in or (phi_value > phi_min or \ 875 phi_value < phi_max) 886 phi_value < phi_max) 876 887 else: 877 888 is_in = is_in or (phi_value >= phi_min and \ … … 879 890 880 891 if not is_in: 881 frac = 0 892 frac = 0 882 893 if frac == 0: 883 894 continue 884 895 # Check which type of averaging we need 885 if run.lower() == 'phi': 896 if run.lower() == 'phi': 886 897 temp_x = (self.nbins) * (phi_value - self.phi_min) 887 898 temp_y = (self.phi_max - self.phi_min) … … 892 903 i_bin = int(math.floor(temp_x / temp_y)) 893 904 894 # Take care of the edge case at phi = 2pi. 895 if i_bin == self.nbins: 896 i_bin = 905 # Take care of the edge case at phi = 2pi. 906 if i_bin == self.nbins: 907 i_bin = self.nbins - 1 897 908 898 ## Get the total y 909 ## Get the total y 899 910 y[i_bin] += frac * data_n 900 911 x[i_bin] += frac * q_value … … 907 918 908 919 if dq_data != None: 909 # To be consistent with dq calculation in 1d reduction, 910 # we need just the averages (not quadratures) because 911 # it should not depend on the number of the q points 920 # To be consistent with dq calculation in 1d reduction, 921 # we need just the averages (not quadratures) because 922 # it should not depend on the number of the q points 912 923 # in the qr bins. 913 924 x_err[i_bin] += frac * dq_data[n] … … 923 934 # The type of averaging: phi,q2, or q 924 935 # Calculate x[i]should be at the center of the bin 925 if run.lower() =='phi':936 if run.lower() == 'phi': 926 937 x[i] = (self.phi_max - self.phi_min) / self.nbins * \ 927 938 (1.0 * i + 0.5) + self.phi_min 928 939 else: 929 # We take the center of ring area, not radius. 940 # We take the center of ring area, not radius. 930 941 # This is more accurate than taking the radial center of ring. 931 942 #delta_r = (self.r_max - self.r_min) / self.nbins … … 934 945 #x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2) 935 946 x[i] = x[i] / y_counts[i] 936 y_err[y_err ==0] = numpy.average(y_err)947 y_err[y_err == 0] = numpy.average(y_err) 937 948 idx = (numpy.isfinite(y) & numpy.isfinite(y_err)) 938 949 if x_err != None: … … 941 952 d_x = None 942 953 if not idx.any(): 943 msg = "Average Error: No points inside sector of ROI to average..." 954 msg = "Average Error: No points inside sector of ROI to average..." 944 955 raise ValueError, msg 945 956 #elif len(y[idx])!= self.nbins: … … 948 959 return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x) 949 960 961 950 962 class SectorPhi(_Sector): 951 963 """ … … 963 975 :return: Data1D object 964 976 """ 965 966 977 return self._agv(data2D, 'phi') 978 967 979 968 980 class SectorQ(_Sector): … … 972 984 973 985 A sector is defined by r_min, r_max, phi_min, phi_max. 974 r_min, r_max, phi_min, phi_max >0. 986 r_min, r_max, phi_min, phi_max >0. 975 987 The number of bin in Q also has to be defined. 976 988 """ … … 984 996 """ 985 997 return self._agv(data2D, 'q2') 998 986 999 987 1000 class Ringcut(object): … … 993 1006 The data returned is the region inside the ring 994 1007 995 Phi_min and phi_max should be defined between 0 and 2*pi 1008 Phi_min and phi_max should be defined between 0 and 2*pi 996 1009 in anti-clockwise starting from the x- axis on the left-hand side 997 1010 """ 998 def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0 1011 def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0): 999 1012 # Minimum radius 1000 1013 self.r_min = r_min … … 1006 1019 self.center_y = center_y 1007 1020 1008 1009 1021 def __call__(self, data2D): 1010 1022 """ … … 1020 1032 1021 1033 # Get data 1022 qx_data = data2D.qx_data 1034 qx_data = data2D.qx_data 1023 1035 qy_data = data2D.qy_data 1024 1036 mask = data2D.mask 1025 1037 q_data = numpy.sqrt(qx_data * qx_data + qy_data * qy_data) 1026 #q_data_max = numpy.max(q_data)1027 1038 1028 1039 # check whether or not the data point is inside ROI … … 1051 1062 1052 1063 :param data2D: Data2D object 1053 :return: mask, 1d array (len = len(data)) 1064 :return: mask, 1d array (len = len(data)) 1054 1065 with Trues where the data points are inside ROI, otherwise False 1055 1066 """ … … 1060 1071 def _find(self, data2D): 1061 1072 """ 1062 Find a rectangular 2D region of interest. 1063 1064 :param data2D: Data2D object 1065 1066 :return: out, 1d array (length = len(data)) 1073 Find a rectangular 2D region of interest. 1074 1075 :param data2D: Data2D object 1076 1077 :return: out, 1d array (length = len(data)) 1067 1078 with Trues where the data points are inside ROI, otherwise Falses 1068 1079 """ 1069 1080 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 1070 1081 raise RuntimeError, "Boxcut take only plottable_2D objects" 1071 # Get qx_ and qy_data 1082 # Get qx_ and qy_data 1072 1083 qx_data = data2D.qx_data 1073 1084 qy_data = data2D.qy_data … … 1080 1091 return (outx & outy) 1081 1092 1093 1082 1094 class Sectorcut(object): 1083 1095 """ 1084 1096 Defines a sector (major + minor) region on a 2D data set. 1085 1097 The sector is defined by phi_min, phi_max, 1086 where phi_min and phi_max are defined by the right 1087 and left lines wrt central line. 1098 where phi_min and phi_max are defined by the right 1099 and left lines wrt central line. 1088 1100 1089 Phi_min and phi_max are given in units of radian 1101 Phi_min and phi_max are given in units of radian 1090 1102 and (phi_max-phi_min) should not be larger than pi 1091 1103 """ … … 1100 1112 :param data2D: Data2D object 1101 1113 1102 :return: mask, 1d array (len = len(data)) 1114 :return: mask, 1d array (len = len(data)) 1103 1115 1104 1116 with Trues where the data points are inside ROI, otherwise False … … 1110 1122 def _find(self, data2D): 1111 1123 """ 1112 Find a rectangular 2D region of interest. 1113 1114 :param data2D: Data2D object 1115 1116 :return: out, 1d array (length = len(data)) 1124 Find a rectangular 2D region of interest. 1125 1126 :param data2D: Data2D object 1127 1128 :return: out, 1d array (length = len(data)) 1117 1129 1118 1130 with Trues where the data points are inside ROI, otherwise Falses 1119 1131 """ 1120 1132 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 1121 raise RuntimeError, "Sectorcut take only plottable_2D objects" 1133 raise RuntimeError, "Sectorcut take only plottable_2D objects" 1122 1134 Pi = math.pi 1123 1135 # Get data 1124 1136 qx_data = data2D.qx_data 1125 qy_data = data2D.qy_data 1137 qy_data = data2D.qy_data 1126 1138 phi_data = numpy.zeros(len(qx_data)) 1127 1139 … … 1131 1143 # Get the min and max into the region: -pi <= phi < Pi 1132 1144 phi_min_major = flip_phi(self.phi_min + Pi) - Pi 1133 phi_max_major = flip_phi(self.phi_max + Pi) - Pi 1145 phi_max_major = flip_phi(self.phi_max + Pi) - Pi 1134 1146 # check for major sector 1135 1147 if phi_min_major > phi_max_major: … … 1146 1158 if phi_min_minor > phi_max_minor: 1147 1159 out_minor = (phi_min_minor <= phi_data) + \ 1148 (phi_max_minor >= phi_data) 1160 (phi_max_minor >= phi_data) 1149 1161 else: 1150 1162 out_minor = (phi_min_minor <= phi_data) & \ 1151 (phi_max_minor >= phi_data) 1163 (phi_max_minor >= phi_data) 1152 1164 out = out_major + out_minor 1153 1165 1154 1166 return out 1155 1156 if __name__ == "__main__":1157 1158 from loader import Loader1159 1160 1161 d = Loader().load('test/MAR07232_rest.ASC')1162 #d = Loader().load('test/MP_New.sans')1163 1164 1165 r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)1166 o = r(d)1167 1168 s = Ring(r_min=.000001, r_max=.01)1169 p = s(d)1170 1171 for i in range(len(o.x)):1172 print o.x[i], o.y[i], o.dy[i]1173 1174 1175
Note: See TracChangeset
for help on using the changeset viewer.
