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

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

more changes and additons to display VSANS data

adding functions for IvsQ plotting

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

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

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

File size: 27.1 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma version=1.0
3#pragma IgorVersion=6.1
4
5
6
7//
8// functions to apply corrections to the detector panels
9//
10// these are meant to be called by the procedures that convert "raw" data to
11// "adjusted" or corrected data sets
12//
13// may be relocated in the future
14//
15
16
17
18//
19// detector dead time
20//
21// input is the data array (N tubes x M pixels)
22// input of N x 1 array of dead time values
23//
24// output is the corrected counts in data, overwriting the input data
25//
26// Note that the equation in Roe (eqn 2.15, p. 63) looks different, but it is really the
27// same old equation, just written in a more complex form.
28//
29// TODO
30// -- verify the direction of the tubes and indexing
31// x- decide on the appropriate functional form for the tubes
32// x- need count time as input
33// -- be sure I'm working in the right data folder
34// -- clean up when done
35// -- calculate + return the error contribution?
36//
37Function DeadTimeCorrectionTubes(dataW,data_errW,dtW,ctTime)
38        Wave dataW,data_errW,dtW
39        Variable ctTime
40       
41        // do I count on the orientation as an input, or do I just figure it out on my own?
42        String orientation
43        Variable dimX,dimY
44        dimX = DimSize(dataW,0)
45        dimY = DimSize(dataw,1)
46        if(dimX > dimY)
47                orientation = "horizontal"
48        else
49                orientation = "vertical"
50        endif
51       
52        // sum the counts in each tube and divide by time for total cr per tube
53        Variable npt
54       
55        if(cmpstr(orientation,"vertical")==0)
56                //      this is data dimensioned as (Ntubes,Npix)
57               
58                MatrixOp/O sumTubes = sumRows(dataW)            // n x 1 result
59                sumTubes /= ctTime              //now count rate per tube
60               
61                dataW[][] = dataW[p][q]/(1-sumTubes[p][0]*dtW[p])               //correct the data
62
63        elseif(cmpstr(orientation,"horizontal")==0)
64        //      this is data (horizontal) dimensioned as (Npix,Ntubes)
65
66                MatrixOp/O sumTubes = sumCols(dataW)            // 1 x m result
67                sumTubes /= ctTime
68               
69                dataW[][] = dataW[p][q]/(1-sumTubes[0][q]*dtW[q])
70       
71        else           
72                DoAlert 0,"Orientation not correctly passed in DeadTimeCorrectionTubes(). No correction done."
73        endif
74       
75        return(0)
76end
77
78//
79// Non-linear data correction
80//
81// input is the data array (N tubes x M pixels)
82// input of N x M array of quadratic coefficients
83//
84// output is wave of corrected real space distance corresponding to each pixel of the data
85//
86//
87// TODO
88// -- verify the direction of the tubes and indexing
89// -- be sure I'm working in the right data folder
90// -- clean up when done
91// -- calculate + return the error contribution?
92// -- do I want this to return a wave?
93// -- do I need to write a separate function that returns the distance wave for later calculations?
94// -- do I want to make the distance array 3D to keep the x and y dims together? Calculate them all right now?
95// -- what else do I need to pass to the function? (fname=folder? detStr?)
96//
97//
98//
99Function NonLinearCorrection(dataW,coefW,tube_width,detStr,destPath)
100        Wave dataW,coefW
101        Variable tube_width
102        String detStr,destPath
103       
104         
105        // do I count on the orientation as an input, or do I just figure it out on my own?
106        String orientation
107        Variable dimX,dimY
108        dimX = DimSize(dataW,0)
109        dimY = DimSize(dataw,1)
110        if(dimX > dimY)
111                orientation = "horizontal"
112        else
113                orientation = "vertical"
114        endif
115
116        // make a wave of the same dimensions, in the same data folder for the distance
117        // ?? or a 3D wave?
118        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistX")
119        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistY")
120        Wave data_realDistX = $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistX")
121        Wave data_realDistY = $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistY")
122       
123        // then per tube, do the quadratic calculation to get the real space distance along the tube
124        // the distance perpendicular to the tube is n*(8.4mm) per tube index
125       
126        if(cmpstr(orientation,"vertical")==0)
127                //      this is data dimensioned as (Ntubes,Npix)
128                data_realDistY[][] = coefW[0][p] + coefW[1][p]*q + coefW[2][p]*q*q
129
130                data_realDistX[][] = tube_width*p
131               
132        elseif(cmpstr(orientation,"horizontal")==0)
133                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
134                data_realDistX[][] = coefW[0][q] + coefW[1][q]*p + coefW[2][q]*p*p
135       
136                data_realDistY[][] = tube_width*q
137
138        else           
139                DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
140        endif
141       
142        return(0)
143end
144
145
146/////
147//
148// non-linear corrections to the tube pixels
149// - returns the distance in mm (although this may change)
150//
151// c0,c1,c2,pix
152// c0-c2 are the fit coefficients
153// pix is the test pixel
154//
155// returns the distance in mm (relative to ctr pixel)
156// ctr is the center pixel, as defined when fitting to quadratic was done
157//
158Function V_TubePixel_to_mm(c0,c1,c2,pix)
159        Variable c0,c1,c2,pix
160       
161        Variable dist
162        dist = c0 + c1*pix + c2*pix*pix
163       
164        return(dist)
165End
166//
167////
168
169// TODO
170// get rid of this in the real data
171//
172// TESTING ONLY
173Macro MakeFakeCalibrationWaves()
174        // make these in the RAW data folder, before converting to a work folder
175        // - then they will be "found" by get()
176        // -- only for the tube, not the Back det
177       
178        DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
179       
180        fMakeFakeCalibrationWaves()
181End
182
183
184
185// TODO
186// get rid of this in the real data
187//
188// TESTING ONLY
189//
190// orientation does not matter, there are 48 tubes in each bank
191// so dimension (3,48) for everything.
192//
193// -- but the orientation does indicate TB vs LR, which has implications for
194//  the (fictional) dimension of the pixel along the tube axis, at least as far
195// as for making the fake coefficients.
196//
197Function fMakeFakeCalibrationWaves()
198
199        Variable ii,pixSize
200        String detStr,fname="RAW",orientation
201       
202        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
203                detStr = StringFromList(ii, ksDetectorListNoB, ";")
204//              Wave w = V_getDetectorDataW(fname,detStr)
205                Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:entry:instrument:detector_"+detStr+":spatial_calibration")
206                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:entry:instrument:detector_"+detStr+":spatial_calibration")
207                // !!!! this overwrites what is there
208
209                orientation = V_getDet_tubeOrientation(fname,detStr)
210                if(cmpstr(orientation,"vertical")==0)
211                //      this is vertical tube data dimensioned as (Ntubes,Npix)
212                        pixSize = 8                     //V_getDet_y_pixel_size(fname,detStr)
213                       
214                elseif(cmpstr(orientation,"horizontal")==0)
215                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
216                        pixSize = 4                     //V_getDet_x_pixel_size(fname,detStr)
217                       
218                else           
219                        DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
220                endif
221               
222                calib[0][] = -(128/2)*pixSize                   //approx (n/2)*pixSixe
223                calib[1][] = pixSize
224                calib[2][] = 2e-4
225               
226        endfor
227       
228        return(0)
229End
230
231
232
233Function V_Detector_CalcQVals(fname,detStr,destPath)
234        String fname,detStr,destPath
235
236        String orientation
237        Variable xCtr,yCtr,lambda,sdd
238       
239// get all of the geometry information 
240        orientation = V_getDet_tubeOrientation(fname,detStr)
241        sdd = V_getDet_distance(fname,detStr)
242        xCtr = V_getDet_beam_center_x(fname,detStr)
243        yCtr = V_getDet_beam_center_y(fname,detStr)
244        lambda = V_getWavelength(fname)
245        Wave data_realDistX = $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistX")
246        Wave data_realDistY = $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistY")
247
248// make the new waves
249        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qTot_"+detStr)
250        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qx_"+detStr)
251        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qy_"+detStr)
252        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qz_"+detStr)
253        Wave qTot = $(destPath + ":entry:entry:instrument:detector_"+detStr+":qTot_"+detStr)
254        Wave qx = $(destPath + ":entry:entry:instrument:detector_"+detStr+":qx_"+detStr)
255        Wave qy = $(destPath + ":entry:entry:instrument:detector_"+detStr+":qy_"+detStr)
256        Wave qz = $(destPath + ":entry:entry:instrument:detector_"+detStr+":qz_"+detStr)
257
258// calculate all of the q-values
259        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
260        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
261        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
262        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
263       
264       
265        return(0)
266End
267
268
269//function to calculate the overall q-value, given all of the necesary trig inputs
270//
271// TODO:
272// -- verify the calculation (accuracy - in all input conditions)
273// -- verify the units of everything here, it's currently all jumbled an wrong...
274// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
275//    to each pixel
276//
277//sdd is in meters
278//wavelength is in Angstroms
279//
280//returned magnitude of Q is in 1/Angstroms
281//
282Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
283        Variable xaxval,yaxval,xctr,yctr,sdd,lam
284        Wave distX,distY
285       
286        Variable dx,dy,qval,two_theta,dist
287               
288        sdd *=100               //convert to cm
289        dx = (distX[xaxval][yaxval] - xctr)             //delta x in cm
290        dy = (distY[xaxval][yaxval] - yctr)             //delta y in cm
291        dist = sqrt(dx^2 + dy^2)
292       
293        two_theta = atan(dist/sdd)
294
295        qval = 4*Pi/lam*sin(two_theta/2)
296       
297        return qval
298End
299
300//calculates just the q-value in the x-direction on the detector
301// TODO:
302// -- verify the calculation (accuracy - in all input conditions)
303// -- verify the units of everything here, it's currently all jumbled an wrong...
304// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
305//    to each pixel
306//
307//
308// this properly accounts for qz
309//
310Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
311        Variable xaxval,yaxval,xctr,yctr,sdd,lam
312        Wave distX,distY
313
314        Variable qx,qval,phi,dx,dy,dist,two_theta
315       
316        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
317       
318        sdd *=100               //convert to cm
319        dx = (distX[xaxval][yaxval] - xctr)             //delta x in cm
320        dy = (distY[xaxval][yaxval] - yctr)             //delta y in cm
321        phi = V_FindPhi(dx,dy)
322       
323        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
324        dist = sqrt(dx^2 + dy^2)
325        two_theta = atan(dist/sdd)
326
327        qx = qval*cos(two_theta/2)*cos(phi)
328       
329        return qx
330End
331
332//calculates just the q-value in the y-direction on the detector
333// TODO:
334// -- verify the calculation (accuracy - in all input conditions)
335// -- verify the units of everything here, it's currently all jumbled an wrong...
336// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
337//    to each pixel
338//
339//
340// this properly accounts for qz
341//
342Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
343        Variable xaxval,yaxval,xctr,yctr,sdd,lam
344        Wave distX,distY
345
346        Variable qy,qval,phi,dx,dy,dist,two_theta
347       
348        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
349       
350        sdd *=100               //convert to cm
351        dx = (distX[xaxval][yaxval] - xctr)             //delta x in cm
352        dy = (distY[xaxval][yaxval] - yctr)             //delta y in cm
353        phi = V_FindPhi(dx,dy)
354       
355        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
356        dist = sqrt(dx^2 + dy^2)
357        two_theta = atan(dist/sdd)
358
359        qy = qval*cos(two_theta/2)*sin(phi)
360       
361        return qy
362End
363
364//calculates just the q-value in the z-direction on the detector
365// TODO:
366// -- verify the calculation (accuracy - in all input conditions)
367// -- verify the units of everything here, it's currently all jumbled an wrong...
368// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
369//    to each pixel
370//
371// not actually used for anything, but here for completeness if anyone asks
372//
373// this properly accounts for qz
374//
375Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
376        Variable xaxval,yaxval,xctr,yctr,sdd,lam
377        Wave distX,distY
378
379        Variable qz,qval,phi,dx,dy,dist,two_theta
380       
381        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
382       
383        sdd *=100               //convert to cm
384        dx = (distX[xaxval][yaxval] - xctr)             //delta x in cm
385        dy = (distY[xaxval][yaxval] - yctr)             //delta y in cm
386       
387        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
388        dist = sqrt(dx^2 + dy^2)
389        two_theta = atan(dist/sdd)
390
391        qz = qval*sin(two_theta/2)
392       
393        return qz
394End
395
396
397
398////////////
399// TODO: all of below is untested code
400//   copied from SANS
401//
402//
403// TODO :
404//   -- DoAlert 0,"This has not yet been updated for VSANS"
405//
406//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
407//function is called by Raw_to_work() and Add_raw_to_work() functions
408//works on the actual data array, assumes that is is already on LINEAR scale
409//
410Function DetCorr(data,data_err,realsread,doEfficiency,doTrans)
411        Wave data,data_err,realsread
412        Variable doEfficiency,doTrans
413
414        DoAlert 0,"This has not yet been updated for VSANS"
415       
416        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
417        Variable ii,jj,dtdist,dtdis2
418        Variable xi,xd,yd,rad,ratio,domega,xy
419        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
420       
421//      Print "...doing jacobian and non-linear corrections"
422
423        NVAR pixelsX = root:myGlobals:gNPixelsX
424        NVAR pixelsY = root:myGlobals:gNPixelsY
425       
426        //set up values to send to auxiliary trig functions
427        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
428        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
429
430        x0 = realsread[16]
431        y0 = realsread[17]
432        sx = realsread[10]
433        sx3 = realsread[11]
434        sy = realsread[13]
435        sy3 = realsread[14]
436       
437        dtdist = 1000*realsread[18]     //sdd in mm
438        dtdis2 = dtdist^2
439       
440        lambda = realsRead[26]
441        trans = RealsRead[4]
442        trans_err = RealsRead[41]               //new, March 2011
443       
444
445        //waves to contain repeated function calls
446        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
447        ii=0
448        do
449                xi = ii
450//              fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
451//              xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
452//              yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
453                ii+=1
454        while(ii<pixelsX)
455       
456        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
457       
458        ii=0
459        do
460                xi = ii
461//              xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
462                jj=0
463                do
464                        yd = fyy[jj]-yy0
465                        //rad is the distance of pixel ij from the sample
466                        //domega is the ratio of the solid angle of pixel ij versus center pixel
467                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
468                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
469                        rad = sqrt(dtdis2 + xd^2 + yd^2)
470                        domega = rad/dtdist
471                        ratio = domega^3
472                        xy = xx[ii]*yy[jj]
473                       
474                        data[ii][jj] *= xy*ratio
475                       
476                        solidAngle[ii][jj] = xy*ratio           //testing only 
477                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
478                       
479                       
480                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
481                        // correction inserted 11/2007 SRK
482                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
483                        // so divide here to get the correct answer (5/22/08 SRK)
484                        if(doEfficiency)
485                                data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
486                                data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
487//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
488                        endif
489                       
490                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
491                        // so divide here to get the correct answer
492                        if(doTrans)
493                       
494                                if(trans<0.1 && ii==0 && jj==0)
495                                        Print "***transmission is less than 0.1*** and is a significant correction"
496                                endif
497                               
498                                if(trans==0)
499                                        if(ii==0 && jj==0)
500                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
501                                        endif
502                                        trans = 1
503                                endif
504                               
505                                // pass in the transmission error, and the error in the correction is returned as the last parameter
506                                lat_corr = LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)             //moved from 1D avg SRK 11/2007
507                                data[ii][jj] /= lat_corr                        //divide by the correction factor
508                                //
509                                //
510                                //
511                                // relative errors add in quadrature
512                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
513                                tmp_err = sqrt(tmp_err)
514                               
515                                data_err[ii][jj] = tmp_err
516                               
517//                              solidAngle[ii][jj] = lat_err
518
519                               
520                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
521                        endif
522                       
523                        jj+=1
524                while(jj<pixelsX)
525                ii+=1
526        while(ii<pixelsX)
527       
528        //clean up waves
529       
530        Return(0)
531End
532
533
534
535//distances passed in are in mm
536// dtdist is SDD
537// xd and yd are distances from the beam center to the current pixel
538//
539// TODO:
540//   --         DoAlert 0,"This has not yet been updated for VSANS"
541//
542Function DetEffCorr(lambda,dtdist,xd,yd)
543        Variable lambda,dtdist,xd,yd
544
545        DoAlert 0,"This has not yet been updated for VSANS"
546       
547        Variable theta,cosT,ff,stAl,stHe
548       
549        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )
550        cosT = cos(theta)
551       
552        stAl = 0.00967*lambda*0.8               //dimensionless, constants from JGB memo
553        stHe = 0.146*lambda*2.5
554       
555        ff = exp(-stAl/cosT)*(1-exp(-stHe/cosT)) / ( exp(-stAl)*(1-exp(-stHe)) )
556               
557        return(ff)
558End
559
560// DIVIDE the intensity by this correction to get the right answer
561// TODO:
562//   --         DoAlert 0,"This has not yet been updated for VSANS"
563//
564//
565Function LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err)
566        Variable trans,dtdist,xd,yd,trans_err,&err
567
568        DoAlert 0,"This has not yet been updated for VSANS"
569       
570        //angle dependent transmission correction
571        Variable uval,arg,cos_th,correction,theta
572       
573        ////this section is the trans_correct() VAX routine
574//      if(trans<0.1)
575//              Print "***transmission is less than 0.1*** and is a significant correction"
576//      endif
577//      if(trans==0)
578//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
579//              trans = 1
580//      endif
581       
582        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )              //theta at the input pixel
583       
584        //optical thickness
585        uval = -ln(trans)               //use natural logarithm
586        cos_th = cos(theta)
587        arg = (1-cos_th)/cos_th
588       
589        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
590        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
591        // OR
592        if((uval<0.01) || (cos_th>0.99))       
593                //small arg, approx correction
594                correction= 1-0.5*uval*arg
595        else
596                //large arg, exact correction
597                correction = (1-exp(-uval*arg))/(uval*arg)
598        endif
599
600        Variable tmp
601       
602        if(trans == 1)
603                err = 0         //no correction, no error
604        else
605                //sigT, calculated from the Taylor expansion
606                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
607                tmp *= tmp
608                tmp *= trans_err^2
609                tmp = sqrt(tmp)         //sigT
610               
611                err = tmp
612        endif
613       
614//      Printf "trans error = %g\r",trans_err
615//      Printf "correction = %g +/- %g\r", correction, err
616       
617        //end of transmission/pathlength correction
618
619        return(correction)
620end
621
622
623//
624// TODO:
625//   --         DoAlert 0,"This has not yet been updated for VSANS"
626//
627//test procedure, not called anymore
628Proc AbsoluteScaling(type,c0,c1,c2,c3,c4,c5)
629        String type
630        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0
631        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
632        Prompt c0, "Sample Transmission"
633        Prompt c1, "Sample Thickness (cm)"
634        Prompt c2, "Standard Transmission"
635        Prompt c3, "Standard Thickness (cm)"
636        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
637        Prompt c5, "Standard Cross-Section (cm-1)"
638
639        Variable err
640        //call the function to do the math
641        //data from "type" will be scaled and deposited in ABS
642        err = Absolute_Scale(type,c0,c1,c2,c3,c4,c5)
643       
644        if(err)
645                Abort "Error in Absolute_Scale()"
646        endif
647       
648        //contents are always dumped to ABS
649        type = "ABS"
650       
651        String newTitle = "WORK_"+type
652        DoWindow/F SANS_Data
653        DoWindow/T SANS_Data, newTitle
654        KillStrings/Z newTitle
655       
656        //need to update the display with "data" from the correct dataFolder
657        //reset the current displaytype to "type"
658        String/G root:myGlobals:gDataDisplayType=Type
659       
660        fRawWindowHook()
661       
662End
663
664//
665// TODO:
666//   --         DoAlert 0,"This has not yet been updated for VSANS"
667//
668//s_ is the standard
669//w_ is the "work" file
670//both are work files and should already be normalized to 10^8 monitor counts
671Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
672        String type
673        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
674
675        DoAlert 0,"This has not yet been updated for VSANS"
676               
677        //convert the "type" data to absolute scale using the given standard information
678        //copying the "type" waves to ABS
679       
680        //check for existence of data, rescale to linear if needed
681        String destPath
682        //check for "type"
683        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
684                Print "There is no work file in "+type+"--Aborting"
685                Return(1)               //error condition
686        Endif
687        //check for log-scaling of the "type" data and adjust if necessary
688        destPath = "root:Packages:NIST:"+Type
689        NVAR gIsLogScale = $(destPath + ":gIsLogScale")
690        if(gIsLogScale)
691                Duplicate/O $(destPath + ":linear_data") $(destPath + ":data")//back to linear scale
692                Variable/G $(destPath + ":gIsLogScale")=0       //the "type" data is not logscale anymore
693        endif
694       
695        //copy "oldtype" information to ABS
696        //overwriting out the old contents of the ABS folder (/O option in Duplicate)
697        //copy over the waves data,vlegend,text,integers,reals(read)
698
699        String oldType= "root:Packages:NIST:"+type              //this is where the data to be absoluted is
700        //copy from current dir (type) to ABS, defined by destPath
701        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data"
702        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data"
703        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error"
704//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend"
705        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread"
706        Duplicate/O $(oldType + ":integersread"),$"root:Packages:NIST:ABS:integersread"
707        Duplicate/O $(oldType + ":realsread"),$"root:Packages:NIST:ABS:realsread"
708        //need to save a copy of filelist string too (from the current type folder)
709        SVAR oldFileList = $(oldType + ":fileList")
710        //need to copy filelist string too
711        String/G $"root:Packages:NIST:ABS:fileList" = oldFileList
712       
713        //now switch to ABS folder
714        //make appropriate wave references
715        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS
716        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS
717        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display
718        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on
719        WAVE integersread=$"root:Packages:NIST:ABS:integersread"
720        WAVE realsread=$"root:Packages:NIST:ABS:realsread"
721        Variable/G $"root:Packages:NIST:ABS:gIsLogscale"=0                      //make new flag in ABS folder, data is linear scale
722       
723        //do the actual absolute scaling here, modifying the data in ABS
724        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
725       
726        w_moncount = realsread[0]               //monitor count in "type"
727        if(w_moncount == 0)
728                //zero monitor counts will give divide by zero ---
729                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
730                Return(1)               //report error
731        Endif
732       
733        //calculate scale factor
734        Variable scale,trans_err
735        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1)
736        s2 = s_thick/w_thick
737        s3 = s_trans/w_trans
738        s4 = s_cross/s_izero
739       
740        // kappa comes in as s_izero, so be sure to use 1/kappa_err
741       
742        data *= s1*s2*s3*s4
743       
744        scale = s1*s2*s3*s4
745        trans_err = realsRead[41]
746       
747//      print scale
748//      print data[0][0]
749       
750        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
751
752//      print data_err[0][0]
753       
754// keep "data" in sync with linear_data
755        data_copy = data
756       
757        //********* 15APR02
758        // DO NOt correct for atenuators here - the COR step already does this, putting all of the data one equal
759        // footing (zero atten) before doing the subtraction.
760        //
761        //Print "ABS data multiplied by  ",s1*s2*s3*s4/attenFactor
762       
763        //update the ABS header information
764        textread[1] = date() + " " + time()             //date + time stamp
765       
766        Return (0) //no error
767End
768
769
770//
771// TODO:
772//   --         DoAlert 0,"This has not yet been updated for VSANS"
773//
774//
775// match the attenuation of the RAW data to the "type" data
776// so that they can be properly added
777//
778// are the attenuator numbers the same? if so exit
779//
780// if not, find the attenuator number for type
781// - find both attenuation factors
782//
783// rescale the raw data to match the ratio of the two attenuation factors
784// -- adjust the detector count (rw)
785// -- the linear data
786//
787//
788Function Adjust_RAW_Attenuation(type)
789        String type
790
791        DoAlert 0,"This has not yet been updated for VSANS"
792       
793        WAVE rw=$("root:Packages:NIST:RAW:realsread")
794        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
795        WAVE data=$("root:Packages:NIST:RAW:data")
796        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
797        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
798       
799        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
800
801        Variable dest_atten,raw_atten,tol
802        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
803        String fileStr
804
805        dest_atten = dest_reals[3]
806        raw_atten = rw[3]
807       
808        tol = 0.1               // within 0.1 atten units is OK
809        if(abs(dest_atten - raw_atten) < tol )
810                return(0)
811        endif
812
813        fileStr = tw[3]
814        lambda = rw[26]
815        // TODO access correct values
816        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
817        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
818               
819        rw[2] *= dest_AttenFactor/raw_AttenFactor
820        linear_data *= dest_AttenFactor/raw_AttenFactor
821       
822        // to keep "data" and linear_data in sync
823        data = linear_data
824       
825        return(0)
826End
827
828//
829// TODO:
830//   --         DoAlert 0,"This has not yet been updated for VSANS"
831//
832//************************
833//unused testing procedure, may not be up-to-date with other procedures
834//check before re-implementing
835//
836Macro DIV_a_Workfile(type)
837        String type
838        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
839       
840        //macro will take whatever is in SELECTED folder and DIVide it by the current
841        //contents of the DIV folder - the function will check for existence
842        //before proceeding
843
844//DoAlert 0,"This has not yet been updated for VSANS"
845       
846        Variable err
847        err = DIVCorrection(type)               //returns err = 1 if data doesn't exist in specified folders
848       
849        if(err)
850                Abort "error in DIVCorrection()"
851        endif
852       
853        //contents are NOT always dumped to CAL, but are in the new type folder
854       
855        String newTitle = "WORK_"+type
856        DoWindow/F VSANS_Data
857        DoWindow/T VSANS_Data, newTitle
858        KillStrings/Z newTitle
859       
860        //need to update the display with "data" from the correct dataFolder
861        //reset the current displaytype to "type"
862        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
863       
864        UpdateDisplayInformation(type)
865       
866End
867
868
869//
870// TODO:
871//   --         DoAlert 0,"This has not yet been updated for VSANS"
872//
873//function will divide the contents of "workType" folder with the contents of
874//the DIV folder + detStr
875// all data is linear scale for the calculation
876//
877Function DIVCorrection(data,data_err,detStr,workType)
878        Wave data,data_err
879        String detStr,workType
880       
881        //check for existence of data in type and DIV
882        // if the desired data doesn't exist, let the user know, and abort
883        String destPath=""
884
885        if(WaveExists(data) == 0)
886                Print "The data wave does not exist in DIVCorrection()"
887                Return(1)               //error condition
888        Endif
889       
890        //check for DIV
891        // if the DIV workfile doesn't exist, let the user know,and abort
892
893        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:entry:instrument:detector_"+detStr+":data")
894        if(WaveExists(div_data) == 0)
895                Print "The DIV wave does not exist in DIVCorrection()"
896                Return(1)               //error condition
897        Endif
898        //files exist, proceed
899
900        data /= div_data
901       
902//      data_err /= div_data
903       
904        Return(0)
905End
906
907
908//////////////////////////
Note: See TracBrowser for help on using the repository browser.