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

Last change on this file since 982 was 982, checked in by srkline, 7 years ago

more changes and additons to display VSANS data

adding functions for IvsQ plotting

coverted much of VCALC to have similar folder structure as HDF to allow re-use of the Q-binning procedures from VCALC with real data in work files.

re-working the beam center finder to get it to work with work file data rather then only VCALC.

new plotting routines for the panels to rescale to the beam center (still in pixels, though)

File size: 28.8 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//
14// - copying workfiles to another folder
15//
16// - absolute scaling
17//
18// - (no) the WorkFile Math panel for simple image math (not done - maybe in the future?)
19// -
20// - (no) adding work.drk data without normalizing to monitor counts (the case not currently handled)
21//***************************
22
23//
24// Functions used for manipulation of the local Igor "WORK" folder
25// structure as raw data is displayed and processed.
26//
27//
28Strconstant ksDetectorListNoB = "FT;FB;FL;FR;MT;MB;ML;MR;"
29Strconstant ksDetectorListAll = "FT;FB;FL;FR;MT;MB;ML;MR;B;"
30
31
32//
33//Entry procedure from main panel
34//
35Proc 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        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//                      -- when do I make the 2D error waves?
50//
51// TODO - decide what exactly I need to copy over. May be best to copy all, and delete
52//       what I know that I don't need
53//
54// TODO !!! DuplicateDataFolder will FAIL - in the base case of RAW data files, the
55//  data is actually in use - so it will fail every time. need an alternate solution. in SANS,
56// there are a limited number of waves to carry over, so Dupliate/O is used for rw, tw, data, etc.
57//
58// TODO : I also need a list of what is generated during processing that may be hanging around - that I need to
59//     be sure to get rid of - like the calibration waves, solidAngle, etc.
60//
61// hdfDF is the name only of the data in storage. May be full file name with extension (clean as needed)
62// type is the destination WORK folder for the copy
63//
64Function CopyHDFToWorkFolder(fromStr,toStr)
65        String fromStr,toStr
66       
67        String fromDF, toDF
68       
69        // make the DF paths - source and destination
70        fromDF = "root:Packages:NIST:VSANS:"+fromStr
71        toDF = "root:Packages:NIST:VSANS:"+toStr
72       
73//      // make a copy of the file name for my own use, since it's not in the file
74//      String/G $(toDF+":file_name") = root:
75       
76        // copy the folders
77        KillDataFolder/Z $toDF                  //DuplicateDataFolder will not overwrite, so Kill
78       
79        if(V_flag == 0)         // kill DF was OK
80                DuplicateDataFolder $("root:Packages:NIST:VSANS:"+fromStr),$("root:Packages:NIST:VSANS:"+toStr)
81               
82                // I can delete these if they came along with RAW
83                //   DAS_logs
84                //   top-level copies of data (duplicate links)
85                KillDataFolder/Z $(toDF+":entry:entry:DAS_logs")
86                KillDataFolder/Z $(toDF+":entry:entry:data")
87                KillDataFolder/Z $(toDF+":entry:entry:data_B")
88                KillDataFolder/Z $(toDF+":entry:entry:data_ML")
89                KillDataFolder/Z $(toDF+":entry:entry:data_MR")
90                KillDataFolder/Z $(toDF+":entry:entry:data_MT")
91                KillDataFolder/Z $(toDF+":entry:entry:data_MB")
92                KillDataFolder/Z $(toDF+":entry:entry:data_FL")
93                KillDataFolder/Z $(toDF+":entry:entry:data_FR")
94                KillDataFolder/Z $(toDF+":entry:entry:data_FT")
95                KillDataFolder/Z $(toDF+":entry:entry:data_FB")
96
97                return(0)
98        else
99                // need to do this the hard way, duplicate/O recursively
100                // see V_CopyToWorkFolder()
101               
102                // everything on the top level
103                V_DuplicateDataFolder($(fromDF+":entry:entry"),fromStr,toStr,0,"",0)    //no recursion here
104                // control
105                V_DuplicateDataFolder($(fromDF+":entry:entry:control"),fromStr,toStr,0,"",1)    //yes recursion here
106                // instrument
107                V_DuplicateDataFolder($(fromDF+":entry:entry:instrument"),fromStr,toStr,0,"",1) //yes recursion here
108                // reduction
109                V_DuplicateDataFolder($(fromDF+":entry:entry:reduction"),fromStr,toStr,0,"",1)  //yes recursion here
110                // sample
111                V_DuplicateDataFolder($(fromDF+":entry:entry:sample"),fromStr,toStr,0,"",1)     //yes recursion here
112
113        endif   
114       
115        return(0)
116end
117
118
119////////
120// see the help entry for IndexedDir for help on (possibly) how to do this faster
121// -- see the function Function ScanDirectories(pathName, printDirNames)
122//
123
124
125// from IgorExchange On July 17th, 2011 jtigor
126// started from "Recursively List Data Folder Contents"
127// Posted July 15th, 2011 by hrodstein
128//
129//
130//
131Proc V_CopyWorkFolderTest(dataFolderStr, fromStr, toStr, level, sNBName, recurse)
132        String dataFolderStr="root:Packages:NIST:VSANS:RAW"
133        String fromStr = "RAW"
134        String toStr="SAM"
135        Variable level=0
136        String sNBName="DataFolderTree"
137        Variable recurse = 1
138       
139        V_DuplicateDataFolder($dataFolderStr, fromStr, toStr, level, sNBName, recurse)
140
141
142end
143
144// ListDataFolder(dfr, level)
145// Recursively lists objects in data folder.
146// Pass data folder path for dfr and 0 for level.
147// pass level == 0 for the first call
148//  sNBName = "" prints nothing. any name will generate a notebook
149//
150// recurse == 0 will do only the specified folder, anything else will recurse all levels
151// toStr is the string name of the top-level folder only, not the full path
152//
153//
154Function V_DuplicateDataFolder(dfr, fromStr, toStr, level, sNBName,recurse)
155        DFREF dfr
156        String fromStr
157        String toStr
158        Variable level                  // Pass 0 to start
159        String sNBName
160        Variable recurse
161 
162        String name
163        String dfName
164        String sString
165       
166        String toDF = ""
167 
168        if (level == 0)         // this is the data folder, generate if needed in the destination
169                name = GetDataFolder(1, dfr)
170//              sPrintf sString, "%s (data folder)\r", name
171                toDF = ReplaceString(fromStr,name,toStr,1)              // case-sensitive replace
172                sprintf sString, "NewDataFolder/O %s\r",toDF
173                NewDataFolder/O $(RemoveEnding(toDF,":"))                       // remove trailing semicolon if it's there
174               
175                V_WriteBrowserInfo(sString, 1, sNBName)
176        endif
177 
178        dfName = GetDataFolder(1, dfr)
179        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
180        Variable i
181 
182        String indentStr = "\t"
183        for(i=0; i<level; i+=1)
184                indentStr += "\t"
185        endfor
186 
187        Variable numWaves = CountObjectsDFR(dfr, 1)
188        for(i=0; i<numWaves; i+=1)
189                name = GetIndexedObjNameDFR(dfr, 1, i)
190                //
191                // wave type does not matter now. Duplicate does not care
192                //
193                sPrintf sString, "Duplicate/O  %s,  %s\r",dfName+name,toDF+name
194                Duplicate/O $(dfName+name),$(toDF+name)
195               
196                V_WriteBrowserInfo(sString, 2, sNBName)
197        endfor 
198 
199        Variable numNumericVariables = CountObjectsDFR(dfr, 2) 
200        for(i=0; i<numNumericVariables; i+=1)
201                name = GetIndexedObjNameDFR(dfr, 2, i)
202                sPrintf sString, "%s%s (numeric variable)\r", indentStr, name
203                V_WriteBrowserInfo(sString, 3, sNBName)
204        endfor 
205 
206        Variable numStringVariables = CountObjectsDFR(dfr, 3)   
207        for(i=0; i<numStringVariables; i+=1)
208                name = GetIndexedObjNameDFR(dfr, 3, i)
209                sPrintf sString, "%s%s (string variable)\r", indentStr, name
210                V_WriteBrowserInfo(sString, 4, sNBName)
211        endfor 
212
213        if(recurse)
214                Variable numDataFolders = CountObjectsDFR(dfr, 4)       
215                for(i=0; i<numDataFolders; i+=1)
216                        name = GetIndexedObjNameDFR(dfr, 4, i)
217//                      sPrintf sString, "%s%s (data folder)\r", indentStr, name
218                         dfName = GetDataFolder(1, dfr)
219                         
220                        toDF = ReplaceString(fromStr,dfName,toStr,1)            // case-sensitive replace
221                        sprintf sString, "NewDataFolder/O %s\r",toDF+name
222                        NewDataFolder/O $(toDF+name)
223                       
224                       
225                        V_WriteBrowserInfo(sString, 1, sNBName)
226                        DFREF childDFR = dfr:$(name)
227                        V_DuplicateDataFolder(childDFR, fromStr, toStr, level+1, sNBName, recurse)
228                endfor 
229        endif
230         
231//when finished walking tree, save as RTF with dialog   
232//      if(level == 0 && strlen(sNBName) != 0)
233//              SaveNotebook /I /S=4  $sNBName
234//      endif
235End
236 
237Function V_WriteBrowserInfo(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 duplicate the error for WM...
304Macro 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 = 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 = 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        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 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        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
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                        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        // -- currently, the "B" detector is skipped
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        // -- still not sure whether to duplicate these calculations as the RAW data is loaded. It would allow the RAW
420        //    data to be properly displayed, but without all of the (complete) set of 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 -other corrections modify the data
423        NVAR gDoNonLinearCor = root:Packages:NIST:VSANS:Globals:gDoNonLinearCor
424        // generate a distance matrix for each of the detectors
425        if (gDoNonLinearCor == 1)
426                Print "Doing Non-linear correction"// for "+ detStr
427                for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
428                        detStr = StringFromList(ii, ksDetectorListNoB, ";")
429                        Wave w = V_getDetectorDataW(fname,detStr)
430//                      Wave w_err = V_getDetectorDataErrW(fname,detStr)
431                        Wave w_calib = V_getDetTube_spatialCalib(fname,detStr)
432                        Variable tube_width = V_getDet_tubeWidth(fname,detStr)
433                        NonLinearCorrection(w,w_calib,tube_width,detStr,destPath)
434                       
435                        // (2.5) Calculate the q-values
436                        // calculating q-values can't be done unless the non-linear corrections are calculated
437                        // so go ahead and put it in this loop.
438                        // TODO :
439                        // -- make sure that everything is present before the calculation
440                        // -- beam center must be properly defined in terms of real distance
441                        // -- distances/zero location/ etc. must be clearly documented for each detector
442                        //      ** this assumes that NonLinearCorrection() has been run to generate data_RealDistX and Y
443                        // ** this routine Makes the waves QTot, qx, qy, qz in each detector folder.
444                        //
445                        V_Detector_CalcQVals(fname,detStr,destPath)
446                       
447                endfor
448        else
449                Print "Non-linear correction not done"
450        endif
451
452        // (3) solid angle correction
453        NVAR gDoSolidAngleCor = root:Packages:NIST:VSANS:Globals:gDoSolidAngleCor
454        if (gDoSolidAngleCor == 1)
455                Print "Doing Solid Angle correction"// for "+ detStr
456                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
457                        detStr = StringFromList(ii, ksDetectorListAll, ";")
458                        Wave w = V_getDetectorDataW(fname,detStr)
459                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
460                        Wave w_dt = V_getDetector_deadtime(fname,detStr)
461//                      SolidAngleCorrection(fill this in)
462                       
463                endfor
464        else
465                Print "Solid Angle correction not done"
466        endif   
467       
468        // (4) dead time correction
469        // TODO: -- remove the hard-wired test
470        // -- test for correct operation
471        // -- loop over all of the detectors
472        // -- B detector is a special case (do separately, then loop over NoB)
473        NVAR gDoDeadTimeCor = root:Packages:NIST:VSANS:Globals:gDoDeadTimeCor
474        ctTime = V_getCount_time(fname)
475        if (gDoDeadTimeCor == 1)
476                Print "Doing DeadTime correction"// for "+ detStr
477                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
478                        detStr = StringFromList(ii, ksDetectorListAll, ";")
479                        Wave w = V_getDetectorDataW(fname,detStr)
480                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
481                        Wave w_dt = V_getDetector_deadtime(fname,detStr)
482//                      DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
483                                //deadtime corrections
484//      itim = integersread[2]
485//      cntrate = sum(data,-inf,inf)/itim               //use sum of detector counts rather than scaler value
486//      //TODO - do correct dead time correction for tubes
487//      deadtime = 1//DetectorDeadtime(textread[3],textread[9],dateAndTimeStr=textRead[1],dtime=realsRead[48])  //pick the correct deadtime
488//      dscale = 1/(1-deadTime*cntrate)
489//     
490       
491// dead time correction
492//      data *= dscale          //deadtime correction for everyone else, including NCNR
493//      data_err *= dscale
494
495                endfor
496        else
497                Print "Dead Time correction not done"
498        endif   
499       
500        // (5) angle-dependent tube shadowing
501        NVAR gDoTubeShadowCor = root:Packages:NIST:VSANS:Globals:gDoTubeShadowCor
502        if (gDoTubeShadowCor == 1)
503       
504        else
505                Print "Tube shadowing correction not done"
506        endif   
507               
508        // (6) angle dependent transmission correction
509        NVAR gDoTrans = root:Packages:NIST:VSANS:Globals:gDoTransmissionCor
510        if (gDoTrans == 1)
511                Print "Doing Large-angle transmission correction"// for "+ detStr
512                for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
513                        detStr = StringFromList(ii, ksDetectorListAll, ";")
514                        Wave w = V_getDetectorDataW(fname,detStr)
515                        Wave w_err = V_getDetectorDataErrW(fname,detStr)
516                        Wave w_dt = V_getDetector_deadtime(fname,detStr)
517//                      TransmissionCorrection(fill this in)
518                       
519                endfor
520        else
521                Print "Sample Transmission correction not done"
522        endif   
523       
524        // (7) normalize to default monitor counts
525        // TODO -- each detector is rescaled separately, but the rescaling factor is global (only one monitor!)
526        // TODO -- but there are TWO monitors - so how to switch?
527        // TODO -- what do I really need to save?
528        Print "Doing monitor normalization"// for "+ detStr
529
530        defmon=1e8                      //default monitor counts
531        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
532                detStr = StringFromList(ii, ksDetectorListAll, ";")
533                Wave w = V_getDetectorDataW(fname,detStr)
534                Wave w_err = V_getDetectorDataErrW(fname,detStr)
535                Variable monCt = V_getBeamMonNormData(fname)
536//                      MonitorNormalization(fill this in)
537        //scale the data to the default montor counts
538       
539        // TODO -- un-comment these three lines once monitor counts are reasonable - currently monCt = 9!!!
540//              scale = defmon/monCt
541//              w *= scale
542//              w_err *= scale          //assumes total monitor count is so large there is essentially no error
543
544// TODO
545// -- to write back to the local value, get the wave reference rather than the value, then I can
546//    re-assign the value directly, rather than this method (which is not terrible)     
547                // V_getBeamMonNormSaved_count()
548                // save the true monitor counts? save the scaling factor?
549                String path = "entry:instrument:beam_monitor_norm:saved_count"
550                Wave/Z savW = $("root:Packages:NIST:VSANS:"+fname+":entry:"+path)
551                savW[0] = scale
552        endfor
553       
554       
555        // (not done) angle dependent efficiency correction
556        NVAR doEfficiency = root:Packages:NIST:VSANS:Globals:gDoDetectorEffCor
557
558       
559// this function, in the past did the non-linear, solid angle, transmission, and efficiency corrections all at once
560//      DetCorr(data,data_err,realsread,doEfficiency,doTrans)           //the parameters are waves, and will be changed by the function
561
562
563       
564        //update totals to put in the work header (at the end of the function)
565//      total_mon += realsread[0]
566//
567//      total_det += dscale*realsread[2]
568//
569//      total_trn += realsread[39]
570//      total_rtime += integersread[2]
571//      total_numruns +=1
572//     
573
574        //all is done, except for the bookkeeping, updating the header information in the work folder
575
576//      integersread[3] = total_numruns                                         //numruns = 1
577//      realsread[1] = total_mon                        //save the true monitor count
578//      realsread[0] = defmon                                   //monitor ct = defmon
579//      realsread[2] = scale*total_det                  //scaled detector counts
580//     
581        //reset the current displaytype to "newtype"
582        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=newType
583       
584        //return to root folder (redundant)
585        SetDataFolder root:
586       
587        Return(0)
588End
589
590
591//will "ADD" the current contents of the RAW folder to the newType work folder
592//and will ADD the RAW contents to the existing content of the newType folder
593// - used when adding multiple runs together
594//(the function Raw_to_work(type) makes a fresh workfile)
595//
596//the current display type is updated to newType (global)
597Function Add_raw_to_work(newType)
598        String newType
599       
600        // NEW OCT 2014
601        // this corrects for adding raw data files with different attenuation   
602        // does nothing if the attenuation of RAW and destination are the same
603        NVAR doAdjustRAW_Atten = root:Packages:NIST:gDoAdjustRAW_Atten
604        if(doAdjustRAW_Atten)
605                Adjust_RAW_Attenuation(newType)
606        endif
607       
608        String destPath=""
609       
610        // if the desired workfile doesn't exist, let the user know, and just make a new one
611        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
612                Print "There is no old work file to add to - a new one will be created"
613                //call Raw_to_work(), then return from this function
614                Raw_to_Work(newType)
615                Return(0)               //does not generate an error - a single file was converted to work.newtype
616        Endif
617       
618        NVAR pixelsX = root:myGlobals:gNPixelsX
619        NVAR pixelsY = root:myGlobals:gNPixelsY
620       
621        //now make references to data in newType folder
622        DestPath="root:Packages:NIST:"+newType 
623        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data
624        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data
625        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data
626        WAVE/T textread=$(destPath + ":textread")
627        WAVE integersread=$(destPath + ":integersread")
628        WAVE realsread=$(destPath + ":realsread")
629       
630        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
631        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
632
633
634        defmon=1e8                      //default monitor counts
635       
636        //Yes, add to previous run(s) in work, that does exist
637        //use the actual monitor count run.savmon rather than the normalized monitor count
638        //in run.moncnt and unscale the work data
639       
640        total_mon = realsread[1]        //saved monitor count
641        uscale = total_mon/defmon               //unscaling factor
642        total_det = uscale*realsread[2]         //unscaled detector count
643        total_trn = uscale*realsread[39]        //unscaled trans det count
644        total_numruns = integersread[3] //number of runs in workfile
645        total_rtime = integersread[2]           //total counting time in workfile
646        //retrieve workfile beamcenter
647        wrk_beamx = realsread[16]
648        wrk_beamy = realsread[17]
649        //unscale the workfile data in "newType"
650        //
651        //check for log-scaling and adjust if necessary
652        // should not be needed now - using display flag instead
653//      ConvertFolderToLinearScale(newType)
654        //
655        //then unscale the data array
656        data *= uscale
657        dest_data_err *= uscale
658       
659        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
660        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
661        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
662        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread"
663        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
664        WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
665       
666        //check for log-scaling of the raw data - make sure it's linear
667        // should not be needed now - using display flag instead
668//      ConvertFolderToLinearScale("RAW")
669       
670        // switches to control what is done, don't do the transmission correction for the BGD measurement
671        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCorr
672        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCorr
673        Variable doTrans = gDoTrans
674        if(cmpstr("BGD",newtype) == 0)
675                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
676        endif   
677       
678        DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans)   //applies correction to raw_data, and overwrites it
679       
680        //deadtime corrections to raw data
681        // TODO - do the tube correction for dead time now
682        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
683        itim = raw_ints[2]
684        cntrate = sum(raw_data,-inf,inf)/itim           //080802 use data sum, rather than scaler value
685        dscale = 1/(1-deadTime*cntrate)
686
687#if (exists("ILL_D22")==6)
688        Variable tubeSum
689        // for D22 detector might need to use cntrate/128 as it is the tube response
690        for(ii=0;ii<pixelsX;ii+=1)
691                //sum the counts in each tube
692                tubeSum = 0
693                for(jj=0;jj<pixelsY;jj+=1)
694                        tubeSum += data[jj][ii]
695                endfor
696                // countrate in tube ii
697                cntrate = tubeSum/itim
698                // deadtime scaling in tube ii
699                dscale = 1/(1-deadTime*cntrate)
700                // multiply data[ii][] by the dead time
701                raw_data[][ii] *= dscale
702                raw_data_err[][ii] *= dscale
703        endfor
704#else
705        // dead time correction on all other RAW data, including NCNR
706        raw_data *= dscale
707        raw_data_err *= dscale
708#endif
709
710        //update totals by adding RAW values to the local ones (write to work header at end of function)
711        total_mon += raw_reals[0]
712
713        total_det += dscale*raw_reals[2]
714
715        total_trn += raw_reals[39]
716        total_rtime += raw_ints[2]
717        total_numruns +=1
718       
719        //do the beamcenter shifting if there is a mismatch
720        //and then add the two data sets together, changing "data" since it is the workfile data
721        xshift = raw_reals[16] - wrk_beamx
722        yshift = raw_reals[17] - wrk_beamy
723       
724        If((xshift != 0) || (yshift != 0))
725                DoAlert 1,"Do you want to ignore the beam center mismatch?"
726                if(V_flag==1)
727                        xshift=0
728                        yshift=0
729                endif
730        endif
731       
732        If((xshift == 0) && (yshift == 0))              //no shift, just add them
733                data += raw_data                //deadtime correction has already been done to the raw data
734                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum
735        Endif
736       
737        //scale the data to the default montor counts
738        scale = defmon/total_mon
739        data *= scale
740        dest_data_err *= scale
741       
742        // keep "data" and linear_data in sync in the destination folder
743        data_copy = data
744       
745        //all is done, except for the bookkeeping of updating the header info in the work folder
746        textread[1] = date() + " " + time()             //date + time stamp
747        integersread[3] = total_numruns                                         //numruns = more than one
748        realsread[1] = total_mon                        //save the true monitor count
749        realsread[0] = defmon                                   //monitor ct = defmon
750        integersread[2] = total_rtime                   // total counting time
751        realsread[2] = scale*total_det                  //scaled detector counts
752        realsread[39] = scale*total_trn                 //scaled transmission counts
753       
754        //Add the added raw filename to the list of files in the workfile
755        String newfile = ";" + raw_text[0]
756        SVAR oldList = $(destPath + ":fileList")
757        String/G $(destPath + ":fileList") = oldList + newfile
758       
759        //reset the current displaytype to "newtype"
760        String/G root:myGlobals:gDataDisplayType=newType
761       
762        //return to root folder (redundant)
763        SetDataFolder root:
764       
765        Return(0)
766End
767
768
769//used for adding DRK (beam shutter CLOSED) data to a workfile
770//force the monitor count to 1, since it's irrelevant
771// run data through normal "add" step, then unscale default monitor counts
772//to get the data back on a simple time basis
773//
774Function Raw_to_Work_NoNorm(type)
775        String type
776       
777        WAVE reals=$("root:Packages:NIST:RAW:realsread")
778        reals[1]=1              //true monitor counts, still in raw
779        Raw_to_work(type)
780        //data is now in "type" folder
781        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
782        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
783        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
784        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
785       
786        Variable norm_mon,tot_mon,scale
787       
788        norm_mon = new_reals[0]         //should be 1e8
789        tot_mon = new_reals[1]          //should be 1
790        scale= norm_mon/tot_mon
791       
792        data /= scale           //unscale the data
793        data_err /= scale
794       
795        // to keep "data" and linear_data in sync
796        data_copy = data
797       
798        return(0)
799End
800
801//used for adding DRK (beam shutter CLOSED) data to a workfile
802//force the monitor count to 1, since it's irrelevant
803// run data through normal "add" step, then unscale default monitor counts
804//to get the data back on a simple time basis
805//
806Function Add_Raw_to_Work_NoNorm(type)
807        String type
808       
809        WAVE reals=$("root:Packages:NIST:RAW:realsread")
810        reals[1]=1              //true monitor counts, still in raw
811        Add_Raw_to_work(type)
812        //data is now in "type" folder
813        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
814        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
815        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
816        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
817       
818        Variable norm_mon,tot_mon,scale
819       
820        norm_mon = new_reals[0]         //should be 1e8
821        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
822        scale= norm_mon/tot_mon
823       
824        data /= scale           //unscale the data
825        data_err /= scale
826       
827        // to keep "data" and linear_data in sync
828        data_copy = data
829       
830        return(0)
831End
832
Note: See TracBrowser for help on using the repository browser.