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

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

corrected issues with transmission calculation where the popup menu of sample files was limited in number. this limitation has been removed.

fixed the re-calculation of transmission when the same value is to be patched to multiple sample files with the same group ID. now the transmission is calculated once, for the first file in the popup and the values are simply written to the remaining files.

when the box for the open beam is defined, the panel where it is located is written to the file in a new field under /reduction. it is later recalled in the transmission panel.

Defined a "Reference" beam center position for each carriage as the RIGHT panel center. then all other panel centers (L, T, B) can be derived from this value. if the beam center is measured on the Left panel, it is converted to "right" coordinates before reporting.

File size: 35.7 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//
386// x- 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// x- be sure the data is either properly written as 2D in the file, or converted to 2D before
389//    duplicating here
390// x- ? do I recast to DP here? No- V_MakeDataWaves_DP() is called directly prior to this call, so data
391//     coming in is already DP
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 - anything else???
501        // -- ?? some of this is done at load time for RAW data. Not a problem 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 each 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                detStr = "B"
596                Wave w = V_getDetectorDataW(fname,detStr)
597                Wave cal_x = V_getDet_cal_x(fname,detStr)
598                Wave cal_y = V_getDet_cal_y(fname,detStr)
599               
600                V_NonLinearCorrection_B(fname,w,cal_x,cal_y,detStr,destPath)
601               
602                if(kBCTR_CM)
603
604                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
605                        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
606                        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
607                        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
608                        x_mm[0] = V_getDet_beam_center_x(fname,detStr) * 10             // convert cm to mm
609                        y_mm[0] = V_getDet_beam_center_y(fname,detStr) * 10             // convert cm to mm
610                       
611                        // now I need to convert the beam center in mm to pixels
612                        // and have some rational place to look for it...
613                        V_ConvertBeamCtr_to_pixB(fname,detStr,destPath)
614                else
615                        // beam center is in pixels, so use the old routine
616                        V_ConvertBeamCtr_to_mmB(fname,detStr,destPath)
617
618                endif
619                V_Detector_CalcQVals(fname,detStr,destPath)
620               
621        else
622                Print "Non-linear correction NOT DONE"
623        endif
624
625
626        // (3) dead time correction
627        // DONE:
628        // x- test for correct operation
629        // x- loop over all of the detectors
630        // x- B detector is a special case (do separately, then loop over NoB)
631        // x- this DOES alter the data
632        // x- verify the error propagation (not done yet)
633        //
634        Variable countRate
635        NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
636        if (gDoDeadTimeCor == 1)
637                Print "Doing DeadTime correction"// for "+ detStr
638                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
639                        detStr = StringFromList(ii, ksDetectorListAll, ";")
640                        Wave w = V_getDetectorDataW(fname,detStr)
641                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
642                        ctTime = V_getCount_time(fname)
643
644                        if(cmpstr(detStr,"B") == 0)
645                                Variable b_dt = V_getDetector_deadtime_B(fname,detStr)
646                                // do the correction for the back panel
647                                countRate = sum(w,-inf,inf)/ctTime              //use sum of detector counts
648
649                                w = w/(1-countRate*b_dt)
650                                w_err = w_err/(1-countRate*b_dt)
651                                                               
652                        else
653                                // do the corrections for 8 tube panels
654                                Wave w_dt = V_getDetector_deadtime(fname,detStr)
655                                V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
656                               
657                        endif
658                endfor
659               
660        else
661                Print "Dead Time correction NOT DONE"
662        endif   
663       
664
665        // (4) solid angle correction
666        //  -- this currently calculates the correction factor AND applys it to the data
667        //  -- as a result, the data values are very large since they are divided by a very small
668        //     solid angle per pixel. But all of the count values are now on the basis of
669        //    counts/(solid angle) --- meaning that they can all be binned together for I(q)
670        //    -and- - this is taken into account for absolute scaling (this part is already done)
671        //
672        // the solid angle correction is calculated for ALL detector panels.
673        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
674        if (gDoSolidAngleCor == 1)
675                Print "Doing Solid Angle correction"// for "+ detStr
676                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
677                        detStr = StringFromList(ii, ksDetectorListAll, ";")
678                        Wave w = V_getDetectorDataW(fname,detStr)
679                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
680                        // any other dimensions to pass in?
681                        V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
682                       
683                endfor
684        else
685                Print "Solid Angle correction NOT DONE"
686        endif   
687       
688               
689        // (5) angle-dependent tube shadowing
690        // TODO:
691        // -- not sure about this correction yet...
692        //
693        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
694        if (gDoTubeShadowCor == 1)
695                Print "(stub)Tube Shadow correction"
696        else
697                Print "Tube shadowing correction NOT DONE"
698        endif   
699               
700        // (6) angle dependent transmission correction
701        // TODO:
702        // x- this still needs to be filled in
703        // -- still some debate of when/where in the corrections that this is best applied
704        //    - do it here, and it's done whether the output is 1D or 2D
705        //    - do it later (where SAMPLE information is used) since this section is ONLY instrument-specific
706        // -- verify that the calculation is correct
707        // -- verify that the error propagation (in 2D) is correct
708        //
709        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
710        if (gDoTrans == 1)
711                Print "Doing Large-angle transmission correction"// for "+ detStr
712                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
713                        detStr = StringFromList(ii, ksDetectorListAll, ";")
714                        Wave w = V_getDetectorDataW(fname,detStr)
715                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
716                       
717                        V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
718                       
719                endfor
720        else
721                Print "Sample Transmission correction NOT DONE"
722        endif   
723       
724        // (7) normalize to default monitor counts
725        // TODO(DONE) x- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
726        // TODO -- but there are TWO monitors - so how to switch?
727        //  --- AND, there is also /entry/control/monitor_counts !!! Which one is the correct value? Which will NICE write
728        // TODO -- what do I really need to save?
729       
730        NVAR gDoMonitorNormalization = root:Packages:NIST:VSANS:Globals:gDoMonitorNormalization
731        if (gDoMonitorNormalization == 1)
732               
733                Variable monCount,savedMonCount
734                defmon=1e8                      //default monitor counts
735                monCount = V_getMonitorCount(fname)                     // TODO -- this is read in since VCALC fakes this on output
736//              monCount = V_getBeamMonNormData(fname)          // TODO -- I think this is the *real* one to read
737                savedMonCount   = monCount
738                scale = defMon/monCount         // scale factor to MULTIPLY data by to rescale to defmon
739
740                // write to newType=fname will put new values in the destination WORK folder
741                V_writeBeamMonNormSaved_count(fname,savedMonCount)                      // save the true count
742                V_writeBeamMonNormData(fname,defMon)            // mon ct is now 10^8
743                                       
744//                      // TODO
745//                      // the low efficiency monitor, expect to use this for white beam mode
746//                              -- need code switch here to determine which monitor to use.
747//
748//                      V_getBeamMonLowData(fname)
749//                      V_getBeamMonLowSaved_count(fname)       
750                       
751                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
752                        detStr = StringFromList(ii, ksDetectorListAll, ";")
753                        Wave w = V_getDetectorDataW(fname,detStr)
754                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
755
756                        // do the calculation right here. It's a simple scaling and not worth sending to another function.     
757                        //scale the data and error to the default monitor counts
758               
759//
760                        w *= scale
761                        w_err *= scale          //assumes total monitor count is so large there is essentially no error
762
763                        // TODO
764                        // -- do I want to update and save the integrated detector count?
765                        Variable integratedCount = V_getDet_IntegratedCount(fname,detStr)
766                        V_writeDet_IntegratedCount(fname,detStr,integratedCount*scale)
767
768                endfor
769        else
770                Print "Monitor Normalization correction NOT DONE"
771        endif
772       
773       
774        // (not done) angle dependent efficiency correction
775        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
776
777//
778///// these lines are if files are added together, not done (yet) for VSANS
779//
780//update totals to put in the work header (at the end of the function)
781//      total_mon += realsread[0]
782//
783//      total_det += dscale*realsread[2]
784//
785//      total_trn += realsread[39]
786//      total_rtime += integersread[2]
787//      total_numruns +=1
788//     
789
790//     
791        //reset the current displaytype to "newtype"
792        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
793       
794        //return to root folder (redundant)
795        SetDataFolder root:
796       
797        Return(0)
798End
799
800
801//will "ADD" the current contents of the RAW folder to the newType work folder
802//and will ADD the RAW contents to the existing content of the newType folder
803// - used when adding multiple runs together
804//(the function Raw_to_work(type) makes a fresh workfile)
805//
806//the current display type is updated to newType (global)
807Function V_Add_raw_to_work(newType)
808        String newType
809       
810        // NEW OCT 2014
811        // this corrects for adding raw data files with different attenuation   
812        // does nothing if the attenuation of RAW and destination are the same
813        NVAR doAdjustRAW_Atten = root:Packages:NIST:gDoAdjustRAW_Atten
814        if(doAdjustRAW_Atten)
815                V_Adjust_RAW_Attenuation(newType)
816        endif
817       
818        String destPath=""
819       
820        // if the desired workfile doesn't exist, let the user know, and just make a new one
821        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
822                Print "There is no old work file to add to - a new one will be created"
823                //call Raw_to_work(), then return from this function
824                V_Raw_to_Work(newType)
825                Return(0)               //does not generate an error - a single file was converted to work.newtype
826        Endif
827       
828        NVAR pixelsX = root:myGlobals:gNPixelsX
829        NVAR pixelsY = root:myGlobals:gNPixelsY
830       
831        //now make references to data in newType folder
832        DestPath="root:Packages:NIST:"+newType 
833        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data
834        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data
835        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data
836        WAVE/T textread=$(destPath + ":textread")
837        WAVE integersread=$(destPath + ":integersread")
838        WAVE realsread=$(destPath + ":realsread")
839       
840        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
841        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
842
843
844        defmon=1e8                      //default monitor counts
845       
846        //Yes, add to previous run(s) in work, that does exist
847        //use the actual monitor count run.savmon rather than the normalized monitor count
848        //in run.moncnt and unscale the work data
849       
850        total_mon = realsread[1]        //saved monitor count
851        uscale = total_mon/defmon               //unscaling factor
852        total_det = uscale*realsread[2]         //unscaled detector count
853        total_trn = uscale*realsread[39]        //unscaled trans det count
854        total_numruns = integersread[3] //number of runs in workfile
855        total_rtime = integersread[2]           //total counting time in workfile
856        //retrieve workfile beamcenter
857        wrk_beamx = realsread[16]
858        wrk_beamy = realsread[17]
859        //unscale the workfile data in "newType"
860        //
861        //check for log-scaling and adjust if necessary
862        // should not be needed now - using display flag instead
863//      ConvertFolderToLinearScale(newType)
864        //
865        //then unscale the data array
866        data *= uscale
867        dest_data_err *= uscale
868       
869        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
870        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
871        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
872        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread"
873        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
874        WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
875       
876        //check for log-scaling of the raw data - make sure it's linear
877        // should not be needed now - using display flag instead
878//      ConvertFolderToLinearScale("RAW")
879       
880        // switches to control what is done, don't do the transmission correction for the BGD measurement
881        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCor
882        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCor
883        Variable doTrans = gDoTrans
884        if(cmpstr("BGD",newtype) == 0)
885                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
886        endif   
887       
888        V_DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans) //applies correction to raw_data, and overwrites it
889       
890        //deadtime corrections to raw data
891        // TODO - do the tube correction for dead time now
892        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
893        itim = raw_ints[2]
894        cntrate = sum(raw_data,-inf,inf)/itim           //080802 use data sum, rather than scaler value
895        dscale = 1/(1-deadTime*cntrate)
896
897#if (exists("ILL_D22")==6)
898        Variable tubeSum
899        // for D22 detector might need to use cntrate/128 as it is the tube response
900        for(ii=0;ii<pixelsX;ii+=1)
901                //sum the counts in each tube
902                tubeSum = 0
903                for(jj=0;jj<pixelsY;jj+=1)
904                        tubeSum += data[jj][ii]
905                endfor
906                // countrate in tube ii
907                cntrate = tubeSum/itim
908                // deadtime scaling in tube ii
909                dscale = 1/(1-deadTime*cntrate)
910                // multiply data[ii][] by the dead time
911                raw_data[][ii] *= dscale
912                raw_data_err[][ii] *= dscale
913        endfor
914#else
915        // dead time correction on all other RAW data, including NCNR
916        raw_data *= dscale
917        raw_data_err *= dscale
918#endif
919
920        //update totals by adding RAW values to the local ones (write to work header at end of function)
921        total_mon += raw_reals[0]
922
923        total_det += dscale*raw_reals[2]
924
925        total_trn += raw_reals[39]
926        total_rtime += raw_ints[2]
927        total_numruns +=1
928       
929        //do the beamcenter shifting if there is a mismatch
930        //and then add the two data sets together, changing "data" since it is the workfile data
931        xshift = raw_reals[16] - wrk_beamx
932        yshift = raw_reals[17] - wrk_beamy
933       
934        If((xshift != 0) || (yshift != 0))
935                DoAlert 1,"Do you want to ignore the beam center mismatch?"
936                if(V_flag==1)
937                        xshift=0
938                        yshift=0
939                endif
940        endif
941       
942        If((xshift == 0) && (yshift == 0))              //no shift, just add them
943                data += raw_data                //deadtime correction has already been done to the raw data
944                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum
945        Endif
946       
947        //scale the data to the default montor counts
948        scale = defmon/total_mon
949        data *= scale
950        dest_data_err *= scale
951       
952        // keep "data" and linear_data in sync in the destination folder
953        data_copy = data
954       
955        //all is done, except for the bookkeeping of updating the header info in the work folder
956        textread[1] = date() + " " + time()             //date + time stamp
957        integersread[3] = total_numruns                                         //numruns = more than one
958        realsread[1] = total_mon                        //save the true monitor count
959        realsread[0] = defmon                                   //monitor ct = defmon
960        integersread[2] = total_rtime                   // total counting time
961        realsread[2] = scale*total_det                  //scaled detector counts
962        realsread[39] = scale*total_trn                 //scaled transmission counts
963       
964        //Add the added raw filename to the list of files in the workfile
965        String newfile = ";" + raw_text[0]
966        SVAR oldList = $(destPath + ":fileList")
967        String/G $(destPath + ":fileList") = oldList + newfile
968       
969        //reset the current display type to "newtype"
970        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
971        gCurDispType = newType
972       
973        //return to root folder (redundant)
974        SetDataFolder root:
975       
976        Return(0)
977End
978
979
980//used for adding DRK (beam shutter CLOSED) data to a workfile
981//force the monitor count to 1, since it's irrelevant
982// run data through normal "add" step, then unscale default monitor counts
983//to get the data back on a simple time basis
984//
985Function V_Raw_to_Work_NoNorm(type)
986        String type
987       
988        WAVE reals=$("root:Packages:NIST:RAW:realsread")
989        reals[1]=1              //true monitor counts, still in raw
990        V_Raw_to_work(type)
991        //data is now in "type" folder
992        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
993        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
994        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
995        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
996       
997        Variable norm_mon,tot_mon,scale
998       
999        norm_mon = new_reals[0]         //should be 1e8
1000        tot_mon = new_reals[1]          //should be 1
1001        scale= norm_mon/tot_mon
1002       
1003        data /= scale           //unscale the data
1004        data_err /= scale
1005       
1006        // to keep "data" and linear_data in sync
1007        data_copy = data
1008       
1009        return(0)
1010End
1011
1012//used for adding DRK (beam shutter CLOSED) data to a workfile
1013//force the monitor count to 1, since it's irrelevant
1014// run data through normal "add" step, then unscale default monitor counts
1015//to get the data back on a simple time basis
1016//
1017Function V_Add_Raw_to_Work_NoNorm(type)
1018        String type
1019       
1020        WAVE reals=$("root:Packages:NIST:RAW:realsread")
1021        reals[1]=1              //true monitor counts, still in raw
1022        V_Add_Raw_to_work(type)
1023        //data is now in "type" folder
1024        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
1025        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
1026        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
1027        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
1028       
1029        Variable norm_mon,tot_mon,scale
1030       
1031        norm_mon = new_reals[0]         //should be 1e8
1032        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
1033        scale= norm_mon/tot_mon
1034       
1035        data /= scale           //unscale the data
1036        data_err /= scale
1037       
1038        // to keep "data" and linear_data in sync
1039        data_copy = data
1040       
1041        return(0)
1042End
1043
Note: See TracBrowser for help on using the repository browser.