Ignore:
Timestamp:
Feb 1, 2019 2:25:12 PM (4 years ago)
Author:
srkline
Message:

added procedures to output QxQy_ASCII data. Each panel is output into its own file. Output format is the same as for 2D SANS data, including the 2D resolution function. However, reading the data back in with the SANS analysis macros currently does not redimension the data back to the matrix correctly as it assumes a square detector.

I will need to add the X and Y dimensions of each panel into the header, and then make use of these values when they are read in. Also, writing the QxQy? data is quick for the M and F panels ( 0.8 s) but is extremely slow for the back, High-res panel (120 s) since there are 1.1.E6 points there vs. 6144 pts. I'll need to find a way to speed this operation up.

File:
1 edited

Legend:

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

    r1108 r1119  
    340340        V_SaveIQ_ButtonProc(ba) 
    341341end 
     342 
     343 
     344///////// QxQy Export  ////////// 
     345// 
     346// (see the similar-named SANS routine for additonal steps - like resolution, etc.) 
     347//ASCII export of data as 8-columns qx-qy-Intensity-err-qz-sigmaQ_parall-sigmaQ_perp-fShad 
     348// + limited header information 
     349// 
     350//      Jan 2019 -- first version, simply exports the basic matrix of data with no resolution information 
     351// 
     352// 
     353Function V_QxQy_Export(type,fullpath,newFileName,dialog) 
     354        String type,fullpath,newFileName 
     355        Variable dialog         //=1 will present dialog for name 
     356         
     357        String typeStr="" 
     358        Variable refnum 
     359        String detStr="",detSavePath 
     360 
     361        SVAR gProtoStr = root:Packages:NIST:VSANS:Globals:Protocols:gProtoStr 
     362         
     363        // declare, or make a fake protocol if needed (if the export type is RAW) 
     364        String rawTag="" 
     365        if(cmpstr(type,"RAW")==0) 
     366                Make/O/T/N=(kNumProtocolSteps) proto 
     367                RawTag = "RAW Data File: "       
     368        else 
     369                Wave/T proto=$("root:Packages:NIST:VSANS:Globals:Protocols:"+gProtoStr)  
     370        endif 
     371         
     372        SVAR samFiles = $("root:Packages:NIST:VSANS:"+type+":gFileList") 
     373         
     374        //check each wave - MUST exist, or will cause a crash 
     375//      If(!(WaveExists(data))) 
     376//              Abort "data DNExist QxQy_Export()" 
     377//      Endif 
     378 
     379        if(dialog) 
     380                PathInfo/S catPathName 
     381                fullPath = DoSaveFileDialog("Save data as") 
     382                If(cmpstr(fullPath,"")==0) 
     383                        //user cancel, don't write out a file 
     384                        Close/A 
     385                        Abort "no data file was written" 
     386                Endif 
     387                //Print "dialog fullpath = ",fullpath 
     388        Endif 
     389 
     390        // data values to populate the file header 
     391        String fileName,fileDate,fileLabel 
     392        Variable monCt,lambda,offset,dist,trans,thick 
     393        Variable bCentX,bCentY,a2,a1a2_dist,deltaLam,bstop 
     394        String a1Str 
     395        Variable pixX,pixY 
     396        Variable numTextLines=19,ii,jj,kk 
     397 
     398        Make/O/T/N=(numTextLines) labelWave 
     399 
     400        // 
     401         
     402        //loop over all of the detector panels 
     403        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB 
     404 
     405        String detList 
     406        if(gIgnoreDetB) 
     407                detList = ksDetectorListNoB 
     408        else 
     409                detList = ksDetectorListAll 
     410        endif 
     411         
     412        for(kk=0;kk<ItemsInList(detList);kk+=1) 
     413 
     414                detStr = StringFromList(kk, detList, ";") 
     415                detSavePath = fullPath + "_" + detStr 
     416                 
     417                pixX = V_getDet_pixel_num_x(type,detStr) 
     418                pixY = V_getDet_pixel_num_y(type,detStr) 
     419                 
     420                fileName = newFileName 
     421                fileDate = V_getDataStartTime(type)             // already a string 
     422                fileLabel = V_getSampleDescription(type) 
     423                 
     424                monCt = V_getBeamMonNormData(type) 
     425                lambda = V_getWavelength(type) 
     426         
     427        // TODO - switch based on panel type 
     428        //      V_getDet_LateralOffset(fname,detStr) 
     429        //      V_getDet_VerticalOffset(fname,detStr) 
     430                offset = V_getDet_LateralOffset(type,detStr) 
     431         
     432                dist = V_getDet_ActualDistance(type,detStr) 
     433                trans = V_getSampleTransmission(type) 
     434                thick = V_getSampleThickness(type) 
     435                 
     436                bCentX = V_getDet_beam_center_x(type,detStr) 
     437                bCentY = V_getDet_beam_center_y(type,detStr) 
     438                a1Str = V_getSourceAp_size(type)                //already a string 
     439                a2 = V_getSampleAp2_size(type) 
     440                a1a2_dist = V_getSourceAp_distance(type) 
     441                deltaLam = V_getWavelength_spread(type) 
     442                // TODO -- decipher which beamstop, if any is actually in place 
     443        // or -- V_getBeamStopC3_size(type) 
     444                bstop = V_getBeamStopC2_size(type) 
     445                 
     446        ///////// 
     447                labelWave[0] = "FILE: "+fileName+"   CREATED: "+fileDate 
     448                labelWave[1] = "LABEL: "+fileLabel 
     449                labelWave[2] = "MON CNT   LAMBDA (A)  DET_OFF(cm)   DET_DIST(cm)   TRANS   THICK(cm)" 
     450                labelWave[3] = num2str(monCt)+"  "+num2str(lambda)+"       "+num2str(offset)+"     "+num2str(dist) 
     451                labelWave[3] += "     "+num2str(trans)+"     "+num2str(thick) 
     452                labelWave[4] = "BCENT(X,Y)(cm)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)" 
     453                labelWave[5] = num2str(bCentX)+"  "+num2str(bCentY)+"  "+a1Str+"    "+num2str(a2)+"    " 
     454                labelWave[5] += num2str(a1a2_dist)+"    "+num2str(deltaLam)+"    "+num2str(bstop) 
     455                labelWave[6] =  "SAM: "+rawTag+samFiles 
     456                labelWave[7] =  "BGD: "+proto[0] 
     457                labelWave[8] =  "EMP: "+proto[1] 
     458                labelWave[9] =  "DIV: "+proto[2] 
     459                labelWave[10] =  "MASK: "+proto[3] 
     460                labelWave[11] =  "ABS Parameters (3-6): "+proto[4] 
     461                labelWave[12] = "Average Choices: "+proto[5] 
     462                labelWave[13] = "Collimation type: "+proto[9] 
     463                labelWave[14] = "" 
     464                labelWave[15] = "*** Data written from "+type+" folder and may not be a fully corrected data file ***" 
     465//              labelWave[16] = "Data columns are Qx - Qy - Qz - I(Qx,Qy) - Err I(Qx,Qy)" 
     466        //      labelWave[16] = "Data columns are Qx - Qy - I(Qx,Qy) - Qz - SigmaQ_parall - SigmaQ_perp - fSubS(beam stop shadow)" 
     467                labelWave[16] = "Data columns are Qx - Qy - I(Qx,Qy) - err(I) - Qz - SigmaQ_parall - SigmaQ_perp - fSubS(beam stop shadow)" 
     468                labelWave[17] = "The error wave may not be properly propagated (1/2019)" 
     469                labelWave[18] = "ASCII data created " +date()+" "+time() 
     470                //strings can be too long to print-- must trim to 255 chars 
     471                for(jj=0;jj<numTextLines;jj+=1) 
     472                        labelWave[jj] = (labelWave[jj])[0,240] 
     473                endfor 
     474         
     475         
     476        // get the data waves for output 
     477        // QxQyQz have already been calculated for VSANS data 
     478                 
     479                WAVE data = V_getDetectorDataW(type,detStr) 
     480                WAVE data_err = V_getDetectorDataErrW(type,detStr) 
     481                 
     482                // TOOD - replace hard wired paths with Read functions 
     483                // hard-wired 
     484                Wave qx_val = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qx_"+detStr) 
     485                Wave qy_val = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qy_"+detStr) 
     486                Wave qz_val = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qz_"+detStr) 
     487                Wave qTot = $("root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_"+detStr+":qTot_"+detStr) 
     488                 
     489        ///// calculation of the resolution function (2D) 
     490 
     491        // 
     492                Variable acc,ssd,lambda0,yg_d,qstar,g,L1,L2,vz_1,sdd 
     493                // L1 = source to sample distance [cm]  
     494                L1 = V_getSourceAp_distance(type) 
     495         
     496        // L2 = sample to detector distance [cm] 
     497                L2 = V_getDet_ActualDistance(type,detStr)               //cm 
     498 
     499        //               
     500                G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2 
     501                vz_1 =  3.956E5 //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC 
     502                acc = vz_1 
     503                SDD = L2                //1317 
     504                SSD = L1                //1627          //cm 
     505                lambda0 = lambda                //              15 
     506                YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2 
     507                Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG_d 
     508        ////            Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd 
     509                qstar = -2*pi/lambda0*2*yg_d/sdd 
     510        //       
     511        // 
     512        //// the gravity center is not the resolution center 
     513        //// gravity center = beam center 
     514        //// resolution center = offset y = dy + (2)*yg_d 
     515        /////************ 
     516        //// do everything to write out the resolution too 
     517        //      // un-comment these if you want to write out qz_val and qval too, then use the proper save command 
     518        //      qval = CalcQval(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     519                Duplicate/O qTot,phi,r_dist 
     520                Variable pixSizeX,pixSizeY 
     521                pixSizeX = V_getDet_x_pixel_size(type,detStr) 
     522                pixSizeY = V_getDet_y_pixel_size(type,detStr) 
     523 
     524                phi = V_FindPhi( pixSizeX*((p+1)-bCentX) , pixSizeY*((q+1)-bCentY)+(2)*yg_d)            //(dx,dy+yg_d) 
     525                r_dist = sqrt(  (pixSizeX*((p+1)-bCentX))^2 +  (pixSizeY*((q+1)-bCentY)+(2)*yg_d)^2 )           //radial distance from ctr to pt 
     526         
     527                //make everything in 1D now 
     528                Duplicate/O qTot SigmaQX,SigmaQY,fsubS,qval      
     529                Redimension/N=(pixX*pixY) SigmaQX,SigmaQY,fsubS,qval,phi,r_dist 
     530 
     531                Variable ret1,ret2,ret3,nq 
     532                String collimationStr 
     533                 
     534                 
     535                collimationStr = proto[9] 
     536                 
     537                nq = pixX*pixY 
     538                ii=0 
     539 
     540s_tic()  
     541// TODO 
     542// this loop is the slow step. it takes Å 0.7 s for F or M panels, and Å 120 s for the Back panel (6144 pts vs. 1.12e6 pts) 
     543// find some way to speed this up! 
     544// MultiThreading will be difficult as it requires all the dependent functions (HDF5 reads, etc.) to be threadsafe as well 
     545// and there are a lot of them... 
     546                //type = work folder 
     547                 
     548//              (this doesn't work...and isn't any faster) 
     549//              Duplicate/O qval dum 
     550//              dum = V_get2DResolution(qval,phi,r_dist,type,detStr,collimationStr,SigmaQX,SigmaQY,fsubS) 
     551 
     552                do 
     553                        V_get2DResolution(qval[ii],phi[ii],r_dist[ii],type,detStr,collimationStr,ret1,ret2,ret3) 
     554                        SigmaQX[ii] = ret1       
     555                        SigmaQY[ii] = ret2       
     556                        fsubs[ii] = ret3         
     557                        ii+=1 
     558                while(ii<nq)     
     559s_toc() 
     560         
     561        ////*********************        
     562                Duplicate/O qx_val,qx_val_s 
     563                Duplicate/O qy_val,qy_val_s 
     564                Duplicate/O qz_val,qz_val_s 
     565                Duplicate/O data,z_val_s 
     566                Duplicate/O SigmaQx,sigmaQx_s 
     567                Duplicate/O SigmaQy,sigmaQy_s 
     568                Duplicate/O fSubS,fSubS_s 
     569                Duplicate/O data_err,sw_s 
     570                 
     571                //so that double precision data is not written out 
     572                Redimension/S qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s 
     573                Redimension/S SigmaQx_s,SigmaQy_s,fSubS_s 
     574         
     575                Redimension/N=(pixX*pixY) qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s 
     576                 
     577                //not demo-compatible, but approx 8x faster!!    
     578#if(strsearch(stringbykey("IGORKIND",IgorInfo(0),":",";"), "demo", 0 ) == -1) 
     579                 
     580//              Save/O/G/M="\r\n" labelWave,qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s as detSavePath      // without resolution 
     581                Save/O/G/M="\r\n" labelWave,qx_val_s,qy_val_s,z_val_s,sw_s,qz_val_s,SigmaQx_s,SigmaQy_s,fSubS_s as detSavePath  // write out the resolution information 
     582#else 
     583                Open refNum as detSavePath 
     584                wfprintf refNum,"%s\r\n",labelWave 
     585                fprintf refnum,"\r\n" 
     586//              wfprintf refNum,"%8g\t%8g\t%8g\t%8g\t%8g\r\n",qx_val_s,qy_val_s,qz_val_s,z_val_s,sw_s 
     587                wfprintf refNum,"%8g\t%8g\t%8g\t%8g\t%8g\t%8g\t%8g\t%8g\r\n",qx_val,qy_val,z_val,sw,qz_val,SigmaQx,SigmaQy,fSubS 
     588                Close refNum 
     589#endif 
     590                 
     591                KillWaves/Z qx_val_s,qy_val_s,z_val_s,qz_val_s,SigmaQx_s,SigmaQy_s,fSubS_s,sw,sw_s 
     592                 
     593                Killwaves/Z qx_val,qy_val,z_val,qval,qz_val,sigmaQx,SigmaQy,fSubS 
     594                 
     595                Print "QxQy_Export File written: ", V_GetFileNameFromPathNoSemi(detSavePath) 
     596         
     597        endfor 
     598         
     599        KillWaves/Z labelWave,dum 
     600        return(0) 
     601End 
     602 
     603// this assumes that: 
     604// --QxQy data was written out in the format specified by the Igor macros, that is the x varies most rapidly 
     605// 
     606// TODO -- this needs to be made generic for reading in different panels with different XY dimensions 
     607// -- add the XY dimensions to the QxQyASCII file header somewhere so that it can be read in and used here 
     608// 
     609// the SANS analysis 2D loader assumes that the matrix is square, mangling the VSANS data. 
     610// the column data (for fitting) is still fine, but the matrix representation is incorrect. 
     611// 
     612Function V_ConvertQxQy2Mat(Qx,Qy,inten,matStr) 
     613        Wave Qx,Qy,inten 
     614        String matStr 
     615         
     616        String folderStr=GetWavesDataFolder(Qx,1) 
     617         
     618        Variable numX,numY 
     619        numX=48 
     620        numY=128 
     621        Make/O/D/N=(numX,numY) $(folderStr + matStr) 
     622        Wave mat=$(folderStr + matStr) 
     623         
     624        WaveStats/Q Qx 
     625        SetScale/I x, V_min, V_max, "", mat 
     626        WaveStats/Q Qy 
     627        SetScale/I y, V_min, V_max, "", mat 
     628         
     629        Variable xrows=numX 
     630         
     631        mat = inten[q*xrows+p] 
     632         
     633        return(0) 
     634End 
     635 
     636 
Note: See TracChangeset for help on using the changeset viewer.