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

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

lots of changes to how VCALC simulations can be written to Nexus files to effectively make test data with different callues before the virtual machine is ready. many changes for visualization to effectively handle zeros in log scaled images

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: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: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.