source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_WorkFolderUtils.ipf @ 1119

Last change on this file since 1119 was 1119, checked in by srkline, 4 years ago

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 size: 36.5 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma version=1.0
3#pragma IgorVersion=6.1
4
5//*******************
6// Vers 1.0 JAN2016
7//
8//*******************
9//  VSANS Utility procedures for handling of workfiles (each is housed in a separate datafolder)
10//
11// - adding RAW data to a workfile
12// -- **this conversion applies the detector corrections**
13// -- Raw_to_work(newType) IS THE MAJOR ROUTINE TO APPLY DETECTOR CORRECTIONS
14//
15//
16// - copying workfiles to another folder
17//
18// - absolute scaling
19//
20// - (no) the WorkFile Math panel for simple image math (not done - maybe in the future?)
21// -
22// - (no) adding work.drk data without normalizing to monitor counts (the case not currently handled)
23//***************************
24
25//
26// Functions used for manipulation of the local Igor "WORK" folder
27// structure as raw data is displayed and processed.
28//
29//
30
31
32//
33//Entry procedure from main panel
34//
35Proc V_CopyWorkFolder(oldType,newType)
36        String oldType,newType
37        Prompt oldType,"Source WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
38        Prompt newType,"Destination WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
39
40        // data folder "old" will be copied to "new" (either kills/copies or will overwrite)
41        V_CopyHDFToWorkFolder(oldtype,newtype)
42End
43
44//
45// copy what is needed for data processing (not the DAS_logs)
46// from the RawVSANS storage folder to the local WORK folder as needed
47//
48// TODO -- at what stage do I make copies of data in linear/log forms for data display?
49//      -- do I need to make copies, if I'm displaying with the lookup wave (no copy needed if this works)
50//                      -- when do I make the 2D error waves?
51//
52// TODO - decide what exactly I need to copy over. May be best to copy all, and delete
53//       what I know that I don't need
54//
55// TODO !!! DuplicateDataFolder will FAIL - in the base case of RAW data files, the
56//  data is actually in use - so it will fail every time. need an alternate solution. in SANS,
57// there are a limited number of waves to carry over, so Dupliate/O is used for rw, tw, data, etc.
58//
59// TODO : I also need a list of what is generated during processing that may be hanging around - that I need to
60//     be sure to get rid of - like the calibration waves, solidAngle, etc.
61//
62// hdfDF is the name only of the data in storage. May be full file name with extension (clean as needed)
63// type is the destination WORK folder for the copy
64//
65Function V_CopyHDFToWorkFolder(fromStr,toStr)
66        String fromStr,toStr
67       
68        String fromDF, toDF
69       
70        // make the DF paths - source and destination
71        fromDF = "root:Packages:NIST:VSANS:"+fromStr
72        toDF = "root:Packages:NIST:VSANS:"+toStr
73       
74//      // make a copy of the file name for my own use, since it's not in the file
75//      String/G $(toDF+":file_name") = root:
76       
77        // copy the folders
78        KillDataFolder/Z $toDF                  //DuplicateDataFolder will not overwrite, so Kill
79       
80        if(V_flag == 0)         // kill DF was OK
81                DuplicateDataFolder $("root:Packages:NIST:VSANS:"+fromStr),$("root:Packages:NIST:VSANS:"+toStr)
82               
83                // I can delete these if they came along with RAW
84                //   DAS_logs
85                //   top-level copies of data (duplicate links, these should not be present in a proper NICE file)
86                KillDataFolder/Z $(toDF+":entry:DAS_logs")
87                KillDataFolder/Z $(toDF+":entry:data")
88                KillDataFolder/Z $(toDF+":entry:data_B")
89                KillDataFolder/Z $(toDF+":entry:data_ML")
90                KillDataFolder/Z $(toDF+":entry:data_MR")
91                KillDataFolder/Z $(toDF+":entry:data_MT")
92                KillDataFolder/Z $(toDF+":entry:data_MB")
93                KillDataFolder/Z $(toDF+":entry:data_FL")
94                KillDataFolder/Z $(toDF+":entry:data_FR")
95                KillDataFolder/Z $(toDF+":entry:data_FT")
96                KillDataFolder/Z $(toDF+":entry:data_FB")
97
98                return(0)
99        else
100       
101                V_KillWavesFullTree($toDF,toStr,0,"",1)                 // this will traverse the whole tree, trying to kill what it can
102
103                if(DataFolderExists("root:Packages:NIST:VSANS:"+toStr) == 0)            //if the data folder (RAW, SAM, etc.) was just killed?
104                        NewDataFolder/O $("root:Packages:NIST:VSANS:"+toStr)
105                endif   
106                       
107                // need to do this the hard way, duplicate/O recursively
108                // see V_CopyToWorkFolder()
109
110                // gFileList is above "entry" which is my additions
111                SVAR fileList_dest = $("root:Packages:NIST:VSANS:"+toStr+":gFileList")
112                SVAR fileList_tmp = $("root:Packages:NIST:VSANS:"+fromStr+":gFileList")
113                fileList_dest = fileList_tmp
114       
115                // everything on the top level
116                V_DuplicateDataFolder($(fromDF+":entry"),fromStr,toStr,0,"",0)  //no recursion here
117                // control
118                V_DuplicateDataFolder($(fromDF+":entry:control"),fromStr,toStr,0,"",1)  //yes recursion here
119                // instrument
120                V_DuplicateDataFolder($(fromDF+":entry:instrument"),fromStr,toStr,0,"",1)       //yes recursion here
121                // reduction
122                V_DuplicateDataFolder($(fromDF+":entry:reduction"),fromStr,toStr,0,"",1)        //yes recursion here
123                // sample
124                V_DuplicateDataFolder($(fromDF+":entry:sample"),fromStr,toStr,0,"",1)   //yes recursion here
125                // user
126                V_DuplicateDataFolder($(fromDF+":entry:user"),fromStr,toStr,0,"",1)     //yes recursion here
127               
128        endif   
129       
130        return(0)
131end
132
133
134
135Function V_KillWavesInFolder(folderStr)
136        String folderStr
137       
138        if(DataFolderExists(folderStr) && strlen(folderStr) != 0)
139                SetDataFolder $folderStr
140                KillWaves/A/Z
141        endif
142       
143        SetDataFolder root:
144        return(0)
145end
146
147////////
148// see the help entry for IndexedDir for help on (possibly) how to do this faster
149// -- see  Function ScanDirectories(pathName, printDirNames)
150//
151
152
153// from IgorExchange On July 17th, 2011 jtigor
154// started from "Recursively List Data Folder Contents"
155// Posted July 15th, 2011 by hrodstein
156//
157//
158//
159Proc V_CopyWorkFolderProc(dataFolderStr, fromStr, toStr, level, sNBName, recurse)
160        String dataFolderStr="root:Packages:NIST:VSANS:RAW"
161        String fromStr = "RAW"
162        String toStr="SAM"
163        Variable level=0
164        String sNBName="DataFolderTree"
165        Variable recurse = 1
166       
167        V_DuplicateDataFolder($dataFolderStr, fromStr, toStr, level, sNBName, recurse)
168
169
170end
171
172// ListDataFolder(dfr, level)
173// Recursively lists objects in data folder.
174// Pass data folder path for dfr and 0 for level.
175// pass level == 0 for the first call
176//  sNBName = "" prints nothing. any name will generate a notebook
177//
178// recurse == 0 will do only the specified folder, anything else will recurse all levels
179// toStr is the string name of the top-level folder only, not the full path
180//
181//
182Function V_DuplicateDataFolder(dfr, fromStr, toStr, level, sNBName,recurse)
183        DFREF dfr
184        String fromStr
185        String toStr
186        Variable level                  // Pass 0 to start
187        String sNBName
188        Variable recurse
189 
190        String name
191        String dfName
192        String sString
193       
194        String toDF = ""
195 
196        if (level == 0)         // this is the data folder, generate if needed in the destination
197                name = GetDataFolder(1, dfr)
198//              sPrintf sString, "%s (data folder)\r", name
199                toDF = ReplaceString(fromStr,name,toStr,1)              // case-sensitive replace
200                sprintf sString, "NewDataFolder/O %s\r",toDF
201                NewDataFolder/O $(RemoveEnding(toDF,":"))                       // remove trailing semicolon if it's there
202               
203                V_WriteBrowserInfo_test(sString, 1, sNBName)
204        endif
205 
206        dfName = GetDataFolder(1, dfr)
207        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
208        Variable i
209 
210        String indentStr = "\t"
211        for(i=0; i<level; i+=1)
212                indentStr += "\t"
213        endfor
214 
215        Variable numWaves = CountObjectsDFR(dfr, 1)
216        for(i=0; i<numWaves; i+=1)
217                name = GetIndexedObjNameDFR(dfr, 1, i)
218                //
219                // wave type does not matter now. Duplicate does not care
220                //
221                sPrintf sString, "Duplicate/O  %s,  %s\r",dfName+name,toDF+name
222                Duplicate/O $(dfName+name),$(toDF+name)
223               
224                V_WriteBrowserInfo_test(sString, 2, sNBName)
225        endfor 
226 
227        Variable numNumericVariables = CountObjectsDFR(dfr, 2) 
228        for(i=0; i<numNumericVariables; i+=1)
229                name = GetIndexedObjNameDFR(dfr, 2, i)
230                sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
231                V_WriteBrowserInfo_test(sString, 3, sNBName)
232        endfor 
233 
234        Variable numStringVariables = CountObjectsDFR(dfr, 3)   
235        for(i=0; i<numStringVariables; i+=1)
236                name = GetIndexedObjNameDFR(dfr, 3, i)
237                sPrintf sString, "%s%s (string variable)\r", indentStr, name
238                V_WriteBrowserInfo_test(sString, 4, sNBName)
239        endfor 
240
241        if(recurse)
242                Variable numDataFolders = CountObjectsDFR(dfr, 4)       
243                for(i=0; i<numDataFolders; i+=1)
244                        name = GetIndexedObjNameDFR(dfr, 4, i)
245//                      sPrintf sString, "%s%s (data folder)\r", indentStr, name
246                         dfName = GetDataFolder(1, dfr)
247                         
248                        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
249                        sprintf sString, "NewDataFolder/O %s\r",toDF+name
250                        NewDataFolder/O $(toDF+name)
251                       
252                       
253                        V_WriteBrowserInfo_test(sString, 1, sNBName)
254                        DFREF childDFR = dfr:$(name)
255                        V_DuplicateDataFolder(childDFR, fromStr, toStr, level+1, sNBName, recurse)
256                endfor 
257        endif
258         
259
260End
261
262
263// ListDataFolder(dfr, level)
264// Recursively lists objects in data folder.
265// Pass data folder path for dfr and 0 for level.
266// pass level == 0 for the first call
267//  sNBName = "" prints nothing. any name will generate a notebook
268//
269// recurse == 0 will do only the specified folder, anything else will recurse all levels
270// toStr is the string name of the top-level folder only, not the full path
271//
272//
273Function V_KillWavesFullTree(dfr, fromStr, level, sNBName,recurse)
274        DFREF dfr
275        String fromStr
276//      String toStr
277        Variable level                  // Pass 0 to start
278        String sNBName
279        Variable recurse
280 
281        String name
282        String dfName
283        String sString
284       
285        String toDF = ""
286 
287        if (level == 0)         // this is the data folder, generate if needed in the destination
288                name = GetDataFolder(1, dfr)
289                sPrintf sString, "%s (data folder)\r", name
290//              toDF = ReplaceString(fromStr,name,toStr,1)              // case-sensitive replace
291//              sprintf sString, "NewDataFolder/O %s\r",toDF
292//              NewDataFolder/O $(RemoveEnding(toDF,":"))                       // remove trailing semicolon if it's there
293               
294                V_WriteBrowserInfo_test(sString, 1, sNBName)
295        endif
296 
297        dfName = GetDataFolder(1, dfr)
298//      toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
299        Variable i
300 
301        String indentStr = "\t"
302        for(i=0; i<level; i+=1)
303                indentStr += "\t"
304        endfor
305 
306        Variable numWaves = CountObjectsDFR(dfr, 1)
307        for(i=0; i<numWaves; i+=1)
308                name = GetIndexedObjNameDFR(dfr, 1, i)
309                //
310                // wave type does not matter now. Kill does not care
311                //
312                sPrintf sString, "Killing  %s\r",dfName+name
313                KillWaves/Z $(dfName+name)
314               
315                V_WriteBrowserInfo_test(sString, 2, sNBName)
316        endfor 
317 
318 // now kill the data folder if possible
319        KillDataFolder/Z $dfName
320       
321       
322        Variable numNumericVariables = CountObjectsDFR(dfr, 2) 
323        for(i=0; i<numNumericVariables; i+=1)
324                name = GetIndexedObjNameDFR(dfr, 2, i)
325                sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
326                V_WriteBrowserInfo_test(sString, 3, sNBName)
327        endfor 
328 
329        Variable numStringVariables = CountObjectsDFR(dfr, 3)   
330        for(i=0; i<numStringVariables; i+=1)
331                name = GetIndexedObjNameDFR(dfr, 3, i)
332                sPrintf sString, "%s%s (string variable)\r", indentStr, name
333                V_WriteBrowserInfo_test(sString, 4, sNBName)
334        endfor 
335
336        if(recurse)
337                Variable numDataFolders = CountObjectsDFR(dfr, 4)       
338                for(i=0; i<numDataFolders; i+=1)
339                        name = GetIndexedObjNameDFR(dfr, 4, i)
340                        sPrintf sString, "%s%s (data folder)\r", indentStr, name
341                         dfName = GetDataFolder(1, dfr)
342                         
343//                      toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
344//                      sprintf sString, "NewDataFolder/O %s\r",toDF+name
345//                      NewDataFolder/O $(toDF+name)
346                       
347                       
348                        V_WriteBrowserInfo_test(sString, 1, sNBName)
349                        DFREF childDFR = dfr:$(name)
350                        V_KillWavesFullTree(childDFR, fromStr, level+1, sNBName, recurse)
351                endfor 
352        endif
353         
354
355End
356 
357
358
359Function V_WriteBrowserInfo_test(sString, vType, sNBName)
360        String sString
361        Variable vType
362        String sNBName
363 
364        if(strlen(sNBName) == 0)
365//              print sString
366                return 0
367        endif
368        DoWindow $sNBName
369        if(V_flag != 1)
370                NewNoteBook/F=0 /N=$sNBName /V=1 as sNBName
371        else
372                DoWindow/F $sNBName
373        endif
374        Notebook $sNBName selection={endOfFile, endOfFile}
375        if(vType == 1)          // a data folder
376//              Notebook $sNBName fstyle=1
377                Notebook $sNBName text=sString
378//              Notebook $sNBName fstyle=-1
379        else
380                Notebook $sNBName text=sString 
381        endif
382 
383End
384
385///////////////////////////////
386
387
388//
389// given the folder, duplicate the data -> linear_data and generate the error
390//
391// x- do I want to use different names here? If it turns out that I don't need to drag a copy of
392//    the data around as "linear_data", then I can eliminate that, and rename the error wave
393// x- be sure the data is either properly written as 2D in the file, or converted to 2D before
394//    duplicating here
395// x- ? do I recast to DP here? No- V_MakeDataWaves_DP() is called directly prior to this call, so data
396//     coming in is already DP
397Function V_MakeDataError(folderStr)
398        String folderStr
399       
400        SetDataFolder $folderStr
401        Wave data=data
402        Duplicate/O data linear_data                    // at this point, the data is still the raw data, and is linear_data
403       
404        // proper error for counting statistics, good for low count values too
405        // rather than just sqrt(n)
406        // see N. Gehrels, Astrophys. J., 303 (1986) 336-346, equation (7)
407        // for S = 1 in eq (7), this corresponds to one sigma error bars
408        Duplicate/O linear_data linear_data_error
409        linear_data_error = 1 + sqrt(linear_data + 0.75)
410       
411        Duplicate/O linear_data_error data_error                       
412        //
413       
414        SetDataFolder root:
415        return(0)
416End
417
418
419/////////////////////
420
421
422
423//testing procedure
424// TODO -- can't duplicate this with another proceudre, but if I change the name of the variable
425//   "newType" to "type", then when Raw_to_work() gets to CopyHDFToWorkFolder(), the KillDataFolder/Z
426//   line fails (but reports no error), then DuplicateDataFolder fails, and reports an error. Trying
427//   to simplify this condition, I can't reproduce the bug for WM...
428Proc V_Convert_to_Workfile(newtype, doadd)
429        String newtype,doadd
430        Prompt newtype,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
431        Prompt doadd,"Add to current WORK contents?",popup,"No;Yes;"
432       
433        //macro will take whatever is in RAW folder and "ADD" it to the folder specified
434        //in the popup menu
435       
436        //"add" = yes/no, don't add to previous runs
437        //switch here - two separate functions to avoid (my) confusion
438        Variable err// = Raw_to_work(newtype)
439        if(cmpstr(doadd,"No")==0)
440                //don't add to prev work contents, copy RAW contents to work and convert
441                err = V_Raw_to_work(newtype)
442        else
443                //yes, add RAW to the current work folder contents
444                //Abort "Adding RAW data files is currently unsupported"
445                err = V_Add_raw_to_work(newtype)
446        endif
447       
448        String newTitle = "WORK_"+newtype
449        DoWindow/F VSANS_Data
450        DoWindow/T VSANS_Data, newTitle
451        KillStrings/Z newTitle
452       
453        //need to update the display with "data" from the correct dataFolder
454        V_UpdateDisplayInformation(newtype)
455       
456End
457
458
459//
460// THIS IS THE MAJOR ROUTINE TO APPLY DATA CORRECTIONS
461//
462//will copy the current contents of the RAW folder to the newType work folder
463//and do the geometric corrections and normalization to monitor counts
464//(the function Add_Raw_to_work(type) adds multiple runs together - and is LOW priority)
465//
466// JUL 2018
467// now removes a constant kReadNoiseLevel from the back detector (this is a constant value, not time
468//  or count dependent). It does not appear to be dependent on gain, and is hoepfully stable over time.
469//
470//the current display type is updated to newType (global)
471//
472Function V_Raw_to_work(newType)
473        String newType
474       
475        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
476        Variable ii,jj,itim,cntrate,dscale,scale,uscale
477        String destPath
478       
479        String fname = newType
480        String detStr
481        Variable ctTime
482
483        //initialize values before normalization
484        total_mon=0
485        total_det=0
486        total_trn=0
487        total_numruns=0
488        total_rtime=0
489       
490        //Not adding multiple runs, so wipe out the old contents of the work folder and
491        // replace with the contents of raw
492
493        destPath = "root:Packages:NIST:VSANS:" + newType
494       
495        //copy from current dir (RAW) to work, defined by newType
496        V_CopyHDFToWorkFolder("RAW",newType)
497       
498        // now work with the waves from the destination folder.
499       
500        // apply corrections ---
501        // switches to control what is done, don't do the transmission correction for the BGD measurement
502        // start with the DIV correction, before conversion to mm
503        // then do all of the other corrections, order doesn't matter.
504        // rescaling to default monitor counts however, must be LAST.
505
506// each correction must loop over each detector. tedious.
507
508        //except for removing the read noise of the back detector
509        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
510
511        if(gIgnoreDetB == 0)
512                Wave w = V_getDetectorDataW(fname,"B")
513               
514                NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
515               
516                switch(gHighResBinning)
517                        case 1:
518                                w -= kReadNoiseLevel_bin1               // a constant value
519                               
520                                MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)
521                               
522                                Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
523                                break
524                        case 4:
525                                w -= kReadNoiseLevel_bin4               // a constant value
526                               
527                                MatrixFilter /N=3 /P=3 median w                 //              /P=n flag sets the number of passes (default is 1 pass)
528                               
529                                Print "*** median noise filter 3x3 applied to the back detector (3 passes) ***"
530                                break
531                        default:
532                                Abort "No binning case matches in V_Raw_to_Work"
533                endswitch       
534
535        endif
536       
537       
538        // (0) Redimension the data waves in the destination folder
539        //     so that they are DP, not integer
540        // TODO
541        // -- currently only redimensioning the data and linear_data_error - anything else???
542        // -- ?? some of this is done at load time for RAW data. Not a problem to re-do the redimension
543        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
544                detStr = StringFromList(ii, ksDetectorListAll, ";")
545                Wave w = V_getDetectorDataW(fname,detStr)
546                Wave w_err = V_getDetectorDataErrW(fname,detStr)
547                Redimension/D w,w_err
548        endfor
549       
550       
551        // (1) DIV correction
552        // do this in terms of pixels.
553        // TODO : This must also exist at the time the first work folder is generated.
554        //   So it must be in the user folder at the start of the experiment, and defined.
555        NVAR gDoDIVCor = root:Packages:NIST:VSANS:Globals:gDoDIVCor
556        if (gDoDIVCor == 1)
557                // need extra check here for file existence
558                // if not in DIV folder, load.
559                // if unable to load, skip correction and report error (Alert?) (Ask to Load?)
560                Print "Doing DIV correction"// for "+ detStr
561                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
562                        detStr = StringFromList(ii, ksDetectorListAll, ";")
563                       
564                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
565                                // do nothing
566                        else
567                                Wave w = V_getDetectorDataW(fname,detStr)
568                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
569                               
570                                V_DIVCorrection(w,w_err,detStr,newType)
571                        endif
572                endfor
573        else
574                Print "DIV correction NOT DONE"         // not an error since correction was unchecked
575        endif
576       
577        // (2) non-linear correction   
578        // TODO:
579        // x-  the "B" detector is calculated in its own routines
580        // -- document what is generated here:
581        //    **in each detector folder: data_realDistX and data_realDistY (2D waves of the [mm] position of each pixel)
582        // x- these spatial calculations ARE DONE as the RAW data is loaded. It allows the RAW
583        //    data to be properly displayed, but without all of the (complete) set of detector corrections
584        // * the corrected distances are calculated into arrays, but nothing is done with them yet
585        // * there is enough information now to calculate the q-arrays, so it is done now
586        // - other corrections may modify the data, this calculation does NOT modify the data
587        NVAR gDoNonLinearCor = root:Packages:NIST:VSANS:Globals:gDoNonLinearCor
588        // generate a distance matrix for each of the detectors
589        if (gDoNonLinearCor == 1)
590                Print "Doing Non-linear correction"// for "+ detStr
591                for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
592                        detStr = StringFromList(ii, ksDetectorListNoB, ";")
593                        Wave w = V_getDetectorDataW(fname,detStr)
594//                      Wave w_err = V_getDetectorDataErrW(fname,detStr)
595                        Wave w_calib = V_getDetTube_spatialCalib(fname,detStr)
596                        Variable tube_width = V_getDet_tubeWidth(fname,detStr)
597                        V_NonLinearCorrection(fname,w,w_calib,tube_width,detStr,destPath)
598
599                        //(2.4) Convert the beam center values from pixels to mm
600                        // TODO -- there needs to be a permanent location for these values??
601                        //
602                                // TODO
603                                // -- the beam center value in mm needs to be present - it is used in calculation of Qvalues
604                                // -- but having both the same is wrong...
605                                // -- the pixel value is needed for display of the panels
606                                if(kBCTR_CM)
607                                        //V_ConvertBeamCtrPix_to_mm(folder,detStr,destPath)
608                                        //
609       
610                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
611                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
612                                        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
613                                        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
614                                        x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
615                                        y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
616                                       
617                                        // TODO:::
618                                // now I need to convert the beam center in mm to pixels
619                                // and have some rational place to look for it...
620                                        V_ConvertBeamCtr_to_pix(fname,detStr,destPath)
621                                else
622                                        // beam center is in pixels, so use the old routine
623                                        V_ConvertBeamCtrPix_to_mm(fname,detStr,destPath)
624                                endif           
625                                                       
626                        // (2.5) Calculate the q-values
627                        // calculating q-values can't be done unless the non-linear corrections are calculated
628                        // so go ahead and put it in this loop.
629                        // TODO :
630                        // -- make sure that everything is present before the calculation
631                        // -- beam center must be properly defined in terms of real distance
632                        // -- distances/zero location/ etc. must be clearly documented for each detector
633                        //      ** this assumes that NonLinearCorrection() has been run to generate data_RealDistX and Y
634                        // ** this routine Makes the waves QTot, qx, qy, qz in each detector folder.
635                        //
636                        V_Detector_CalcQVals(fname,detStr,destPath)
637                       
638                endfor
639
640                if(gIgnoreDetB==1)             
641                        //"B" is separate
642                        detStr = "B"
643                        Wave w = V_getDetectorDataW(fname,detStr)
644                        Wave cal_x = V_getDet_cal_x(fname,detStr)
645                        Wave cal_y = V_getDet_cal_y(fname,detStr)
646                       
647                        V_NonLinearCorrection_B(fname,w,cal_x,cal_y,detStr,destPath)
648       
649        // "B" is always naturally defined in terms of a pixel center. This can be converted to mm,
650        // but the experiment will measure pixel x,y - just like ordela detectors.
651                       
652        //              if(kBCTR_CM)
653        //
654        //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
655        //                      Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
656        //                      WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
657        //                      WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
658        //                      x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
659        //                      y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
660        //                     
661        //                      // now I need to convert the beam center in mm to pixels
662        //                      // and have some rational place to look for it...
663        //                      V_ConvertBeamCtr_to_pixB(fname,detStr,destPath)
664        //              else
665                                // beam center is in pixels, so use the old routine
666                                V_ConvertBeamCtrPix_to_mmB(fname,detStr,destPath)
667       
668        //              endif
669                        V_Detector_CalcQVals(fname,detStr,destPath)
670               
671                endif
672               
673        else
674                Print "Non-linear correction NOT DONE"
675        endif
676
677
678        // (3) dead time correction
679        // DONE:
680        // x- test for correct operation
681        // x- loop over all of the detectors
682        // x- B detector is a special case (do separately, then loop over NoB)
683        // x- this DOES alter the data
684        // x- verify the error propagation (not done yet)
685        //
686        Variable countRate
687        NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
688        if (gDoDeadTimeCor == 1)
689                Print "Doing DeadTime correction"// for "+ detStr
690                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
691                        detStr = StringFromList(ii, ksDetectorListAll, ";")
692                        Wave w = V_getDetectorDataW(fname,detStr)
693                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
694                        ctTime = V_getCount_time(fname)
695
696                        if(cmpstr(detStr,"B") == 0 )
697                                if(gIgnoreDetB == 0)
698                                        Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
699                                        // do the correction for the back panel
700                                        countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts
701       
702                                        w = w/(1-countRate*b_dt)
703                                        w_err = w_err/(1-countRate*b_dt)
704                                endif                   
705                        else
706                                // do the corrections for 8 tube panels
707                                Wave w_dt = V_getDetector_deadtime(fname,detStr)
708                                V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
709                               
710                        endif
711                endfor
712               
713        else
714                Print "Dead Time correction NOT DONE"
715        endif   
716       
717
718        // (4) solid angle correction
719        //  -- this currently calculates the correction factor AND applys it to the data
720        //  -- as a result, the data values are very large since they are divided by a very small
721        //     solid angle per pixel. But all of the count values are now on the basis of
722        //    counts/(solid angle) --- meaning that they can all be binned together for I(q)
723        //    -and- - this is taken into account for absolute scaling (this part is already done)
724        //
725        // the solid angle correction is calculated for ALL detector panels.
726        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
727        if (gDoSolidAngleCor == 1)
728                Print "Doing Solid Angle correction"// for "+ detStr
729                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
730                        detStr = StringFromList(ii, ksDetectorListAll, ";")
731                       
732                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
733                                // do nothing
734                        else
735                                Wave w = V_getDetectorDataW(fname,detStr)
736                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
737                                // any other dimensions to pass in?
738                                V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
739                        endif
740                       
741                endfor
742        else
743                Print "Solid Angle correction NOT DONE"
744        endif   
745       
746               
747        // (5) angle-dependent tube shadowing
748        // TODO:
749        // -- not sure about this correction yet...
750        //
751        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
752        if (gDoTubeShadowCor == 1)
753                Print "(stub)Tube Shadow correction"
754        else
755                Print "Tube shadowing correction NOT DONE"
756        endif   
757               
758        // (6) angle dependent transmission correction
759        // TODO:
760        // x- this still needs to be filled in
761        // -- still some debate of when/where in the corrections that this is best applied
762        //    - do it here, and it's done whether the output is 1D or 2D
763        //    - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
764        // -- verify that the calculation is correct
765        // -- verify that the error propagation (in 2D) is correct
766        //
767        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
768        if (gDoTrans == 1)
769                Print "Doing Large-angle transmission correction"// for "+ detStr
770                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
771                        detStr = StringFromList(ii, ksDetectorListAll, ";")
772                       
773                        if(cmpstr(detStr,"B") == 0  && gIgnoreDetB == 1)
774                                // do nothing
775                        else
776                                Wave w = V_getDetectorDataW(fname,detStr)
777                                Wave w_err = V_getDetectorDataErrW(fname,detStr)
778                               
779                                V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
780                        endif
781                endfor
782        else
783                Print "Sample Transmission correction NOT DONE"
784        endif   
785       
786        // (7) normalize to default monitor counts
787        // TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
788        // TODO -- but there are TWO monitors - so how to switch?
789        //  --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
790        // TODO -- what do I really need to save?
791       
792        NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
793        if (gDoMonitorNormalization == 1)
794               
795                Variable monCount,savedMonCount
796                defmon=1e8                      //default monitor counts
797//              monCount = V_getControlMonitorCount(fname)                      // TODO -- this is read in since VCALC fakes this on output
798                monCount = V_getBeamMonNormData(fname)          // TODO -- I think this is the *real* one to read
799                savedMonCount   = monCount
800                scale = defMon/monCount         // scale factor to MULTIPLY data by to rescale to defmon
801
802                // write to newType=fname will put new values in the destination WORK folder
803                V_writeBeamMonNormSaved_count(fname,savedMonCount)                      // save the true count
804                V_writeBeamMonNormData(fname,defMon)            // mon ct is now 10^8
805                                       
806//                      // TODO
807//                      // the low efficiency monitor, expect to use this for white beam mode
808//                              -- need code switch here to determine which monitor to use.
809//
810//                      V_getBeamMonLowData(fname)
811//                      V_getBeamMonLowSaved_count(fname)       
812                       
813                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
814                        detStr = StringFromList(ii, ksDetectorListAll, ";")
815                        Wave w = V_getDetectorDataW(fname,detStr)
816                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
817
818                        // do the calculation right here. It's a simple scaling and not worth sending to another function.     
819                        //scale the data and error to the default monitor counts
820               
821//
822                        w *= scale
823                        w_err *= scale          //assumes total monitor count is so large there is essentially no error
824
825                        // TODO
826                        // -- do I want to update and save the integrated detector count?
827                        Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
828                        V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
829
830                endfor
831        else
832                Print "Monitor Normalization correction NOT DONE"
833        endif
834       
835       
836        // (not done) angle dependent efficiency correction
837        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
838
839//     
840        //reset the current displaytype to "newtype"
841        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
842       
843        //return to root folder (redundant)
844        SetDataFolder root:
845       
846        Return(0)
847End
848
849//
850//will "ADD" the current contents of the RAW folder to the newType work folder
851//
852// - used when adding multiple runs together
853//(the function V_Raw_to_work(type) makes a fresh workfile)
854//
855//the current display type is updated to newType (global)
856Function V_Add_raw_to_work(newType)
857        String newType
858       
859        String destPath=""
860        // if the desired workfile doesn't exist, let the user know, and just make a new one
861        if(WaveExists($("root:Packages:NIST:VSANS:"+newType + ":entry:instrument:detector_FL:data")) == 0)
862                Print "There is no old work file to add to - a new one will be created"
863                //call V_Raw_to_work(), then return from this function
864                V_Raw_to_Work(newType)
865                Return(0)               //does not generate an error - a single file was converted to work.newtype
866        Endif
867
868
869        // convert the RAW data to a WORK file.
870        // this will do all of the necessary corrections to the data
871        // put this in some separate work folder that can be cleaned out at the end (ADJ)
872        String tmpType="ADJ"
873       
874        //this step removes the read noise from the back so that neither added file will have this constant
875        V_Raw_to_Work(tmpType)         
876                       
877        //now make references to data in newType folder
878        DestPath="root:Packages:NIST:VSANS:"+newType   
879
880
881/////////////////
882//fields that need to be added together
883// entry block
884        // collection_time              V_getCollectionTime(fname)              V_putCollectionTime(fname,val)
885
886// instrument block
887        // beam_monitor_norm
888                // data (this will be 1e8)                              V_getBeamMonNormData(fname)             V_putBeamMonNormData(fname,val)
889                // saved_count (this is the original monitor count)  V_getBeamMonNormSaved_count(fname)         V_putBeamMonNormSaved_count(fname,val)
890
891        // for each detector
892        // data         V_getDetectorDataW(fname,detStr)
893        // integrated_count             V_getDet_IntegratedCount(fname,detStr)   V_putDet_IntegratedCount(fname,detStr,val)
894        // linear_data           V_getDetectorLinearDataW(fname,detStr)
895        // RECALCULATE (or add properly) linear_data_error              V_getDetectorDataErrW(fname,detStr)
896
897
898// control block (these may not actually be used?)
899        // count_time                           V_getCount_time(fname)                                          V_putCount_time(fname,val)
900        // detector_counts              V_getDetector_counts(fname)                             V_putDetector_counts(fname,val)
901        // monitor_counts               V_getControlMonitorCount(fname)         V_putControlMonitorCount(fname,val)
902
903// sample block - nothing
904// reduction block - nothing
905// user block - nothing
906
907// ?? need to add the file name to a list of what was actually added - so it will be saved with I(q)
908//
909////////////////////
910
911//      total_mon = realsread[1]        //saved monitor count
912//      uscale = total_mon/defmon               //unscaling factor
913//      total_det = uscale*realsread[2]         //unscaled detector count
914//      total_trn = uscale*realsread[39]        //unscaled trans det count
915//      total_numruns = integersread[3] //number of runs in workfile
916//      total_rtime = integersread[2]           //total counting time in workfile
917       
918
919        String detStr
920       
921        Variable saved_mon_dest,scale_dest,saved_mon_tmp,scale_tmp
922        Variable collection_time_dest,collection_time_tmp,count_time_dest,count_time_tmp
923        Variable detCount_dest,detCount_tmp,det_integrated_ct_dest,det_integrated_ct_tmp
924        Variable ii,new_scale,defMon
925       
926        defMon=1e8                      //default monitor counts
927
928        // find the scaling factors, one for each folder
929        saved_mon_dest = V_getBeamMonNormSaved_count(newType)
930        scale_dest = saved_mon_dest/defMon              //un-scaling factor
931       
932        saved_mon_tmp = V_getBeamMonNormSaved_count(tmpType)
933        scale_tmp = saved_mon_tmp/defMon                        //un-scaling factor
934
935        new_scale = defMon / (saved_mon_dest+saved_mon_tmp)
936       
937       
938        // get the count time for each (two locations)
939        collection_time_dest = V_getCollectionTime(newType)
940        collection_time_tmp = V_getCollectionTime(tmpType)
941       
942        count_time_dest = V_getCount_time(newType)
943        count_time_tmp = V_getCount_time(tmpType)
944       
945        detCount_dest = V_getDetector_counts(newType)
946        detCount_tmp = V_getDetector_counts(tmpType)
947
948// update the fields that are not in the detector blocks
949// in entry
950        V_putCollectionTime(newType,collection_time_dest+collection_time_tmp)
951
952// in control block
953        V_putCount_time(newType,count_time_dest+count_time_tmp)
954        V_putDetector_counts(newType,detCount_dest+detCount_tmp)
955        V_putControlMonitorCount(newType,saved_mon_dest+saved_mon_tmp)
956
957
958// now loop over all of the detector panels
959        // data
960        // data_err
961        // integrated count
962        // linear_data
963        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
964                detStr = StringFromList(ii, ksDetectorListAll, ";")
965               
966                Wave data_dest = V_getDetectorDataW(newType,detStr)
967                Wave data_err_dest = V_getDetectorDataErrW(newType,detStr)
968                Wave linear_data_dest = V_getDetectorLinearDataW(newType,detStr)
969                det_integrated_ct_dest = V_getDet_IntegratedCount(newType,detStr)
970
971                Wave data_tmp = V_getDetectorDataW(tmpType,detStr)
972                Wave data_err_tmp = V_getDetectorDataErrW(tmpType,detStr)
973                Wave linear_data_tmp = V_getDetectorLinearDataW(tmpType,detStr)
974                det_integrated_ct_tmp = V_getDet_IntegratedCount(tmpType,detStr)
975               
976                // unscale the data arrays
977                data_dest *= scale_dest
978                data_err_dest *= scale_dest
979                linear_data_dest *= scale_dest
980               
981                data_tmp *= scale_tmp
982                data_err_tmp *= scale_tmp
983                linear_data_tmp *= scale_tmp
984                               
985                // add them together, the dest is a wave so it is automatically changed in the "dest" folder
986                V_putDet_IntegratedCount(tmpType,detStr,det_integrated_ct_dest+det_integrated_ct_tmp)
987                data_dest += data_tmp
988                data_err_dest = sqrt(data_err_dest^2 + data_err_tmp^2)          // add in quadrature
989                linear_data_dest += linear_data_tmp
990               
991                // now rescale the data_dest to the monitor counts
992                data_dest *= new_scale
993                data_err_dest *= new_scale
994                linear_data_dest *= new_scale
995               
996        endfor
997
998       
999        //Add the added raw filename to the list of files in the workfile
1000        SVAR fileList_dest = $("root:Packages:NIST:VSANS:"+newType+":gFileList")
1001        SVAR fileList_tmp = $("root:Packages:NIST:VSANS:"+tmpType+":gFileList")
1002       
1003        fileList_dest += ";" + fileList_tmp
1004       
1005        //reset the current display type to "newtype"
1006        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1007        gCurDispType = newType
1008       
1009        //return to root folder (redundant)
1010        SetDataFolder root:
1011       
1012        Return(0)
1013End
1014
1015
1016//used for adding DRK (beam shutter CLOSED) data to a workfile
1017//force the monitor count to 1, since it's irrelevant
1018// run data through normal "add" step, then unscale default monitor counts
1019//to get the data back on a simple time basis
1020//
1021Function V_Raw_to_Work_NoNorm(type)
1022        String type
1023       
1024        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1025        reals[1]=1              //true monitor counts, still in raw
1026        V_Raw_to_work(type)
1027        //data is now in "type" folder
1028        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1029        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1030        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1031        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1032       
1033        Variable norm_mon,tot_mon,scale
1034       
1035        norm_mon = new_reals[0]         //should be 1e8
1036        tot_mon = new_reals[1]          //should be 1
1037        scale= norm_mon/tot_mon
1038       
1039        data /= scale           //unscale the data
1040        data_err /= scale
1041       
1042        // to keep "data" and linear_data in sync
1043        data_copy = data
1044       
1045        return(0)
1046End
1047
1048//used for adding DRK (beam shutter CLOSED) data to a workfile
1049//force the monitor count to 1, since it's irrelevant
1050// run data through normal "add" step, then unscale default monitor counts
1051//to get the data back on a simple time basis
1052//
1053Function V_Add_Raw_to_Work_NoNorm(type)
1054        String type
1055       
1056        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1057        reals[1]=1              //true monitor counts, still in raw
1058        V_Add_Raw_to_work(type)
1059        //data is now in "type" folder
1060        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1061        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1062        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1063        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1064       
1065        Variable norm_mon,tot_mon,scale
1066       
1067        norm_mon = new_reals[0]         //should be 1e8
1068        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
1069        scale= norm_mon/tot_mon
1070       
1071        data /= scale           //unscale the data
1072        data_err /= scale
1073       
1074        // to keep "data" and linear_data in sync
1075        data_copy = data
1076       
1077        return(0)
1078End
1079
Note: See TracBrowser for help on using the repository browser.