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

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

Added the angle dependent transmission correction to the data correction in the raw_to_work step, in 2D

added a testing file that can generate fake event data, read, write, and decode it. Read is based on GBLoadWave. Hoepfully I'll not need to write an XOP. manipulation of the 64 bit words are done with simple bit shifts and logic.

also added are a number of error checking routines to improve behavior when wave, folders, etc. are missing.

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