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

Last change on this file since 540 was 540, checked in by srkline, 13 years ago

Changed T-Wide calculation to be normalized to the monitor counts, hopefully to avoid oddities of MCR fluctuations.

Also handle the case where data greater than 2 degrees is not present. in this case, a warning is presented, and the last 4 points are used, whatever the angles are.

  • Property eol-style set to native
  • Property svn:executable set to *
File size: 7.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma Version=2.20
3#pragma IgorVersion=6.0
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=200,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        Wave Angle = $(USANSFolder+":"+type+":Angle")
40        Wave DetCts = $(USANSFolder+":"+type+":DetCts")
41        Wave ErrDetCts = $(USANSFolder+":"+type+":ErrDetCts")
42        Wave MonCts = $(USANSFolder+":"+type+":MonCts")
43        Wave TransCts = $(USANSFolder+":"+type+":TransCts")
44       
45        Open/R refNum as fname          //if fname is "", a dialog will be presented
46        if(refnum==0)
47                return(1)               //user cancelled
48        endif
49        //read in the ASCII data line-by-line
50        Variable numLinesLoaded = 0,firstchar,countTime
51        Variable v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ii,valuesRead
52        String buffer ="",s1,s2,s3,s4,s5,s6,s7,s8,s9,s10
53        String filen="",fileLabel=""
54       
55        //parse the first line
56        FReadLine refnum,buffer
57        sscanf buffer, "%s%s%s%s%s%s%g%g%s%g%s",s1,s2,s3,s4,s5,s6,v1,v2,s7,v3,s8
58        ii=strlen(s1)
59        filen=s1[1,(ii-2)]                      //remove the ' from beginning and end
60        //Print "filen = ",filen,strlen(filen)
61        //v2 is the monitor prefactor. multiply monitor time by prefactor. AJJ 5 March 07
62        countTime=v1*v2
63        //Print "time (sec) = ",countTime
64       
65        //skip line 2
66        FReadLine refnum,buffer
67        //the next line is the sample label, use it all, minus the terminator
68        FReadLine refnum,filelabel
69       
70        //skip the next 10 lines
71        For(ii=0;ii<10;ii+=1)
72                FReadLine refnum,buffer
73        EndFor
74       
75        //read the data until EOF - assuming always a pair or lines
76        do
77                FReadLine refNum, buffer
78                if (strlen(buffer) == 0)
79                        break                                                   // End of file
80                endif
81                firstChar = char2num(buffer[0])
82                if ( (firstChar==10) || (firstChar==13) )
83                        break                                                   // Hit blank line. End of data in the file.
84                endif
85                sscanf buffer,"%g%g%g",v1,v2,v3         //v2,v3 not used
86                angle[numlinesloaded] = v1              //[0] is the ANGLE
87               
88                FReadLine refNum,buffer //assume a 2nd line is there, w/16 values
89                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
90                //valuesRead = V_flag
91                //print valuesRead
92                MonCts[numlinesloaded] = v1             //monitor
93                DetCts[numlinesloaded] = v2 + v3 + v5 + v6 + v7         //5 detectors
94                TransCts[numlinesloaded] = v4           //trans detector
95                ErrDetCts[numlinesloaded] = sqrt(detCts[numlinesloaded])
96                //values 8-16 are always zero
97                numlinesloaded += 1             //2 lines in file read, one "real" line of data
98        while(1)
99               
100        Close refNum            // Close the file.
101        //Print "numlines = ",numlinesloaded
102        //trim the waves to the correct number of points
103        Redimension/N=(numlinesloaded) Angle,MonCts,DetCts,TransCts,ErrDetCts
104       
105        //remove LF from end of filelabel
106        filelabel=fileLabel[0,(strlen(fileLabel)-2)]
107       
108        //set the wave note for the DetCts
109        String str=""
110        str = "FILE:"+filen+";"
111        str += "TIMEPT:"+num2str(counttime)+";"
112        str += "PEAKANG:;"              //no value yet
113        str += "STATUS:;"               //no value yet
114        str += "LABEL:"+fileLabel+";"
115        str += "PEAKVAL:;"              //no value yet
116        str += "TWIDE:0;"               //no value yet
117        Note DetCts,str
118       
119        String/G fileLbl = filelabel
120        return err                      // Zero signifies no error.     
121End
122
123
124//will search for peak in root:type:DetCts,  find and returns the angle
125// of the peak poistion (= q=0)
126//returns -9999 if error, appends value to wavenote
127// !can't return -1,0,1, since these may be the peak angle!
128//
129Function FindZeroAngle(type)
130        String type
131       
132        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
133       
134        Variable pkNotFound,pkPt,pkAngle,pkVal,temp
135        Wave angle = $(USANSFolder+":"+type+":Angle")
136        Wave detCts = $(USANSFolder+":"+type+":DetCts")
137
138
139        WaveStats/Q detcts
140        temp=V_Maxloc           //in points
141        FindPeak/P/Q/M=10/R=[(temp-10),(temp+10)] DetCts                //+/- 10 pts from maximum
142       
143        //FindPeak/P/Q/M=(1000) detCts          // /M=min ht -- peak must be at least 100x higher than 1st pt
144        pkNotFound=V_Flag               //V_Flag==1 if no pk found
145        pkPt = V_PeakLoc
146        pkVal = V_PeakVal               //for calc of T_rock
147        if(pkNotFound)
148                //DoAlert 0, "Peak not found"
149                Return (-9999)          //fatal error
150        Endif
151        pkAngle = Angle[pkPt]
152        //Print "Peak Angle = ",pkAngle
153        //update the note
154        String str=""
155        str=note(DetCts)
156        str = ReplaceNumberByKey("PEAKANG",str,pkangle,":",";")
157        str = ReplaceNumberByKey("PEAKVAL",str,pkVal,":",";")
158        Note/K DetCts
159        Note detCts,str
160       
161        return(pkAngle)
162End
163
164// given the peakAngle (q=0) and the "type" folder of scattering angle,
165// convert the angle to Q-values [degrees to (1/A)]
166// makes a new Qvals wave, duplicating the Angle wave, which is assumed to exist
167// Uses a conversion constant supplied by John Barker, and is hard-wired in
168//
169Function ConvertAngle2Qvals(type,pkAngle)
170        String type
171        Variable pkAngle
172
173        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
174       
175        Wave angle = $(USANSFolder+":"+type+":Angle")
176        Variable num=numpnts(angle)
177        Variable deg2QConv=5.55e-5              //JGB -- 2/24/01
178       
179        Make/O/N=(num) $(USANSFolder+":"+type+":Qvals")
180        Wave qvals = $(USANSFolder+":"+type+":Qvals")   
181        Qvals = deg2QConv*(angle[p] - pkAngle)
182       
183        return(0)       //no error
184End
185
186//updates the wavenote with the average Trans det cts for calculation of T_Wide
187//  finds the average number of counts on the transmission detector at angles
188// greater than 2 deg. (per John, practical experience where the trans data levels off)
189//
190// error value of 1 is returned and wavenote not updated if level is not found
191//
192// now normalizes to the monitor counts
193//
194Function FindTWideCts(type)
195        String type
196       
197        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
198       
199        Variable levNotFound,levPt,Cts,num,ii,sumMonCts
200       
201        Wave angle = $(USANSFolder+":"+type+":Angle")
202        Wave detCts = $(USANSFolder+":"+type+":DetCts")
203        Wave TransCts = $(USANSFolder+":"+type+":TransCts")
204        Wave monCts = $(USANSFolder+":"+type+":MonCts")
205       
206        num=numpnts(TransCts)
207
208
209        FindLevel/Q/P angle,2           //use angles greater than 2 deg
210        levNotFound=V_Flag              //V_Flag==1 if no pk found
211       
212        if(levNotFound)
213                //post a warning, and use just the last 4 points...
214                DoAlert 0,"You don't have "+type+" data past 2 degrees - so Twide may not be reliable"
215                levPt = num-4
216        else
217                levPt = trunc(V_LevelX)         // in points, force to integer
218        endif
219       
220        //average the trans cts from levPt to the end of the dataset
221        Cts=0
222        sumMonCts=0
223       
224        for(ii=levPt;ii<num;ii+=1)
225                Cts += transCts[ii]
226                sumMonCts += monCts[ii]
227        endfor
228        // get the average over that number of data points
229        sumMonCts /= (num-levPt-1)
230        Cts /= (num-levPt-1)
231       
232//      Print "cts = ",cts
233//      Print "sumMoncts = ",sumMoncts
234       
235        //normalize to the average monitor counts
236        Cts /= sumMonCts
237       
238//      Print "normalized counts = ",cts
239
240
241        //update the note
242        Wave DetCts = $(USANSFolder+":"+type+":DetCts")
243        String str,strVal
244        str=note(DetCts)
245        str = ReplaceNumberByKey("TWIDE",str,Cts,":",";")
246        Note/K DetCts
247        Note detCts,str
248
249        return(0)
250End
Note: See TracBrowser for help on using the repository browser.