/*******************************************************************************
* Simple simulation of PRISMA2 with RITA-style analyser backend.
*
* Written by Kristian Nielsen and Mark Hagen August 1998.
*
* Demonstrates how the standard components from the component library
* may be easily modified for special purposes; in this case to have
* the individual analyser blades paint a "color" on the neutrons to
* differentiate them in the detector.
*
* Output is in the file "sim/prisma2.tof". The format is ASCII; each
* line consists of the time-of-flight in microseconds followed by seven
* intensities of neutrons from each individual analyser blade.
*
* Examples:
*
* prisma2 --ncount=2e6 TT=-30 PHA=22 PHA1=-3 PHA2=-2 PHA3=-1 PHA4=0 PHA5=1 PHA6=2 PHA7=3 TTA=44
* prisma2 --ncount=2e6 TT=-30 PHA=22 PHA1=3 PHA2=2 PHA3=1 PHA4=0 PHA5=-1 PHA6=-2 PHA7=-3 TTA=44
*******************************************************************************/
/* Modified from Monochromator.comp to paint a "color" on the neutron
if it is scattered. */
DEFINE COMPONENT Monochromator_color
DEFINITION PARAMETERS (zmin, zmax, ymin, ymax, mosaich, mosaicv, r0, Q, color)
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
%}
TRACE
%{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
%}
END
/* Modified from TOF_monitor.comp to bin neutrons according to their
"color". */
DEFINE COMPONENT TOF_monitor_color
DEFINITION PARAMETERS (xmin, xmax, ymin, ymax, nchan, dt, filename, maxcolor)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (TOF_N, TOF_p, TOF_p2)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
int TOF_N[maxcolor+1][nchan];
double TOF_p[maxcolor+1][nchan];
double TOF_p2[maxcolor+1][nchan];
%}
INITIALIZE
%{
int i,c;
for (i=0; i<nchan; i++)
for (c=0; c<=maxcolor; c++)
{
TOF_N[c][i] = 0;
TOF_p[c][i] = 0;
TOF_p2[c][i] = 0;
}
%}
TRACE
%{
int i;
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
i = floor(1E6*t/dt); /* Bin number */
if(i >= nchan) i = nchan;
if(i < 0)
{
printf("FATAL ERROR: negative time-of-flight.\n");
exit(1);
}
if(neu_color < 0 || neu_color > maxcolor)
{
printf("FATAL ERROR: wrong color neutron.\n");
exit(1);
}
TOF_N[neu_color][i]++;
TOF_p[neu_color][i] += p;
TOF_p2[neu_color][i] += p*p;
}
%}
FINALLY
%{
int i,c;
FILE *outfile;
outfile=fopen(filename,"w");
for (i=0; i<nchan; i++)
{
int empty = 1;
for(c=0; c<=maxcolor; c++)
if(TOF_p[c][i] != 0.0)
{
empty = 0;
break;
}
if(empty)
continue;
fprintf(outfile,"%g", (double)i*dt);
for(c=0; c<=maxcolor; c++)
fprintf(outfile," %g", TOF_p[c][i]);
fprintf(outfile,"\n");
}
fclose(outfile);
%}
END
DEFINE INSTRUMENT
prisma2(TT,PHA,PHA1,PHA2,PHA3,PHA4,PHA5,PHA6,PHA7,TTA)
DECLARE
%{
int neu_color; /* "Color" of current neutron */
/* 30' mosaicity used on analysator */
double prisma_ana_mosaic = 30;
/* Q vector for bragg scattering with monochromator and analysator */
double prisma_ana_q = 1.87325;
double prisma_ana_r0 = 0.6;
double focus_x,focus_z;
double apos1, apos2, apos3, apos4, apos5, apos6, apos7;
%}
INITIALIZE
%{
focus_x = 0.52 * sin(TT*DEG2RAD);
focus_z = 0.52 * cos(TT*DEG2RAD);
/* Rita-style analyser. */
{
double l = 0.0125;
apos1 = -3*l;
apos2 = -2*l;
apos3 = -1*l;
apos4 = 0*l;
apos5 = 1*l;
apos6 = 2*l;
apos7 = 3*l;
}
%}
TRACE
COMPONENT mod = Moderator(
radius = 0.0707,
dist = 9.035,
xw = 0.021,
yh = 0.021,
E0 = 10, E1 = 15,
Ec = 9.0, t0 = 37.15, gam = 39.1)
AT (0,0,0) ABSOLUTE
COMPONENT modslit = Slit(xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05)
AT(0,0,0) RELATIVE mod ROTATED (0,0,0) RELATIVE mod
COMPONENT tof_test = TOF_monitor(xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05,
nchan = 10000, dt = 10,
filename = "sim/prisma2.mon")
AT (0,0,0.005) RELATIVE mod ROTATED (0,0,0) RELATIVE mod
COMPONENT mon1 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
AT(0,0,0.01) RELATIVE mod ROTATED (0,0,0) RELATIVE mod
COMPONENT slit1 = Slit(xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05)
AT(0,0,1.7) RELATIVE mod ROTATED (0,0,0) RELATIVE mod
COMPONENT slit2 = Slit(xmin = -0.02, xmax = 0.02,
ymin = -0.03, ymax = 0.03)
AT(0,0,7) RELATIVE slit1 ROTATED (0,0,0) RELATIVE mod
COMPONENT mon2 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
AT(0,0,9) RELATIVE mod ROTATED (0,0,0) RELATIVE mod
COMPONENT sample = V_sample(
radius_i = 0.00001, radius_o = 0.01,
h = 0.02,
focus_r = 0.03,
pack = 1,
target_x = focus_x, target_y = 0, target_z = focus_z)
AT (0, 0, 9.035) RELATIVE mod ROTATED (0,0,0) RELATIVE mod
COMPONENT a2 = Arm() AT (0,0,0) RELATIVE sample ROTATED (0,TT,0) RELATIVE sample
COMPONENT mon3 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
AT(0,0,0.39) RELATIVE a2 ROTATED (0,0,0) RELATIVE a2
COMPONENT coll2 = Soller(xmin = -0.015, xmax = 0.015,
ymin = -0.025, ymax = 0.025,
len = 0.12, divergence = 120)
AT(0,0,0.40) RELATIVE a2 ROTATED (0,0,0) RELATIVE a2
COMPONENT mon4 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
AT(0,0,0.521) RELATIVE a2 ROTATED (0,0,0) RELATIVE a2
COMPONENT rita_ana = Arm()
AT(0, 0, 0.58) relative a2 ROTATED (0, PHA, 0) RELATIVE a2
COMPONENT ana1 = Monochromator_color(
ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
r0=prisma_ana_r0, Q=prisma_ana_q, color = 0)
AT (0, 0, apos1) RELATIVE rita_ana
ROTATED (0, PHA1, 0) RELATIVE rita_ana
COMPONENT ana2 = Monochromator_color(
ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
r0=prisma_ana_r0, Q=prisma_ana_q, color = 1)
AT (0, 0, apos2) RELATIVE rita_ana
ROTATED (0, PHA2, 0) RELATIVE rita_ana
COMPONENT ana3 = Monochromator_color(
ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
r0=prisma_ana_r0, Q=prisma_ana_q, color = 2)
AT (0, 0, apos3) RELATIVE rita_ana
ROTATED (0, PHA3, 0) RELATIVE rita_ana
COMPONENT ana4 = Monochromator_color(
ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
r0=prisma_ana_r0, Q=prisma_ana_q, color = 3)
AT (0, 0, apos4) RELATIVE rita_ana
ROTATED (0, PHA4, 0) RELATIVE rita_ana
COMPONENT ana5 = Monochromator_color(
ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
r0=prisma_ana_r0, Q=prisma_ana_q, color = 4)
AT (0, 0, apos5) RELATIVE rita_ana
ROTATED (0, PHA5, 0) RELATIVE rita_ana
COMPONENT ana6 = Monochromator_color(
ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
r0=prisma_ana_r0, Q=prisma_ana_q, color = 5)
AT (0, 0, apos6) RELATIVE rita_ana
ROTATED (0, PHA6, 0) RELATIVE rita_ana
COMPONENT ana7 = Monochromator_color(
ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
r0=prisma_ana_r0, Q=prisma_ana_q, color = 6)
AT (0, 0, apos7) RELATIVE rita_ana
ROTATED (0, PHA7, 0) RELATIVE rita_ana
COMPONENT a3 = Arm()
AT (0,0,0) relative rita_ana ROTATED (0,TTA,0) RELATIVE a2
COMPONENT mon5 = Monitor(xmin = -0.05, xmax = 0.05, ymin = -0.05, ymax = 0.05)
AT(0,0,0.06) RELATIVE a3 ROTATED (0,0,0) RELATIVE a3
COMPONENT mon6 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
AT(0,0,0.161) RELATIVE a3 ROTATED (0,0,0) RELATIVE a3
COMPONENT psd = PSD_monitor(xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05,
nx = 100, ny = 100,
filename = "sim/prisma2.psd")
AT(0,0,0.20) RELATIVE a3 ROTATED (0,0,0) RELATIVE a3
COMPONENT detector = TOF_monitor_color(xmin = -0.05, xmax = 0.05,
ymin = -0.05, ymax = 0.05,
nchan = 10000, dt = 10, maxcolor = 6,
filename = "sim/prisma2.tof")
AT (0,0,0.20) RELATIVE a3 ROTATED (0,0,0) RELATIVE a3
COMPONENT mon9 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
AT(0,0,0.01) RELATIVE detector ROTATED (0,0,0) RELATIVE a3
END
Legal notice