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

Last change on this file since 1074 was 1073, checked in by srkline, 5 years ago

lots of changes here:
many little fixes to clean up TODO items and marke them DONE

changed the handling of the panel "gap" to split the gap evenly. Q-calculations have been re-verified with this change.

re-named the list of "bin Type" values, and added a few more choices. Streamlined how the averaging and plotting works with this list so that it can be more easily modified as different combinations of binning are envisioned. This resulted in a lot of excess code being cut out and replaced with cleaner logic. This change has also been verified to work as intended.

Attenuation is now always calculated from the table. The table also by (NEW) definition has values for the white beam (one waelength) and graphite (multiple possible wavelengths) where the wavelengths are artificially scaled (*1000) or *1e6) so that the interpolations can be done internally without the need for multiple attenuator tables.

File size: 34.9 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                // everything on the top level
111                V_DuplicateDataFolder($(fromDF+":entry"),fromStr,toStr,0,"",0)  //no recursion here
112                // control
113                V_DuplicateDataFolder($(fromDF+":entry:control"),fromStr,toStr,0,"",1)  //yes recursion here
114                // instrument
115                V_DuplicateDataFolder($(fromDF+":entry:instrument"),fromStr,toStr,0,"",1)       //yes recursion here
116                // reduction
117                V_DuplicateDataFolder($(fromDF+":entry:reduction"),fromStr,toStr,0,"",1)        //yes recursion here
118                // sample
119                V_DuplicateDataFolder($(fromDF+":entry:sample"),fromStr,toStr,0,"",1)   //yes recursion here
120                // user
121                V_DuplicateDataFolder($(fromDF+":entry:user"),fromStr,toStr,0,"",1)     //yes recursion here
122               
123        endif   
124       
125        return(0)
126end
127
128
129
130Function V_KillWavesInFolder(folderStr)
131        String folderStr
132       
133        if(DataFolderExists(folderStr) && strlen(folderStr) != 0)
134                SetDataFolder $folderStr
135                KillWaves/A/Z
136        endif
137       
138        SetDataFolder root:
139        return(0)
140end
141
142////////
143// see the help entry for IndexedDir for help on (possibly) how to do this faster
144// -- see  Function ScanDirectories(pathName, printDirNames)
145//
146
147
148// from IgorExchange On July 17th, 2011 jtigor
149// started from "Recursively List Data Folder Contents"
150// Posted July 15th, 2011 by hrodstein
151//
152//
153//
154Proc V_CopyWorkFolderProc(dataFolderStr, fromStr, toStr, level, sNBName, recurse)
155        String dataFolderStr="root:Packages:NIST:VSANS:RAW"
156        String fromStr = "RAW"
157        String toStr="SAM"
158        Variable level=0
159        String sNBName="DataFolderTree"
160        Variable recurse = 1
161       
162        V_DuplicateDataFolder($dataFolderStr, fromStr, toStr, level, sNBName, recurse)
163
164
165end
166
167// ListDataFolder(dfr, level)
168// Recursively lists objects in data folder.
169// Pass data folder path for dfr and 0 for level.
170// pass level == 0 for the first call
171//  sNBName = "" prints nothing. any name will generate a notebook
172//
173// recurse == 0 will do only the specified folder, anything else will recurse all levels
174// toStr is the string name of the top-level folder only, not the full path
175//
176//
177Function V_DuplicateDataFolder(dfr, fromStr, toStr, level, sNBName,recurse)
178        DFREF dfr
179        String fromStr
180        String toStr
181        Variable level                  // Pass 0 to start
182        String sNBName
183        Variable recurse
184 
185        String name
186        String dfName
187        String sString
188       
189        String toDF = ""
190 
191        if (level == 0)         // this is the data folder, generate if needed in the destination
192                name = GetDataFolder(1, dfr)
193//              sPrintf sString, "%s (data folder)\r", name
194                toDF = ReplaceString(fromStr,name,toStr,1)              // case-sensitive replace
195                sprintf sString, "NewDataFolder/O %s\r",toDF
196                NewDataFolder/O $(RemoveEnding(toDF,":"))                       // remove trailing semicolon if it's there
197               
198                V_WriteBrowserInfo_test(sString, 1, sNBName)
199        endif
200 
201        dfName = GetDataFolder(1, dfr)
202        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
203        Variable i
204 
205        String indentStr = "\t"
206        for(i=0; i<level; i+=1)
207                indentStr += "\t"
208        endfor
209 
210        Variable numWaves = CountObjectsDFR(dfr, 1)
211        for(i=0; i<numWaves; i+=1)
212                name = GetIndexedObjNameDFR(dfr, 1, i)
213                //
214                // wave type does not matter now. Duplicate does not care
215                //
216                sPrintf sString, "Duplicate/O  %s,  %s\r",dfName+name,toDF+name
217                Duplicate/O $(dfName+name),$(toDF+name)
218               
219                V_WriteBrowserInfo_test(sString, 2, sNBName)
220        endfor 
221 
222        Variable numNumericVariables = CountObjectsDFR(dfr, 2) 
223        for(i=0; i<numNumericVariables; i+=1)
224                name = GetIndexedObjNameDFR(dfr, 2, i)
225                sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
226                V_WriteBrowserInfo_test(sString, 3, sNBName)
227        endfor 
228 
229        Variable numStringVariables = CountObjectsDFR(dfr, 3)   
230        for(i=0; i<numStringVariables; i+=1)
231                name = GetIndexedObjNameDFR(dfr, 3, i)
232                sPrintf sString, "%s%s (string variable)\r", indentStr, name
233                V_WriteBrowserInfo_test(sString, 4, sNBName)
234        endfor 
235
236        if(recurse)
237                Variable numDataFolders = CountObjectsDFR(dfr, 4)       
238                for(i=0; i<numDataFolders; i+=1)
239                        name = GetIndexedObjNameDFR(dfr, 4, i)
240//                      sPrintf sString, "%s%s (data folder)\r", indentStr, name
241                         dfName = GetDataFolder(1, dfr)
242                         
243                        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
244                        sprintf sString, "NewDataFolder/O %s\r",toDF+name
245                        NewDataFolder/O $(toDF+name)
246                       
247                       
248                        V_WriteBrowserInfo_test(sString, 1, sNBName)
249                        DFREF childDFR = dfr:$(name)
250                        V_DuplicateDataFolder(childDFR, fromStr, toStr, level+1, sNBName, recurse)
251                endfor 
252        endif
253         
254
255End
256
257
258// ListDataFolder(dfr, level)
259// Recursively lists objects in data folder.
260// Pass data folder path for dfr and 0 for level.
261// pass level == 0 for the first call
262//  sNBName = "" prints nothing. any name will generate a notebook
263//
264// recurse == 0 will do only the specified folder, anything else will recurse all levels
265// toStr is the string name of the top-level folder only, not the full path
266//
267//
268Function V_KillWavesFullTree(dfr, fromStr, level, sNBName,recurse)
269        DFREF dfr
270        String fromStr
271//      String toStr
272        Variable level                  // Pass 0 to start
273        String sNBName
274        Variable recurse
275 
276        String name
277        String dfName
278        String sString
279       
280        String toDF = ""
281 
282        if (level == 0)         // this is the data folder, generate if needed in the destination
283                name = GetDataFolder(1, dfr)
284                sPrintf sString, "%s (data folder)\r", name
285//              toDF = ReplaceString(fromStr,name,toStr,1)              // case-sensitive replace
286//              sprintf sString, "NewDataFolder/O %s\r",toDF
287//              NewDataFolder/O $(RemoveEnding(toDF,":"))                       // remove trailing semicolon if it's there
288               
289                V_WriteBrowserInfo_test(sString, 1, sNBName)
290        endif
291 
292        dfName = GetDataFolder(1, dfr)
293//      toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
294        Variable i
295 
296        String indentStr = "\t"
297        for(i=0; i<level; i+=1)
298                indentStr += "\t"
299        endfor
300 
301        Variable numWaves = CountObjectsDFR(dfr, 1)
302        for(i=0; i<numWaves; i+=1)
303                name = GetIndexedObjNameDFR(dfr, 1, i)
304                //
305                // wave type does not matter now. Kill does not care
306                //
307                sPrintf sString, "Killing  %s\r",dfName+name
308                KillWaves/Z $(dfName+name)
309               
310                V_WriteBrowserInfo_test(sString, 2, sNBName)
311        endfor 
312 
313 // now kill the data folder if possible
314        KillDataFolder/Z $dfName
315       
316       
317        Variable numNumericVariables = CountObjectsDFR(dfr, 2) 
318        for(i=0; i<numNumericVariables; i+=1)
319                name = GetIndexedObjNameDFR(dfr, 2, i)
320                sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
321                V_WriteBrowserInfo_test(sString, 3, sNBName)
322        endfor 
323 
324        Variable numStringVariables = CountObjectsDFR(dfr, 3)   
325        for(i=0; i<numStringVariables; i+=1)
326                name = GetIndexedObjNameDFR(dfr, 3, i)
327                sPrintf sString, "%s%s (string variable)\r", indentStr, name
328                V_WriteBrowserInfo_test(sString, 4, sNBName)
329        endfor 
330
331        if(recurse)
332                Variable numDataFolders = CountObjectsDFR(dfr, 4)       
333                for(i=0; i<numDataFolders; i+=1)
334                        name = GetIndexedObjNameDFR(dfr, 4, i)
335                        sPrintf sString, "%s%s (data folder)\r", indentStr, name
336                         dfName = GetDataFolder(1, dfr)
337                         
338//                      toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
339//                      sprintf sString, "NewDataFolder/O %s\r",toDF+name
340//                      NewDataFolder/O $(toDF+name)
341                       
342                       
343                        V_WriteBrowserInfo_test(sString, 1, sNBName)
344                        DFREF childDFR = dfr:$(name)
345                        V_KillWavesFullTree(childDFR, fromStr, level+1, sNBName, recurse)
346                endfor 
347        endif
348         
349
350End
351 
352
353
354Function V_WriteBrowserInfo_test(sString, vType, sNBName)
355        String sString
356        Variable vType
357        String sNBName
358 
359        if(strlen(sNBName) == 0)
360//              print sString
361                return 0
362        endif
363        DoWindow $sNBName
364        if(V_flag != 1)
365                NewNoteBook/F=0 /N=$sNBName /V=1 as sNBName
366        else
367                DoWindow/F $sNBName
368        endif
369        Notebook $sNBName selection={endOfFile, endOfFile}
370        if(vType == 1)          // a data folder
371//              Notebook $sNBName fstyle=1
372                Notebook $sNBName text=sString
373//              Notebook $sNBName fstyle=-1
374        else
375                Notebook $sNBName text=sString 
376        endif
377 
378End
379
380///////////////////////////////
381
382
383//
384// given the folder, duplicate the data -> linear_data and generate the error
385// TODO
386// -- do I want to use different names here? If it turns out that I don't need to drag a copy of
387//    the data around as "linear_data", then I can eliminate that, and rename the error wave
388// -- be sure the data is either properly written as 2D in the file, or converted to 2D before
389//    duplicating here
390// -- ? do I recast to DP here. Probably necessary since I'm doing a DP calculation, but Redimension
391//    is done in the Raw_to_Work step too. very confusing.
392Function V_MakeDataError(folderStr)
393        String folderStr
394       
395        SetDataFolder $folderStr
396        Wave data=data
397        Duplicate/O data linear_data                    // at this point, the data is still the raw data, and is linear_data
398       
399        // proper error for counting statistics, good for low count values too
400        // rather than just sqrt(n)
401        // see N. Gehrels, Astrophys. J., 303 (1986) 336-346, equation (7)
402        // for S = 1 in eq (7), this corresponds to one sigma error bars
403        Duplicate/O linear_data linear_data_error
404        linear_data_error = 1 + sqrt(linear_data + 0.75)                               
405        //
406       
407        SetDataFolder root:
408        return(0)
409End
410
411
412/////////////////////
413
414
415
416//testing procedure
417// TODO -- can't duplicate this with another proceudre, but if I change the name of the variable
418//   "newType" to "type", then when Raw_to_work() gets to CopyHDFToWorkFolder(), the KillDataFolder/Z
419//   line fails (but reports no error), then DuplicateDataFolder fails, and reports an error. Trying
420//   to simplify this condition, I can't reproduce the bug for WM...
421Proc V_Convert_to_Workfile(newtype, doadd)
422        String newtype,doadd
423        Prompt newtype,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
424        Prompt doadd,"Add to current WORK contents?",popup,"No;Yes;"
425       
426        //macro will take whatever is in RAW folder and "ADD" it to the folder specified
427        //in the popup menu
428       
429        //"add" = yes/no, don't add to previous runs
430        //switch here - two separate functions to avoid (my) confusion
431        Variable err// = Raw_to_work(newtype)
432        if(cmpstr(doadd,"No")==0)
433                //don't add to prev work contents, copy RAW contents to work and convert
434                err = V_Raw_to_work(newtype)
435        else
436                //yes, add RAW to the current work folder contents
437                Abort "Adding RAW data files is currently unsupported"
438                err = V_Add_raw_to_work(newtype)
439        endif
440       
441        String newTitle = "WORK_"+newtype
442        DoWindow/F VSANS_Data
443        DoWindow/T VSANS_Data, newTitle
444        KillStrings/Z newTitle
445       
446        //need to update the display with "data" from the correct dataFolder
447        V_UpdateDisplayInformation(newtype)
448       
449End
450
451
452//
453// THIS IS THE MAJOR ROUTINE TO APPLY DATA CORRECTIONS
454//
455//will copy the current contents of the RAW folder to the newType work folder
456//and do the geometric corrections and normalization to monitor counts
457//(the function Add_Raw_to_work(type) adds multiple runs together - and is LOW priority)
458//
459//the current display type is updated to newType (global)
460//
461Function V_Raw_to_work(newType)
462        String newType
463       
464        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
465        Variable ii,jj,itim,cntrate,dscale,scale,uscale
466        String destPath
467       
468        String fname = newType
469        String detStr
470        Variable ctTime
471
472        //initialize values before normalization
473        total_mon=0
474        total_det=0
475        total_trn=0
476        total_numruns=0
477        total_rtime=0
478       
479        //Not adding multiple runs, so wipe out the old contents of the work folder and
480        // replace with the contents of raw
481
482        destPath = "root:Packages:NIST:VSANS:" + newType
483       
484        //copy from current dir (RAW) to work, defined by newType
485        V_CopyHDFToWorkFolder("RAW",newType)
486       
487        // now work with the waves from the destination folder.
488       
489        // apply corrections ---
490        // switches to control what is done, don't do the transmission correction for the BGD measurement
491        // start with the DIV correction, before conversion to mm
492        // then do all of the other corrections, order doesn't matter.
493        // rescaling to default monitor counts however, must be LAST.
494
495// each correction must loop over each detector. tedious.
496
497        // (0) Redimension the data waves in the destination folder
498        //     so that they are DP, not integer
499        // TODO
500        // -- currently only redimensioning the data and linear_data_error - What else???
501        // -- ?? some of this is done at load time for RAW data. shouldn't be an issue to re-do the redimension
502        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
503                detStr = StringFromList(ii, ksDetectorListAll, ";")
504                Wave w = V_getDetectorDataW(fname,detStr)
505                Wave w_err = V_getDetectorDataErrW(fname,detStr)
506                Redimension/D w,w_err
507        endfor
508       
509       
510        // (1) DIV correction
511        // do this in terms of pixels.
512        // TODO : This must also exist at the time the first work folder is generated.
513        //   So it must be in the user folder at the start of the experiment, and defined.
514        NVAR gDoDIVCor = root:Packages:NIST:VSANS:Globals:gDoDIVCor
515        if (gDoDIVCor == 1)
516                // need extra check here for file existence
517                // if not in DIV folder, load.
518                // if unable to load, skip correction and report error (Alert?) (Ask to Load?)
519                Print "Doing DIV correction"// for "+ detStr
520                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
521                        detStr = StringFromList(ii, ksDetectorListAll, ";")
522                        Wave w = V_getDetectorDataW(fname,detStr)
523                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
524                       
525                        V_DIVCorrection(w,w_err,detStr,newType)
526                endfor
527        else
528                Print "DIV correction NOT DONE"         // not an error since correction was unchecked
529        endif
530       
531        // (2) non-linear correction   
532        // TODO:
533        // x-  the "B" detector is calculated in its own routines
534        // -- document what is generated here:
535        //    **in each detector folder: data_realDistX and data_realDistY (2D waves of the mm? position of the pixel)
536        // x- these spatial calculations ARE DONE as the RAW data is loaded. It allows the RAW
537        //    data to be properly displayed, but without all of the (complete) set of detector corrections
538        // * the corrected distances are calculated into arrays, but nothing is done with them yet
539        // * there is enough information now to calculate the q-arrays, so it is done now
540        // - other corrections may modify the data, this calculation does NOT modify the data
541        NVAR gDoNonLinearCor = root:Packages:NIST:VSANS:Globals:gDoNonLinearCor
542        // generate a distance matrix for each of the detectors
543        if (gDoNonLinearCor == 1)
544                Print "Doing Non-linear correction"// for "+ detStr
545                for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
546                        detStr = StringFromList(ii, ksDetectorListNoB, ";")
547                        Wave w = V_getDetectorDataW(fname,detStr)
548//                      Wave w_err = V_getDetectorDataErrW(fname,detStr)
549                        Wave w_calib = V_getDetTube_spatialCalib(fname,detStr)
550                        Variable tube_width = V_getDet_tubeWidth(fname,detStr)
551                        V_NonLinearCorrection(fname,w,w_calib,tube_width,detStr,destPath)
552
553                        //(2.4) Convert the beam center values from pixels to mm
554                        // TODO -- there needs to be a permanent location for these values??
555                        //
556                                // TODO
557                                // -- the beam center value in mm needs to be present - it is used in calculation of Qvalues
558                                // -- but having both the same is wrong...
559                                // -- the pixel value is needed for display of the panels
560                                if(kBCTR_CM)
561                                        //V_ConvertBeamCtr_to_mm(folder,detStr,destPath)
562                                        //
563       
564                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
565                                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
566                                        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
567                                        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
568                                        x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
569                                        y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
570                                       
571                                        // TODO:::
572                                // now I need to convert the beam center in mm to pixels
573                                // and have some rational place to look for it...
574                                        V_ConvertBeamCtr_to_pix(fname,detStr,destPath)
575                                else
576                                        // beam center is in pixels, so use the old routine
577                                        V_ConvertBeamCtr_to_mm(fname,detStr,destPath)
578                                endif           
579                                                       
580                        // (2.5) Calculate the q-values
581                        // calculating q-values can't be done unless the non-linear corrections are calculated
582                        // so go ahead and put it in this loop.
583                        // TODO :
584                        // -- make sure that everything is present before the calculation
585                        // -- beam center must be properly defined in terms of real distance
586                        // -- distances/zero location/ etc. must be clearly documented for each detector
587                        //      ** this assumes that NonLinearCorrection() has been run to generate data_RealDistX and Y
588                        // ** this routine Makes the waves QTot, qx, qy, qz in each detector folder.
589                        //
590                        V_Detector_CalcQVals(fname,detStr,destPath)
591                       
592                endfor
593               
594                //"B" is separate
595                Wave w = V_getDetectorDataW(fname,"B")
596                Wave cal_x = V_getDet_cal_x(fname,"B")
597                Wave cal_y = V_getDet_cal_y(fname,"B")
598                V_NonLinearCorrection_B(fname,w,cal_x,cal_y,"B",destPath)
599                V_ConvertBeamCtr_to_mmB(fname,"B",destPath)
600                V_Detector_CalcQVals(fname,"B",destPath)
601               
602        else
603                Print "Non-linear correction NOT DONE"
604        endif
605
606
607        // (3) dead time correction
608        // TODO:
609        // -- test for correct operation
610        // x- loop over all of the detectors
611        // x- B detector is a special case (do separately, then loop over NoB)
612        // -- this DOES alter the data
613        // -- verify the error propagation (not done yet)
614        //
615        Variable countRate
616        NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
617        if (gDoDeadTimeCor == 1)
618                Print "Doing DeadTime correction"// for "+ detStr
619                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
620                        detStr = StringFromList(ii, ksDetectorListAll, ";")
621                        Wave w = V_getDetectorDataW(fname,detStr)
622                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
623                        ctTime = V_getCount_time(fname)
624
625                        if(cmpstr(detStr,"B") == 0)
626                                Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
627                                // do the correction for the back panel
628                                countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts
629
630                                w = w/(1-countRate*b_dt)
631                                w_err = w_err/(1-countRate*b_dt)
632                                                               
633                        else
634                                // do the corrections for 8 tube panels
635                                Wave w_dt = V_getDetector_deadtime(fname,detStr)
636                                V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
637                               
638                        endif
639                endfor
640               
641        else
642                Print "Dead Time correction NOT DONE"
643        endif   
644       
645
646        // (4) solid angle correction
647        //  -- this currently calculates the correction factor AND applys it to the data
648        //  -- as a result, the data values are very large since they are divided by a very small
649        //     solid angle per pixel. But all of the count values are now on the basis of
650        //    counts/(solid angle) --- meaning that they can all be binned together for I(q)
651        //    -and- - this is taken into account for absolute scaling (this part is already done)
652        //
653        // the solid angle correction is calculated for ALL detector panels.
654        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
655        if (gDoSolidAngleCor == 1)
656                Print "Doing Solid Angle correction"// for "+ detStr
657                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
658                        detStr = StringFromList(ii, ksDetectorListAll, ";")
659                        Wave w = V_getDetectorDataW(fname,detStr)
660                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
661                        // any other dimensions to pass in?
662                        V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
663                       
664                endfor
665        else
666                Print "Solid Angle correction NOT DONE"
667        endif   
668       
669               
670        // (5) angle-dependent tube shadowing
671        // TODO:
672        // -- not sure about this correction yet...
673        //
674        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
675        if (gDoTubeShadowCor == 1)
676                Print "(stub)Tube Shadow correction"
677        else
678                Print "Tube shadowing correction NOT DONE"
679        endif   
680               
681        // (6) angle dependent transmission correction
682        // TODO:
683        // x- this still needs to be filled in
684        // -- still some debate of when/where in the corrections that this is best applied
685        //    - do it here, and it's done whether the output is 1D or 2D
686        //    - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
687        // -- verify that the calculation is correct
688        // -- verify that the error propagation (in 2D) is correct
689        //
690        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
691        if (gDoTrans == 1)
692                Print "Doing Large-angle transmission correction"// for "+ detStr
693                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
694                        detStr = StringFromList(ii, ksDetectorListAll, ";")
695                        Wave w = V_getDetectorDataW(fname,detStr)
696                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
697                       
698                        V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
699                       
700                endfor
701        else
702                Print "Sample Transmission correction NOT DONE"
703        endif   
704       
705        // (7) normalize to default monitor counts
706        // TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
707        // TODO -- but there are TWO monitors - so how to switch?
708        //  --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
709        // TODO -- what do I really need to save?
710       
711        NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
712        if (gDoMonitorNormalization == 1)
713               
714                Variable monCount,savedMonCount
715                defmon=1e8                      //default monitor counts
716                monCount = V_getMonitorCount(fname)                     // TODO -- this is read in since VCALC fakes this on output
717//              monCount = V_getBeamMonNormData(fname)          // TODO -- I think this is the *real* one to read
718                savedMonCount   = monCount
719                scale = defMon/monCount         // scale factor to MULTIPLY data by to rescale to defmon
720
721                // write to newType=fname will put new values in the destination WORK folder
722                V_writeBeamMonNormSaved_count(fname,savedMonCount)                      // save the true count
723                V_writeBeamMonNormData(fname,defMon)            // mon ct is now 10^8
724                                       
725//                      // TODO
726//                      // the low efficiency monitor, expect to use this for white beam mode
727//                              -- need code switch here to determine which monitor to use.
728//
729//                      V_getBeamMonLowData(fname)
730//                      V_getBeamMonLowSaved_count(fname)       
731                       
732                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
733                        detStr = StringFromList(ii, ksDetectorListAll, ";")
734                        Wave w = V_getDetectorDataW(fname,detStr)
735                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
736
737                        // do the calculation right here. It's a simple scaling and not worth sending to another function.     
738                        //scale the data and error to the default montor counts
739               
740//
741                        w *= scale
742                        w_err *= scale          //assumes total monitor count is so large there is essentially no error
743
744                        // TODO
745                        // -- do I want to update and save the integrated detector count?
746                        Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
747                        V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
748
749                endfor
750        else
751                Print "Monitor Normalization correction NOT DONE"
752        endif
753       
754       
755        // (not done) angle dependent efficiency correction
756        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
757
758//
759///// these lines are if files are added together, not done (yet) for VSANS
760//
761//update totals to put in the work header (at the end of the function)
762//      total_mon += realsread[0]
763//
764//      total_det += dscale*realsread[2]
765//
766//      total_trn += realsread[39]
767//      total_rtime += integersread[2]
768//      total_numruns +=1
769//     
770
771//     
772        //reset the current displaytype to "newtype"
773        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
774       
775        //return to root folder (redundant)
776        SetDataFolder root:
777       
778        Return(0)
779End
780
781
782//will "ADD" the current contents of the RAW folder to the newType work folder
783//and will ADD the RAW contents to the existing content of the newType folder
784// - used when adding multiple runs together
785//(the function Raw_to_work(type) makes a fresh workfile)
786//
787//the current display type is updated to newType (global)
788Function V_Add_raw_to_work(newType)
789        String newType
790       
791        // NEW OCT 2014
792        // this corrects for adding raw data files with different attenuation   
793        // does nothing if the attenuation of RAW and destination are the same
794        NVAR doAdjustRAW_Atten = root:Packages:NIST:gDoAdjustRAW_Atten
795        if(doAdjustRAW_Atten)
796                V_Adjust_RAW_Attenuation(newType)
797        endif
798       
799        String destPath=""
800       
801        // if the desired workfile doesn't exist, let the user know, and just make a new one
802        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
803                Print "There is no old work file to add to - a new one will be created"
804                //call Raw_to_work(), then return from this function
805                V_Raw_to_Work(newType)
806                Return(0)               //does not generate an error - a single file was converted to work.newtype
807        Endif
808       
809        NVAR pixelsX = root:myGlobals:gNPixelsX
810        NVAR pixelsY = root:myGlobals:gNPixelsY
811       
812        //now make references to data in newType folder
813        DestPath="root:Packages:NIST:"+newType 
814        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data
815        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data
816        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data
817        WAVE/T textread=$(destPath + ":textread")
818        WAVE integersread=$(destPath + ":integersread")
819        WAVE realsread=$(destPath + ":realsread")
820       
821        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
822        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
823
824
825        defmon=1e8                      //default monitor counts
826       
827        //Yes, add to previous run(s) in work, that does exist
828        //use the actual monitor count run.savmon rather than the normalized monitor count
829        //in run.moncnt and unscale the work data
830       
831        total_mon = realsread[1]        //saved monitor count
832        uscale = total_mon/defmon               //unscaling factor
833        total_det = uscale*realsread[2]         //unscaled detector count
834        total_trn = uscale*realsread[39]        //unscaled trans det count
835        total_numruns = integersread[3] //number of runs in workfile
836        total_rtime = integersread[2]           //total counting time in workfile
837        //retrieve workfile beamcenter
838        wrk_beamx = realsread[16]
839        wrk_beamy = realsread[17]
840        //unscale the workfile data in "newType"
841        //
842        //check for log-scaling and adjust if necessary
843        // should not be needed now - using display flag instead
844//      ConvertFolderToLinearScale(newType)
845        //
846        //then unscale the data array
847        data *= uscale
848        dest_data_err *= uscale
849       
850        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
851        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
852        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
853        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread"
854        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
855        WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
856       
857        //check for log-scaling of the raw data - make sure it's linear
858        // should not be needed now - using display flag instead
859//      ConvertFolderToLinearScale("RAW")
860       
861        // switches to control what is done, don't do the transmission correction for the BGD measurement
862        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCor
863        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCor
864        Variable doTrans = gDoTrans
865        if(cmpstr("BGD",newtype) == 0)
866                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
867        endif   
868       
869        V_DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans) //applies correction to raw_data, and overwrites it
870       
871        //deadtime corrections to raw data
872        // TODO - do the tube correction for dead time now
873        deadTime = 1//DetectorDeadtime(raw_text[3],raw_text[9],dateAndTimeStr=raw_text[1],dtime=raw_reals[48])          //pick the correct detector deadtime, switch on date too
874        itim = raw_ints[2]
875        cntrate = sum(raw_data,-inf,inf)/itim           //080802 use data sum, rather than scaler value
876        dscale = 1/(1-deadTime*cntrate)
877
878#if (exists("ILL_D22")==6)
879        Variable tubeSum
880        // for D22 detector might need to use cntrate/128 as it is the tube response
881        for(ii=0;ii<pixelsX;ii+=1)
882                //sum the counts in each tube
883                tubeSum = 0
884                for(jj=0;jj<pixelsY;jj+=1)
885                        tubeSum += data[jj][ii]
886                endfor
887                // countrate in tube ii
888                cntrate = tubeSum/itim
889                // deadtime scaling in tube ii
890                dscale = 1/(1-deadTime*cntrate)
891                // multiply data[ii][] by the dead time
892                raw_data[][ii] *= dscale
893                raw_data_err[][ii] *= dscale
894        endfor
895#else
896        // dead time correction on all other RAW data, including NCNR
897        raw_data *= dscale
898        raw_data_err *= dscale
899#endif
900
901        //update totals by adding RAW values to the local ones (write to work header at end of function)
902        total_mon += raw_reals[0]
903
904        total_det += dscale*raw_reals[2]
905
906        total_trn += raw_reals[39]
907        total_rtime += raw_ints[2]
908        total_numruns +=1
909       
910        //do the beamcenter shifting if there is a mismatch
911        //and then add the two data sets together, changing "data" since it is the workfile data
912        xshift = raw_reals[16] - wrk_beamx
913        yshift = raw_reals[17] - wrk_beamy
914       
915        If((xshift != 0) || (yshift != 0))
916                DoAlert 1,"Do you want to ignore the beam center mismatch?"
917                if(V_flag==1)
918                        xshift=0
919                        yshift=0
920                endif
921        endif
922       
923        If((xshift == 0) && (yshift == 0))              //no shift, just add them
924                data += raw_data                //deadtime correction has already been done to the raw data
925                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum
926        Endif
927       
928        //scale the data to the default montor counts
929        scale = defmon/total_mon
930        data *= scale
931        dest_data_err *= scale
932       
933        // keep "data" and linear_data in sync in the destination folder
934        data_copy = data
935       
936        //all is done, except for the bookkeeping of updating the header info in the work folder
937        textread[1] = date() + " " + time()             //date + time stamp
938        integersread[3] = total_numruns                                         //numruns = more than one
939        realsread[1] = total_mon                        //save the true monitor count
940        realsread[0] = defmon                                   //monitor ct = defmon
941        integersread[2] = total_rtime                   // total counting time
942        realsread[2] = scale*total_det                  //scaled detector counts
943        realsread[39] = scale*total_trn                 //scaled transmission counts
944       
945        //Add the added raw filename to the list of files in the workfile
946        String newfile = ";" + raw_text[0]
947        SVAR oldList = $(destPath + ":fileList")
948        String/G $(destPath + ":fileList") = oldList + newfile
949       
950        //reset the current display type to "newtype"
951        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
952        gCurDispType = newType
953       
954        //return to root folder (redundant)
955        SetDataFolder root:
956       
957        Return(0)
958End
959
960
961//used for adding DRK (beam shutter CLOSED) data to a workfile
962//force the monitor count to 1, since it's irrelevant
963// run data through normal "add" step, then unscale default monitor counts
964//to get the data back on a simple time basis
965//
966Function V_Raw_to_Work_NoNorm(type)
967        String type
968       
969        WAVE reals=$("root:Packages:NIST:RAW:realsread")
970        reals[1]=1              //true monitor counts, still in raw
971        V_Raw_to_work(type)
972        //data is now in "type" folder
973        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
974        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
975        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
976        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
977       
978        Variable norm_mon,tot_mon,scale
979       
980        norm_mon = new_reals[0]         //should be 1e8
981        tot_mon = new_reals[1]          //should be 1
982        scale= norm_mon/tot_mon
983       
984        data /= scale           //unscale the data
985        data_err /= scale
986       
987        // to keep "data" and linear_data in sync
988        data_copy = data
989       
990        return(0)
991End
992
993//used for adding DRK (beam shutter CLOSED) data to a workfile
994//force the monitor count to 1, since it's irrelevant
995// run data through normal "add" step, then unscale default monitor counts
996//to get the data back on a simple time basis
997//
998Function V_Add_Raw_to_Work_NoNorm(type)
999        String type
1000       
1001        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1002        reals[1]=1              //true monitor counts, still in raw
1003        V_Add_Raw_to_work(type)
1004        //data is now in "type" folder
1005        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1006        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1007        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1008        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1009       
1010        Variable norm_mon,tot_mon,scale
1011       
1012        norm_mon = new_reals[0]         //should be 1e8
1013        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
1014        scale= norm_mon/tot_mon
1015       
1016        data /= scale           //unscale the data
1017        data_err /= scale
1018       
1019        // to keep "data" and linear_data in sync
1020        data_copy = data
1021       
1022        return(0)
1023End
1024
Note: See TracBrowser for help on using the repository browser.