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

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

LOTS of changes to accommodate the beam center being reported in cm rather than pixels. Required a lot of changes to VCALC (to fill in simulated data), and in the reading and interpreting of data for display, and most importantly, the calculation of q.

There may still be a few residual bugs with this. I am still re-testing with new sample data sets.

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