Ignore:
Timestamp:
May 18, 2009 5:15:08 PM (14 years ago)
Author:
srkline
Message:

(1) fix in ProtocolAsPanel? that correctly calculates Kappa for data taken with the Cerca detectors. In the Kappa calculation, the data is manually scaled, bypassing the normal "ADD" step where the 4x detector data is correctly renormalized. Added a flag to divide by 4 if ILL detector type found.

(2) lots of bits to try to accomodate 2D resolution smearing. Changes to the QxQy? output to have everything needed: Qz, sigX, sigY, and shadowing. So presumably, all of the information is in the reduced data, where it should be.

  • added a 2D resolution calculation based on David Mildner's notes.
  • added 2D quadrature calculation for smearing, and new structure definitions

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

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Sphere_2D_v40.ipf

    r253 r509  
    3232        Duplicate/O $(str+"_qy") ywave_sf2D,zwave_sf2D                   
    3333                 
    34         Variable/G gs_sf2D=0 
    35         gs_sf2D := Sphere2D(coef_sf2D,zwave_sf2D,xwave_sf2D,ywave_sf2D) //AAO 2D calculation 
     34        Variable/G g_sf2D=0 
     35        g_sf2D := Sphere2D(coef_sf2D,zwave_sf2D,xwave_sf2D,ywave_sf2D)  //AAO 2D calculation 
    3636         
    3737        Display ywave_sf2D vs xwave_sf2D 
     
    5151 
    5252        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation 
    53         Variable/G gs_sf2Dmat=0 
    54         gs_sf2Dmat := UpdateQxQy2Mat(xwave_sf2D,ywave_sf2D,zwave_sf2D,sf2D_lin,sf2D_mat) 
     53        Variable/G g_sf2Dmat=0 
     54        g_sf2Dmat := UpdateQxQy2Mat(xwave_sf2D,ywave_sf2D,zwave_sf2D,sf2D_lin,sf2D_mat) 
    5555         
    5656         
    5757        SetDataFolder root: 
    5858        AddModelToStrings("Sphere2D","coef_sf2D","sf2D") 
     59End 
     60 
     61// - sets up a dependency to a wrapper, not the actual SmearedModelFunction 
     62Proc PlotSmearedSphere2D(str)                                                            
     63        String str 
     64        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4) 
     65         
     66        // if any of the resolution waves are missing => abort 
     67//      if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils) 
     68//              Abort 
     69//      endif 
     70         
     71        SetDataFolder $("root:"+str) 
     72         
     73        // Setup parameter table for model function 
     74        Make/O/D smear_coef_sf2D = {1.,60,1e-6,6.3e-6,0.01}                                      
     75        make/o/t smear_parameters_sf2D = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)","bkgd (cm-1)"} 
     76        Edit smear_parameters_sf2D,smear_coef_sf2D                                       
     77         
     78        Duplicate/O $(str+"_qx") smeared_sf2D   //1d place for the smeared model 
     79        SetScale d,0,0,"1/cm",smeared_sf2D                                       
     80                 
     81        Variable/G gs_sf2D=0 
     82        gs_sf2D := fSmearedSphere2D(smear_coef_sf2D,smeared_sf2D)       //this wrapper fills the STRUCT 
     83 
     84        Display $(str+"_qy") vs $(str+"_qx") 
     85        modifygraph log=0 
     86        ModifyGraph mode=3,marker=16,zColor($(str+"_qy"))={smeared_sf2D,*,*,YellowHot,0} 
     87        ModifyGraph standoff=0 
     88        ModifyGraph width={Aspect,1} 
     89        ModifyGraph lowTrip=0.001 
     90        Label bottom "qx (A\\S-1\\M)" 
     91        Label left "qy (A\\S-1\\M)" 
     92        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) 
     93         
     94        SetDataFolder root: 
     95        AddModelToStrings("SmearedSphere2D","smear_coef_sf2D","sf2D") 
    5996End 
    6097 
     
    106143        return(0) 
    107144End 
     145 
     146 
     147//threaded version of the function 
     148Function Sphere2D_noThread(cw,zw,xw,yw) 
     149        WAVE cw,zw, xw,yw 
     150         
     151#if exists("Sphere_2DX")                        //to hide the function if XOP not installed 
     152         
     153        zw= Sphere_2DX(cw,xw,yw) 
     154 
     155#endif 
     156 
     157        return 0 
     158End 
     159 
     160 
     161Function SmearedSphere2D(s) 
     162        Struct ResSmear_2D_AAOStruct &s 
     163         
     164        Smear_2DModel_5(Sphere2D_noThread,s) 
     165        return(0) 
     166end 
     167 
     168 
     169Function fSmearedSphere2D(coefW,resultW) 
     170        Wave coefW,resultW 
     171         
     172        String str = getWavesDataFolder(resultW,0) 
     173        String DF="root:"+str+":" 
     174         
     175        WAVE qx = $(DF+str+"_qx") 
     176        WAVE qy = $(DF+str+"_qy") 
     177        WAVE qz = $(DF+str+"_qz") 
     178        WAVE sigQx = $(DF+str+"_sigQx") 
     179        WAVE sigQy = $(DF+str+"_sigQy") 
     180        WAVE shad = $(DF+str+"_fs") 
     181         
     182        STRUCT ResSmear_2D_AAOStruct s 
     183        WAVE s.coefW = coefW     
     184        WAVE s.zw = resultW      
     185        WAVE s.qx = qx 
     186        WAVE s.qy = qy 
     187        WAVE s.qz = qz 
     188        WAVE s.sigQx = sigQx 
     189        WAVE s.sigQy = sigQy 
     190        WAVE s.fs = shad 
     191         
     192        Variable err 
     193        err = SmearedSphere2D(s) 
     194         
     195        return (0) 
     196End 
  • sans/Dev/trunk/NCNR_User_Procedures/Common/GaussUtils_v40.ipf

    r480 r509  
    6161EndStructure 
    6262 
     63 
     64// tentative pass at 2D resolution smearing 
     65// 
     66Structure ResSmear_2D_AAOStruct 
     67        Wave coefW 
     68        Wave zw                 //answer 
     69        Wave qy                 // q-value 
     70        Wave qx 
     71        Wave qz 
     72        Wave sigQx              //resolution 
     73        Wave sigQy 
     74        Wave fs 
     75        String info 
     76EndStructure 
    6377 
    6478 
     
    10621076end 
    10631077 
     1078// prototype function for 2D smearing routine 
     1079Function SANS_2D_ModelAAO_proto(w,zw,xw,yw) 
     1080        Wave w,zw,xw,yw 
     1081         
     1082        Print "in SANSModelAAO_proto function" 
     1083        return(1) 
     1084end 
     1085 
     1086// prototype function for fit wrapper using 2D smearing 
     1087// not used (yet) 
     1088Function SANS_2D_ModelSTRUCT_proto(s) 
     1089        Struct ResSmear_2D_AAOStruct &s  
     1090 
     1091        Print "in SANSModelSTRUCT_proto function" 
     1092        return(1) 
     1093end 
     1094 
    10641095//Numerical Recipes routine to calculate the nn(th) stage 
    10651096//refinement of a trapezoid integration 
  • sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtils2D_v40.ipf

    r490 r509  
    1515 
    1616// 
     17// changed May 2009 to read in the resolution information (now 7 columns) 
     18// -- subject to change -- 
    1719// 
    1820Proc LoadQxQy() 
     
    2426 
    2527        String w0,w1,w2,n0,n1,n2 
     28        String w3,w4,w5,w6,n3,n4,n5,n6 
    2629                 
    2730        // put the names of the three loaded waves into local names 
     
    2932        n1 = StringFromList(1, S_waveNames ,";" ) 
    3033        n2 = StringFromList(2, S_waveNames ,";" ) 
     34        n3 = StringFromList(3, S_waveNames ,";" ) 
     35        n4 = StringFromList(4, S_waveNames ,";" ) 
     36        n5 = StringFromList(5, S_waveNames ,";" ) 
     37        n6 = StringFromList(6, S_waveNames ,";" ) 
    3138         
    3239        //remove the semicolon AND period from files from the VAX 
     
    3441        w1 = CleanupName((S_fileName + "_qy"),0) 
    3542        w2 = CleanupName((S_fileName + "_i"),0) 
    36          
     43        w3 = CleanupName((S_fileName + "_qz"),0) 
     44        w4 = CleanupName((S_fileName + "_sigQx"),0) 
     45        w5 = CleanupName((S_fileName + "_sigQy"),0) 
     46        w6 = CleanupName((S_fileName + "_fs"),0) 
     47 
    3748        String baseStr=w1[0,strlen(w1)-4] 
    3849        if(DataFolderExists("root:"+baseStr)) 
     
    4051                        if(V_flag==2)   //user selected No, don't load the data 
    4152                                SetDataFolder root: 
    42                                 KillWaves $n0,$n1,$n2           // kill the default waveX that were loaded 
     53                                KillWaves $n0,$n1,$n2,$n3,$n4,$n5,$n6           // kill the default waveX that were loaded 
    4354                                return          //quits the macro 
    4455                        endif 
     
    6980        Duplicate/O $("root:"+n1), $w1 
    7081        Duplicate/O $("root:"+n2), $w2 
     82        Duplicate/O $("root:"+n3), $w3 
     83        Duplicate/O $("root:"+n4), $w4 
     84        Duplicate/O $("root:"+n5), $w5 
     85        Duplicate/O $("root:"+n6), $w6 
    7186         
    7287        Variable/G gIsLogScale = 0 
     
    88103        //clean up               
    89104        SetDataFolder root: 
    90         KillWaves/Z $n0,$n1,$n2 
     105        KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6 
    91106EndMacro 
    92107 
     
    413428                                //make log scale 
    414429                                w=log(linW) 
     430                                w = w[p][q] == -inf ? NaN : w[p][q]             //remove the -inf for display 
    415431                        endif 
    416432        endfor   
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NCNR_Utils.ipf

    r412 r509  
    152152        Return results 
    153153End 
     154 
     155// 
     156// 
     157//      this should end up in NCNR_Utils once it's fully functional 
     158// 
     159// lots of unnecessary junk... 
     160// 
     161//********************** 
     162// 2D resolution function calculation - in terms of X and Y 
     163// 
     164// based on notes from David Mildner, 2008 
     165// 
     166// - 21 MAR 07 uses projected BS diameter on the detector 
     167// - APR 07 still need to add resolution with lenses. currently there is no flag in the  
     168//          raw data header to indicate the presence of lenses. 
     169// 
     170// - Aug 07 - added input to switch calculation based on lenses (==1 if in) 
     171// 
     172// passed values are read from RealsRead 
     173// except DDet and apOff, which are set from globals before passing 
     174// 
     175// phi is the azimuthal angle, CCW from +x axis 
     176// r_dist is the real-space distance from ctr of detector to QxQy pixel location 
     177// 
     178Function/S get2DResolution(inQ,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,SigmaQX,SigmaQY,fSubS) 
     179        Variable inQ, phi,lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r,usingLenses,r_dist 
     180        Variable &SigmaQX,&SigmaQY,&fSubS               //these are the output quantities at the input Q value 
     181         
     182        //lots of calculation variables 
     183        Variable a2, lp, v_lambda, v_b, v_d, vz, yg, v_g 
     184        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r 
     185 
     186        //Constants 
     187        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron 
     188        Variable g = 981.0                              //gravity acceleration [cm/s^2] 
     189        Variable h_m = 3995                             // h/m [=] A*m/s 
     190 
     191        String results 
     192        results ="Failure" 
     193 
     194        S1 *= 0.5*0.1                   //convert to radius and [cm] 
     195        S2 *= 0.5*0.1 
     196 
     197        L1 *= 100.0                     // [cm] 
     198        L1 -= apOff                             //correct the distance 
     199 
     200        L2 *= 100.0 
     201        L2 += apOff 
     202        del_r *= 0.1                            //width of annulus, convert mm to [cm] 
     203         
     204        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm] 
     205        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture 
     206        Variable LB 
     207        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical) 
     208        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax 
     209         
     210        //Start resolution calculation 
     211        a2 = S1*L2/L1 + S2*(L1+L2)/L1 
     212        lp = 1.0/( 1.0/L1 + 1.0/L2) 
     213 
     214        v_lambda = lambdaWidth^2/6.0 
     215         
     216//      if(usingLenses==1)                      //SRK 2007 
     217        if(usingLenses != 0)                    //SRK 2008 allows for the possibility of different numbers of lenses in header 
     218                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term 
     219        else 
     220                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form 
     221        endif 
     222         
     223        v_d = (DDet/2.3548)^2 + del_r^2/12.0 
     224        vz = vz_1 / lambda 
     225        yg = 0.5*g*L2*(L1+L2)/vz^2 
     226        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007 
     227 
     228        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) )) 
     229        delta = 0.5*(BS - r0)^2/v_d 
     230 
     231        if (r0 < BS)  
     232                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta)) 
     233        else 
     234                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta)) 
     235        endif 
     236 
     237        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) ) 
     238        if (fSubS <= 0.0)  
     239                fSubS = 1.e-10 
     240        endif 
     241//      fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi)) 
     242//      fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d 
     243// 
     244//      rmd = fr*r0 
     245//      v_r1 = v_b + fv*v_d +v_g 
     246// 
     247//      rm = rmd + 0.5*v_r1/rmd 
     248//      v_r = v_r1 - 0.5*(v_r1/rmd)^2 
     249//      if (v_r < 0.0)  
     250//              v_r = 0.0 
     251//      endif 
     252 
     253        Variable kap,a_val 
     254         
     255        kap = 2*pi/lambda 
     256        a_val = (L1+L2)*g/2/(h_m)^2 
     257         
     258        SigmaQX = 3*(S1/L1)^2 + 3*(S2/LP)^2 + 2*(DDet/L2)^2 + 2*(r_dist/L2)^2*(lambdaWidth)^2*(cos(phi))^2 
     259 
     260        SigmaQY = 3*(S1/L1)^2 + 3*(S2/LP)^2 + 2*(DDet/L2)^2 + 2*(r_dist/L2)^2*(lambdaWidth)^2*(sin(phi))^2 + 8*(a_val/L2)^2*lambda^4*lambdaWidth^2 
     261 
     262        SigmaQX = sqrt(kap*kap/12*SigmaQX) 
     263        SigmaQy = sqrt(kap*kap/12*SigmaQY) 
     264 
     265        results = "success" 
     266        Return results 
     267End 
     268 
     269 
    154270 
    155271 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/ProtocolAsPanel.ipf

    r418 r509  
    20732073//              Print "box is now ",x1,x2,y1,y2 
    20742074                detCnt = SumCountsInBox(x1,x2,y1,y2,"RAW") 
     2075                if(cmpstr(tw[9],"ILL   ")==0) 
     2076                        detCnt /= 4             // for cerca detector, header is right, sum(data) is 4x too large 
     2077                                                                // this is usually corrected in the Add step 
     2078                endif 
    20752079                //               
    20762080                kappa = detCnt/countTime/attenTrans*1.0e8/(monCnt/countTime)*(pixel/sdd)^2 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/WriteQIS.ipf

    r474 r509  
    630630 
    631631 
    632 //ASCII export of data as 3-columns qx-qy-Intensity 
     632// NEW additions - May 2009 
     633//ASCII export of data as 7-columns qx-qy-Intensity-qz-sigmaQx-sigmaQy-fShad 
    633634//limited header information? 
    634635// 
     
    718719        labelWave[13] = "" 
    719720        labelWave[14] = "*** Data written from "+type+" folder and may not be a fully corrected data file ***" 
    720         labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy)" 
     721        labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy) - Qz - SigmaQx - SigmaQy - fSubS(beam stop shadow)" 
    721722        labelWave[16] = "" 
    722723        labelWave[17] = "ASCII data created " +date()+" "+time() 
     
    736737//      Endif 
    737738         
    738         Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val 
     739        Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist 
    739740         
    740741//      Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val 
    741742//      MyMat2XYZ(data,qx_val,qy_val,z_val)             //x and y are [p][q] indexes, not q-vals yet 
    742743         
     744        Variable xctr,yctr,sdd,lambda,pixSize 
     745        xctr = rw[16] 
     746        yctr = rw[17] 
     747        sdd = rw[18] 
     748        lambda = rw[26] 
     749        pixSize = rw[13]/10             //convert mm to cm (x and y are the same size pixels) 
     750         
    743751        qx_val = CalcQx(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)          //+1 converts to detector coordinate system 
    744752        qy_val = CalcQy(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    745753         
    746754        Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val 
    747          
     755 
     756///************ 
     757// do everything to write out the resolution too 
    748758        // un-comment these if you want to write out qz_val and qval too, then use the proper save command 
    749 //      qval = CalcQval(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    750 //      qz_val = CalcQz(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    751 //      Redimension/N=(pixelsX*pixelsY) qz_val,qval 
     759        qval = CalcQval(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     760        qz_val = CalcQz(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     761        phi = FindPhi( pixSize*((p+1)-xctr) , pixSize*((q+1)-yctr))             //(dx,dy) 
     762        r_dist = sqrt(  (pixSize*((p+1)-xctr))^2 +  (pixSize*((q+1)-yctr))^2 )          //radial distance from ctr to pt 
     763        Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist 
     764        //everything in 1D now 
     765        Duplicate/O qval SigmaQX,SigmaQY,fsubS 
     766 
     767        Variable L2 = rw[18] 
     768        Variable BS = rw[21] 
     769        Variable S1 = rw[23] 
     770        Variable S2 = rw[24] 
     771        Variable L1 = rw[25] 
     772        Variable lambdaWidth = rw[27]    
     773        Variable usingLenses = rw[28]           //new 2007 
     774 
     775        //Two parameters DDET and APOFF are instrument dependent.  Determine 
     776        //these from the instrument name in the header. 
     777        //From conversation with JB on 01.06.99 these are the current good values 
     778        Variable DDet 
     779        NVAR apOff = root:myGlobals:apOff               //in cm 
     780        DDet = rw[10]/10                        // header value (X) is in mm, want cm here 
     781 
     782        Variable ret1,ret2,ret3,nq 
     783        nq = pixelsX*pixelsY 
     784        ii=0 
     785         
     786        do 
     787                get2DResolution(qval[ii],phi[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,pixSize,usingLenses,r_dist[ii],ret1,ret2,ret3) 
     788                SigmaQX[ii] = ret1       
     789                SigmaQY[ii] = ret2       
     790                fsubs[ii] = ret3         
     791                ii+=1 
     792        while(ii<nq)     
     793 
     794//********************* 
    752795 
    753796        //not demo-compatible, but approx 8x faster!!    
    754797#if(cmpstr(stringbykey("IGORKIND",IgorInfo(0),":",";"),"pro") == 0)      
    755         Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val as fullpath       // /M=termStr specifies terminator       
    756 //      Save/G/M="\r\n" labelWave,qx_val,qy_val,qz_val,qval,z_val as fullpath   // for debugging, write out everything 
     798//      Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val as fullpath       // /M=termStr specifies terminator       
     799        Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val,qz_val,SigmaQx,SigmaQy,fSubS as fullpath  // for debugging, write out everything 
    757800#else 
    758801        Open refNum as fullpath 
     
    763806#endif 
    764807         
    765         Killwaves/Z spWave,labelWave,qx_val,qy_val,z_val,qval,qz_val 
     808        Killwaves/Z spWave,labelWave,qx_val,qy_val,z_val,qval,qz_val,sigmaQx,SigmaQy,fSubS,phi,r_dist 
    766809         
    767810        Print "QxQy_Export File written: ", GetFileNameFromPathNoSemi(fullPath) 
     
    812855        return(0) 
    813856End 
    814  
Note: See TracChangeset for help on using the changeset viewer.