Changeset 811 for sans/Dev/trunk/NCNR_User_Procedures/Reduction
- Timestamp:
- Jun 21, 2011 1:07:46 PM (12 years ago)
- Location:
- sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/CatVSTable.ipf
r707 r811 151 151 //write to notebook that file was not found 152 152 //if string is not a number, report the error 153 if( str2num(partialName) == NaN)153 if(numtype(str2num(partialName)) == 2) 154 154 str = "this file was not found: "+partialName+"\r\r" 155 155 //Notebook CatWin,font="Times",fsize=12,text=str … … 508 508 //write to notebook that file was not found 509 509 //if string is not a number, report the error 510 if( str2num(partialName) == NaN)510 if(numtype(str2num(partialName)) == 2) 511 511 str = "this file was not found: "+partialName+"\r\r" 512 512 Notebook CatWin,font="Times",fsize=12,text=str … … 677 677 //write to notebook that file was not found 678 678 //if string is not a number, report the error 679 if( str2num(partialName) == NaN)679 if(numtype(str2num(partialName)) == 2) 680 680 str = "this file was not found: "+partialName+"\r\r" 681 681 Notebook CatWin,font="Times",fsize=12,text=str -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r800 r811 20 20 // whatever XOP crashing was happining during threading is really unlikely to be from the RNG 21 21 // 22 // -- June 2010 - calls from different threads to the same RNG really seems tocause a crash. Probably as soon22 // -- June 2010 - calls from different threads to the same RNG will absolutely cause a crash. Probably as soon 23 23 // as the different threads try to call at the same time. Found this out by accident doing the 24 24 // wavelength spread. Each thread called ran3 at that point, and the crash came quickly. Went … … 85 85 // --- TO ADD --- 86 86 // X- wavelength distribution = another RNG to select the wavelength 87 // ---- doneJun 2010, approximating the wavelength distribution as a Gaussian, based on the triangular87 // DONE Jun 2010, approximating the wavelength distribution as a Gaussian, based on the triangular 88 88 // FWHM. Wavelength distribution added to XOP too, and now very accurately matches the shape of the 1D 89 89 // simulation. … … 96 96 // data can be plotted, but only by hand right now. EC and blocked beam are combined. 97 97 // 98 // -- divergence / size of the incoming beam. Currently everything is parallel, and anything that is transmitted 98 99 // X- divergence / size of the incoming beam. Currently everything is parallel, and anything that is transmitted 99 100 // simply ends up in (xCtr,yCtr), and the "real" profile of the beam is not captured. 100 101 // DONE, 21 JUN 11 SRK 102 // point is picked in source aperture, then sample aperture. This sets the initial direction of the neutron. 103 // Gravity is included, lowering every neutron by yg_d. This affects qy only.Simulated beam center makes sense 104 // now both in size and xy positon. Gravity effects can be seen in the beam spot itself (fall + spread) and in the scattering 105 // pattern at extreme cases (sharp peaks, wide spread, long SDD, long wavelength, small source aperture) 101 106 102 107 … … 138 143 nthreads = 1 139 144 #endif 140 141 // nthreads = 1 145 146 147 // nthreads = 2 148 149 142 150 NVAR mt=root:myGlobals:gThreadGroupID 143 151 mt = ThreadGroupCreate(nthreads) … … 382 390 Variable NMultipleCoherent,NCoherentEvents 383 391 Variable deltaLam,v1,v2,currWavelength,rsq,fac //for simulating wavelength distribution 392 Variable SSD, sourAp, souXX, souYY, magn //source-to-sample, and source Ap radius for initlal trajectory 393 394 Variable vz_1 = 3.956e5 //velocity [cm/s] of 1 A neutron 395 Variable g = 981.0 //gravity acceleration [cm/s^2] 396 Variable yg_d //fall due to gravity 397 384 398 385 399 // don't set to other than one here. Detector efficiency is handled outside, only passing the number of … … 399 413 sig_sas = inputWave[10] 400 414 deltaLam = inputWave[11] 415 SSD = inputWave[12] // in cm, like SDD 416 sourAp = inputWave[13] // radius, in cm, like r1 and r2 417 418 // print SSD, sourAp 401 419 402 420 // SetRandomSeed 0.1 //to get a reproduceable sequence … … 462 480 // NOW, start the loop, throwing neutrons at the sample. 463 481 do 464 Vx = 0.0 // Initialize direction vector.465 Vy = 0.0466 Vz = 1.0482 // Vx = 0.0 // Initialize direction vector. 483 // Vy = 0.0 484 // Vz = 1.0 467 485 468 486 Theta = 0.0 // Initialize scattering angle. … … 474 492 incoherentEvent = 0 475 493 coherentEvent = 0 476 494 495 // pick point in source aperture area 496 do // Makes sure position is within circle. 497 ran = abs(enoise(1)) //[0,1] 498 souxx = 2.0*sourAp*(Ran-0.5) //X beam position of neutron entering sample. 499 ran = abs(enoise(1)) //[0,1] 500 souyy = 2.0*sourAp*(Ran-0.5) //Y beam position ... 501 RR = SQRT(souxx*souxx+souyy*souyy) //Radial position of neutron in incident beam. 502 while(rr>sourAp) 503 504 // pick point in sample aperture 477 505 do // Makes sure position is within circle. 478 506 ran = abs(enoise(1)) //[0,1] … … 497 525 currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength 498 526 527 // if(n1 == 1) 528 // Print "lambda ", currWavelength 529 // endif 530 // 531 magn = sqrt((souXX - xx)^2 + (souYY - yy)^2 + ssd^2) 532 Vx = (souXX - xx)/magn // Initialize direction vector. 533 Vy = (souYY - yy)/magn 534 Vz = (ssd - 0)/magn 535 536 // Vx = 0.0 // Initialize direction vector. 537 // Vy = 0.0 538 // Vz = 1.0 539 540 541 // 542 // if(n1 == 1) 543 // Print "vx, vy, vz, mag",vx,vy,vz,sqrt(vx^2+vy^2+vz^2) 544 // endif 545 499 546 500 547 do //Scattering Loop, will exit when "done" == 1 … … 505 552 err = NewDirection(vx,vy,vz,Theta,Phi) //vx,vy,vz is updated, theta, phi unchanged by function 506 553 554 // if(vx != 0) 555 // Print "vx, vy, vz, mag" = vx,vy,vz,sqrt(vx^2+vy^2+vz^2) 556 // endif 557 558 507 559 //X,Y,Z-POSITION OF SCATTERING EVENT. 508 560 xx += ll*vx … … 578 630 Endif 579 631 632 // calculate fall due to gravity (in cm) (note that it is negative) 633 YG_d = -0.5*g*SDD*(SSD+SDD)*(currWavelength/vz_1)^2 634 635 // yg_d=0 636 637 if(n1 == 1) 638 // Print "gravity fall (cm) = ",Yg_d 639 Print "gravity fall at mean lam (cm) = ",-0.5*g*SDD*(SSD+SDD)*(wavelength/vz_1)^2 640 endif 641 580 642 if(index != 0) //the neutron interacted at least once, figure out where it ends up 581 643 … … 590 652 591 653 // is it on the detector? 592 FindPixel(testQ,testPhi,currWavelength, sdd,pixSize,xCtr,yCtr,xPixel,yPixel)654 FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 593 655 594 656 if(xPixel != -1 && yPixel != -1) 595 657 //if(index==1) // only the single scattering events 658 if( xPixel > 127 || yPixel > 127) 659 print "error XY=",xPixel,yPixel 660 endif 596 661 MC_linear_data[xPixel][yPixel] += 1 //this is the total scattering, including multiple scattering 597 662 //endif … … 641 706 642 707 // then it must be a transmitted neutron 643 // don't need to calculate, just increment the proper counters708 // Now, calculate where it lands 644 709 645 MC_linear_data[xCtr+xx/pixsize][yCtr+yy/pixsize] += 1 646 isOn += 1 647 nt[0] += 1 710 Theta_z = acos(Vz) // Angle WITH respect to z axis. 711 testQ = 2*pi*sin(theta_z)/currWavelength 712 713 testPhi = MC_FindPhi(Vx,Vy) //use the exiting phi value as defined by Vx and Vy 714 715 // is it on the detector? 716 FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 717 718 if(xPixel != -1 && yPixel != -1) 719 //if(index==1) // only the single scattering events 720 if( xPixel > 127 || yPixel > 127) 721 print "error XY=",xPixel,yPixel 722 endif 723 MC_linear_data[xPixel][yPixel] += 1 //this is the total scattering, including multiple scattering 724 //endif 725 isOn += 1 // neutron that lands on detector 726 endif 727 728 // MC_linear_data[xCtr+xx/pixsize][yCtr+yy/pixsize] += 1 729 // isOn += 1 730 If(theta_z < theta_max) 731 //Choose index for scattering angle array. 732 //IND = NINT(THETA_z/DTH + 0.4999999) 733 ind = round(THETA_z/DTH + 0.4999999) //round is eqivalent to nint() 734 NT[ind] += 1 //Increment bin for angle. 735 endif 736 //nt[0] += 1 // may not be at zero angle anmore 648 737 649 738 endif //if interacted … … 652 741 while(n1 < imon) 653 742 654 // Print "Monte Carlo Done" 655 results[0] = n1 656 results[1] = n2 657 results[2] = isOn 658 results[3] = NScatterEvents //sum of # of times that neutrons scattered (coh+incoh) 659 results[4] = NSingleCoherent //# of events that are single, coherent 660 results[5] = NMultipleCoherent //# of scattered neutrons that are coherently scattered more than once 661 results[6] = NMultipleScatter //# of scattered neutrons that are scattered more than once (coh and/or incoh) 662 results[7] = NCoherentEvents //# of scattered neutrons that are scattered coherently one or more times 663 743 Print "Non-XOP Monte Carlo Done" 744 // results[0] = n1 745 // results[1] = n2 746 // results[2] = isOn 747 // results[3] = NScatterEvents //sum of # of times that neutrons scattered (coh+incoh) 748 // results[4] = NSingleCoherent //# of events that are single, coherent 749 // results[5] = NMultipleCoherent //# of scattered neutrons that are coherently scattered more than once 750 // results[6] = NMultipleScatter //# of scattered neutrons that are scattered more than once (coh and/or incoh) 751 // results[7] = NCoherentEvents //# of scattered neutrons that are scattered coherently one or more times 752 753 754 Variable xc,yc 755 xc=inputWave[3] 756 yc=inputWave[4] 757 results[0] = inputWave[9]+inputWave[10] //total XS 758 results[1] = inputWave[10] //SAS XS 759 results[2] = n2 //number that interact n2 760 results[3] = isOn - MC_linear_data[xc][yc] //# reaching detector minus Q(0) 761 results[4] = NScatterEvents/n2 //avg# times scattered 762 results[5] = NSingleCoherent/NCoherentEvents //single coherent fraction 763 results[6] = NMultipleCoherent/NCoherentEvents //multiple coherent fraction 764 results[7] = NMultipleScatter/n2 //multiple scatter fraction 765 results[8] = (n1-n2)/n1 //transmitted fraction 766 767 664 768 // Print "# absorbed = ",n3 665 769 … … 790 894 End 791 895 792 ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 793 Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel 896 897 898 // xCtr and yCtr here are the "optical" center of the detector ~(65,65) and the full fall due to 899 // gravity is calculated from this horizontal axis 900 // 901 ThreadSafe Function FindPixel(testQ,testPhi,lam,yg_d,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 902 Variable testQ,testPhi,lam,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel 794 903 795 904 Variable theta,dy,dx,qx,qy … … 798 907 qy = testQ*sin(testPhi) 799 908 909 // correct qy for gravity 910 // qy = 4*pi/lam * (theta/2) 911 qy += 4*pi/lam*(yg_d/sdd/2) 912 913 800 914 //convert qx,qy to pixel locations relative to # of pixels x, y from center 801 915 theta = 2*asin(qy*lam/4/pi) … … 811 925 812 926 //if on detector, return xPix and yPix values, otherwise -1 813 if(yPixel > pixelsY || yPixel < 0)927 if(yPixel >= pixelsY || yPixel < 0) 814 928 yPixel = -1 815 929 endif 816 if(xPixel > pixelsX || xPixel < 0)930 if(xPixel >= pixelsX || xPixel < 0) 817 931 xPixel = -1 818 932 endif … … 1332 1446 inputWave[10] = sig_sas 1333 1447 inputWave[11] = deltaLam 1334 // inputWave[] 12-14 are currently unused 1448 1449 inputWave[12] = sourceToSampleDist() // function returns dist in cm 1450 inputWave[13] = sourceApertureDiam()/2 // function returns diam in cm, this is radius 1451 // inputWave[] 14 are currently unused 1335 1452 1336 1453 linear_data = 0 //initialize … … 1402 1519 1403 1520 Variable rad=beamstopDiam()/2 //beamstop radius in cm 1521 // the function detectorOffset() does not take gravity into account, and resets the beam center to (64.5,64.5) 1522 // after the MC simulation is done (displayConfigurationText is always called) 1523 // 1524 // so this (locally) changes the beam center, just for the correct placement of the beam stop 1525 // the data files will have the wrong 1526 Variable vz_1 = 3.956e5 //velocity [cm/s] of 1 A neutron 1527 Variable g = 981.0 //gravity acceleration [cm/s^2] 1528 Variable yg_d,ssd //fall due to gravity 1529 1530 ssd = sourceToSampleDist() 1531 YG_d = -0.5*g*SDD*(SSD+SDD)*(wavelength/vz_1)^2 // fall in cm (negative value) 1532 1533 xCtr = 64.5 + round(2*rw[19]) // I'm always off by one for some reason, so start at 65.5? 1534 yCtr = 65 + yg_d/0.5 // this will lower the beam center 1535 // rw[16] = xCtr 1536 // rw[17] = yCtr 1537 1538 Print "Gravity q* = ",-2*pi/wavelength*2*yg_d/sdd 1539 1540 1404 1541 if(MC_BS_in) 1405 1542 rad /= 0.5 //convert cm to pixels … … 1407 1544 Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data 1408 1545 WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask 1546 1409 1547 tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1 //behind beamstop = 0, away = 1 1410 1548 … … 1458 1596 // re-average the 2D data 1459 1597 S_CircularAverageTo1D("SAS") 1460 1598 1461 1599 // put the new result into the simulation folder 1462 1600 Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation") … … 1474 1612 endif 1475 1613 endif 1476 1477 1614 1478 1615 return(0) -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultipleReduce.ipf
r764 r811 607 607 //write to notebook that file was not found 608 608 //if string is not a number, report the error 609 if( str2num(partialName) == NaN)609 if(numtype(str2num(partialName)) == 2) 610 610 str = "this file was not found: "+partialName+"\r\r" 611 611 //Notebook CatWin,font="Times",fsize=12,text=str … … 664 664 loc = strsearch(labels[ii], findThisStr, 0 ,2) //2==case insensitive, but Igor 5 specific 665 665 if(loc != -1) 666 Print "Remove w[ii] = ", num," ",labels[ii]666 Print "Remove w[ii] = ",runnum[ii]," ",labels[ii] 667 667 DeletePoints ii, 1, filenames,suffix,labels,sdd,runnum,isTrans 668 668 endif -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NCNR_DataReadWrite.ipf
r806 r811 2004 2004 WAVE w=tmp_data 2005 2005 2006 // check for data values that are too large. the maximum VAX compressed data value is 2767000 2007 // 2008 WaveStats/Q w 2009 if(V_max > 2767000) 2010 Abort "Some individual pixel values are > 2767000 and the data can't be saved in VAX format" 2011 Endif 2006 2012 2007 2013 //check each wave … … 2044 2050 tmpFile=0 2045 2051 2046 Make/O/W/N=16401 dataWRecMarkers 2052 // Make/O/W/N=16401 dataWRecMarkers // don't convert to 16 bit here, rather write to file as 16 bit later 2053 Make/O/I/N=16401 dataWRecMarkers 2047 2054 AddRecordMarkers(w,dataWRecMarkers) 2048 2055 2049 2056 // need to re-compress?? maybe never a problem, but should be done for the odd case 2050 2057 dataWRecMarkers = CompressI4toI2(dataWRecMarkers) //unless a pixel value is > 32767, the same values are returned -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NCNR_Utils.ipf
r801 r811 320 320 321 321 ///////////////////////////////////////////////// 322 /// 323 ////// this is all experimental, inclusion of gravity effect into the parallel component324 //325 Variable yg_d,acc,sdd,ssd,lambda0,DL_L,sig_l326 Variable var_qlx,var_qly,var_ql,qx,qy,sig_perp,sig_para322 ///// 323 // ////// this is all experimental, inclusion of gravity effect into the parallel component 324 //// // 325 // Variable yg_d,acc,sdd,ssd,lambda0,DL_L,sig_l 326 // Variable var_qlx,var_qly,var_ql,qx,qy,sig_perp,sig_para 327 327 // G = 981. //! ACCELERATION OF GRAVITY, CM/SEC^2 328 acc = vz_1 // 3.956E5 //! CONVERT WAVELENGTH TO VELOCITY CM/SEC 329 SDD = L2 //1317 330 SSD = L1 //1627 //cm 331 lambda0 = lambda // 15 332 DL_L = lambdaWidth //0.236 333 SIG_L = DL_L/sqrt(6) 334 YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2 335 // Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG 336 337 sig_perp = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2) 338 sig_perp = sqrt(sig_perp) 339 340 341 FindQxQy(inQ,phi,qx,qy) 342 343 VAR_QLY = SIG_L^2 * (QY+4*PI*YG_d/(SDD*LAMBDA0))^2 344 VAR_QLX = (SIG_L*QX)^2 345 VAR_QL = VAR_QLY + VAR_QLX //! WAVELENGTH CONTRIBUTION TO VARIANCE 346 sig_para = (sig_perp^2 + VAR_QL)^0.5 347 348 //return values PBR 349 SigmaQX = sig_para 350 SigmaQy = sig_perp 351 328 // acc = vz_1 // 3.956E5 //! CONVERT WAVELENGTH TO VELOCITY CM/SEC 329 // SDD = L2 //1317 330 // SSD = L1 //1627 //cm 331 // lambda0 = lambda // 15 332 // DL_L = lambdaWidth //0.236 333 // SIG_L = DL_L/sqrt(6) 334 // YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2 335 /////// Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG 352 336 // 337 // sig_perp = kap*kap/12 * (3*(S1/L1)^2 + 3*(S2/LP)^2 + (proj_DDet/L2)^2) 338 // sig_perp = sqrt(sig_perp) 339 // 340 // 341 // FindQxQy(inQ,phi,qx,qy) 342 // 343 // VAR_QLY = SIG_L^2 * (QY+4*PI*YG_d/(SDD*LAMBDA0))^2 344 // VAR_QLX = (SIG_L*QX)^2 345 // VAR_QL = VAR_QLY + VAR_QLX //! WAVELENGTH CONTRIBUTION TO VARIANCE 346 // sig_para = (sig_perp^2 + VAR_QL)^0.5 347 // 348 /////// return values PBR 349 // SigmaQX = sig_para 350 // SigmaQy = sig_perp 351 // 352 ////// 353 353 354 354 results = "success" -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NSORT.ipf
r741 r811 1535 1535 //write to notebook that file was not found 1536 1536 //if string is not a number, report the error 1537 if( str2num(partialName) == NaN)1537 if(numtype(str2num(partialName)) == 2) 1538 1538 str = "this file was not found: "+partialName+"\r\r" 1539 1539 //Notebook CatWin,font="Times",fsize=12,text=str -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf
r800 r811 294 294 // real wave 295 295 rw = 0 296 rw[16] = 6 4// beamcenter X (pixels)296 rw[16] = 65 // beamcenter X (pixels) 297 297 rw[17] = 64 // beamcenter Y 298 298 … … 754 754 NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo 755 755 756 if( CB_Struct.eventMod == 2) //if the shift key is down- go to 2D mode756 if((CB_Struct.eventMod & 2^1) != 0 ) //if the shift key is down (bit 1)- go to 2D mode 757 757 do2D = 1 758 758 endif … … 946 946 if(doMonteCarlo == 1) 947 947 //2D simulation (in MultiScatter_MonteCarlo_2D.ipf) 948 949 // may want to take this back out--- 20 JUN 2011 SRK 950 // this call to detectorOffset() sets the (xCtr,yCtr) in RealsRead to default values of (64.5,64.5) 951 // -- if user changes the values in RealsRead to be able to do a proper average, the next sim starts wtih the 952 // wrong values... 953 954 detectorOffset() 948 955 949 956 Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs) … … 1799 1806 rw[16] = 64 + round(2*rw[19]) + 0.5 //approximate beam X is 64 w/no offset, 114 w/25 cm offset 1800 1807 rw[17] = 64 + 0.5 //typical value 1801 1808 1802 1809 return(val) 1803 1810 end
Note: See TracChangeset
for help on using the changeset viewer.