Changeset 1076


Ignore:
Timestamp:
Dec 6, 2017 1:36:11 PM (5 years ago)
Author:
srkline
Message:

woking on routines to correctly reduce data colloected with slit apertures. trying to get the data to absolute scale, with proper errro bars and resolution information.

Still not there yet.

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

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Slit.ipf

    r992 r1076  
    2424 
    2525 
     26// TODO: 
     27// -- verify the error calculation 
     28// -- add in functionality to handle FLR and MLR cases (2 panels of data) 
    2629// 
    2730// seems backwards to call this "byRows", but this is the way that Igor indexes 
     
    3639//      SetDataFolder root:Packages:NIST:VSANS:VCALC     
    3740         
    38         Variable pixSizeX,pixSizeY,delQx, delQy 
    39         Variable isVCALC=0 
     41        Variable pixSizeX,pixSizeY,delQy 
     42        Variable isVCALC=0,maskMissing,nsets=0 
     43         
     44        maskMissing = 1         // set to zero if a mask is actually present 
    4045         
    4146        if(cmpstr(folderStr,"VCALC") == 0) 
     
    4550        String folderPath = "root:Packages:NIST:VSANS:"+folderStr 
    4651        String instPath = ":entry:instrument:detector_"  
    47          
    48         if(isVCALC) 
    49                 WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data 
    50                 WAVE/Z iErr = $("asdf_iErr_"+detStr)                    // TODO: 2D errors -- may not exist, especially for simulation 
    51         else 
    52                 Wave inten = V_getDetectorDataW(folderStr,detStr) 
    53                 Wave iErr = V_getDetectorDataErrW(folderStr,detStr) 
     52 
     53        strswitch(detStr)       // string switch         
     54// only one panel, simply pick that panel and move on out of the switch 
     55                case "FL": 
     56                case "FR": 
     57                case "ML": 
     58                case "MR": 
     59                case "B": 
     60                        if(isVCALC) 
     61                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data 
     62                                WAVE/Z iErr = $("asdf_iErr_"+detStr)                    // TODO: 2D errors -- may not exist, especially for simulation 
     63                        else 
     64                                Wave inten = V_getDetectorDataW(folderStr,detStr) 
     65                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr) 
     66                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data") 
     67                                if(WaveExists(mask) == 1) 
     68                                        maskMissing = 0 
     69                                endif 
     70                        endif 
     71                         
     72                        nsets = 1 
     73                        break 
     74 
     75//              case "FLR": 
     76//              // detStr has multiple values now, so unfortuntely, I'm hard-wiring things... 
     77//                      if(isVCALC) 
     78//                              WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL") 
     79//                              WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation                 
     80//                              WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR") 
     81//                              WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation         
     82//                      else 
     83//                              Wave inten = V_getDetectorDataW(folderStr,"FL") 
     84//                              Wave iErr = V_getDetectorDataErrW(folderStr,"FL") 
     85//                              Wave inten2 = V_getDetectorDataW(folderStr,"FR") 
     86//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR") 
     87//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data") 
     88//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data") 
     89//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) 
     90//                                      maskMissing = 0 
     91//                              endif 
     92//                      endif    
     93// 
     94//                      nsets = 2 
     95//                      break 
     96//                       
     97//              case "MLR": 
     98//                      if(isVCALC) 
     99//                              WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML") 
     100//                              WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation                 
     101//                              WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR") 
     102//                              WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation         
     103//                      else 
     104//                              Wave inten = V_getDetectorDataW(folderStr,"ML") 
     105//                              Wave iErr = V_getDetectorDataErrW(folderStr,"ML") 
     106//                              Wave inten2 = V_getDetectorDataW(folderStr,"MR") 
     107//                              Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR") 
     108//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data") 
     109//                              Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data") 
     110//                              if(WaveExists(mask) == 1 && WaveExists(mask2) == 1) 
     111//                                      maskMissing = 0 
     112//                              endif 
     113//                      endif    
     114//               
     115//                      nSets = 2 
     116//                      break                    
     117 
     118                default: 
     119                        nSets = 0                                                        
     120                        Print "ERROR   ---- type is not recognized " 
     121        endswitch 
     122 
     123        if(nSets == 0) 
     124                SetDataFolder root: 
     125                return(0) 
    54126        endif 
     127 
    55128 
    56129        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
     
    58131        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr) 
    59132 
    60 // ?? TODO not needed here?      
    61 //      pixSizeX = VCALC_getPixSizeX(detStr) 
    62 //      pixSizeY = VCALC_getPixSizeY(detStr) 
    63          
    64         delQx = abs(qx[0][0] - qx[1][0]) 
    65         delQy = abs(qy[0][1] - qy[0][0]) 
    66          
    67         // delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution 
    68  
     133 
     134// delta Qx is set by the pixel X dimension of the detector, which is the limiting resolution 
     135//      delQx = abs(qx[0][0] - qx[1][0]) 
     136         
    69137        Variable nq,val 
    70138        nq = DimSize(inten,0)           //nq == the number of columns (x dimension) 
     
    86154        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy_"+detStr) 
    87155 
     156 
     157        iBin_qxqy = 0 
     158        iBin2_qxqy = 0 
     159        eBin_qxqy = 0 
     160        eBin2D_qxqy = 0 
     161        nBin_qxqy = 0    
     162 
     163 
    88164// sum the rows  
    89         MatrixOp/O iBin_qxqy = sumRows(inten)   //automatically generates the destination 
    90  
     165 
     166// MatrixOp would be fast, but I can't figure out how to apply the mask with sumRows??? 
     167//      MatrixOp/O iBin_qxqy = sumRows(inten)   //automatically generates the destination 
     168//               
     169//// how to properly calculate the error? 
     170//// This MatrixOp gives values way too large -- larger than the intensity (be sure to correct *delQy below... 
     171// 
     172//      MatrixOp/O tmp = sqrt(varCols(inten^t))         // variance: no varRows operation, so use the transpose of the matrix 
     173//      eBin_qxqy = tmp[0][p] 
     174 
     175        Variable ii,jj,ntube,npix,sum_inten, sum_n, sum_inten2,avesq,aveisq,var,mask_val 
     176        ntube = DimSize(inten,0) 
     177        npix = DimSize(inten,1) 
     178         
     179        for(ii=0;ii<ntube;ii+=1) 
     180                sum_inten = 0                   // initialize the sum 
     181                sum_n = 0 
     182                sum_inten2 = 0 
     183                 
     184                for(jj=0;jj<npix;jj+=1) 
     185                                val = inten[ii][jj] 
     186                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point 
     187                                        mask_val = 0 
     188                                else 
     189                                        mask_val = mask[ii][jj] 
     190                                endif 
     191                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf 
     192                                        sum_inten += val 
     193                                        sum_n += 1 
     194                                        sum_inten2 += val*val 
     195                                endif 
     196                endfor 
     197                iBin_qxqy[ii] = sum_inten 
     198                 
     199                avesq = sum_inten 
     200                aveisq = sum_inten2/sum_n 
     201                var = aveisq-avesq 
     202                if(var<=0) 
     203                        eBin_qxqy[ii] = 1e-6 
     204                else 
     205                        eBin_qxqy[ii] = sqrt(var/(sum_n - 1)) 
     206                endif 
     207 
     208        endfor 
     209 
     210        qBin_qxqy =  qx[p][npix/2]              //TODO:  use only the Qx component in the y-center of the detector, not Qtotal 
     211                 
    91212// if the detectors are "L", then the values are all negative... 
    92 // if the detectors are T/B, then half is negative, and there's a very nearly zero point in the middle...        
    93 // and it may make no sense to use T/B anyways... 
    94213        qBin_qxqy = abs(qx[p][0]) 
    95214 
     215// for the L panels, sort the q-values (and data) after the abs() step, otherwise the data is reversed 
     216// won't hurt to sort the R data 
     217        Sort qBin_qxqy, qBin_qxqy,iBin_qxqy,eBin_qxqy 
    96218         
    97219        //now get the scaling correct 
    98220        // q-integration (rectangular), matrixOp simply summed, so I need to multiply by dy (pixelSizeY -> as Qy?) 
    99          
     221 
     222//      delQy = abs(qy[0][1] - qy[0][0])        // this is only one pixel 
     223        delQy = abs(qy[0][npix-1] - qy[0][0])   // TODO: do I use dQ for the height of the panel? 
     224                 
    100225        iBin_qxqy *= delQy 
    101          
     226        eBin_qxqy *= delQy 
     227         
     228 
     229/// TODO -- this is not necessary, but just for getting the I(Q) display to look "pretty" 
     230// clean out the near-zero Q point in the T/B  and Back detectors 
     231        qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]                        
     232 
     233 
     234// clear out zero data values before exiting... 
     235        // find the last non-zero point, working backwards 
     236        val = numpnts(iBin_qxqy) 
     237        do 
     238                val -= 1 
     239        while((iBin_qxqy[val] == 0) && val > 0) 
     240         
     241        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,eBin_qxqy 
     242 
     243// work forwards? this doesn't work... 
     244//      val = -1 
     245//      do 
     246//              val += 1 
     247//      while(nBin_qxqy[val] == 0 && val < numpnts(iBin_qxqy)-1)         
     248//      DeletePoints 0, val, iBin_qxqy,qBin_qxqy,eBin_qxqy 
     249 
     250 
     251// TODO: 
     252// -- calculate the slit resolution here. don't really have any idea how to represent this in VSANS 
     253//    since it's not an infinite slit, and not anything that I expect could be represented as a Gaussian 
     254 
     255        // TODO: 
     256        // -- This is where I calculate the resolution in SANS (see CircSectAve) 
     257        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary 
     258        // -- from the top of the function, folderStr = work folder, detStr = "FLRTB" or other type of averaging 
     259        // 
     260 
     261        nq = numpnts(qBin_qxqy) 
     262        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+detStr) 
     263        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+detStr) 
     264        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+detStr) 
     265        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+detStr) 
     266        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+detStr) 
     267        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+detStr) 
     268 
    102269// TODO 
    103 //      iBin_qxqy *= 4          //why the factor of 4??? -- this is what I needed to do with FFT->USANS. Do I need it here? 
    104  
    105  
    106 /// TODO -- this is not correct, but just for getting the I(Q) display to look "pretty" 
    107         qBin_qxqy = (abs(qBin_qxqy[p][q]) < 1e-5) ? NaN : qBin_qxqy[p][q]                       // clean out the near-zero Q point in the T/B  and Back detectors 
     270// -- these are DUMMY VALUES!!! 
     271        sigmaq = -delQy 
     272        qbar = -delQy 
     273        fsubs = -delQy 
    108274 
    109275 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Utils.ipf

    r1075 r1076  
    799799        endif 
    800800         
     801        detStr = type 
     802 
    801803// now switch on the type to determine which waves to declare and create 
    802804// since there may be more than one panel to step through. There may be two, there may be four 
     
    902904                case "MB":                       
    903905                case "B":        
    904                         detStr = type 
    905906                        if(isVCALC) 
    906907                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr) 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VSANS_Includes.ipf

    r1073 r1076  
    1212// These files are COMMON NCNR FILES 
    1313// the first three are necessary for loading and plotting of 1D data sets 
    14 // using the PlotManager 
     14// using the PlotManager, including loading of slit-smeared VSANS data 
    1515#include "PlotUtilsMacro_v40" 
    1616#include "PlotManager_v40" 
     17#include "GaussUtils_v40" 
    1718#include "NIST_XML_v40" 
    1819// 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_DetectorCorrections.ipf

    r1075 r1076  
    158158// TODO -- replace all of this with V_getDet_panel_gap(fname,detStr) once it is added to the file 
    159159// these hard-wired values were determined from 6A and WB beam centers. LR values were exactly the same for 
    160 // both beam considitions (+/- 0.0 mm). FTB was +/- 0.8 mm, MTB +/- 2 mm 
     160// both beam conditions (+/- 0.0 mm). FTB was +/- 0.8 mm, MTB +/- 2 mm 
    161161        if(cmpstr(detStr,"FL") == 0 || cmpstr(detStr,"FR") == 0) 
    162162                gap = 3.2               //mm 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_FileCatalog.ipf

    r1073 r1076  
    1616// -- can I make the choice of columns customizable? There are "sets" of columns that are not used for  
    1717//    some experiments (magnetic, rotation, temperature scans, etc.) but are necessary for others. 
     18// -- SortColumns operation may be of help in managing the long list of files to sort 
    1819// 
    1920// (DONE): 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_IQ_Utilities.ipf

    r1073 r1076  
    193193        // remove the q=0 point from the back detector, if it's there 
    194194        // does not need to know binType 
    195         V_RemoveQ0_B(type) 
     195        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB 
     196        if(!gIgnoreDetB) 
     197                V_RemoveQ0_B(type) 
     198        endif 
    196199 
    197200// concatenate the data sets 
     
    644647 
    645648 
    646 // 
    647 // returns 1 if the val is non-negative, other value 
    648 // indicates that the resoution data is USANS data. 
    649 // 
    650 // TODO: 
    651 // -- this DUPLICATES a same-named SANS procedure, so there could be a clash at some point 
    652 // -- bigger issue - I'll need a better way to identify and load the different resolution  
    653 //              conditions with VSANS 
    654 // 
    655 // 
    656 Function isSANSResolution(val) 
    657         Variable val 
    658          
    659         if(val >= 0) 
    660                 return(1) 
    661         else 
    662                 return(0) 
    663         endif 
    664 End 
    665  
    666  
     649//// 
     650//// returns 1 if the val is non-negative, other value 
     651//// indicates that the resoution data is USANS data. 
     652//// 
     653//// TODO: 
     654//// -- this DUPLICATES a same-named SANS procedure, so there could be a clash at some point 
     655//// -- bigger issue - I'll need a better way to identify and load the different resolution  
     656////            conditions with VSANS 
     657//// 
     658//// 
     659//xFunction isSANSResolution(val) 
     660//      Variable val 
     661//       
     662//      if(val >= 0) 
     663//              return(1) 
     664//      else 
     665//              return(0) 
     666//      endif 
     667//End 
     668 
     669 
Note: See TracChangeset for help on using the changeset viewer.