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

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

Many changes:

Made the VCALC panel aware of all of the binning options
Corrected the behavior of the VCALC preset conditions
Adjusted how the Slit data is binned so that there are not duplicated q-values in the output

Made Absolute scaling aware of the back detector. Now the ABS String in the protocol has a second
set of scaling constants tagged with "_B" for the back detector. There is an added button
on the protocol panel to set the second set of constants. For the back detector, the read noise
is subtracted by reading it from the empty beam file (shifting over to the right by one box width)
All of the associated abs procedures are now aware of this.
More error checking needs to be added.

Back detector image is now shifted upon loading of the data. the default mask takes this into account
and masks out the padded (zero) regions.

in the protocol, DIV and MSK do not use grep any longer. it was just way too slow. Now it depends on

the file name having DIV or MASK respectively.



Raw data files can now be added together, in the usual way from the protocol panel.



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