source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/USANS/BT5_Loader.ipf

Last change on this file was 1238, checked in by srkline, 3 years ago

1) bug fix in the BT5 loader for USANS where an incorrect detector (of the 5) was specified (from B. Maranville). No impact on reduced data.

2) Cleaned up steps of generating DIV files and added file generation time into the DIV file that is saved.

  • Property eol-style set to native
  • Property svn:executable set to *
File size: 13.8 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma Version=2.20
3#pragma IgorVersion=6.1
4
5/////////////////////////
6// 101001 Vers. 1
7//
8// functions to load and parse the ASCII ICP files generated at BT5 USANS
9// Nearly all of the important information about the data is dragged around
10// with the DetCts wave as a wave note
11//
12
13        // as of Feb 6 2019, USANS is using NICE for data collection rather than ICP
14        // as a consequence, the data file is somewhat different. Specifically, the
15        // order of the detectors in the comma-delimited list is different, with the
16        // 5 main detectors all listed in order, rather than separated.
17        // --- I can either flag on the date, or set a global to switch
18        // I also do not know if there will be any other changes in the data file format
19       
20        // as of March 2019, data collected from NICE is collected in a single file, in terms of q-values
21        // not A2 angle. An additional consequence is that each data point can be collected for a different
22        // length of time. This is now read in per line and stored in a wave - and is used for the dead time
23        // correction and the normalization to monitor count.
24       
25// - thes wave note is a string of KEY:value items
26//
27//      str = "FILE:"+filen+";"
28//      str += "TIMEPT:"+num2str(counttime)+";"
29//      str += "PEAKANG:;"              //no value yet
30//      str += "STATUS:;"               //no value yet
31//      str += "LABEL:"+fileLabel+";"
32//      str += "PEAKVAL:;"              //no value yet
33//      str += "TWIDE:0;"               //no value yet
34//
35/////////////////////////
36
37
38//given the filename, loads a single data file into a "TYPE" folder, into a
39//series of 1D waves. counting time is associated with the DetCts wave as a note
40//"blank" note entries are set up here as well, so later they can be polled/updated
41Function LoadBT5File(fname,type)
42        String fname,type
43
44        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
45       
46        Variable num=500,err=0,refnum
47        Make/O/D/N=(num) $(USANSFolder+":"+type+":Angle")
48        Make/O/D/N=(num) $(USANSFolder+":"+type+":DetCts")
49        Make/O/D/N=(num) $(USANSFolder+":"+type+":ErrDetCts")
50        Make/O/D/N=(num) $(USANSFolder+":"+type+":MonCts")
51        Make/O/D/N=(num) $(USANSFolder+":"+type+":TransCts")
52        Make/O/D/N=(num) $(USANSFolder+":"+type+":CtTime")
53        Wave Angle = $(USANSFolder+":"+type+":Angle")
54        Wave DetCts = $(USANSFolder+":"+type+":DetCts")
55        Wave ErrDetCts = $(USANSFolder+":"+type+":ErrDetCts")
56        Wave MonCts = $(USANSFolder+":"+type+":MonCts")
57        Wave TransCts = $(USANSFolder+":"+type+":TransCts")
58        Wave CtTime = $(USANSFolder+":"+type+":CtTime")
59       
60        Open/R refNum as fname          //if fname is "", a dialog will be presented
61        if(refnum==0)
62                return(1)               //user cancelled
63        endif
64        //read in the ASCII data line-by-line
65        Variable numLinesLoaded = 0,firstchar,countTime
66        Variable v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ii,valuesRead
67        String buffer ="",s1,s2,s3,s4,s5,s6,s7,s8,s9,s10
68        String filen="",fileLabel="",filedt=""
69       
70        Variable MainDeadTime,TransDeadTime
71       
72        //parse the first line
73        FReadLine refnum,buffer
74        sscanf buffer, "%s%s%s%s%s%s%g%g%s%g%s",s1,s2,s3,s4,s5,s6,v1,v2,s7,v3,s8
75        ii=strlen(s1)
76        filen=s1[1,(ii-2)]                      //remove the ' from beginning and end
77        //Print "filen = ",filen,strlen(filen)
78        //v2 is the monitor prefactor. multiply monitor time by prefactor. AJJ 5 March 07
79        countTime=v1*v2
80        //Print "time (sec) = ",countTime
81       
82        //Deadtime correction
83        //Discovered significant deadtime on detectors in Oct 2010
84        //JGB and AJJ changed pre-amps during shutdown Oct 27 - Nov 7 2010
85        //Need different deadtime before and after 8th November 2010
86        filedt = s2+" "+s3+" "+s4+" "+s5
87        ii=strlen(filedt)
88        filedt = filedt[1,(ii-2)]
89               
90//      print filedt
91//      print BT5Date2Secs(filedt)
92//      print date2secs(2010,11,7)
93       
94        // as of Feb 6 2019, USANS is using NICE for data collection rather than ICP
95        // as a consequence, the data file is somewhat different. Specifically, the
96        // order of the detectors in the comma-delimited list is different, with the
97        // 5 main detectors all listed in order, rather than separated.
98        // --- I can either flag on the date, or set a global to switch
99        // I also do not know if there will be any other changes in the data file format
100       
101        Variable useNewDataFormat
102        NVAR    gRawUSANSisQvalues = root:Packages:NIST:gRawUSANSisQvalues
103       
104       
105// test by date
106        Variable thisFileSecs
107        NVAR switchSecs = root:Packages:NIST:USANS:Globals:MainPanel:gFileSwitchSecs
108        thisFileSecs = BT5DateTime2Secs(filedt)         // could use BT5Date2Secs() to exclude HR:MIN
109        if(thisFileSecs >= switchSecs)
110                useNewDataFormat = 1
111        else
112                useNewDataFormat = 0
113        endif
114
115       
116        USANS_DetectorDeadtime(filedt,MainDeadTime,TransDeadTime)
117       
118        //skip line 2
119        FReadLine refnum,buffer
120        //the next line is the sample label, use it all, minus the terminator
121        FReadLine refnum,filelabel
122       
123        //skip the next 10 lines
124        For(ii=0;ii<10;ii+=1)
125                FReadLine refnum,buffer
126        EndFor
127       
128        //read the data until EOF - assuming always a pair or lines
129        do
130                FReadLine refNum, buffer
131                if (strlen(buffer) == 0)
132                        break                                                   // End of file
133                endif
134                firstChar = char2num(buffer[0])
135                if ( (firstChar==10) || (firstChar==13) )
136                        break                                                   // Hit blank line. End of data in the file.
137                endif
138                //1st line of pair
139                sscanf buffer,"%g%g%g%g%g",v1,v2,v3,v4,v5               // 5 values here now
140                angle[numlinesloaded] = v1              //[0] is the ANGLE
141                if(gRawUSANSisQvalues==1)
142                        // in this mode, each data point is collected for a different time
143                        countTime = v2 * 60             // convert MIN to seconds
144                endif
145               
146               
147                FReadLine refNum,buffer //assume a 2nd line is there, w/16 values
148                sscanf buffer,"%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g",v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16
149                //valuesRead = V_flag
150                //print valuesRead
151                if(useNewDataFormat == 1)
152                // new order is MonCt, det1,det2,det3,det4,det5, Trans
153                        MonCts[numlinesloaded] = v1             //monitor
154                        v2 = v2/(1.0-v2*MainDeadTime/countTime)   // Deadtime correction
155                        v3 = v3/(1.0-v3*MainDeadTime/countTime)   // Deadtime correction
156                        v4 = v4/(1.0-v4*MainDeadTime/countTime)   // Deadtime correction        // bug fix 1/31/20 from B. Maranville
157                        v5 = v5/(1.0-v5*MainDeadTime/countTime)   // Deadtime correction
158                        v6 = v6/(1.0-v6*MainDeadTime/countTime)   // Deadtime correction
159                        DetCts[numlinesloaded] = v2 + v3 + v4 + v5 + v6         //5 detectors
160                        TransCts[numlinesloaded] = v7/(1.0-v7*TransDeadTime/countTime)          //trans detector+deadtime correction
161                        ErrDetCts[numlinesloaded] = sqrt(detCts[numlinesloaded])
162                        //values 8-16 are always zero
163                else
164                // this is the original format from ICP
165                        MonCts[numlinesloaded] = v1             //monitor
166                        v2 = v2/(1.0-v2*MainDeadTime/countTime)   // Deadtime correction
167                        v3 = v3/(1.0-v3*MainDeadTime/countTime)   // Deadtime correction
168                        v5 = v5/(1.0-v5*MainDeadTime/countTime)   // Deadtime correction
169                        v6 = v6/(1.0-v6*MainDeadTime/countTime)   // Deadtime correction
170                        v7 = v7/(1.0-v7*MainDeadTime/countTime)   // Deadtime correction
171                        DetCts[numlinesloaded] = v2 + v3 + v5 + v6 + v7         //5 detectors
172                        TransCts[numlinesloaded] = v4/(1.0-v4*TransDeadTime/countTime)          //trans detector+deadtime correction
173                        ErrDetCts[numlinesloaded] = sqrt(detCts[numlinesloaded])
174                        //values 8-16 are always zero
175                endif
176               
177                CtTime[numlinesloaded] = countTime
178               
179                numlinesloaded += 1             //2 lines in file read, one "real" line of data
180        while(1)
181               
182        Close refNum            // Close the file.
183        //Print "numlines = ",numlinesloaded
184        //trim the waves to the correct number of points
185        Redimension/N=(numlinesloaded) Angle,MonCts,DetCts,TransCts,ErrDetCts,CtTime
186       
187        //remove LF from end of filelabel
188        filelabel=fileLabel[0,(strlen(fileLabel)-2)]
189       
190        //set the wave note for the DetCts
191        String str=""
192        str = "FILE:"+filen+";"
193        str += "TIMEPT:"+num2str(counttime)+";"
194        str += "PEAKANG:;"              //no value yet
195        str += "STATUS:;"               //no value yet
196        str += "LABEL:"+fileLabel+";"
197        str += "PEAKVAL:;"              //no value yet
198        str += "TWIDE:0;"               //no value yet
199        Note DetCts,str
200       
201        String/G fileLbl = filelabel
202        return err                      // Zero signifies no error.     
203End
204
205
206//will search for peak in root:type:DetCts,  find and returns the angle
207// of the peak poistion (= q=0)
208//returns -9999 if error, appends value to wavenote
209// !can't return -1,0,1, since these may be the peak angle!
210//
211// MAR 2019
212// if the data is colleced with NICE, the "angle" is really q-values.
213// in this case the returned value is a q-value, and when this shift is applied to the
214// nominal (raw) q-values, the peak value is not converted from angle to q
215//
216//
217Function FindZeroAngle(type)
218        String type
219       
220        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
221       
222        Variable pkNotFound,pkPt,pkAngle,pkVal,temp
223        Wave angle = $(USANSFolder+":"+type+":Angle")
224        Wave detCts = $(USANSFolder+":"+type+":DetCts")
225
226
227        WaveStats/Q detcts
228        temp=V_Maxloc           //in points
229        FindPeak/P/Q/M=10/R=[(temp-10),(temp+10)] DetCts                //+/- 10 pts from maximum
230       
231        //FindPeak/P/Q/M=(1000) detCts          // /M=min ht -- peak must be at least 100x higher than 1st pt
232        pkNotFound=V_Flag               //V_Flag==1 if no pk found
233        pkPt = V_PeakLoc
234        pkVal = V_PeakVal               //for calc of T_rock
235        if(pkNotFound)
236                //DoAlert 0, "Peak not found"
237                Return (-9999)          //fatal error
238        Endif
239        pkAngle = Angle[pkPt]
240       
241        Print "Peak Angle = ",pkAngle
242        //update the note
243        String str=""
244        str=note(DetCts)
245        str = ReplaceNumberByKey("PEAKANG",str,pkangle,":",";")
246        str = ReplaceNumberByKey("PEAKVAL",str,pkVal,":",";")
247        Note/K DetCts
248        Note detCts,str
249       
250        return(pkAngle)
251End
252
253// given the peakAngle (q=0) and the "type" folder of scattering angle,
254// convert the angle to Q-values [degrees to (1/A)]
255// makes a new Qvals wave, duplicating the Angle wave, which is assumed to exist
256// Uses a conversion constant supplied by John Barker, and is hard-wired in
257//
258//      Mar 2019
259// if the data is collected from NICE, the data is in terms of Q, not angle since the
260// conversion factor has already been applied. In this case, the "angle" and pkAngle
261// are actually in terms of Q, and the deg2QConv value is set to == 1
262// -- then nothing needs to be changed
263//
264Function ConvertAngle2Qvals(type,pkAngle)
265        String type
266        Variable pkAngle
267
268        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
269       
270        Wave angle = $(USANSFolder+":"+type+":Angle")
271        Variable num=numpnts(angle)
272//      Variable deg2QConv=5.55e-5              //JGB -- 2/24/01
273        NVAR deg2QConv=root:Packages:NIST:USANS:Globals:MainPanel:deg2QConv
274
275       
276        Make/O/N=(num) $(USANSFolder+":"+type+":Qvals")
277        Wave qvals = $(USANSFolder+":"+type+":Qvals")   
278        Qvals = deg2QConv*(angle[p] - pkAngle)
279       
280        return(0)       //no error
281End
282
283//updates the wavenote with the average Trans det cts for calculation of T_Wide
284//  finds the average number of counts on the transmission detector at angles
285// greater than 2 deg. (per John, practical experience where the trans data levels off)
286//
287// error value of 1 is returned and wavenote not updated if level is not found
288//
289// now normalizes to the monitor counts
290//
291// 17 APR 2013 SRK - coverted the 2 degree "cutoff" to use q-values instead. Then the
292//  location is generic and can be used for either NCNR or KUSANS
293//
294Function FindTWideCts(type)
295        String type
296       
297        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
298       
299        Variable levNotFound,levPt,Cts,num,ii,sumMonCts
300       
301        Wave angle = $(USANSFolder+":"+type+":Angle")
302        Wave detCts = $(USANSFolder+":"+type+":DetCts")
303        Wave TransCts = $(USANSFolder+":"+type+":TransCts")
304        Wave monCts = $(USANSFolder+":"+type+":MonCts")
305        Wave Qvals = $(USANSFolder+":"+type+":Qvals")
306       
307        num=numpnts(TransCts)
308
309//      FindLevel/Q/P angle,2           //use angles greater than 2 deg
310//      Print "Using angle, pt = ",V_levelX
311       
312        FindLevel/Q/P Qvals,1e-4                //use angles greater than 2 deg = 2*5.55e-5 = 1e-4 (1/A)
313//      Print "Using Qval, pt = ",V_levelX
314
315        levNotFound=V_Flag              //V_Flag==1 if no pk found
316       
317        if(levNotFound)
318                //post a warning, and use just the last 4 points...
319//              DoAlert 0,"You don't have "+type+" data past 2 degrees - so Twide may not be reliable"
320                DoAlert 0,"You don't have "+type+" data past 1e-4 A-1 (~2 degrees) - so Twide may not be reliable"
321                levPt = num-4
322        else
323                levPt = trunc(V_LevelX)         // in points, force to integer
324        endif
325       
326        //average the trans cts from levPt to the end of the dataset
327        Cts=0
328        sumMonCts=0
329       
330        for(ii=levPt;ii<num;ii+=1)
331                Cts += transCts[ii]
332                sumMonCts += monCts[ii]
333        endfor
334        // get the average over that number of data points
335        sumMonCts /= (num-levPt-1)
336        Cts /= (num-levPt-1)
337       
338//      Print "cts = ",cts
339//      Print "sumMoncts = ",sumMoncts
340       
341        //normalize to the average monitor counts
342        Cts /= sumMonCts
343       
344//      Print "normalized counts = ",cts
345
346
347        //update the note
348        Wave DetCts = $(USANSFolder+":"+type+":DetCts")
349        String str,strVal
350        str=note(DetCts)
351        str = ReplaceNumberByKey("TWIDE",str,Cts,":",";")
352        Note/K DetCts
353        Note detCts,str
354
355        return(0)
356End
357
358
359Function BT5DateTime2Secs(datestring)
360        String datestring
361       
362        Variable bt5secs
363       
364        String monthnums = "Jan:1;Feb:2;Mar:3;Apr:4;May:5;Jun:6;Jul:7;Aug:8;Sep:9;Oct:10;Nov:11;Dec:12"
365               
366        Variable bt5month = str2num(StringByKey(stringfromlist(0,datestring, " "),monthnums))
367        Variable bt5day = str2num(stringfromlist(1,datestring," "))
368        Variable bt5year = str2num(stringfromlist(2,datestring," "))
369        Variable bt5hours = str2num(stringfromlist(0,stringfromlist(3,datestring," "),":"))
370        Variable bt5mins = str2num(stringfromlist(1,stringfromlist(3,datestring," "),":"))
371       
372        bt5secs = date2secs(bt5year,bt5month,bt5day) + 3600*bt5hours + 60*bt5mins
373
374        return bt5secs
375End
376
377Function BT5Date2Secs(datestring)
378        String datestring
379       
380        Variable bt5secs
381       
382        String monthnums = "Jan:1;Feb:2;Mar:3;Apr:4;May:5;Jun:6;Jul:7;Aug:8;Sep:9;Oct:10;Nov:11;Dec:12"
383               
384        Variable bt5month = str2num(StringByKey(stringfromlist(0,datestring, " "),monthnums))
385        Variable bt5day = str2num(stringfromlist(1,datestring," "))
386        Variable bt5year = str2num(stringfromlist(2,datestring," "))
387       
388        bt5secs = date2secs(bt5year,bt5month,bt5day)
389
390        return bt5secs
391End
Note: See TracBrowser for help on using the repository browser.