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

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

more changes and additons to display VSANS data

adding functions for IvsQ plotting

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

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

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

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