Ignore:
Timestamp:
Mar 18, 2013 10:10:23 AM (10 years ago)
Author:
srkline
Message:

-adjusted 2D simulation to get a proper beam center into the file. Previously was off by 1 pixel in x.

  • conditions to keep locaed pixel in simulation to range [0,127]
  • in 1D sim, -ve values from noise are replaced w/ zero data value and error of one

AutoFit? allows epsilon wave again, checkbox is un-hidden.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/CircSectAve.ipf

    r886 r902  
    252252                                var = aveisq-avesq 
    253253                                if(var<=0) 
    254                                         sigave[kk-1] = small_num 
     254//                                      sigave[kk-1] = small_num 
     255                                        sigave[kk-1] = large_num                //if there are zero counts, make error large_num to be a warning flag 
    255256                                else 
    256257                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1)) 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MC_SimulationScripting.ipf

    r901 r902  
    296296        cw = {1e-18,3.346,0,9,0} 
    297297        Sim_SetIncohXS(0)                                                                       // incoh XS = 0 for empty beam 
     298        Sim_SetThickness(0.1)                                                   // thickness (cm) to give proper level of scattering 
     299 
    298300                 
    299301        totalTime = Sim_RunTrans_2D(confList,ctTimeList,titleStr,runIndex) 
    … …  
    10101012End 
    10111013 
    1012 Proc Sim_SetupRunPanel() : Panel 
    1013         PauseUpdate; Silent 1           // building window... 
    1014         NewPanel /W=(399,44,764,386) 
    1015         ModifyPanel cbRGB=(56639,28443,117) 
    1016         ShowTools/A 
    1017 EndMacro 
     1014//Proc Sim_SetupRunPanel() : Panel 
     1015//      PauseUpdate; Silent 1           // building window... 
     1016//      NewPanel /W=(399,44,764,386) 
     1017//      ModifyPanel cbRGB=(56639,28443,117) 
     1018//      ShowTools/A 
     1019//EndMacro 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r901 r902  
    1515// - not for the casual user at all. 
    1616 
    17  
     17// 
     18// - Mar 2013 - modified the 1D simulation to properly account for very low counts. In a real experiment, 2D data 
     19//    is collected, and there can be lots of zeros on the detector. Then averaging out to 1D, they remain zero. Arbitrarily 
     20//    they are assigned an error of 1. Previously, a (low) but non-zero, non-integer count value was assigned to the 
     21//    1D intensity. Then when noise was added, the resulting I(q) could be negative. Now the negative intensities are 
     22//    replaced with zero, with an error of 1 (1e-10 was previously used as the error, but this put a ridiculous 
     23//    weighting on that point during fitting, making any fit impossible) -- But is this correct? 
     24// 
     25// -- so take it back out for now --?? Search for MAR 2013 in Simulate_1D() - two places there 
     26// 
     27// 
    1828 
    1929// the RNG issue is really not worth the effort. multiple copies with different RNG is as good as I need. Plus, 
    … …  
    716726                                        if(xPixel != -1 && yPixel != -1) 
    717727                                                //if(index==1)  // only the single scattering events 
    718                                                         if( xPixel > 127 || yPixel > 127) 
    719                                                                 print "error XY=",xPixel,yPixel 
     728                                                        if( xPixel > 127 || yPixel > 127 || xPixel < 0 || yPixel < 0) 
     729                                                                print "error XY loc2 =",xPixel,yPixel 
    720730                                                        endif 
    721731                                                        MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering 
    … …  
    900910//    passed in, xCtr and yCtr are in detector coordinates 
    901911// 
     912// everything passed in is in cm, excepty lam, in A 
     913// 
    902914ThreadSafe Function FindPixel(testQ,testPhi,lam,yg_d,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
    903915        Variable testQ,testPhi,lam,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel 
    904916 
    905         Variable theta,dy,dx,qx,qy,qz 
    906          
    907         theta = 2*asin(testQ*lam/4/pi) 
     917        Variable theta,dy,dx,qx,qy,qz,theta_qy,theta_qx,two_theta,dist 
     918         
     919        two_theta = 2*asin(testQ*lam/4/pi) 
    908920         
    909921        //decompose to qx,qy 
    910         qx = testQ*cos(theta/2)*cos(testPhi) 
    911         qy = testQ*cos(theta/2)*sin(testPhi) 
    912         qz = testQ*sin(theta/2) 
     922//      qx = testQ*cos(theta/2)*cos(testPhi) 
     923//      qy = testQ*cos(theta/2)*sin(testPhi) 
     924//      qz = testQ*sin(theta/2) 
    913925         
    914926// correct qy for gravity 
    915927// qy = 4*pi/lam * (theta/2) 
    916         qy += 4*pi/lam*(yg_d/sdd/2) 
    917          
    918          
    919         //convert qx,qy to pixel locations relative to # of pixels x, y from center 
    920         theta = 2*asin(qy*lam/4/pi) 
    921         dy = sdd*tan(theta) 
    922 //      dy +=  
    923         yPixel = round(yCtr + dy/pixSize) 
    924          
    925         theta = 2*asin(qx*lam/4/pi) 
    926         dx = sdd*tan(theta) 
    927         xPixel = round(xCtr + dx/pixSize) 
     928//      qy += 4*pi/lam*(yg_d/sdd/2) 
     929 
     930         
     931        dist = sdd*tan(two_theta)               //hypot in xy plane 
     932 
     933        dx = dist*cos(testPhi) 
     934        dy = dist*sin(testPhi) 
     935         
     936        xPixel = dx/pixSize + xCtr 
     937        yPixel = dy/pixSize + yCtr + yg_d/pixSize               //shift down due to gravity 
     938         
     939        xPixel = round(xPixel) 
     940        yPixel = round(yPixel) 
     941         
     942//      //convert qx,qy to pixel locations relative to # of pixels x, y from center 
     943//      theta = 2*asin(qy*lam/4/pi) 
     944//      dy = sdd*tan(theta) 
     945//      yPixel = abs(round(yCtr + dy/pixSize)) 
     946//       
     947//      theta = 2*asin(qx*lam/4/pi) 
     948//      dx = sdd*tan(theta) 
     949//      xPixel = abs(round(xCtr + dx/pixSize)) 
    928950 
    929951        NVAR pixelsX = root:myGlobals:gNPixelsX 
    930952        NVAR pixelsY = root:myGlobals:gNPixelsY 
    931          
     953 
     954// convert the detector pixel coordinates to [0,127] before passing back...      
    932955        xPixel -= 1 
    933956        yPixel -= 1 
    … …  
    935958         
    936959        //if on detector, return xPix and yPix values, otherwise -1 
    937         if(yPixel >= pixelsY || yPixel < 0) 
     960        // > 127 or < 0 == off detector, no good 
     961        if(yPixel > pixelsY-1 || yPixel < 0) 
    938962                yPixel = -1 
    939963        endif 
    940         if(xPixel >= pixelsX || xPixel < 0) 
     964        if(xPixel > pixelsX-1 || xPixel < 0) 
    941965                xPixel = -1 
    942966        endif 
    … …  
    945969End 
    946970 
     971 
     972// 
     973// this is the original version before I messed with it. 
     974// 
    947975//// xCtr and yCtr here are the "optical" center of the detector ~(64,64) and the full fall due to  
    948976//// gravity is calculated from this horizontal axis 
    … …  
    16001628        YG_d = -0.5*g*SDD*(SSD+SDD)*(wavelength/vz_1)^2         // fall in cm (negative value) 
    16011629         
    1602         xCtr = 64.5      + round(2*rw[19])              // I'm always off by one for some reason, so start at 65.5? 
    1603         yCtr = 65 + yg_d/0.5                                    // this will lower the beam center 
     1630        xCtr = 64        + round(2*rw[19])              // I'm always off by one for some reason, so start at 65.5? 
     1631        yCtr = 64 + yg_d/0.5                                    // this will lower the beam center 
    16041632//      rw[16] = xCtr 
    16051633//      rw[17] = yCtr 
    … …  
    16141642                WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask 
    16151643                 
    1616                 tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1 
     1644                tmp_mask = (sqrt((p-xCtr+1)^2+(q-yCtr+1)^2) < rad) ? 0 : 1              //behind beamstop = 0, away = 1 
    16171645                 
    16181646                linear_data *= tmp_mask 
    … …  
    22072235        sigave = sqrt(aveint/nCells)            // corrected based on John's memo, from 8/9/99 
    22082236         
     2237         
    22092238        // add in random error in aveint based on the sigave 
     2239         
     2240        /// ?? MAR 2013 do I replace the negative intensities with zero? 
     2241        // if the resulting value is less than zero, replace with zero. This is what would happen if the 
     2242        // data came from 2D. Then if any of the aveint points are zero, set the associated error to == 1. This is 
     2243        // done after all of the other scaling... 
    22102244        if(addNoise) 
    22112245                aveint += gnoise(sigave) 
     2246                aveint = (aveint[p] < 0) ? 0 : aveint[p]                        // MAR 2013 -- is this the right thing to do 
    22122247        endif 
    22132248 
    … …  
    22252260                sigave /= detectorEff 
    22262261        endif 
    2227                                  
     2262         
     2263/// ?? MAR 2013 do I replace the "zero" intensity error with 1? 
     2264        if(addNoise) 
     2265                sigave = (aveint[p] == 0) ? 1 : sigave[p] 
     2266        endif                    
    22282267         
    22292268//      Simulate_1D_EmptyCell("TwoLevel_EC",aveint,qval,sigave,sigmaq,qbar,fsubs) 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NSORT.ipf

    r898 r902  
    19581958        String path=S_Path 
    19591959         
     1960////////        Make/O/D/N=(numpnts(lowW)) scale_4m 
     1961        NVAR scale12 = root:myGlobals:NSORT:gScale1_2 
     1962 
    19601963         
    19611964// this variable must exist and be set to 1 to be able to automatically name files 
    … …  
    19941997                WriteNSORTFileButton("") 
    19951998                 
     1999//////          scale_4m[ii] = scale12 
     2000                 
    19962001                Print "wrote file : ",path+saveW[ii] 
    19972002                ii+=1 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/RectAnnulAvg.ipf

    r886 r902  
    252252                                var = aveisq-avesq 
    253253                                if(var<=0) 
    254                                         sigave[kk-1] = small_num 
     254//                                      sigave[kk-1] = small_num 
     255                                        sigave[kk-1] = large_num                //if there are zero counts, make error large_num to be a warning flag 
    255256                                else 
    256257                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1)) 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf

    r901 r902  
    15971597                                var = aveisq-avesq 
    15981598                                if(var<=0) 
    1599                                         sigave[kk-1] = small_num 
     1599//                                      sigave[kk-1] = small_num 
     1600                                        sigave[kk-1] = large_num                //if there are zero counts, make error large_num to be a warning flag 
    16001601                                else 
    16011602                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1)) 
Note: See TracChangeset for help on using the changeset viewer.