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

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

1) speed improvement (2x + faster) in recalculation of VCALC configuration by using multithreading where possible.

2) Important bug fix for USANS to account for data collected in q-values, where each point is collected for a different lenght of time. This change had unintended consequences for normalization and dead time corrections.

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