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

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

removed the doubled "entry" field from the VSANS file load.

appears now to work fine with R/W routines and with VCALC.

File size: 32.5 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3/////////////////////////
4//
5// Utility functions to:
6//              calculate Q, Qx, Qy, Qz
7//              fill the detector panels with simulated data (the model functions are here)
8//              bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width)
9//
10/////////////////////////
11
12
13
14
15// TODO: hard wired for a sphere - change this to allow minimal selections and altering of coefficients
16// TODO: add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
17//
18// NOTE - this is a VCALC only routine, so it's not been made completely generic
19//
20Function FillPanel_wModelData(det,qTot,type)
21        Wave det,qTot
22        String type
23
24//      SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
25
26        // q-values and detector arrays already allocated and calculated
27        Duplicate/O det tmpInten,tmpSig,prob_i
28               
29        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
30
31        //imon = VC_BeamIntensity()*CountTime
32        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
33        trans = 0.8
34        thick = 0.1
35       
36        // need SDD
37        // need pixel dimensions
38        // nominal sdd in meters, offset in mm, want result in cm !
39        sdd = VCALC_getSDD(type)*100    +  VCALC_getTopBottomSDDOffset(type) / 10               // result is sdd in [cm]
40
41        pixSizeX = VCALC_getPixSizeX(type)              // cm
42        pixSizeY = VCALC_getPixSizeY(type)
43       
44       
45        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
46        //  keep the panel to keep people used to using it.)
47        // peak @ 0.1 ~ AgBeh
48        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
49        //
50        // peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
51        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
52        String funcStr = VCALC_getModelFunctionStr()
53        strswitch(funcStr)
54                case "Big Debye":
55                        tmpInten = VC_Debye(100,3000,0.0001,qTot[p][q])
56                        break
57                case "Big Sphere":
58                        tmpInten = VC_SphereForm(1,900,1e-6,0.01,qTot[p][q])   
59                        break
60                case "Debye":
61                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
62                        break
63                case "Sphere":
64                        tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])   
65                        break
66                case "AgBeh":
67                        tmpInten = VC_BroadPeak(1e-9,3,20,100.0,0.1,3,0.1,qTot[p][q])
68                        break
69                case "Vycor":
70                        tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
71                        break   
72                case "Empty Cell":
73                        tmpInten = VC_EC_Empirical(2.2e-8,3.346,0.0065,9.0,0.016,qTot[p][q])
74                        break
75                case "Blocked Beam":
76                        tmpInten = VC_BlockedBeam(1,qTot[p][q])
77                        break
78                default:
79                        tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
80        endswitch
81
82// TODO: this is faked to get around the singularity at the center of the back detector
83//
84//
85        if(cmpstr(type,"B") == 0)
86                Variable nx,ny,px,py
87                nx = VCALC_get_nPix_X(type)
88                ny = VCALC_get_nPix_Y(type)
89                px = trunc(nx/2)
90                py = trunc(ny/2)
91               
92                tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
93        endif
94
95
96
97///////////////
98//      // calculate the scattering cross section simply to be able to estimate the transmission
99//      Variable sig_sas=0
100//     
101//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
102//      // subtracted before the calculation.
103//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
104//
105//      if(sig_sas > 100)
106//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
107//      endif           
108//
109//      // calculate the multiple scattering fraction for display (10/2009)
110//      Variable ii,nMax=10,tau
111//      mScat=0
112//      tau = thick*sig_sas
113//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
114//      // neutrons out of those that were scattered
115//      for(ii=2;ii<nMax;ii+=1)
116//              mScat += tau^(ii)/factorial(ii)
117////            print tau^(ii)/factorial(ii)
118//      endfor
119//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
120//      mscat *= (estTrans)/(1-estTrans)
121//
122////    if(mScat > 0.1)         //  Display warning
123//
124//      Print "Sig_sas = ",sig_sas
125////////////////////
126       
127        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
128               
129        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
130
131
132/// **** can I safely assume a Gaussian error in the count rate??
133        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
134
135        tmpInten += gnoise(tmpSig)
136        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
137        tmpInten = trunc(tmpInten)
138               
139       
140        det = tmpInten
141
142// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
143//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
144
145       
146        KillWaves/Z tmpInten,tmpSig,prob_i     
147        SetDataFolder root:
148
149        return(0)
150End
151
152
153// For a given detector panel, calculate the q-values
154// -work with everything as arrays
155// Input needed:
156// detector data
157// detector type (LRTB?)
158// beam center (may be off the detector)
159// SDD
160// lambda
161//
162// pixel dimensions for detector type (global constants)
163// - data dimensions read directly from array
164//
165// --What is calculated:
166// array of Q
167// array of qx,qy,qz
168// array of error already exists
169//
170//
171// -- sdd in meters
172// -- lambda in Angstroms
173Function VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
174        Wave data,qTot,qx,qy,qz
175        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
176               
177        // loop over the array and calculate the values - this is done as a wave assignment
178// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
179        qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
180        qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
181        qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
182        qz = VC_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
183       
184        return(0)
185End
186
187
188//////////////////////
189// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
190// accout the different x/y pixel sizes and the beam center not being on the detector -
191// off a different edge for each LRTB type
192/////////////////////
193
194//function to calculate the overall q-value, given all of the necesary trig inputs
195//and are in detector coordinates (1,128) rather than axis values
196//the pixel locations need not be integers, reals are ok inputs
197//sdd is in meters
198//wavelength is in Angstroms
199//
200//returned magnitude of Q is in 1/Angstroms
201//
202Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
203        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
204       
205        Variable dx,dy,qval,two_theta,dist
206               
207        sdd *=100               //convert to cm
208        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
209        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
210        dist = sqrt(dx^2 + dy^2)
211       
212        two_theta = atan(dist/sdd)
213
214        qval = 4*Pi/lam*sin(two_theta/2)
215       
216        return qval
217End
218
219//calculates just the q-value in the x-direction on the detector
220//input/output is the same as CalcQval()
221//ALL inputs are in detector coordinates
222//
223//sdd is in meters
224//wavelength is in Angstroms
225//
226// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
227// now properly accounts for qz
228//
229Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
230        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
231
232        Variable qx,qval,phi,dx,dy,dist,two_theta
233       
234        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
235       
236        sdd *=100               //convert to cm
237        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
238        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
239        phi = V_FindPhi(dx,dy)
240       
241        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
242        dist = sqrt(dx^2 + dy^2)
243        two_theta = atan(dist/sdd)
244
245        qx = qval*cos(two_theta/2)*cos(phi)
246       
247        return qx
248End
249
250//calculates just the q-value in the y-direction on the detector
251//input/output is the same as CalcQval()
252//ALL inputs are in detector coordinates
253//sdd is in meters
254//wavelength is in Angstroms
255//
256// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
257// now properly accounts for qz
258//
259Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
260        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
261       
262        Variable dy,qval,dx,phi,qy,dist,two_theta
263       
264        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
265       
266        sdd *=100               //convert to cm
267        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
268        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
269        phi = V_FindPhi(dx,dy)
270       
271        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
272        dist = sqrt(dx^2 + dy^2)
273        two_theta = atan(dist/sdd)
274       
275        qy = qval*cos(two_theta/2)*sin(phi)
276       
277        return qy
278End
279
280//calculates just the z-component of the q-vector, not measured on the detector
281//input/output is the same as CalcQval()
282//ALL inputs are in detector coordinates
283//sdd is in meters
284//wavelength is in Angstroms
285//
286// not actually used, but here for completeness if anyone asks
287//
288Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
289        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
290       
291        Variable dy,qval,dx,phi,qz,dist,two_theta
292       
293        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
294       
295        sdd *=100               //convert to cm
296       
297        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
298        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
299        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
300        dist = sqrt(dx^2 + dy^2)
301        two_theta = atan(dist/sdd)
302       
303        qz = qval*sin(two_theta/2)
304       
305        return qz
306End
307
308//phi is defined from +x axis, proceeding CCW around [0,2Pi]
309Threadsafe Function V_FindPhi(vx,vy)
310        variable vx,vy
311       
312        variable phi
313       
314        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
315       
316        // special cases
317        if(vx==0 && vy > 0)
318                return(pi/2)
319        endif
320        if(vx==0 && vy < 0)
321                return(3*pi/2)
322        endif
323        if(vx >= 0 && vy == 0)
324                return(0)
325        endif
326        if(vx < 0 && vy == 0)
327                return(pi)
328        endif
329       
330       
331        if(vx > 0 && vy > 0)
332                return(phi)
333        endif
334        if(vx < 0 && vy > 0)
335                return(phi + pi)
336        endif
337        if(vx < 0 && vy < 0)
338                return(phi + pi)
339        endif
340        if( vx > 0 && vy < 0)
341                return(phi + 2*pi)
342        endif
343       
344        return(phi)
345end
346
347Function VC_SphereForm(scale,radius,delrho,bkg,x)                               
348        Variable scale,radius,delrho,bkg
349        Variable x
350       
351        // variables are:                                                       
352        //[0] scale
353        //[1] radius (A)
354        //[2] delrho (A-2)
355        //[3] background (cm-1)
356       
357//      Variable scale,radius,delrho,bkg                               
358//      scale = w[0]
359//      radius = w[1]
360//      delrho = w[2]
361//      bkg = w[3]
362       
363       
364        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
365        // and is rescaled to give [=] cm^-1
366       
367        Variable bes,f,vol,f2
368        //
369        //handle q==0 separately
370        If(x==0)
371                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
372                return(f)
373        Endif
374       
375//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
376       
377        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
378       
379        vol = 4*pi/3*radius^3
380        f = vol*bes*delrho              // [=] A
381        // normalize to single particle volume, convert to 1/cm
382        f2 = f * f / vol * 1.0e8                // [=] 1/cm
383       
384        return (scale*f2+bkg)   // Scale, then add in the background
385       
386End
387
388Function VC_Debye(scale,rg,bkg,x)
389        Variable scale,rg,bkg
390        Variable x
391       
392        // variables are:
393        //[0] scale factor
394        //[1] radius of gyration [A]
395        //[2] background        [cm-1]
396       
397        // calculates (scale*debye)+bkg
398        Variable Pq,qr2
399       
400        qr2=(x*rg)^2
401        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
402       
403        //scale
404        Pq *= scale
405        // then add in the background
406        return (Pq+bkg)
407End
408
409// a sum of a power law and debye to approximate the scattering from a real empty cell
410//
411//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
412//
413Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
414        Variable aa,mm,scale,rg,bkg
415        Variable x
416       
417        // variables are:
418        //[0] = A
419        //[1] = power m
420        //[2] scale factor
421        //[3] radius of gyration [A]
422        //[4] background        [cm-1]
423       
424        Variable Iq
425       
426        // calculates (scale*debye)+bkg
427        Variable Pq,qr2
428       
429//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
430//              return(scale+bkg)
431//      endif
432       
433        Iq = aa*x^-mm
434       
435        qr2=(x*rg)^2
436        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
437       
438        //scale
439        Pq *= scale
440        // then add the terms up
441        return (Iq + Pq + bkg)
442End
443
444// blocked beam
445//
446Function VC_BlockedBeam(bkg,x)
447        Variable bkg
448        Variable x
449       
450        return (bkg)
451End
452
453
454//
455// a broad peak to simulate silver behenate or vycor
456//
457// peak @ 0.1 ~ AgBeh
458//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
459//
460//
461// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
462//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
463//
464//
465Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
466        Variable aa,nn,cc,LL,Qzero,mm,bgd
467        Variable x
468       
469        // variables are:                                                       
470        //[0] Porod term scaling
471        //[1] Porod exponent
472        //[2] Lorentzian term scaling
473        //[3] Lorentzian screening length [A]
474        //[4] peak location [1/A]
475        //[5] Lorentzian exponent
476        //[6] background
477       
478//      local variables
479        Variable inten, qval
480//      x is the q-value for the calculation
481        qval = x
482//      do the calculation and return the function value
483       
484        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
485
486        Return (inten)
487       
488End
489
490//
491// updated to new folder structure Feb 2016
492// folderStr = RAW,SAM, VCALC or other
493// detStr is the panel identifer "ML", etc.
494//
495Function SetDeltaQ(folderStr,detStr)
496        String folderStr,detStr
497
498        Variable isVCALC
499        if(cmpstr(folderStr,"VCALC") == 0)
500                isVCALC = 1
501        endif
502       
503        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
504        String instPath = ":entry:instrument:detector_"
505               
506        if(isVCALC)
507                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
508        else
509                Wave inten = V_getDetectorDataW(folderStr,detStr)
510        endif
511
512        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
513        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
514       
515        Variable xDim,yDim,delQ
516       
517        xDim=DimSize(inten,0)
518        yDim=DimSize(inten,1)
519       
520        if(xDim<yDim)
521                delQ = abs(qx[0][0] - qx[1][0])/2
522        else
523                delQ = abs(qy[0][1] - qy[0][0])/2
524        endif
525       
526        // set the global
527        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
528//      Print "SET delQ = ",delQ," for ",type
529       
530        return(delQ)
531end
532
533
534//TODO -- need a switch here to dispatch to the averaging type
535Proc VC_BinQxQy_to_1D(folderStr,type)
536        String folderStr
537        String type
538//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
539//      Prompt type,"detector identifier"
540
541
542        VC_fDoBinning_QxQy2D(folderStr, type)
543
544
545/// this is for a tall, narrow slit mode       
546//      VC_fBinDetector_byRows(folderStr,type)
547       
548End
549
550
551// folderStr is RAW, VCALC, SAM, etc.
552// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
553//
554Proc VC_Graph_1D_detType(folderStr,type)
555        String folderStr,type
556       
557        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
558       
559        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
560        ModifyGraph mirror=2,grid=1,log=1
561        ModifyGraph mode=4,marker=19,msize=2
562//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
563        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
564        legend
565       
566        SetDataFolder root:
567
568End
569
570
571
572//////////
573//
574//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
575//              - each pixel QxQyQz has been calculated beforehand
576//              - if multiple panels are selected to be combined, it is done here during the binning
577//              - the setting of deltaQ step is still a little suspect (TODO)
578//
579//
580// see the equivalent function in PlotUtils2D_v40.ipf
581//
582//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
583//
584// this has been modified to accept different detector panels and to take arrays
585// -- type = FL or FR or...other panel identifiers
586//
587// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
588//
589//
590// updated Feb2016 to take new folder structure
591// TODO
592// -- VERIFY
593// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
594//    but this is a lousy choice.
595//
596Function VC_fDoBinning_QxQy2D(folderStr,type)
597        String folderStr,type
598       
599        Variable nSets = 0
600        Variable xDim,yDim
601        Variable ii,jj
602        Variable qVal,nq,var,avesq,aveisq
603        Variable binIndex,val,isVCALC=0
604
605        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
606        String instPath = ":entry:instrument:detector_"
607        String detStr
608               
609        if(cmpstr(folderStr,"VCALC") == 0)
610                isVCALC = 1
611        endif
612       
613// now switch on the type to determine which waves to declare and create
614// since there may be more than one panel to step through. There may be two, there may be four
615//
616
617        strswitch(type) // string switch
618                case "FL":              // execute if case matches expression
619                case "FR":
620                        detStr = type
621                        if(isVCALC)
622                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
623                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
624                        else
625                                Wave inten = V_getDetectorDataW(folderStr,detStr)
626                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
627                        endif
628                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
629                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
630                        nSets = 1
631                        break   
632                                                               
633                case "FT":             
634                case "FB":
635                        detStr = type
636                        if(isVCALC)
637                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
638                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
639                        else
640                                Wave inten = V_getDetectorDataW(folderStr,detStr)
641                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
642                        endif
643                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
644                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
645                        nSets = 1
646                        break
647                       
648                case "ML":             
649                case "MR":
650                        detStr = type
651                        if(isVCALC)
652                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
653                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
654                        else
655                                Wave inten = V_getDetectorDataW(folderStr,detStr)
656                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
657                        endif   
658                        //TODO:
659                        // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
660                        //    determined specifically for that panel (currently using one tube width as deltaQ)
661                        // -- this is repeated multiple times in this switch
662                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
663//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
664                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
665                        nSets = 1
666                        break   
667                                       
668                case "MT":             
669                case "MB":
670                        detStr = type
671                        if(isVCALC)
672                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
673                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
674                        else
675                                Wave inten = V_getDetectorDataW(folderStr,detStr)
676                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
677                        endif   
678                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
679                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
680                        nSets = 1
681                        break   
682                                       
683                case "B":       
684                        detStr = type
685                        if(isVCALC)
686                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
687                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
688                        else
689                                Wave inten = V_getDetectorDataW(folderStr,detStr)
690                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
691                        endif   
692                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_B")
693                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
694                        nSets = 1
695                        break   
696                       
697                case "FLR":
698                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
699                // TODO
700                // -- see if I can un-hard-wire some of this below when more than one panel is combined
701                        if(isVCALC)
702                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
703                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
704                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
705                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
706                        else
707                                Wave inten = V_getDetectorDataW(folderStr,"FL")
708                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
709                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
710                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
711                        endif   
712                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
713                       
714                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
715                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
716               
717                        nSets = 2
718                        break                   
719               
720                case "FTB":
721                        if(isVCALC)
722                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
723                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
724                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
725                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
726                        else
727                                Wave inten = V_getDetectorDataW(folderStr,"FT")
728                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
729                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
730                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
731                        endif   
732                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
733                       
734                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
735                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
736       
737                        nSets = 2
738                        break           
739               
740                case "FLRTB":
741                        if(isVCALC)
742                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
743                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
744                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
745                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
746                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
747                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
748                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
749                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
750                        else
751                                Wave inten = V_getDetectorDataW(folderStr,"FL")
752                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
753                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
754                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
755                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
756                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
757                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
758                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
759                        endif   
760                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
761                       
762                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
763                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
764                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
765                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
766               
767                        nSets = 4
768                        break           
769                       
770                case "MLR":
771                        if(isVCALC)
772                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
773                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
774                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
775                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
776                        else
777                                Wave inten = V_getDetectorDataW(folderStr,"ML")
778                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
779                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
780                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
781                        endif   
782                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
783                       
784                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
785                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
786               
787                        nSets = 2
788                        break                   
789               
790                case "MTB":
791                        if(isVCALC)
792                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
793                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
794                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
795                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
796                        else
797                                Wave inten = V_getDetectorDataW(folderStr,"MT")
798                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
799                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
800                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
801                        endif   
802                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
803                       
804                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
805                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
806               
807                        nSets = 2
808                        break                           
809               
810                case "MLRTB":
811                        if(isVCALC)
812                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
813                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
814                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
815                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
816                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
817                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
818                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
819                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
820                        else
821                                Wave inten = V_getDetectorDataW(folderStr,"ML")
822                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
823                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
824                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
825                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
826                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
827                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
828                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
829                        endif   
830                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
831                       
832                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
833                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
834                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
835                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
836               
837                        nSets = 4
838                        break                                                                   
839                                       
840                default:
841                        nSets = 0                                                       // optional default expression executed
842                        Print "ERROR   ---- type is not recognized "
843        endswitch
844
845//      Print "delQ = ",delQ," for ",type
846
847        if(nSets == 0)
848                SetDataFolder root:
849                return(0)
850        endif
851
852
853//TODO: properly define the errors here - I'll have this if I do the simulation
854        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
855                Duplicate/O inten,iErr
856                Wave iErr=iErr
857//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
858                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
859        endif
860        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
861                Duplicate/O inten2,iErr2
862                Wave iErr2=iErr2
863//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
864                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
865        endif
866        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
867                Duplicate/O inten3,iErr3
868                Wave iErr3=iErr3
869//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
870                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
871        endif
872        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
873                Duplicate/O inten4,iErr4
874                Wave iErr4=iErr4
875//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
876                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
877        endif
878
879        nq = 600
880
881        // note that the back panel of 320x320 (1mm res) results in 447 data points!
882        // - so I upped nq to 600
883
884//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
885//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
886        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
887        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
888        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
889        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
890        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
891        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
892       
893        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
894        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
895        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
896        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
897        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
898        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
899       
900       
901//      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
902// TODO: not sure if I want to set dQ in x or y direction...
903        // the short dimension is the 8mm tubes, use this direction as dQ?
904        // 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]
905        // WRT the beam center. use qx or qy directly. Still not happy with this way...
906
907
908        qBin_qxqy[] =  p*delQ   
909        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
910
911        iBin_qxqy = 0
912        iBin2_qxqy = 0
913        eBin_qxqy = 0
914        eBin2D_qxqy = 0
915        nBin_qxqy = 0   //number of intensities added to each bin
916
917// now there are situations of:
918// 1 panel
919// 2 panels
920// 4 panels
921//
922// this needs to be a double loop now...
923
924// use set 1 (no number) only
925        if(nSets >= 1)
926                xDim=DimSize(inten,0)
927                yDim=DimSize(inten,1)
928       
929                for(ii=0;ii<xDim;ii+=1)
930                        for(jj=0;jj<yDim;jj+=1)
931                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
932                                qVal = qTotal[ii][jj]
933                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
934                                val = inten[ii][jj]
935                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
936                                        iBin_qxqy[binIndex] += val
937                                        iBin2_qxqy[binIndex] += val*val
938                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
939                                        nBin_qxqy[binIndex] += 1
940                                endif
941                        endfor
942                endfor
943               
944        endif
945
946// add in set 2 (set 1 already done)
947        if(nSets >= 2)
948                xDim=DimSize(inten2,0)
949                yDim=DimSize(inten2,1)
950       
951                for(ii=0;ii<xDim;ii+=1)
952                        for(jj=0;jj<yDim;jj+=1)
953                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
954                                qVal = qTotal2[ii][jj]
955                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
956                                val = inten2[ii][jj]
957                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
958                                        iBin_qxqy[binIndex] += val
959                                        iBin2_qxqy[binIndex] += val*val
960                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
961                                        nBin_qxqy[binIndex] += 1
962                                endif
963                        endfor
964                endfor
965               
966        endif
967
968// add in set 3 and 4 (set 1 and 2 already done)
969        if(nSets == 4)
970                xDim=DimSize(inten3,0)
971                yDim=DimSize(inten3,1)
972       
973                for(ii=0;ii<xDim;ii+=1)
974                        for(jj=0;jj<yDim;jj+=1)
975                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
976                                qVal = qTotal3[ii][jj]
977                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
978                                val = inten3[ii][jj]
979                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
980                                        iBin_qxqy[binIndex] += val
981                                        iBin2_qxqy[binIndex] += val*val
982                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
983                                        nBin_qxqy[binIndex] += 1
984                                endif
985                        endfor
986                endfor
987               
988               
989                xDim=DimSize(inten4,0)
990                yDim=DimSize(inten4,1)
991       
992                for(ii=0;ii<xDim;ii+=1)
993                        for(jj=0;jj<yDim;jj+=1)
994                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
995                                qVal = qTotal4[ii][jj]
996                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
997                                val = inten4[ii][jj]
998                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
999                                        iBin_qxqy[binIndex] += val
1000                                        iBin2_qxqy[binIndex] += val*val
1001                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1002                                        nBin_qxqy[binIndex] += 1
1003                                endif
1004                        endfor
1005                endfor
1006               
1007        endif
1008
1009
1010// after looping through all of the data on the panels, calculate errors on I(q),
1011// just like in CircSectAve.ipf
1012        for(ii=0;ii<nq;ii+=1)
1013                if(nBin_qxqy[ii] == 0)
1014                        //no pixels in annuli, data unknown
1015                        iBin_qxqy[ii] = 0
1016                        eBin_qxqy[ii] = 1
1017                        eBin2D_qxqy[ii] = NaN
1018                else
1019                        if(nBin_qxqy[ii] <= 1)
1020                                //need more than one pixel to determine error
1021                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1022                                eBin_qxqy[ii] = 1
1023                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1024                        else
1025                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1026                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1027                                avesq = iBin_qxqy[ii]^2
1028                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1029                                var = aveisq-avesq
1030                                if(var<=0)
1031                                        eBin_qxqy[ii] = 1e-6
1032                                else
1033                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1034                                endif
1035                                // and calculate as it is propagated pixel-by-pixel
1036                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1037                        endif
1038                endif
1039        endfor
1040       
1041        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1042       
1043        // find the last non-zero point, working backwards
1044        val=nq
1045        do
1046                val -= 1
1047        while((nBin_qxqy[val] == 0) && val > 0)
1048       
1049//      print val, nBin_qxqy[val]
1050        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1051
1052        if(val == 0)
1053                // all the points were deleted
1054                return(0)
1055        endif
1056       
1057       
1058        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1059        // find the first non-zero point, working forwards
1060        val = -1
1061        do
1062                val += 1
1063        while(nBin_qxqy[val] == 0)     
1064        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1065
1066        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1067        val = numpnts(nBin_qxqy)-1
1068        do
1069                if(nBin_qxqy[val] == 0)
1070                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1071                endif
1072                val -= 1
1073        while(val>0)
1074       
1075        SetDataFolder root:
1076       
1077        return(0)
1078End
1079
1080
Note: See TracBrowser for help on using the repository browser.