source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Utils.ipf @ 1242

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

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

File size: 47.6 KB
RevLine 
[954]1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
[1242]2#pragma IgorVersion = 7.00
[954]3
[955]4/////////////////////////
5//
6// Utility functions to:
7//              calculate Q, Qx, Qy, Qz
8//              fill the detector panels with simulated data (the model functions are here)
9//              bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width)
10//
11/////////////////////////
[954]12
[955]13
14
15
[1051]16// x- hard wired for a sphere - change this to allow minimal selections and altering of coefficients
17// x- add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
[982]18//
19// NOTE - this is a VCALC only routine, so it's not been made completely generic
20//
[954]21Function FillPanel_wModelData(det,qTot,type)
22        Wave det,qTot
23        String type
24
[982]25//      SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
[954]26
27        // q-values and detector arrays already allocated and calculated
28        Duplicate/O det tmpInten,tmpSig,prob_i
29               
30        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
31
[982]32        //imon = VC_BeamIntensity()*CountTime
[954]33        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
34        trans = 0.8
35        thick = 0.1
36       
37        // need SDD
38        // need pixel dimensions
[1062]39        // nominal sdd in cm, offset in cm, want result in cm !
[954]40
[1128]41        sdd = VC_getSDD(type)   // setback is already included  VCALC_getTopBottomSDDSetback(type)              // result is sdd in [cm]
42
[954]43        pixSizeX = VCALC_getPixSizeX(type)              // cm
44        pixSizeY = VCALC_getPixSizeY(type)
45       
46       
47        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
48        //  keep the panel to keep people used to using it.)
[955]49        // peak @ 0.1 ~ AgBeh
50        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
51        //
52        // peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
[1022]53        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}     
54        //
55        //
56        Variable addEmpBgd=0
57       
[1117]58// TODOHIGHRES
59// this is a slow step - try to figure out how to multithread this efficiently. simply adding the
60//  keyword has little effect. maybe only do this for "B", maybe rewrite the calculation to not use pq indexing
61//             
[954]62        String funcStr = VCALC_getModelFunctionStr()
63        strswitch(funcStr)
64                case "Big Debye":
[987]65                        tmpInten = VC_Debye(100,3000,0.0001,qTot[p][q])
[954]66                        break
67                case "Big Sphere":
[1100]68                        tmpInten = VC_SphereForm(1,2000,1e-6,0.01,qTot[p][q])   
[954]69                        break
70                case "Debye":
[1117]71                        MultiThread tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
[954]72                        break
73                case "Sphere":
[982]74                        tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])   
[954]75                        break
[955]76                case "AgBeh":
[1023]77                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
[955]78                        break
79                case "Vycor":
[982]80                        tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
[955]81                        break   
82                case "Empty Cell":
[1023]83                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
[955]84                        break
85                case "Blocked Beam":
[1023]86                        tmpInten = VC_BlockedBeam(0.01,qTot[p][q])
[955]87                        break
[1022]88                case "Debye +":
89                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
90                        addEmpBgd = 1
91                        break
92                case "AgBeh +":
[1023]93                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
[1022]94                        addEmpBgd = 1
95                        break
96                case "Empty Cell +":
[1023]97                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
98                        tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
[1022]99                        break
[954]100                default:
[982]101                        tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
[954]102        endswitch
103
[1022]104
105        if(addEmpBgd == 1)
[1023]106                tmpInten += VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
107                tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
[1022]108        endif
109
110       
[1051]111// x- this is faked to get around the singularity at the center of the back detector
[987]112//
113//
114        if(cmpstr(type,"B") == 0)
115                Variable nx,ny,px,py
116                nx = VCALC_get_nPix_X(type)
117                ny = VCALC_get_nPix_Y(type)
118                px = trunc(nx/2)
119                py = trunc(ny/2)
120               
121                tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
122        endif
[954]123
[987]124
125
[954]126///////////////
127//      // calculate the scattering cross section simply to be able to estimate the transmission
128//      Variable sig_sas=0
129//     
130//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
131//      // subtracted before the calculation.
132//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
133//
134//      if(sig_sas > 100)
135//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
136//      endif           
137//
138//      // calculate the multiple scattering fraction for display (10/2009)
139//      Variable ii,nMax=10,tau
140//      mScat=0
141//      tau = thick*sig_sas
142//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
143//      // neutrons out of those that were scattered
144//      for(ii=2;ii<nMax;ii+=1)
145//              mScat += tau^(ii)/factorial(ii)
146////            print tau^(ii)/factorial(ii)
147//      endfor
148//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
149//      mscat *= (estTrans)/(1-estTrans)
150//
151////    if(mScat > 0.1)         //  Display warning
152//
153//      Print "Sig_sas = ",sig_sas
154////////////////////
155       
[1137]156        MultiThread prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                     //probability of a neutron in q-bin(i)
[954]157               
158        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
159
160
161/// **** can I safely assume a Gaussian error in the count rate??
[1137]162        MultiThread tmpSig = sqrt(tmpInten)             // corrected based on John's memo, from 8/9/99
[954]163
[1137]164        MultiThread tmpInten += gnoise(tmpSig)
165        MultiThread tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                        // MAR 2013 -- is this the right thing to do
[954]166        tmpInten = trunc(tmpInten)
167               
168       
169        det = tmpInten
170
171// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
172//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
173
174       
175        KillWaves/Z tmpInten,tmpSig,prob_i     
176        SetDataFolder root:
177
178        return(0)
179End
180
181
182// For a given detector panel, calculate the q-values
183// -work with everything as arrays
184// Input needed:
185// detector data
186// detector type (LRTB?)
187// beam center (may be off the detector)
188// SDD
189// lambda
190//
191// pixel dimensions for detector type (global constants)
192// - data dimensions read directly from array
193//
194// --What is calculated:
195// array of Q
196// array of qx,qy,qz
197// array of error already exists
198//
199//
[1062]200// -- sdd in cm
[954]201// -- lambda in Angstroms
[982]202Function VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
[954]203        Wave data,qTot,qx,qy,qz
204        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
205               
206        // loop over the array and calculate the values - this is done as a wave assignment
207// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
[1137]208        MultiThread qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
209        MultiThread     qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
210        MultiThread     qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
211        MultiThread     qz = VC_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]212       
213        return(0)
214End
215
216
[993]217// for testing, a version that will calculate the q-arrays for VCALC based on whatever nonlinear coefficients
218// exist in the RAW data folder
219//
220// reverts to the "regular" linear detector if waves not found or a flag is set
221//
222// need to call the VSANS V_CalcQval routines (these use the real-space distance, not pixel dims)
223//
[1062]224// ***** everything passed in is [cm], except for wavelength [A]
225//
226// ****  TODO :: calibration constants are still in [mm]
227//
228//
[993]229// TODO:
230// -- tube width is hard-wired in
231//
232//
233Function VC_Detector_2Q_NonLin(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY,detStr)
234        Wave data,qTot,qx,qy,qz
235        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
236        String detStr
237       
[1055]238        String destPath = "root:Packages:NIST:VSANS:VCALC"
239       
240        // be sure that the real distance waves exist
241        // TODO -- this may not be the best location?
242
[1134]243        Variable tube_width = 8.4                       // TODO: UNITS!!! Hard-wired value in [mm]
244
[1055]245        String orientation
246        Variable dimX,dimY
247        dimX = DimSize(data,0)
248        dimY = DimSize(data,1)
249        if(dimX > dimY)
250                orientation = "horizontal"
251        else
252                orientation = "vertical"
253        endif
254       
[1134]255
[993]256        Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
257        Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
258        NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
[1055]259
[1108]260        if(kBCTR_CM && cmpstr(detStr,"B") != 0)
[1055]261                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
[1062]262                        // beam ctr is in cm already
263
[1055]264                        // calculate all of the q-values
265                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
266                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
267                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
268                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
269               
270        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
271        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
272                else
273                        // do the q-calculation using linear detector
274                        //VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
275                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
276                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
277                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
278                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
279                endif   
[993]280       
[1055]281       
282        else
283        // using the old calculation with beam center in pixels
284                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
285                        // convert the beam centers to mm
286//                      String orientation
287                        Variable newX,newY
288                        dimX = DimSize(data_realDistX,0)
289                        dimY = DimSize(data_realDistX,1)
290                        if(dimX > dimY)
291                                orientation = "horizontal"
292                        else
293                                orientation = "vertical"
294                        endif
295                       
[993]296               
[1055]297                //
298                        if(cmpstr(orientation,"vertical")==0)
299                                //      this is data dimensioned as (Ntubes,Npix)
300                                newX = tube_width*xCtr
301                                newY = data_realDistY[0][yCtr]
302                        else
303                                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
304                                newX = data_realDistX[xCtr][0]
305                                newY = tube_width*yCtr
306                        endif   
[993]307       
[1055]308                        //if detector "B", different calculation for the centers (not tubes)
309                        if(cmpstr(detStr,"B")==0)
310                                newX = data_realDistX[xCtr][0]
311                                newY = data_realDistY[0][yCtr]
312                                //newX = xCtr
313                                //newY = yCtr
314                        endif           
315                                       
316                        // calculate all of the q-values
317                        qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
318                        qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
319                        qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
320                        qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
321               
322        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
323        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
[993]324                else
[1055]325                        // do the q-calculation using linear detector
326                        VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
[993]327                endif   
328       
329        endif
330       
[1093]331        KillWaves/Z tmpCalib,tmpCalibX,tmpCalibY
[1055]332       
[993]333        return(0)
334End
335
336
[1134]337
338// make the data_realDistX,Y Waves that are needed for the calculation of q
339Function VC_MakeRealDistXYWaves(data,detStr)
340        Wave data
341        String detStr
342       
343        String destPath = "root:Packages:NIST:VSANS:VCALC"
344
345        // calibration waves do not exist yet, so make some fake ones   '
346        // do I count on the orientation as an input, or do I just figure it out on my own?
347        String orientation
348        Variable dimX,dimY
349        dimX = DimSize(data,0)
350        dimY = DimSize(data,1)
351        if(dimX > dimY)
352                orientation = "horizontal"
353        else
354                orientation = "vertical"
355        endif
356       
357        if(cmpstr(orientation,"vertical")==0)
358                Make/O/D/N=(3,48) tmpCalib
359                // for the "tall" L/R banks
360                tmpCalib[0][] = -512
361                tmpCalib[1][] = 8
362                tmpCalib[2][] = 0
363        else
364                Make/O/D/N=(3,48) tmpCalib
365                // for the "short" T/B banks
366                tmpCalib[0][] = -256
367                tmpCalib[1][] = 4
368                tmpCalib[2][] = 0
369        endif
370        // override if back panel
371        if(cmpstr(detStr,"B") == 0)
372                // and for the back detector "B"
373                Make/O/D/N=3 tmpCalibX,tmpCalibY
374                tmpCalibX[0] = VCALC_getPixSizeX(detStr)                        // pixel size in [cm]  VCALC_getPixSizeX(detStr) is [cm]
375                tmpCalibX[1] = 1
376                tmpcalibX[2] = 10000
377                tmpCalibY[0] = VCALC_getPixSizeY(detStr)                        // pixel size in [cm]  VCALC_getPixSizeX(detStr) is [cm]
378                tmpCalibY[1] = 1
379                tmpcalibY[2] = 10000
380        endif
381
382
383//      Wave w_calib = V_getDetTube_spatialCalib("VCALC",detStr)
384        Variable tube_width = 8.4                       // TODO: UNITS!!! Hard-wired value in [mm]
385        if(cmpstr(detStr,"B") == 0)
386                V_NonLinearCorrection_B("VCALC",data,tmpCalibX,tmpCalibY,detStr,destPath)
387                // beam center is in pixels, so use the old routine
388                V_ConvertBeamCtrPix_to_mmB("VCALC","B",destPath)
389        else
390                V_NonLinearCorrection("VCALC",data,tmpCalib,tube_width,detStr,destPath)
391        endif
392       
393       
394        return(0)
395End
396
397
398
[954]399//////////////////////
400// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
401// accout the different x/y pixel sizes and the beam center not being on the detector -
402// off a different edge for each LRTB type
403/////////////////////
404
405//function to calculate the overall q-value, given all of the necesary trig inputs
406//and are in detector coordinates (1,128) rather than axis values
407//the pixel locations need not be integers, reals are ok inputs
[1062]408//sdd is in [cm]
[954]409//wavelength is in Angstroms
410//
411//returned magnitude of Q is in 1/Angstroms
412//
[1055]413//
[1137]414Threadsafe Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]415        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
416       
417        Variable dx,dy,qval,two_theta,dist
418               
[1062]419
[954]420        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
421        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
422        dist = sqrt(dx^2 + dy^2)
423       
424        two_theta = atan(dist/sdd)
425
426        qval = 4*Pi/lam*sin(two_theta/2)
427       
428        return qval
429End
430
431//calculates just the q-value in the x-direction on the detector
432//input/output is the same as CalcQval()
433//ALL inputs are in detector coordinates
434//
[1062]435//sdd is in [cm]
[954]436//wavelength is in Angstroms
437//
438// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
439// now properly accounts for qz
440//
[1137]441Threadsafe Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]442        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
443
444        Variable qx,qval,phi,dx,dy,dist,two_theta
445       
[982]446        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]447       
[1062]448//      sdd *=100               //convert to cm
[954]449        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
450        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
451        phi = V_FindPhi(dx,dy)
452       
453        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
454        dist = sqrt(dx^2 + dy^2)
455        two_theta = atan(dist/sdd)
456
457        qx = qval*cos(two_theta/2)*cos(phi)
458       
459        return qx
460End
461
462//calculates just the q-value in the y-direction on the detector
463//input/output is the same as CalcQval()
464//ALL inputs are in detector coordinates
[1062]465//sdd is in [cm]
[954]466//wavelength is in Angstroms
467//
468// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
469// now properly accounts for qz
470//
[1137]471Threadsafe Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]472        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
473       
474        Variable dy,qval,dx,phi,qy,dist,two_theta
475       
[982]476        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]477       
[1062]478//      sdd *=100               //convert to cm
[954]479        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
480        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
481        phi = V_FindPhi(dx,dy)
482       
483        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
484        dist = sqrt(dx^2 + dy^2)
485        two_theta = atan(dist/sdd)
486       
487        qy = qval*cos(two_theta/2)*sin(phi)
488       
489        return qy
490End
491
492//calculates just the z-component of the q-vector, not measured on the detector
493//input/output is the same as CalcQval()
494//ALL inputs are in detector coordinates
[1062]495//sdd is in [cm]
[954]496//wavelength is in Angstroms
497//
498// not actually used, but here for completeness if anyone asks
499//
[1137]500Threadsafe Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]501        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
502       
503        Variable dy,qval,dx,phi,qz,dist,two_theta
504       
[982]505        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
[954]506       
[1062]507//      sdd *=100               //convert to cm
[954]508       
509        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
510        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
511        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
512        dist = sqrt(dx^2 + dy^2)
513        two_theta = atan(dist/sdd)
514       
515        qz = qval*sin(two_theta/2)
516       
517        return qz
518End
519
520//phi is defined from +x axis, proceeding CCW around [0,2Pi]
521Threadsafe Function V_FindPhi(vx,vy)
522        variable vx,vy
523       
524        variable phi
525       
526        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
527       
528        // special cases
529        if(vx==0 && vy > 0)
530                return(pi/2)
531        endif
532        if(vx==0 && vy < 0)
533                return(3*pi/2)
534        endif
535        if(vx >= 0 && vy == 0)
536                return(0)
537        endif
538        if(vx < 0 && vy == 0)
539                return(pi)
540        endif
541       
542       
543        if(vx > 0 && vy > 0)
544                return(phi)
545        endif
546        if(vx < 0 && vy > 0)
547                return(phi + pi)
548        endif
549        if(vx < 0 && vy < 0)
550                return(phi + pi)
551        endif
552        if( vx > 0 && vy < 0)
553                return(phi + 2*pi)
554        endif
555       
556        return(phi)
557end
558
[982]559Function VC_SphereForm(scale,radius,delrho,bkg,x)                               
[954]560        Variable scale,radius,delrho,bkg
561        Variable x
562       
563        // variables are:                                                       
564        //[0] scale
565        //[1] radius (A)
566        //[2] delrho (A-2)
567        //[3] background (cm-1)
568       
569//      Variable scale,radius,delrho,bkg                               
570//      scale = w[0]
571//      radius = w[1]
572//      delrho = w[2]
573//      bkg = w[3]
574       
575       
576        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
577        // and is rescaled to give [=] cm^-1
578       
579        Variable bes,f,vol,f2
580        //
581        //handle q==0 separately
582        If(x==0)
583                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
584                return(f)
585        Endif
586       
587//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
588       
589        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
590       
591        vol = 4*pi/3*radius^3
592        f = vol*bes*delrho              // [=] A
593        // normalize to single particle volume, convert to 1/cm
594        f2 = f * f / vol * 1.0e8                // [=] 1/cm
595       
596        return (scale*f2+bkg)   // Scale, then add in the background
597       
598End
599
[1117]600ThreadSafe Function VC_Debye(scale,rg,bkg,x)
[954]601        Variable scale,rg,bkg
602        Variable x
603       
604        // variables are:
605        //[0] scale factor
606        //[1] radius of gyration [A]
607        //[2] background        [cm-1]
608       
609        // calculates (scale*debye)+bkg
610        Variable Pq,qr2
611       
612        qr2=(x*rg)^2
613        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
614       
615        //scale
616        Pq *= scale
617        // then add in the background
618        return (Pq+bkg)
619End
620
[955]621// a sum of a power law and debye to approximate the scattering from a real empty cell
622//
623//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
624//
[982]625Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
[955]626        Variable aa,mm,scale,rg,bkg
627        Variable x
628       
629        // variables are:
630        //[0] = A
631        //[1] = power m
632        //[2] scale factor
633        //[3] radius of gyration [A]
634        //[4] background        [cm-1]
635       
636        Variable Iq
637       
638        // calculates (scale*debye)+bkg
639        Variable Pq,qr2
640       
641//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
642//              return(scale+bkg)
643//      endif
644       
645        Iq = aa*x^-mm
646       
647        qr2=(x*rg)^2
648        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
649       
650        //scale
651        Pq *= scale
652        // then add the terms up
653        return (Iq + Pq + bkg)
654End
[954]655
[955]656// blocked beam
657//
[982]658Function VC_BlockedBeam(bkg,x)
[955]659        Variable bkg
660        Variable x
661       
662        return (bkg)
663End
[954]664
665
[955]666//
667// a broad peak to simulate silver behenate or vycor
668//
669// peak @ 0.1 ~ AgBeh
670//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
671//
672//
673// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
674//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
675//
676//
[982]677Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
[955]678        Variable aa,nn,cc,LL,Qzero,mm,bgd
679        Variable x
680       
681        // variables are:                                                       
682        //[0] Porod term scaling
683        //[1] Porod exponent
684        //[2] Lorentzian term scaling
685        //[3] Lorentzian screening length [A]
686        //[4] peak location [1/A]
687        //[5] Lorentzian exponent
688        //[6] background
689       
690//      local variables
691        Variable inten, qval
692//      x is the q-value for the calculation
693        qval = x
694//      do the calculation and return the function value
695       
696        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
697
698        Return (inten)
699       
700End
701
[982]702//
703// updated to new folder structure Feb 2016
704// folderStr = RAW,SAM, VCALC or other
705// detStr is the panel identifer "ML", etc.
706//
[1095]707Function V_SetDeltaQ(folderStr,detStr)
[982]708        String folderStr,detStr
[955]709
[1192]710        NVAR binWidth = root:Packages:NIST:VSANS:Globals:gBinWidth
711
[982]712        Variable isVCALC
713        if(cmpstr(folderStr,"VCALC") == 0)
714                isVCALC = 1
715        endif
716       
717        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
[992]718        String instPath = ":entry:instrument:detector_"
[982]719               
720        if(isVCALC)
721                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
722        else
723                Wave inten = V_getDetectorDataW(folderStr,detStr)
724        endif
[955]725
[982]726        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
727        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
[954]728       
729        Variable xDim,yDim,delQ
730       
731        xDim=DimSize(inten,0)
732        yDim=DimSize(inten,1)
733       
734        if(xDim<yDim)
735                delQ = abs(qx[0][0] - qx[1][0])/2
736        else
737                delQ = abs(qy[0][1] - qy[0][0])/2
738        endif
[1117]739
740// TODOHIGHRES
741// -- is this how I want to handle the too-fine resolution of 1x1 binning?
742        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
743
744        if(cmpstr(detStr,"B") == 0 && gHighResBinning == 1)
745                delQ = 4*delQ
746                Print "Reset delta Q for binning the back detector to 4x pix = ",delQ
747        endif
[954]748       
[1192]749        // multiply the deltaQ by the binWidth (=multiple of pixels)
750        // this defaults to 1, and is set in VSANS preferences
751        delQ *= binWidth
752       
[954]753        // set the global
[982]754        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
[954]755//      Print "SET delQ = ",delQ," for ",type
756       
[982]757        return(delQ)
[954]758end
759
760
761//TODO -- need a switch here to dispatch to the averaging type
[982]762Proc VC_BinQxQy_to_1D(folderStr,type)
[954]763        String folderStr
764        String type
765//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
766//      Prompt type,"detector identifier"
767
768
[982]769        VC_fDoBinning_QxQy2D(folderStr, type)
[954]770
771
772/// this is for a tall, narrow slit mode       
[982]773//      VC_fBinDetector_byRows(folderStr,type)
[954]774       
775End
776
777
[982]778// folderStr is RAW, VCALC, SAM, etc.
779// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
780//
781Proc VC_Graph_1D_detType(folderStr,type)
[954]782        String folderStr,type
783       
[982]784        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
[954]785       
786        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
787        ModifyGraph mirror=2,grid=1,log=1
788        ModifyGraph mode=4,marker=19,msize=2
789//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
790        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
791        legend
792       
793        SetDataFolder root:
794
795End
796
797
[955]798
799//////////
800//
801//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
802//              - each pixel QxQyQz has been calculated beforehand
803//              - if multiple panels are selected to be combined, it is done here during the binning
804//              - the setting of deltaQ step is still a little suspect (TODO)
805//
806//
[954]807// see the equivalent function in PlotUtils2D_v40.ipf
808//
809//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
810//
811// this has been modified to accept different detector panels and to take arrays
812// -- type = FL or FR or...other panel identifiers
813//
[955]814// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
[954]815//
[982]816//
817// updated Feb2016 to take new folder structure
818// TODO
819// -- VERIFY
820// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
821//    but this is a lousy choice.
[993]822// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
[1044]823// x- Where do I put the solid angle correction? In here as a weight for each point, or later on as
824//    a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the
825//    step where data is added to a WORK file (see Raw_to_Work())
[982]826//
[999]827//
[1070]828// TODO:
829// -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently
830//    hard-wired. these need to be corrected before even the pinhole resolution is correct
831// x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE),
832//    and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated
833//    is separate functions that are resolution-aware.
[999]834//
[1070]835//
[993]836// folderStr = WORK folder, type = the binning type (may include multiple detectors)
[1097]837Function VC_fDoBinning_QxQy2D(folderStr,type,collimationStr)
838        String folderStr,type,collimationStr
[954]839       
840        Variable nSets = 0
841        Variable xDim,yDim
842        Variable ii,jj
843        Variable qVal,nq,var,avesq,aveisq
[993]844        Variable binIndex,val,isVCALC=0,maskMissing
[954]845
[982]846        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
[992]847        String instPath = ":entry:instrument:detector_"
[982]848        String detStr
849               
850        if(cmpstr(folderStr,"VCALC") == 0)
851                isVCALC = 1
852        endif
[954]853       
[1076]854        detStr = type
855
[954]856// now switch on the type to determine which waves to declare and create
857// since there may be more than one panel to step through. There may be two, there may be four
858//
859
[999]860// TODO:
861// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
862//
[993]863// assume that the mask files are missing unless we can find them. If VCALC data,
864//  then the Mask is missing by definition
865        maskMissing = 1
866
[954]867        strswitch(type) // string switch
[1075]868//              case "FL":              // execute if case matches expression
869//              case "FR":
870//                      detStr = type
871//                      if(isVCALC)
872//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
873//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
874//                      else
875//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
876//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
877//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
878//                              if(WaveExists(mask) == 1)
879//                                      maskMissing = 0
880//                              endif
881//                             
882//                      endif
883//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
884//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
885//                      nSets = 1
886//                      break   
887                                                               
888//              case "FT":             
889//              case "FB":
890//                      detStr = type
891//                      if(isVCALC)
892//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
893//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
894//                      else
895//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
896//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
897//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
898//                              if(WaveExists(mask) == 1)
899//                                      maskMissing = 0
900//                              endif
901//                      endif
902//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
903//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
904//                      nSets = 1
905//                      break
906                       
907//              case "ML":             
908//              case "MR":
909//                      detStr = type
910//                      if(isVCALC)
911//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
912//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
913//                      else
914//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
915//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
916//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
917//                              if(WaveExists(mask) == 1)
918//                                      maskMissing = 0
919//                              endif
920//                      endif   
921//                      //TODO:
922//                      // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
923//                      //    determined specifically for that panel (currently using one tube width as deltaQ)
924//                      // -- this is repeated multiple times in this switch
925//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
926//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
927//                      nSets = 1
928//                      break   
929                                       
930//              case "MT":             
931//              case "MB":
932//                      detStr = type
933//                      if(isVCALC)
934//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
935//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
936//                      else
937//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
938//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
939//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
940//                              if(WaveExists(mask) == 1)
941//                                      maskMissing = 0
942//                              endif
943//                      endif   
944//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
945//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
946//                      nSets = 1
947//                      break   
948
949// only one panel, simply pick that panel and move on out of the switch
950                case "FL":
[954]951                case "FR":
[1075]952                case "FT":
[954]953                case "FB":
[1075]954                case "ML":
[982]955                case "MR":
[1075]956                case "MT":
957                case "MB":                     
958                case "B":       
[982]959                        if(isVCALC)
960                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
961                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
962                        else
963                                Wave inten = V_getDetectorDataW(folderStr,detStr)
964                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
965                        endif   
[1133]966                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
967                        if(WaveExists(mask) == 1)
968                                maskMissing = 0
969                        endif
[982]970                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
971                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
[954]972                        nSets = 1
973                        break   
974                       
975                case "FLR":
[982]976                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
977                // TODO
978                // -- see if I can un-hard-wire some of this below when more than one panel is combined
979                        if(isVCALC)
980                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
981                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
982                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
983                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
984                        else
985                                Wave inten = V_getDetectorDataW(folderStr,"FL")
986                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
987                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
988                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
989                        endif   
[1133]990                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
991                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
992                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
993                                maskMissing = 0
994                        endif
995                       
[982]996                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
997                       
998                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
999                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
1000               
[954]1001                        nSets = 2
1002                        break                   
1003               
1004                case "FTB":
[982]1005                        if(isVCALC)
1006                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
1007                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
1008                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
1009                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
1010                        else
1011                                Wave inten = V_getDetectorDataW(folderStr,"FT")
1012                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
1013                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
1014                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
1015                        endif   
[1133]1016                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1017                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1018                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1019                                maskMissing = 0
1020                        endif
1021                       
[982]1022                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
1023                       
1024                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
1025                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
1026       
[954]1027                        nSets = 2
1028                        break           
1029               
1030                case "FLRTB":
[982]1031                        if(isVCALC)
1032                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
1033                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
1034                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
1035                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
1036                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
1037                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
1038                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
1039                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
1040                        else
1041                                Wave inten = V_getDetectorDataW(folderStr,"FL")
1042                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
1043                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
1044                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
1045                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
1046                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
1047                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
1048                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
[1133]1049                               
[982]1050                        endif   
[1133]1051                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
1052                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
1053                        Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1054                        Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1055                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1056                                maskMissing = 0
1057                        endif
1058                       
[982]1059                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
1060                       
1061                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
1062                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
1063                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
1064                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
1065               
[954]1066                        nSets = 4
1067                        break           
1068                       
1069                case "MLR":
[982]1070                        if(isVCALC)
1071                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1072                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1073                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1074                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1075                        else
1076                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1077                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1078                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1079                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
[1133]1080                               
[982]1081                        endif   
[1133]1082                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1083                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1084                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1085                                maskMissing = 0
1086                        endif
1087                       
[982]1088                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1089                       
1090                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1091                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1092               
[954]1093                        nSets = 2
1094                        break                   
1095               
1096                case "MTB":
[982]1097                        if(isVCALC)
1098                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1099                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
1100                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1101                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1102                        else
1103                                Wave inten = V_getDetectorDataW(folderStr,"MT")
1104                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1105                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1106                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
[1133]1107                               
[982]1108                        endif   
[1133]1109                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1110                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1111                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1112                                maskMissing = 0
1113                        endif
1114                       
[982]1115                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1116                       
1117                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
1118                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1119               
[954]1120                        nSets = 2
1121                        break                           
1122               
1123                case "MLRTB":
[982]1124                        if(isVCALC)
1125                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1126                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1127                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1128                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1129                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1130                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
1131                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1132                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1133                        else
1134                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1135                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1136                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1137                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1138                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1139                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1140                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1141                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
[1133]1142                               
1143                        endif
1144                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1145                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1146                        Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1147                        Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1148                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1149                                maskMissing = 0
1150                        endif
1151                               
[982]1152                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1153                       
1154                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1155                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1156                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
1157                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1158               
[954]1159                        nSets = 4
1160                        break                                                                   
1161                                       
1162                default:
[1075]1163                        nSets = 0                                                       
[954]1164                        Print "ERROR   ---- type is not recognized "
1165        endswitch
1166
1167//      Print "delQ = ",delQ," for ",type
1168
1169        if(nSets == 0)
[982]1170                SetDataFolder root:
[954]1171                return(0)
1172        endif
1173
1174
[1075]1175// RAW data is currently read in and the 2D error wave is correctly generated
1176// 2D error is propagated through all reduction steps, but I have not
1177// verified that it is an exact duplication of the 1D error
[999]1178//
[1075]1179//
1180//
1181// IF ther is no 2D error wave present for some reason, make a fake one
[955]1182        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
[954]1183                Duplicate/O inten,iErr
1184                Wave iErr=iErr
1185//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1186                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1187        endif
[955]1188        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1189                Duplicate/O inten2,iErr2
1190                Wave iErr2=iErr2
1191//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1192                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1193        endif
1194        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1195                Duplicate/O inten3,iErr3
1196                Wave iErr3=iErr3
1197//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1198                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1199        endif
1200        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1201                Duplicate/O inten4,iErr4
1202                Wave iErr4=iErr4
1203//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1204                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1205        endif
[954]1206
[1035]1207        // TODO -- nq will need to be larger, once the back detector is installed
1208        //
[954]1209        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1210        // - so I upped nq to 600
1211
[1093]1212        if(cmpstr(type,"B") == 0)
1213                nq = 8000
1214        else
1215                nq=600
1216        endif
[1035]1217
[954]1218//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1219//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
[982]1220        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1221        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1222        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1223        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1224        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1225        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
[954]1226       
[982]1227        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1228        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1229        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1230        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1231        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1232        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
[954]1233       
1234       
1235//      delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2))               //use bins of 1 pixel width
[955]1236// TODO: not sure if I want to set dQ in x or y direction...
[954]1237        // the short dimension is the 8mm tubes, use this direction as dQ?
1238        // but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0]
1239        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1240
1241
1242        qBin_qxqy[] =  p*delQ   
1243        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1244
1245        iBin_qxqy = 0
1246        iBin2_qxqy = 0
1247        eBin_qxqy = 0
1248        eBin2D_qxqy = 0
1249        nBin_qxqy = 0   //number of intensities added to each bin
1250
1251// now there are situations of:
1252// 1 panel
1253// 2 panels
1254// 4 panels
1255//
1256// this needs to be a double loop now...
[999]1257// TODO:
[1051]1258// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
[1075]1259// -- verify the 2D error propagation by reducing it to 1D error
[1070]1260//
[1075]1261//
[1070]1262// The 1D error does not use iErr, and IS CALCULATED CORRECTLY
1263//
1264// x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW
[954]1265
[1070]1266//
[993]1267// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1268        if(maskMissing == 1)
1269                Print "Mask file not found for at least one detector - so all data is used"
1270        endif
[1095]1271       
1272        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1273        if(gIgnoreDetB && cmpstr(type,"B") == 0)
1274                maskMissing = 1
1275                Print "Mask skipped for B due to possible mismatch (Panel B ignored in preferences)"
1276        endif
[993]1277
1278        Variable mask_val
[954]1279// use set 1 (no number) only
1280        if(nSets >= 1)
1281                xDim=DimSize(inten,0)
1282                yDim=DimSize(inten,1)
1283       
1284                for(ii=0;ii<xDim;ii+=1)
1285                        for(jj=0;jj<yDim;jj+=1)
1286                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1287                                qVal = qTotal[ii][jj]
1288                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1289                                val = inten[ii][jj]
[993]1290                               
[1133]1291//                              if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1292                                if(maskMissing)         // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
[999]1293                                        mask_val = 0
[993]1294                                else
1295                                        mask_val = mask[ii][jj]
1296                                endif
[999]1297                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
[954]1298                                        iBin_qxqy[binIndex] += val
1299                                        iBin2_qxqy[binIndex] += val*val
1300                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1301                                        nBin_qxqy[binIndex] += 1
1302                                endif
1303                        endfor
1304                endfor
1305               
1306        endif
1307
1308// add in set 2 (set 1 already done)
1309        if(nSets >= 2)
1310                xDim=DimSize(inten2,0)
1311                yDim=DimSize(inten2,1)
1312       
1313                for(ii=0;ii<xDim;ii+=1)
1314                        for(jj=0;jj<yDim;jj+=1)
1315                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1316                                qVal = qTotal2[ii][jj]
1317                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1318                                val = inten2[ii][jj]
[993]1319                               
[1133]1320//                              if(isVCALC || maskMissing)
1321                                if(maskMissing)
[999]1322                                        mask_val = 0
[993]1323                                else
1324                                        mask_val = mask2[ii][jj]
1325                                endif
[999]1326                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
[954]1327                                        iBin_qxqy[binIndex] += val
1328                                        iBin2_qxqy[binIndex] += val*val
[955]1329                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
[954]1330                                        nBin_qxqy[binIndex] += 1
1331                                endif
1332                        endfor
1333                endfor
1334               
1335        endif
1336
[982]1337// add in set 3 and 4 (set 1 and 2 already done)
[954]1338        if(nSets == 4)
1339                xDim=DimSize(inten3,0)
1340                yDim=DimSize(inten3,1)
1341       
1342                for(ii=0;ii<xDim;ii+=1)
1343                        for(jj=0;jj<yDim;jj+=1)
1344                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1345                                qVal = qTotal3[ii][jj]
1346                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1347                                val = inten3[ii][jj]
[993]1348                               
[1133]1349//                              if(isVCALC || maskMissing)
1350                                if(maskMissing)
[999]1351                                        mask_val = 0
[993]1352                                else
1353                                        mask_val = mask3[ii][jj]
1354                                endif
[999]1355                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
[954]1356                                        iBin_qxqy[binIndex] += val
1357                                        iBin2_qxqy[binIndex] += val*val
[955]1358                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
[954]1359                                        nBin_qxqy[binIndex] += 1
1360                                endif
1361                        endfor
1362                endfor
1363               
1364               
1365                xDim=DimSize(inten4,0)
1366                yDim=DimSize(inten4,1)
1367       
1368                for(ii=0;ii<xDim;ii+=1)
1369                        for(jj=0;jj<yDim;jj+=1)
1370                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1371                                qVal = qTotal4[ii][jj]
1372                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1373                                val = inten4[ii][jj]
[993]1374                               
[1133]1375//                              if(isVCALC || maskMissing)
1376                                if(maskMissing)
[999]1377                                        mask_val = 0
[993]1378                                else
1379                                        mask_val = mask4[ii][jj]
1380                                endif
[999]1381                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
[954]1382                                        iBin_qxqy[binIndex] += val
1383                                        iBin2_qxqy[binIndex] += val*val
[955]1384                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
[954]1385                                        nBin_qxqy[binIndex] += 1
1386                                endif
1387                        endfor
1388                endfor
1389               
1390        endif
1391
1392
1393// after looping through all of the data on the panels, calculate errors on I(q),
1394// just like in CircSectAve.ipf
[999]1395// TODO:
[1075]1396// -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet)
[1019]1397// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
[954]1398        for(ii=0;ii<nq;ii+=1)
1399                if(nBin_qxqy[ii] == 0)
1400                        //no pixels in annuli, data unknown
1401                        iBin_qxqy[ii] = 0
1402                        eBin_qxqy[ii] = 1
1403                        eBin2D_qxqy[ii] = NaN
1404                else
1405                        if(nBin_qxqy[ii] <= 1)
1406                                //need more than one pixel to determine error
1407                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1408                                eBin_qxqy[ii] = 1
1409                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1410                        else
1411                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
[1019]1412                                //  -- this is correctly calculating the error as the standard error of the mean, as
1413                                //    was always done for SANS as well.
[954]1414                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1415                                avesq = iBin_qxqy[ii]^2
1416                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1417                                var = aveisq-avesq
1418                                if(var<=0)
1419                                        eBin_qxqy[ii] = 1e-6
1420                                else
1421                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1422                                endif
1423                                // and calculate as it is propagated pixel-by-pixel
1424                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1425                        endif
1426                endif
1427        endfor
1428       
1429        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1430       
1431        // find the last non-zero point, working backwards
1432        val=nq
1433        do
1434                val -= 1
1435        while((nBin_qxqy[val] == 0) && val > 0)
1436       
1437//      print val, nBin_qxqy[val]
1438        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1439
1440        if(val == 0)
1441                // all the points were deleted
1442                return(0)
1443        endif
1444       
1445       
1446        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1447        // find the first non-zero point, working forwards
1448        val = -1
1449        do
1450                val += 1
1451        while(nBin_qxqy[val] == 0)     
1452        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1453
1454        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1455        val = numpnts(nBin_qxqy)-1
1456        do
1457                if(nBin_qxqy[val] == 0)
1458                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1459                endif
1460                val -= 1
1461        while(val>0)
[1099]1462
1463// utility function to remove NaN values from the waves
1464        V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
1465
[954]1466       
[1044]1467        // TODO:
[1064]1468        // -- This is where I calculate the resolution in SANS (see CircSectAve)
[1044]1469        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
[1064]1470        // -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
[1044]1471        //
[1064]1472        nq = numpnts(qBin_qxqy)
1473        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1474        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+type)
1475        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1476        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1477        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1478        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
[1119]1479        Variable ret1,ret2,ret3                 
[1064]1480
[1119]1481// all of the different collimation conditions are handled within the V_getResolution function
1482// which is responsible for switching based on the different collimation types (white beam, slit, Xtal, etc)
[1090]1483// to calculate the correct resolution, or fill the waves with the correct "flags"
1484//
1485
1486// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
[1095]1487//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
[1090]1488//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
[1095]1489//  with the inner(lambda) calculated first, then the outer(geometry).
[1090]1490//
1491
[1097]1492// possible values are:
1493//
1494// pinhole
1495// pinhole_whiteBeam
[1098]1496// convergingPinholes
1497//
1498// *slit data should be reduced using the slit routine, not here, proceed but warn
[1097]1499// narrowSlit
1500// narrowSlit_whiteBeam
1501//
[1119]1502        ii=0
1503        do
1504                V_getResolution(qBin_qxqy[ii],folderStr,type,collimationStr,ret1,ret2,ret3)
1505                sigmaq[ii] = ret1       
1506                qbar[ii] = ret2
1507                fsubs[ii] = ret3       
1508                ii+=1
1509        while(ii<nq)
[1097]1510
1511
[954]1512        SetDataFolder root:
1513       
1514        return(0)
1515End
1516
1517
Note: See TracBrowser for help on using the repository browser.