source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/WorkFileUtils.ipf

Last change on this file was 948, checked in by srkline, 8 years ago

Correcting the search in the patch panel for certain cases where multiple run numbers could be returned if they matched the year.

Adding routines and preferences to be able to add together raw data files with different attenuation. Not sure why this was not possible in the past. There must be a good reason for this. This is toggled on/off with a SANS preference checkbox. Default is OFF

Adding to Transmission "guessing" to show buttons to bump the guess to one more or one fewer character, rather than forcing an exit if the match is not correct.

File size: 52.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5//*******************
6// Vers 1.2 100901
7//
8//*******************
9// Utility procedures for handling of workfiles (each is housed in a separate datafolder)
10//
11// - adding data to a workfile
12// - copying workfiles to another folder
13//
14// - absolute scaling
15// - DIV detector sensitivity corrections
16//
17// - the WorkFile Math panel for simple image math
18// -
19// - adding work.drk data without normalizing to monitor counts
20//***************************
21
22//
23//
24// MARCH 2011 - changed the references that manipulate the data to explcitly work on the linear_data wave
25//                                      at the end of routines that maniplulate linear_data, it is copied back to "data" which is displayed
26//
27//
28
29//testing procedure, not called anymore
30Proc Add_to_Workfile(type, add)
31        String type,add
32        Prompt type,"WORK data type",popup,"SAM;EMP;BGD"
33        Prompt add,"Add to current WORK contents?",popup,"No;Yes"
34       
35        //macro will take whatever is in RAW folder and "ADD" it to the folder specified
36        //in the popup menu
37       
38        //"add" = yes/no, don't add to previous runs
39        //switch here - two separate functions to avoid (my) confusion
40        Variable err
41        if(cmpstr(add,"No")==0)
42                //don't add to prev work contents, copy RAW contents to work and convert
43                err = Raw_to_work(type)
44        else
45                //yes, add RAW to the current work folder contents
46                err = Add_raw_to_work(type)
47        endif
48       
49        String newTitle = "WORK_"+type
50        DoWindow/F SANS_Data
51        DoWindow/T SANS_Data, newTitle
52        KillStrings/Z newTitle
53       
54        //need to update the display with "data" from the correct dataFolder
55        fRawWindowHook()
56       
57End
58
59//will "ADD" the current contents of the RAW folder to the newType work folder
60//and will ADD the RAW contents to the existing content of the newType folder
61// - used when adding multiple runs together
62//(the function Raw_to_work(type) makes a fresh workfile)
63//
64//the current display type is updated to newType (global)
65Function Add_raw_to_work(newType)
66        String newType
67       
68        // NEW OCT 2014
69        // this corrects for adding raw data files with different attenuation   
70        // does nothing if the attenuation of RAW and destination are the same
71        NVAR doAdjustRAW_Atten = root:Packages:NIST:gDoAdjustRAW_Atten
72        if(doAdjustRAW_Atten)
73                Adjust_RAW_Attenuation(newType)
74        endif
75       
76        String destPath=""
77       
78        // if the desired workfile doesn't exist, let the user know, and just make a new one
79        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
80                Print "There is no old work file to add to - a new one will be created"
81                //call Raw_to_work(), then return from this function
82                Raw_to_Work(newType)
83                Return(0)               //does not generate an error - a single file was converted to work.newtype
84        Endif
85       
86        NVAR pixelsX = root:myGlobals:gNPixelsX
87        NVAR pixelsY = root:myGlobals:gNPixelsY
88       
89        //now make references to data in newType folder
90        DestPath="root:Packages:NIST:"+newType 
91        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data
92        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data
93        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data
94        WAVE/T textread=$(destPath + ":textread")
95        WAVE integersread=$(destPath + ":integersread")
96        WAVE realsread=$(destPath + ":realsread")
97       
98        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
99        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
100
101// 08/01 detector constants are now returned from a function, based on the detector type and beamline
102//      dt_ornl = 3.4e-6                //deadtime of Ordella detectors as of 30-AUG-99
103//      dt_ill=3.0e-6                   //Cerca detector deadtime constant as of 30-AUG-99
104
105        defmon=1e8                      //default monitor counts
106       
107        //Yes, add to previous run(s) in work, that does exist
108        //use the actual monitor count run.savmon rather than the normalized monitor count
109        //in run.moncnt and unscale the work data
110       
111        total_mon = realsread[1]        //saved monitor count
112        uscale = total_mon/defmon               //unscaling factor
113        total_det = uscale*realsread[2]         //unscaled detector count
114        total_trn = uscale*realsread[39]        //unscaled trans det count
115        total_numruns = integersread[3] //number of runs in workfile
116        total_rtime = integersread[2]           //total counting time in workfile
117        //retrieve workfile beamcenter
118        wrk_beamx = realsread[16]
119        wrk_beamy = realsread[17]
120        //unscale the workfile data in "newType"
121        //
122        //check for log-scaling and adjust if necessary
123        ConvertFolderToLinearScale(newType)
124        //
125        //then unscale the data array
126        data *= uscale
127        dest_data_err *= uscale
128       
129        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
130        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
131        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
132        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread"
133        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
134        WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
135       
136        //check for log-scaling of the raw data - make sure it's linear
137        ConvertFolderToLinearScale("RAW")
138       
139        // switches to control what is done, don't do the transmission correction for the BGD measurement
140        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCorr
141        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCorr
142        Variable doTrans = gDoTrans
143        if(cmpstr("BGD",newtype) == 0)
144                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
145        endif   
146       
147        DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans)   //applies correction to raw_data, and overwrites it
148       
149        //if RAW data is ILL type detector, correct raw_data for same counts being written to 4 pixels
150        if(cmpstr(raw_text[9], "ILL   ") == 0 )         //text field in header is 6 characters "ILL---"
151                raw_data /= 4
152                raw_data_err /= 4
153        endif
154       
155        //deadtime corrections to raw data
156        deadTime = DetectorDeadtime(raw_text[3],raw_text[9],dateAndTimeStr=raw_text[1],dtime=raw_reals[48])             //pick the correct detector deadtime, switch on date too
157        itim = raw_ints[2]
158        cntrate = sum(raw_data,-inf,inf)/itim           //080802 use data sum, rather than scaler value
159        dscale = 1/(1-deadTime*cntrate)
160
161#if (exists("ILL_D22")==6)
162        Variable tubeSum
163        // for D22 detector might need to use cntrate/128 as it is the tube response
164        for(ii=0;ii<pixelsX;ii+=1)
165                //sum the counts in each tube
166                tubeSum = 0
167                for(jj=0;jj<pixelsY;jj+=1)
168                        tubeSum += data[jj][ii]
169                endfor
170                // countrate in tube ii
171                cntrate = tubeSum/itim
172                // deadtime scaling in tube ii
173                dscale = 1/(1-deadTime*cntrate)
174                // multiply data[ii][] by the dead time
175                raw_data[][ii] *= dscale
176                raw_data_err[][ii] *= dscale
177        endfor
178#else
179        // dead time correction on all other RAW data, including NCNR
180        raw_data *= dscale
181        raw_data_err *= dscale
182#endif
183
184        //update totals by adding RAW values to the local ones (write to work header at end of function)
185        total_mon += raw_reals[0]
186#if (exists("ILL_D22")==6)
187        total_det += sum(raw_data,-inf,inf)                     //add the newly scaled detector array
188#else
189        total_det += dscale*raw_reals[2]
190#endif
191        total_trn += raw_reals[39]
192        total_rtime += raw_ints[2]
193        total_numruns +=1
194       
195        //do the beamcenter shifting if there is a mismatch
196        //and then add the two data sets together, changing "data" since it is the workfile data
197        xshift = raw_reals[16] - wrk_beamx
198        yshift = raw_reals[17] - wrk_beamy
199       
200        If((xshift != 0) || (yshift != 0))
201                DoAlert 1,"Do you want to ignore the beam center mismatch?"
202                if(V_flag==1)
203                        xshift=0
204                        yshift=0
205                endif
206        endif
207       
208        If((xshift == 0) && (yshift == 0))              //no shift, just add them
209                data += raw_data                //deadtime correction has already been done to the raw data
210                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum
211        else
212                //shift the beamcenter, then add
213                Make/O/N=1 $(destPath + ":noadd")               //needed to get noadd condition back from ShiftSum()
214                WAVE noadd = $(destPath + ":noadd")
215                Variable sh_sum                 //returned value
216                Print "BEAM CENTER MISMATCH - - BEAM CENTER WILL BE SHIFTED TO THIS FILE'S VALUE"
217                //ii,jj are just indices here, not physical locations - so [0,127] is fine
218                ii=0
219                do
220                        jj=0
221                        do
222                                //get the contribution of shifted data
223                                sh_sum = ShiftSum(data,ii,jj,xshift,yshift,noadd)
224                                if(noadd[0])
225                                        //don't do anything to data[][]
226                                else
227                                        //add the raw_data + shifted sum (and do the deadtime correction on both)
228                                        data[ii][jj] += (raw_data[ii][jj]+sh_sum)               //do the deadtime correction on RAW here
229                                        dest_data_err[ii][jj] = sqrt(dest_data_err[ii][jj]^2 + raw_data_err[ii][jj]^2)                  // error of the sum
230                                Endif
231                                jj+=1
232                        while(jj<pixelsY)
233                        ii+=1
234                while(ii<pixelsX)
235        Endif
236       
237        //scale the data to the default montor counts
238        scale = defmon/total_mon
239        data *= scale
240        dest_data_err *= scale
241       
242        // keep "data" and linear_data in sync in the destination folder
243        data_copy = data
244       
245        //all is done, except for the bookkeeping of updating the header info in the work folder
246        textread[1] = date() + " " + time()             //date + time stamp
247        integersread[3] = total_numruns                                         //numruns = more than one
248        realsread[1] = total_mon                        //save the true monitor count
249        realsread[0] = defmon                                   //monitor ct = defmon
250        integersread[2] = total_rtime                   // total counting time
251        realsread[2] = scale*total_det                  //scaled detector counts
252        realsread[39] = scale*total_trn                 //scaled transmission counts
253       
254        //Add the added raw filename to the list of files in the workfile
255        String newfile = ";" + raw_text[0]
256        SVAR oldList = $(destPath + ":fileList")
257        String/G $(destPath + ":fileList") = oldList + newfile
258       
259        //reset the current displaytype to "newtype"
260        String/G root:myGlobals:gDataDisplayType=newType
261       
262        //return to root folder (redundant)
263        SetDataFolder root:
264       
265        Return(0)
266End
267
268//will copy the current contents of the RAW folder to the newType work folder
269//and do the geometric corrections and normalization to monitor counts
270//(the function Add_Raw_to_work(type) adds multiple runs together)
271//
272//the current display type is updated to newType (global)
273//
274Function Raw_to_work(newType)
275        String newType
276       
277        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
278        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy
279        String destPath
280       
281// 08/01 detector constants are now returned from a function, based on the detector type and beamline
282//      dt_ornl = 3.4e-6                //deadtime of Ordella detectors as of 30-AUG-99
283//      dt_ill=3.0e-6                   //Cerca detector deadtime constant as of 30-AUG-99
284        defmon=1e8                      //default monitor counts
285       
286        //initialize values before normalization
287        total_mon=0
288        total_det=0
289        total_trn=0
290        total_numruns=0
291        total_rtime=0
292       
293        //Not adding multiple runs, so wipe out the old contents of the work folder and
294        // replace with the contents of raw
295
296        destPath = "root:Packages:NIST:" + newType
297       
298        //check for log-scaling of the RAW data and adjust if necessary
299        ConvertFolderToLinearScale("RAW")
300        //then continue
301       
302        NVAR pixelsX = root:myGlobals:gNPixelsX
303        NVAR pixelsY = root:myGlobals:gNPixelsY
304       
305        //copy from current dir (RAW) to work, defined by destpath
306        DestPath = "root:Packages:NIST:"+newType
307        Duplicate/O $"root:Packages:NIST:RAW:data",$(destPath + ":data")
308        Duplicate/O $"root:Packages:NIST:RAW:linear_data_error",$(destPath + ":linear_data_error")
309        Duplicate/O $"root:Packages:NIST:RAW:linear_data",$(destPath + ":linear_data")
310//      Duplicate/O $"root:Packages:NIST:RAW:vlegend",$(destPath + ":vlegend")
311        Duplicate/O $"root:Packages:NIST:RAW:textread",$(destPath + ":textread")
312        Duplicate/O $"root:Packages:NIST:RAW:integersread",$(destPath + ":integersread")
313        Duplicate/O $"root:Packages:NIST:RAW:realsread",$(destPath + ":realsread")
314        Variable/G $(destPath + ":gIsLogscale")=0                       //overwite flag in newType folder, data converted (above) to linear scale
315
316        WAVE data=$(destPath + ":linear_data")                                                  // these wave references point to the data in work
317        WAVE data_copy=$(destPath + ":data")                                                    // these wave references point to the data in work
318        WAVE data_err=$(destPath + ":linear_data_error")
319        WAVE/T textread=$(destPath + ":textread")                               //that are to be directly operated on
320        WAVE integersread=$(destPath + ":integersread")
321        WAVE realsread=$(destPath + ":realsread")
322        String/G $(destPath + ":fileList") = textread[0]                        //a list of names of the files in the work file (1)             //02JUL13
323       
324        //apply nonlinear, Jacobian corrections ---
325        // switches to control what is done, don't do the transmission correction for the BGD measurement
326        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCorr
327        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCorr
328        Variable doTrans = gDoTrans
329        if(cmpstr("BGD",newtype) == 0)
330                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
331        endif
332       
333        DetCorr(data,data_err,realsread,doEfficiency,doTrans)           //the parameters are waves, and will be changed by the function
334
335        //if ILL type detector, correct for same counts being written to 4 pixels
336        if(cmpstr(textread[9], "ILL   ") == 0 )         //text field in header is 6 characters "ILL---"
337                data /= 4
338                data_err /= 4           //rescale error
339        endif
340       
341        //deadtime corrections
342        itim = integersread[2]
343        cntrate = sum(data,-inf,inf)/itim               //use sum of detector counts rather than scaler value
344        deadtime = DetectorDeadtime(textread[3],textread[9],dateAndTimeStr=textRead[1],dtime=realsRead[48])     //pick the correct deadtime
345        dscale = 1/(1-deadTime*cntrate)
346       
347#if (exists("ILL_D22")==6)
348        Variable tubeSum
349        // for D22 detector might need to use cntrate/128 as it is the tube response
350        for(ii=0;ii<pixelsX;ii+=1)
351                //sum the counts in each tube
352                tubeSum = 0
353                for(jj=0;jj<pixelsY;jj+=1)
354                        tubeSum += data[jj][ii]
355                endfor
356                // countrate in tube ii
357                cntrate = tubeSum/itim
358                // deadtime scaling in tube ii
359                dscale = 1/(1-deadTime*cntrate)
360                // multiply data[ii][] by the dead time
361                data[][ii] *= dscale
362                data_err[][ii] *= dscale
363        endfor
364#endif
365
366        // NO xcenter,ycenter shifting is done - this is the first (and only) file in the work folder
367       
368        //only ONE data file- no addition of multiple runs in this function, so data is
369        //just simply corrected for deadtime.
370#if (exists("ILL_D22")==0)              //ILL correction done tube-by-tube above
371        data *= dscale          //deadtime correction for everyone else, including NCNR
372        data_err *= dscale
373#endif
374       
375        //update totals to put in the work header (at the end of the function)
376        total_mon += realsread[0]
377#if (exists("ILL_D22")==6)
378        total_det += sum(data,-inf,inf)                 //add the newly scaled detector array
379#else
380        total_det += dscale*realsread[2]
381#endif
382        total_trn += realsread[39]
383        total_rtime += integersread[2]
384        total_numruns +=1
385       
386       
387        //scale the data to the default montor counts
388        scale = defmon/total_mon
389        data *= scale
390        data_err *= scale               //assumes total monitor count is so large there is essentially no error
391       
392        // to keep "data" and linear_data in sync
393        data_copy = data
394       
395        //all is done, except for the bookkeeping, updating the header information in the work folder
396        textread[1] = date() + " " + time()             //date + time stamp
397        integersread[3] = total_numruns                                         //numruns = 1
398        realsread[1] = total_mon                        //save the true monitor count
399        realsread[0] = defmon                                   //monitor ct = defmon
400        integersread[2] = total_rtime                   // total counting time
401        realsread[2] = scale*total_det                  //scaled detector counts
402        realsread[39] = scale*total_trn                 //scaled transmission counts
403       
404        //reset the current displaytype to "newtype"
405        String/G root:myGlobals:gDataDisplayType=newType
406       
407        //return to root folder (redundant)
408        SetDataFolder root:
409       
410        Return(0)
411End
412
413//used for adding DRK (beam shutter CLOSED) data to a workfile
414//force the monitor count to 1, since it's irrelevant
415// run data through normal "add" step, then unscale default monitor counts
416//to get the data back on a simple time basis
417//
418Function Raw_to_Work_NoNorm(type)
419        String type
420       
421        WAVE reals=$("root:Packages:NIST:RAW:realsread")
422        reals[1]=1              //true monitor counts, still in raw
423        Raw_to_work(type)
424        //data is now in "type" folder
425        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
426        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
427        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
428        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
429       
430        Variable norm_mon,tot_mon,scale
431       
432        norm_mon = new_reals[0]         //should be 1e8
433        tot_mon = new_reals[1]          //should be 1
434        scale= norm_mon/tot_mon
435       
436        data /= scale           //unscale the data
437        data_err /= scale
438       
439        // to keep "data" and linear_data in sync
440        data_copy = data
441       
442        return(0)
443End
444
445//used for adding DRK (beam shutter CLOSED) data to a workfile
446//force the monitor count to 1, since it's irrelevant
447// run data through normal "add" step, then unscale default monitor counts
448//to get the data back on a simple time basis
449//
450Function Add_Raw_to_Work_NoNorm(type)
451        String type
452       
453        WAVE reals=$("root:Packages:NIST:RAW:realsread")
454        reals[1]=1              //true monitor counts, still in raw
455        Add_Raw_to_work(type)
456        //data is now in "type" folder
457        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
458        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
459        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
460        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
461       
462        Variable norm_mon,tot_mon,scale
463       
464        norm_mon = new_reals[0]         //should be 1e8
465        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
466        scale= norm_mon/tot_mon
467       
468        data /= scale           //unscale the data
469        data_err /= scale
470       
471        // to keep "data" and linear_data in sync
472        data_copy = data
473       
474        return(0)
475End
476
477//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
478//function is called by Raw_to_work() and Add_raw_to_work() functions
479//works on the actual data array, assumes that is is already on LINEAR scale
480//
481Function DetCorr(data,data_err,realsread,doEfficiency,doTrans)
482        Wave data,data_err,realsread
483        Variable doEfficiency,doTrans
484       
485        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
486        Variable ii,jj,dtdist,dtdis2
487        Variable xi,xd,yd,rad,ratio,domega,xy
488        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
489       
490//      Print "...doing jacobian and non-linear corrections"
491
492        NVAR pixelsX = root:myGlobals:gNPixelsX
493        NVAR pixelsY = root:myGlobals:gNPixelsY
494       
495        //set up values to send to auxiliary trig functions
496        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
497        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
498
499        x0 = realsread[16]
500        y0 = realsread[17]
501        sx = realsread[10]
502        sx3 = realsread[11]
503        sy = realsread[13]
504        sy3 = realsread[14]
505       
506        dtdist = 1000*realsread[18]     //sdd in mm
507        dtdis2 = dtdist^2
508       
509        lambda = realsRead[26]
510        trans = RealsRead[4]
511        trans_err = RealsRead[41]               //new, March 2011
512       
513        xx0 = dc_fx(x0,sx,sx3,xcenter)
514        yy0 = dc_fy(y0,sy,sy3,ycenter)
515       
516
517        //waves to contain repeated function calls
518        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
519        ii=0
520        do
521                xi = ii
522                fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
523                xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
524                yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
525                ii+=1
526        while(ii<pixelsX)
527       
528        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
529       
530        ii=0
531        do
532                xi = ii
533                xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
534                jj=0
535                do
536                        yd = fyy[jj]-yy0
537                        //rad is the distance of pixel ij from the sample
538                        //domega is the ratio of the solid angle of pixel ij versus center pixel
539                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
540                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
541                        rad = sqrt(dtdis2 + xd^2 + yd^2)
542                        domega = rad/dtdist
543                        ratio = domega^3
544                        xy = xx[ii]*yy[jj]
545                       
546                        data[ii][jj] *= xy*ratio
547                       
548                        solidAngle[ii][jj] = xy*ratio           //testing only 
549                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
550                       
551                       
552                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
553                        // correction inserted 11/2007 SRK
554                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
555                        // so divide here to get the correct answer (5/22/08 SRK)
556                        if(doEfficiency)
557#if (exists("ILL_D22")==6)
558                                data[ii][jj]  /= DetEffCorrILL(lambda,dtdist,xd)                //tube-by-tube corrections
559                                data_err[ii][jj] /= DetEffCorrILL(lambda,dtdist,xd)                     //assumes correction is exact
560//                solidAngle[ii][jj] = DetEffCorrILL(lambda,dtdist,xd)
561#else
562                                data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
563                                data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
564//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
565#endif
566                        endif
567                       
568                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
569                        // so divide here to get the correct answer
570                        if(doTrans)
571                       
572                                if(trans<0.1 && ii==0 && jj==0)
573                                        Print "***transmission is less than 0.1*** and is a significant correction"
574                                endif
575                               
576                                if(trans==0)
577                                        if(ii==0 && jj==0)
578                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
579                                        endif
580                                        trans = 1
581                                endif
582                               
583                                // pass in the transmission error, and the error in the correction is returned as the last parameter
584                                lat_corr = LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)             //moved from 1D avg SRK 11/2007
585                                data[ii][jj] /= lat_corr                        //divide by the correction factor
586                                //
587                                //
588                                //
589                                // relative errors add in quadrature
590                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
591                                tmp_err = sqrt(tmp_err)
592                               
593                                data_err[ii][jj] = tmp_err
594                               
595//                              solidAngle[ii][jj] = lat_err
596
597                               
598                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
599                        endif
600                       
601                        jj+=1
602                while(jj<pixelsX)
603                ii+=1
604        while(ii<pixelsX)
605       
606        //clean up waves
607//      Killwaves/Z fyy,xx,yy
608       
609        Return(0)
610End
611
612//trig function used by DetCorr()
613Function dc_fx(x,sx,sx3,xcenter)
614        Variable x,sx,sx3,xcenter
615       
616        Variable result
617       
618        result = sx3*tan((x-xcenter)*sx/sx3)
619        Return(result)
620End
621
622//trig function used by DetCorr()
623Function dc_fy(y,sy,sy3,ycenter)
624        Variable y,sy,sy3,ycenter
625       
626        Variable result
627       
628        result = sy3*tan((y-ycenter)*sy/sy3)
629        Return(result)
630End
631
632//trig function used by DetCorr()
633Function dc_fxn(x,sx,sx3,xcenter)
634        Variable x,sx,sx3,xcenter
635       
636        Variable result
637       
638        result = (cos((x-xcenter)*sx/sx3))^2
639        Return(result)
640End
641
642//trig function used by DetCorr()
643Function dc_fym(y,sy,sy3,ycenter)
644        Variable y,sy,sy3,ycenter
645       
646        Variable result
647       
648        result = (cos((y-ycenter)*sy/sy3))^2
649        Return(result)
650End
651
652//distances passed in are in mm
653// dtdist is SDD
654// xd and yd are distances from the beam center to the current pixel
655//
656Function DetEffCorr(lambda,dtdist,xd,yd)
657        Variable lambda,dtdist,xd,yd
658       
659        Variable theta,cosT,ff,stAl,stHe
660       
661        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )
662        cosT = cos(theta)
663       
664        stAl = 0.00967*lambda*0.8               //dimensionless, constants from JGB memo
665        stHe = 0.146*lambda*2.5
666       
667        ff = exp(-stAl/cosT)*(1-exp(-stHe/cosT)) / ( exp(-stAl)*(1-exp(-stHe)) )
668               
669        return(ff)
670End
671
672// DIVIDE the intensity by this correction to get the right answer
673Function LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err)
674        Variable trans,dtdist,xd,yd,trans_err,&err
675
676        //angle dependent transmission correction
677        Variable uval,arg,cos_th,correction,theta
678       
679        ////this section is the trans_correct() VAX routine
680//      if(trans<0.1)
681//              Print "***transmission is less than 0.1*** and is a significant correction"
682//      endif
683//      if(trans==0)
684//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
685//              trans = 1
686//      endif
687       
688        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )              //theta at the input pixel
689       
690        //optical thickness
691        uval = -ln(trans)               //use natural logarithm
692        cos_th = cos(theta)
693        arg = (1-cos_th)/cos_th
694       
695        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
696        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
697        // OR
698        if((uval<0.01) || (cos_th>0.99))       
699                //small arg, approx correction
700                correction= 1-0.5*uval*arg
701        else
702                //large arg, exact correction
703                correction = (1-exp(-uval*arg))/(uval*arg)
704        endif
705
706        Variable tmp
707       
708        if(trans == 1)
709                err = 0         //no correction, no error
710        else
711                //sigT, calculated from the Taylor expansion
712                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
713                tmp *= tmp
714                tmp *= trans_err^2
715                tmp = sqrt(tmp)         //sigT
716               
717                err = tmp
718        endif
719       
720//      Printf "trans error = %g\r",trans_err
721//      Printf "correction = %g +/- %g\r", correction, err
722       
723        //end of transmission/pathlength correction
724
725        return(correction)
726end
727
728//******************
729//direct port of the FORTRAN code for calculating the weighted
730//shifted element to add when beam centers in data headers do not match
731//(indices updated to [0,n-1] indexing rather than (1,n) of fortran
732//
733// as of IGOR 4.0, could be rewritten to pass-by-reference noadd, rather than wave, but the function
734// is so little used, it's not worth the time
735Function ShiftSum(DATA,ip,jp,XSHIFT,YSHIFT,noadd)
736        Wave data
737        Variable ip,jp,xshift,yshift
738        Wave noadd
739//
740//       COMPUTE WEIGHTED OFFSET ELEMENT SUM FOR USE IN SANS DATA
741//       ANALYSIS MODULES.
742//
743// "data" wave passed in is the current contents of the work file
744// sum_val is the return value of the function
745// "noadd" is passed back to the calling function as a one-point wave
746
747        Variable XDELTA,YDELTA,kk,II,JJ,ISHIFT,JSHIFT,sum_val
748        Make/O/N=4 iii,jjj,a
749
750//       -----------------------------------------------------------------
751
752        ISHIFT = trunc(XSHIFT)          // INTEGER PART, trunc gives int closest in dierction of zero
753        XDELTA = XSHIFT - ISHIFT        //FRACTIONAL PART.
754        JSHIFT = trunc(YSHIFT)
755        YDELTA = YSHIFT - JSHIFT
756        II = ip + ISHIFT
757        JJ = jp + JSHIFT
758
759//       SHIFT IS DEFINED AS A VECTOR ANCHORED AT THE STATIONARY CENTER
760//       AND POINTING TO THE MOVABLE CENTER.  THE MOVABLE FIELD IS THUS
761//       ACTUALLY MOVED BY -SHIFT.
762//
763        IF ((XDELTA>= 0) && (YDELTA >= 0))              // CASE I ---- "%&" is "and"
764                III[0] = II
765                JJJ[0] = JJ
766                III[1] = II + 1
767                JJJ[1] = JJ
768                III[2] = II + 1
769                JJJ[2] = JJ + 1
770                III[3] = II
771                JJJ[3] = JJ + 1
772                A[0] = (1. - XDELTA)*(1. - YDELTA)
773                A[1] = XDELTA*(1. - YDELTA)
774                A[2] = XDELTA*YDELTA
775                A[3] = (1. - XDELTA)*YDELTA
776        Endif
777        IF ((XDELTA >= 0) && (YDELTA < 0))              // CASE II.
778                III[0] = II
779                JJJ[0] = JJ
780                III[1] = II
781                JJJ[1] = JJ - 1
782                III[2] = II + 1
783                JJJ[2] = JJ - 1
784                III[3] = II + 1
785                JJJ[3] = JJ
786                A[0] = (1. - XDELTA)*(1. + YDELTA)
787                A[1] = (1. - XDELTA)*(-YDELTA)
788                A[2] = XDELTA*(-YDELTA)
789                A[3] = XDELTA*(1. + YDELTA)
790        Endif
791        IF ((XDELTA < 0) && (YDELTA >= 0))              // CASE III.
792                III[0] = II
793                JJJ[0] = JJ
794                III[1] = II
795                JJJ[1] = JJ + 1
796                III[2] = II - 1
797                JJJ[2] = JJ + 1
798                III[3] = II - 1
799                JJJ[3] = JJ
800                A[0] = (1. + XDELTA)*(1 - YDELTA)
801                A[1] = (1. + XDELTA)*YDELTA
802                A[2] = -XDELTA*YDELTA
803                A[3] = -XDELTA*(1. - YDELTA)
804        Endif
805        IF ((XDELTA < 0) && (YDELTA < 0))               //CASE IV.
806                III[0] = II
807                JJJ[0] = JJ
808                III[1] = II - 1
809                JJJ[1] = JJ
810                III[2] = II - 1
811                JJJ[2] = JJ - 1
812                III[3] = II
813                JJJ[3] = JJ - 1
814                A[0] = (1. + XDELTA)*(1. + YDELTA)
815                A[1] = -XDELTA*(1. + YDELTA)
816                A[2] = (-XDELTA)*(-YDELTA)
817                A[3] = (1. + XDELTA)*(-YDELTA)
818        Endif
819
820        NVAR pixelsX = root:myGlobals:gNPixelsX
821        NVAR pixelsY = root:myGlobals:gNPixelsY
822//check to see if iii[0],jjj[0] are valid detector elements, in [0,127]
823//if not set noadd[0] to 1, to let calling routine know NOT to add
824//        CALL TESTIJ(III(1),JJJ(1),OKIJ)
825        NOADD[0] = 0
826        if( (iii[0]<0) || (iii[0]>(pixelsX-1)) )
827                noadd[0] = 1
828        endif
829        if((jjj[0]<0) || (jjj[0]>(pixelsY-1)) )
830                noadd[0] = 1
831        endif
832       
833
834       
835        sum_val = 0.
836        kk = 0
837        Do
838                IF(JJJ[kk] == pixelsX)
839                        //do nothing
840                else
841                        sum_val += A[kk]*DATA[III[kk]][JJJ[kk]]
842                endif
843                kk+=1
844        while(kk<4)
845       
846        //clean up waves
847        KillWaves/z iii,jjj,a
848       
849        RETURN (sum_val)
850       
851End             //function ShiftSum
852
853//************************
854//unused testing procedure, may not be up-to-date with other procedures
855//check before re-implementing
856//
857Proc DIV_a_Workfile(type)
858        String type
859        Prompt type,"WORK data type",popup,"COR;SAM;EMP;BGD"
860       
861        //macro will take whatever is in SELECTED folder and DIVide it by the current
862        //contents of the DIV folder - the function will check for existence
863        //before proceeding
864       
865        Variable err
866        err = Divide_work(type)         //returns err = 1 if data doesn't exist in specified folders
867       
868        if(err)
869                Abort "error in Divide_work"
870        endif
871       
872        //contents are always dumped to CAL
873        type = "CAL"
874       
875        String newTitle = "WORK_"+type
876        DoWindow/F SANS_Data
877        DoWindow/T SANS_Data, newTitle
878        KillStrings/Z newTitle
879       
880        //need to update the display with "data" from the correct dataFolder
881        //reset the current displaytype to "type"
882        String/G root:myGlobals:gDataDisplayType=Type
883       
884        fRawWindowHook()
885       
886End
887
888//function will divide the contents of "type" folder with the contents of
889//the DIV folder
890// all data is converted to linear scale for the calculation
891//
892Function Divide_work(type)
893        String type
894       
895        //check for existence of data in type and DIV
896        // if the desired workfile doesn't exist, let the user know, and abort
897        String destPath=""
898
899        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
900                Print "There is no work file in "+type+"--Aborting"
901                Return(1)               //error condition
902        Endif
903        //check for DIV
904        // if the DIV workfile doesn't exist, let the user know,and abort
905
906        if(WaveExists($"root:Packages:NIST:DIV:data") == 0)
907                Print "There is no work file in DIV --Aborting"
908                Return(1)               //error condition
909        Endif
910        //files exist, proceed
911       
912        //check for log-scaling of the "DIV" data and adjust if necessary
913        ConvertFolderToLinearScale("DIV")
914       
915        //copy type information to CAL, wiping out the old contents of the CAL folder first
916       
917        //destPath = "root:Packages:NIST:CAL"
918        //SetDataFolder destPath
919        //KillWaves/A/Z                 //get rid of the old data in CAL folder
920
921        //check for log-scaling of the "type" data and adjust if necessary
922        ConvertFolderToLinearScale(type)
923        //then continue
924
925        //copy from current dir (type)=destPath to CAL, overwriting CAL contents
926        destPath = "root:Packages:NIST:" + type
927        Duplicate/O $(destPath + ":data"),$"root:Packages:NIST:CAL:data"
928        Duplicate/O $(destPath + ":linear_data"),$"root:Packages:NIST:CAL:linear_data"
929        Duplicate/O $(destPath + ":linear_data_error"),$"root:Packages:NIST:CAL:linear_data_error"
930//      Duplicate/O $(destPath + ":vlegend"),$"root:Packages:NIST:CAL:vlegend"
931        Duplicate/O $(destPath + ":textread"),$"root:Packages:NIST:CAL:textread"
932        Duplicate/O $(destPath + ":integersread"),$"root:Packages:NIST:CAL:integersread"
933        Duplicate/O $(destPath + ":realsread"),$"root:Packages:NIST:CAL:realsread"
934        //need to save a copy of filelist string too (from the current type folder)
935        SVAR oldFileList = $(destPath + ":fileList")
936
937        //now switch to reference waves in CAL folder
938        destPath = "root:Packages:NIST:CAL"
939        //make appropriate wave references
940        Wave data=$(destPath + ":linear_data")                                  // these wave references point to the data in CAL
941//      Wave data_err=$(destPath + ":linear_data_err")                                  // these wave references point to the data in CAL
942        Wave data_copy=$(destPath + ":data")                                    // these wave references point to the data in CAL
943        Wave/t textread=$(destPath + ":textread")                       //that are to be directly operated on
944        Wave integersread=$(destPath + ":integersread")
945        Wave realsread=$(destPath + ":realsread")
946        Variable/G $(destPath + ":gIsLogScale")=0                       //make new flag in CAL folder, data is linear scale
947        //need to copy filelist string too
948        String/G $(destPath + ":fileList") = oldFileList
949
950        Wave div_data = $"root:Packages:NIST:DIV:data"          //hard-wired in....
951        //do the division, changing data in CAL
952        data /= div_data
953       
954//      data_err /= div_data
955       
956        // keep "data" in sync with linear_data
957        data_copy = data
958       
959        //update CAL header
960        textread[1] = date() + " " + time()             //date + time stamp
961       
962        Return(0)
963End
964
965//test procedure, not called anymore
966Proc AbsoluteScaling(type,c0,c1,c2,c3,c4,c5)
967        String type
968        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0
969        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
970        Prompt c0, "Sample Transmission"
971        Prompt c1, "Sample Thickness (cm)"
972        Prompt c2, "Standard Transmission"
973        Prompt c3, "Standard Thickness (cm)"
974        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
975        Prompt c5, "Standard Cross-Section (cm-1)"
976
977        Variable err
978        //call the function to do the math
979        //data from "type" will be scaled and deposited in ABS
980        err = Absolute_Scale(type,c0,c1,c2,c3,c4,c5)
981       
982        if(err)
983                Abort "Error in Absolute_Scale()"
984        endif
985       
986        //contents are always dumped to ABS
987        type = "ABS"
988       
989        String newTitle = "WORK_"+type
990        DoWindow/F SANS_Data
991        DoWindow/T SANS_Data, newTitle
992        KillStrings/Z newTitle
993       
994        //need to update the display with "data" from the correct dataFolder
995        //reset the current displaytype to "type"
996        String/G root:myGlobals:gDataDisplayType=Type
997       
998        fRawWindowHook()
999       
1000End
1001
1002//s_ is the standard
1003//w_ is the "work" file
1004//both are work files and should already be normalized to 10^8 monitor counts
1005Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
1006        String type
1007        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1008       
1009        //convert the "type" data to absolute scale using the given standard information
1010        //copying the "type" waves to ABS
1011       
1012        //check for existence of data, rescale to linear if needed
1013        String destPath
1014        //check for "type"
1015        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
1016                Print "There is no work file in "+type+"--Aborting"
1017                Return(1)               //error condition
1018        Endif
1019        //check for log-scaling of the "type" data and adjust if necessary
1020        destPath = "root:Packages:NIST:"+Type
1021        NVAR gIsLogScale = $(destPath + ":gIsLogScale")
1022        if(gIsLogScale)
1023                Duplicate/O $(destPath + ":linear_data") $(destPath + ":data")//back to linear scale
1024                Variable/G $(destPath + ":gIsLogScale")=0       //the "type" data is not logscale anymore
1025        endif
1026       
1027        //copy "oldtype" information to ABS
1028        //overwriting out the old contents of the ABS folder (/O option in Duplicate)
1029        //copy over the waves data,vlegend,text,integers,reals(read)
1030
1031        String oldType= "root:Packages:NIST:"+type              //this is where the data to be absoluted is
1032        //copy from current dir (type) to ABS, defined by destPath
1033        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data"
1034        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data"
1035        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error"
1036//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend"
1037        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread"
1038        Duplicate/O $(oldType + ":integersread"),$"root:Packages:NIST:ABS:integersread"
1039        Duplicate/O $(oldType + ":realsread"),$"root:Packages:NIST:ABS:realsread"
1040        //need to save a copy of filelist string too (from the current type folder)
1041        SVAR oldFileList = $(oldType + ":fileList")
1042        //need to copy filelist string too
1043        String/G $"root:Packages:NIST:ABS:fileList" = oldFileList
1044       
1045        //now switch to ABS folder
1046        //make appropriate wave references
1047        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS
1048        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS
1049        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display
1050        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on
1051        WAVE integersread=$"root:Packages:NIST:ABS:integersread"
1052        WAVE realsread=$"root:Packages:NIST:ABS:realsread"
1053        Variable/G $"root:Packages:NIST:ABS:gIsLogscale"=0                      //make new flag in ABS folder, data is linear scale
1054       
1055        //do the actual absolute scaling here, modifying the data in ABS
1056        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1057       
1058        w_moncount = realsread[0]               //monitor count in "type"
1059        if(w_moncount == 0)
1060                //zero monitor counts will give divide by zero ---
1061                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1062                Return(1)               //report error
1063        Endif
1064       
1065        //calculate scale factor
1066        Variable scale,trans_err
1067        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1)
1068        s2 = s_thick/w_thick
1069        s3 = s_trans/w_trans
1070        s4 = s_cross/s_izero
1071       
1072        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1073       
1074        data *= s1*s2*s3*s4
1075       
1076        scale = s1*s2*s3*s4
1077        trans_err = realsRead[41]
1078       
1079//      print scale
1080//      print data[0][0]
1081       
1082        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1083
1084//      print data_err[0][0]
1085       
1086// keep "data" in sync with linear_data
1087        data_copy = data
1088       
1089        //********* 15APR02
1090        // DO NOt correct for atenuators here - the COR step already does this, putting all of the data one equal
1091        // footing (zero atten) before doing the subtraction.
1092        //
1093        //Print "ABS data multiplied by  ",s1*s2*s3*s4/attenFactor
1094       
1095        //update the ABS header information
1096        textread[1] = date() + " " + time()             //date + time stamp
1097       
1098        Return (0) //no error
1099End
1100
1101//*************
1102// start of section of functions used for workfile math panel
1103//*************
1104
1105
1106//function will copy the contents of oldtype folder to newtype folder
1107//converted to linear scale before copying
1108//******data in newtype is overwritten********
1109//
1110Function CopyWorkContents(oldtype,newtype)
1111        String oldtype,newtype
1112       
1113        //check for existence of data in oldtype
1114        // if the desired workfile doesn't exist, let the user know, and abort
1115        String destPath=""
1116        if(WaveExists($("root:Packages:NIST:"+oldType + ":data")) == 0)
1117                Print "There is no work file in "+oldtype+"--Aborting"
1118                Return(1)               //error condition
1119        Endif
1120       
1121        //check for log-scaling of the "type" data and adjust if necessary
1122        ConvertFolderToLinearScale(oldtype)
1123        Fix_LogLinButtonState(0)                //make sure the button reflects the new linear scaling
1124        //then continue
1125
1126        //copy from current dir (type)=destPath to newtype, overwriting newtype contents
1127        destPath = "root:Packages:NIST:" + oldtype
1128        Duplicate/O $(destPath + ":data"),$("root:Packages:NIST:"+newtype+":data")
1129        Duplicate/O $(destPath + ":linear_data"),$("root:Packages:NIST:"+newtype+":linear_data")
1130        Duplicate/O $(destPath + ":linear_data_error"),$("root:Packages:NIST:"+newtype+":linear_data_error")
1131        Duplicate/O $(destPath + ":textread"),$("root:Packages:NIST:"+newtype+":textread")
1132        Duplicate/O $(destPath + ":integersread"),$("root:Packages:NIST:"+newtype+":integersread")
1133        Duplicate/O $(destPath + ":realsread"),$("root:Packages:NIST:"+newtype+":realsread")
1134        //
1135        // be sure to get rid of the linear_data if it exists in the destination folder
1136//      KillWaves/Z $("root:Packages:NIST:"+newtype+":linear_data")
1137
1138        //need to save a copy of filelist string too (from the current type folder)
1139        SVAR oldFileList = $(destPath + ":fileList")
1140
1141        //now switch to reference waves in newtype folder
1142        destPath = "root:Packages:NIST:"+newtype
1143        Variable/G $(destPath + ":gIsLogScale")=0                       //make new flag in newtype folder, data is linear scale
1144        //need to copy filelist string too
1145        String/G $(destPath + ":fileList") = oldFileList
1146
1147        Return(0)
1148End
1149
1150//Entry procedure from main panel
1151//
1152Proc CopyWorkFolder(oldType,newType)
1153        String oldType,newType
1154        Prompt oldType,"Source WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
1155        Prompt newType,"Destination WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
1156//      Prompt oldType,"Source WORK data type",popup,"AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1157//      Prompt newType,"Destination WORK data type",popup,"AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1158
1159        // data folder "old" will be copied to "new" (and will overwrite)
1160        CopyWorkContents(oldtype,newtype)
1161End
1162
1163//initializes datafolder and constants necessary for the workfilemath panel
1164//
1165Proc Init_WorkMath()
1166        //create the data folder
1167        //String str = "AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1168        String str = "File_1;File_2;Result;"
1169        NewDataFolder/O/S root:Packages:NIST:WorkMath
1170        String/G gFolderList=str
1171        Variable ii=0,num=itemsinlist(str)
1172        do
1173                Execute "NewDataFolder/O "+StringFromList(ii, str ,";")
1174                ii+=1
1175        while(ii<num)
1176        Variable/G const1=1,const2=1
1177       
1178        SetDataFolder root:
1179End
1180
1181//entry procedure to invoke the workfilemath panel, initializes if needed
1182//
1183Proc Show_WorkMath_Panel()
1184        DoWindow/F WorkFileMath
1185        if(V_flag==0)
1186                Init_WorkMath()
1187                WorkFileMath()
1188        Endif
1189End
1190
1191//attempts to perform the selected math operation based on the selections in the panel
1192// aborts with an error condition if operation can't be completed
1193// or puts the final answer in the Result folder, and displays the selected data
1194//
1195Function WorkMath_DoIt_ButtonProc(ctrlName) : ButtonControl
1196        String ctrlName
1197
1198        String str1,str2,oper,dest = "Result"
1199        String pathStr,workMathStr="WorkMath:"
1200       
1201        //get the panel selections (these are the names of the files on disk)
1202        PathInfo catPathName
1203        pathStr=S_Path
1204        ControlInfo popup0
1205        str1=S_Value
1206        ControlInfo popup1
1207        str2=S_Value
1208        ControlInfo popup2
1209        oper=S_Value
1210       
1211        //check that something has been selected for operation and destination
1212        if(cmpstr(oper,"Operation")==0)
1213                Abort "Select a math operand from the popup"
1214        Endif
1215
1216        //constants from globals
1217        NVAR const1=root:Packages:NIST:WorkMath:const1
1218        NVAR const2=root:Packages:NIST:WorkMath:const2
1219        Printf "(%g)%s %s (%g)%s = %s\r", const1,str1,oper,const2,str2,dest
1220        //check for proper folders (all 3 must be different)
1221       
1222        //load the data in here...
1223        //set #1
1224        Load_NamedASC_File(pathStr+str1,workMathStr+"File_1")
1225       
1226        NVAR pixelsX = root:myGlobals:gNPixelsX         //OK, location is correct
1227        NVAR pixelsY = root:myGlobals:gNPixelsY
1228       
1229        WAVE/Z data1=$("root:Packages:NIST:"+workMathStr+"File_1:linear_data")
1230        WAVE/Z err1=$("root:Packages:NIST:"+workMathStr+"File_1:linear_data_error")
1231       
1232        // set # 2
1233        If(cmpstr(str2,"UNIT MATRIX")==0)
1234                Make/O/N=(pixelsX,pixelsY) root:Packages:NIST:WorkMath:data             //don't put in File_2 folder
1235                Wave/Z data2 =  root:Packages:NIST:WorkMath:data                        //it's not real data!
1236                data2=1
1237                Duplicate/O data2 err2
1238                err2 = 0
1239        else
1240                //Load set #2
1241                Load_NamedASC_File(pathStr+str2,workMathStr+"File_2")
1242                WAVE/Z data2=$("root:Packages:NIST:"+workMathStr+"File_2:linear_data")
1243                WAVE/Z err2=$("root:Packages:NIST:"+workMathStr+"File_2:linear_data_error")
1244        Endif
1245
1246        ///////
1247       
1248        //now that we know that data exists, convert each of the operands to linear scale
1249//      ConvertFolderToLinearScale(workMathStr+"File_1")
1250//      If(cmpstr(str2,"UNIT MATRIX")!=0)
1251//              ConvertFolderToLinearScale(workMathStr+"File_2")                //don't need to convert unit matrix to linear
1252//      endif
1253
1254        //copy contents of str1 folder to dest and create the wave ref (it will exist)
1255        CopyWorkContents(workMathStr+"File_1",workMathStr+dest)
1256        WAVE/Z destData=$("root:Packages:NIST:"+workMathStr+dest+":linear_data")
1257        WAVE/Z destData_log=$("root:Packages:NIST:"+workMathStr+dest+":data")
1258        WAVE/Z destErr=$("root:Packages:NIST:"+workMathStr+dest+":linear_data_error")
1259       
1260        //dispatch
1261        strswitch(oper)
1262                case "*":               //multiplication
1263                        destData = const1*data1 * const2*data2
1264                        destErr = const1^2*const2^2*(err1^2*data2^2 + err2^2*data1^2)
1265                        destErr = sqrt(destErr)
1266                        break   
1267                case "_":               //subtraction
1268                        destData = const1*data1 - const2*data2
1269                        destErr = const1^2*err1^2 + const2^2*err2^2
1270                        destErr = sqrt(destErr)
1271                        break
1272                case "/":               //division
1273                        destData = (const1*data1) / (const2*data2)
1274                        destErr = const1^2/const2^2*(err1^2/data2^2 + err2^2*data1^2/data2^4)
1275                        destErr = sqrt(destErr)
1276                        break
1277                case "+":               //addition
1278                        destData = const1*data1 + const2*data2
1279                        destErr = const1^2*err1^2 + const2^2*err2^2
1280                        destErr = sqrt(destErr)
1281                        break                   
1282        endswitch
1283       
1284        destData_log = log(destData)            //for display
1285        //show the result
1286        WorkMath_Display_PopMenuProc("",0,"Result")
1287       
1288        PopupMenu popup4 win=WorkFileMath,mode=3                //3rd item selected == Result
1289End
1290
1291// closes the panel and kills the data folder when done
1292//
1293Function WorkMath_Done_ButtonProc(ctrlName) : ButtonControl
1294        String ctrlName
1295       
1296        DoAlert 1,"Closing the panel will kill all of the data you have loaded into memory. Do you want to continue?"
1297        If(V_Flag==2)
1298                return(0)               //don't do anything
1299        Endif
1300        //kill the panel
1301        DoWindow/K WorkFileMath
1302        //wipe out the data folder of globals
1303        SVAR dataType=root:myGlobals:gDataDisplayType
1304        if(strsearch(dataType, "WorkMath", 0 ) != -1)           //kill the SANS_Data graph if needed
1305                DoWindow/K SANS_Data
1306        Endif
1307        KillDataFolder root:Packages:NIST:WorkMath
1308End
1309
1310// loads data into the specified folder
1311//
1312// currently unused since button has been removed from panel
1313//
1314Function WorkMath_Load_ButtonProc(ctrlName) : ButtonControl
1315        String ctrlName
1316       
1317        String destStr=""
1318        SVAR folderList=root:Packages:NIST:WorkMath:gFolderList
1319        Prompt destStr,"Select the destination folder",popup,folderList
1320        DoPrompt "Folder for ASC Load",destStr
1321       
1322        if(V_flag==1)
1323                return(1)               //user abort, do nothing
1324        Endif
1325       
1326        String destFolder = "WorkMath:"+destStr
1327       
1328        Load_ASC_File("Pick the ASC file",destFolder)
1329End
1330
1331// changes the display of the SANS_Data window based on popped data type
1332// first loads in the data from the File1 or File2 popup as needed
1333// then displays the selcted dataset, if it exists
1334// makes use of procedure from DisplayUtils
1335//
1336// - Always replaces File1 or File2 with fresh data from disk
1337//
1338Function WorkMath_Display_PopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
1339        String ctrlName
1340        Variable popNum
1341        String popStr
1342       
1343        String folder="WorkMath:",pathStr,str1
1344
1345        PathInfo catPathName
1346        pathStr=S_Path
1347       
1348        //if display result, just do it and return
1349        if(cmpstr(popStr,"Result")==0)
1350                Execute "ChangeDisplay(\""+folder+popstr+"\")"
1351                return(0)
1352        endif
1353        // if file1 or file2, load in the data and display
1354        if(cmpstr(popStr,"File_1")==0)
1355                ControlInfo/W=WorkFileMath popup0
1356                str1 = S_Value
1357        Endif
1358        if(cmpstr(popStr,"File_2")==0)
1359                ControlInfo/W=WorkFileMath popup1
1360                str1 = S_Value
1361        Endif
1362        //don't load or display the unit matrix
1363        Print str1
1364        if(cmpstr(str1,"UNIT MATRIX")!=0)
1365                Load_NamedASC_File(pathStr+str1,folder+popStr)
1366                //change the display
1367                Execute "ChangeDisplay(\""+folder+popstr+"\")"
1368        endif
1369        return(0)       
1370End
1371
1372//simple panel to do workfile arithmetic
1373//
1374Proc WorkFileMath()
1375        PauseUpdate; Silent 1           // building window...
1376        NewPanel /W=(610,211,880,490)/K=2 as "Work File Math"           // replace /K=2 flag
1377        DoWindow/C WorkFileMath
1378        ModifyPanel cbRGB=(47802,54484,6682)
1379        ModifyPanel fixedSize=1
1380        SetDrawLayer UserBack
1381        DrawLine 6,166,214,166
1382        SetDrawEnv fstyle= 4
1383        DrawText 10,34,"File #1:"
1384        SetDrawEnv fstyle= 4
1385        DrawText 13,129,"File #2:"
1386        DrawText 78,186,"= Result"
1387        Button button0,pos={28,245},size={50,20},proc=WorkMath_DoIt_ButtonProc,title="Do It"
1388        Button button0,help={"Performs the specified arithmetic"}
1389        Button button1,pos={183,245},size={50,20},proc=WorkMath_Done_ButtonProc,title="Done"
1390        Button button1,help={"Closes the panel"}
1391//      Button button2,pos={30,8},size={110,20},proc=WorkMath_Load_ButtonProc,title="Load ASC Data"
1392//      Button button2,help={"Loads ASC data files into the specified folder"}
1393        Button button3,pos={205,8},size={25,20},proc=ShowWorkMathHelp,title="?"
1394        Button button3,help={"Show help file for math operations on 2-D data sets"}
1395        SetVariable setvar0,pos={9,46},size={70,15},title=" "
1396        SetVariable setvar0,limits={-Inf,Inf,0},value= root:Packages:NIST:WorkMath:const1
1397        SetVariable setvar0,help={"Multiplicative constant for the first dataset"}
1398        PopupMenu popup0,pos={89,44},size={84,20},title="X  "
1399        PopupMenu popup0,mode=1,popvalue="1st Set",value= ASC_FileList()
1400        PopupMenu popup0,help={"Selects the first dataset for the operation"}
1401        PopupMenu popup1,pos={93,136},size={89,20},title="X  "
1402        PopupMenu popup1,mode=1,popvalue="2nd Set",value= "UNIT MATRIX;"+ASC_FileList()
1403        PopupMenu popup1,help={"Selects the second dataset for the operation"}
1404        PopupMenu popup2,pos={50,82},size={70,20},title="Operator  "
1405        PopupMenu popup2,mode=3,popvalue="Operation",value= #"\"+;_;*;/;\""
1406        PopupMenu popup2,help={"Selects the mathematical operator"}
1407        SetVariable setvar1,pos={13,139},size={70,15},title=" "
1408        SetVariable setvar1,limits={-Inf,Inf,0},value= root:Packages:NIST:WorkMath:const2
1409        SetVariable setvar1,help={"Multiplicative constant for the second dataset"}
1410//      PopupMenu popup3,pos={27,167},size={124,20},title=" = Destination"
1411//      PopupMenu popup3,mode=1,popvalue="Destination",value= root:Packages:NIST:WorkMath:gFolderList
1412//      PopupMenu popup3,help={"Selects the destination folder"}
1413        PopupMenu popup4,pos={55,204},size={103,20},proc=WorkMath_Display_PopMenuProc,title="Display"
1414        PopupMenu popup4,mode=3,value= "File_1;File_2;Result;"
1415        PopupMenu popup4,help={"Displays the data in the specified folder"}
1416EndMacro
1417
1418//jump to the help topic
1419Function ShowWorkMathHelp(ctrlName) : ButtonControl
1420        String ctrlName
1421        DisplayHelpTopic/Z/K=1 "SANS Data Reduction Tutorial[2-D Work File Arithmetic]"
1422        if(V_flag !=0)
1423                DoAlert 0,"The SANS Data Reduction Tutorial Help file could not be found"
1424        endif
1425End
1426
1427//utility function to clear the contents of a data folder
1428//won't clear data that is in use -
1429//
1430Function ClearDataFolder(type)
1431        String type
1432       
1433        SetDataFolder $("root:Packages:NIST:"+type)
1434        KillWaves/a/z
1435        KillStrings/a/z
1436        KillVariables/a/z
1437        SetDataFolder root:
1438End
1439
1440
1441
1442//fileStr must be the FULL path and filename on disk
1443//destFolder is the path to the Igor folder where the data is to be deposited
1444// - "Packages:NIST:WorkMath:File_1" for example, compatible with SANS_Data display type
1445//
1446Function Load_NamedASC_File(fileStr,destFolder)
1447        String fileStr,destFolder
1448
1449        Variable refnum
1450       
1451        //read in the data
1452        ReadASCData(fileStr,destFolder)
1453
1454        //the calling macro must change the display type
1455        String/G root:myGlobals:gDataDisplayType=destFolder
1456       
1457        FillFakeHeader_ASC(destFolder)          //uses info on the panel, if available
1458
1459        //data is displayed here, and needs header info
1460       
1461        fRawWindowHook()
1462       
1463        return(0)
1464End
1465
1466
1467//function called by the main entry procedure (the load button)
1468//loads data into the specified folder, and wipes out what was there
1469//
1470//aborts if no file was selected
1471//
1472// reads the data in if all is OK
1473//
1474// currently unused, as load button has been replaced
1475//
1476Function Load_ASC_File(msgStr,destFolder)
1477        String msgStr,destFolder
1478
1479        String filename="",pathStr=""
1480        Variable refnum
1481
1482        //check for the path
1483        PathInfo catPathName
1484        If(V_Flag==1)           //      /D does not open the file
1485                Open/D/R/T="????"/M=(msgStr)/P=catPathName refNum
1486        else
1487                Open/D/R/T="????"/M=(msgStr) refNum
1488        endif
1489        filename = S_FileName           //get the filename, or null if canceled from dialog
1490        if(strlen(filename)==0)
1491                //user cancelled, abort
1492                SetDataFolder root:
1493                Abort "No file selected, action aborted"
1494        Endif
1495       
1496        //read in the data
1497        ReadASCData(filename,destFolder)
1498
1499        //the calling macro must change the display type
1500        String/G root:myGlobals:gDataDisplayType=destFolder
1501       
1502        FillFakeHeader_ASC(destFolder)          //uses info on the panel, if available
1503
1504        //data is displayed here, and needs header info
1505       
1506        fRawWindowHook()
1507       
1508        return(0)
1509End
1510
1511
1512
1513//function called by the popups to get a file list of data that can be sorted
1514// this procedure simply removes the raw data files from the string - there
1515//can be lots of other junk present, but this is very fast...
1516//
1517// -- simplified to do in a single call -- extension must be .ASC
1518// - so far, only used in WorkFileMath popups, which require .ASC
1519//
1520Function/S ASC_FileList()
1521
1522        String list="",newList="",item=""
1523        Variable num,ii
1524       
1525        //check for the path
1526        PathInfo catPathName
1527        if(V_Flag==0)
1528                DoAlert 0, "Data path does not exist - pick the data path from the button on the main panel"
1529                Return("")
1530        Endif
1531
1532        list = IndexedFile(catpathName,-1,".ASC")
1533
1534
1535//      list = IndexedFile(catpathName,-1,"????")
1536//     
1537//      list = RemoveFromList(ListMatch(list,"*.SA1*",";"), list, ";", 0)
1538//      list = RemoveFromList(ListMatch(list,"*.SA2*",";"), list, ";", 0)
1539//      list = RemoveFromList(ListMatch(list,"*.SA3*",";"), list, ";", 0)
1540//      list = RemoveFromList(ListMatch(list,".*",";"), list, ";", 0)
1541//      list = RemoveFromList(ListMatch(list,"*.pxp",";"), list, ";", 0)
1542//      list = RemoveFromList(ListMatch(list,"*.DIV",";"), list, ";", 0)
1543//      list = RemoveFromList(ListMatch(list,"*.GSP",";"), list, ";", 0)
1544//      list = RemoveFromList(ListMatch(list,"*.MASK",";"), list, ";", 0)
1545//
1546//      //remove VAX version numbers
1547//      list = RemoveVersNumsFromList(List)
1548        //sort
1549        newList = SortList(List,";",0)
1550
1551        return newlist
1552End
1553
1554//
1555// tests if two values are close enough to each other
1556// very useful since ICE came to be
1557//
1558// tol is an absolute value (since input v1 or v2 may be zero, can't reliably
1559// use a percentage
1560Function CloseEnough(v1,v2,tol)
1561        Variable v1, v2, tol
1562
1563        if(abs(v1-v2) < tol)
1564                return(1)
1565        else
1566                return(0)
1567        endif
1568End
1569
1570//
1571//
1572// match the attenuation of the RAW data to the "type" data
1573// so that they can be properly added
1574//
1575// are the attenuator numbers the same? if so exit
1576//
1577// if not, find the attenuator number for type
1578// - find both attenuation factors
1579//
1580// rescale the raw data to match the ratio of the two attenuation factors
1581// -- adjust the detector count (rw)
1582// -- the linear data
1583//
1584//
1585Function Adjust_RAW_Attenuation(type)
1586        String type
1587       
1588        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1589        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1590        WAVE data=$("root:Packages:NIST:RAW:data")
1591        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1592        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1593       
1594        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1595
1596        Variable dest_atten,raw_atten,tol
1597        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1598        String fileStr
1599
1600        dest_atten = dest_reals[3]
1601        raw_atten = rw[3]
1602       
1603        tol = 0.1               // within 0.1 atten units is OK
1604        if(abs(dest_atten - raw_atten) < tol )
1605                return(0)
1606        endif
1607
1608        fileStr = tw[3]
1609        lambda = rw[26]
1610        raw_AttenFactor = AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1611        dest_AttenFactor = AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1612               
1613        rw[2] *= dest_AttenFactor/raw_AttenFactor
1614        linear_data *= dest_AttenFactor/raw_AttenFactor
1615       
1616        // to keep "data" and linear_data in sync
1617        data = linear_data
1618       
1619        return(0)
1620End
Note: See TracBrowser for help on using the repository browser.