9 Commits

19 changed files with 1447 additions and 443 deletions
+93
View File
@@ -0,0 +1,93 @@
#!/usr/bin/env python
from numpy import *
import mcstas_reader as mr
l=arange(3,11,0.01)
lp=(l[1:]+l[:-1])/2.
t=arange(0.0, 0.18, 0.0001)
tp=(t[1:]+t[:-1])/2.
def get_pulses(fn, freq, offset=0.002):
det=mr.McSim(fn)['tof_detector']
pulses=[]
for i in range(int(freq/7)):
li, Ii=det.project1d('L', bins=l, newcols=[('tc', 't-L*0.00406-%f'%offset)],
fltr='(tc>%f)&(tc<%f)'%(i*1./freq, (i+1)*1./freq))
ti, Iti=det.project1d('t', bins=t, newcols=[('tc', 't-L*0.00406-%f'%offset)],
fltr='(tc>%f)&(tc<%f)'%(i*1./freq, (i+1)*1./freq))
pulses.append(((li[1:]+li[:-1])/2., Ii, (ti[1:]+ti[:-1])/2., Iti))
return det, pulses
def get_filtered_pulses(det, pulses, freq, offset=0.0001, threashold=0.01):
Is=array([pi[3] for pi in pulses])
fpulses=[]
for i in range(1,int(freq/7)-1):
fltr=(Is[:i].sum(axis=0)+Is[i+1:].sum(axis=0))<(threashold*Is[i])
tf=tp[fltr]
if len(tf)<2:
continue
li, Ii=det.project1d('L', bins=l, fltr='(t>%f)&(t<%f)'%(tf[0], tf[-1]))
ti, Iti=det.project1d('t', bins=t, fltr='(t>%f)&(t<%f)'%(tf[0], tf[-1]))
fpulses.append(((li[1:]+li[:-1])/2., Ii, (ti[1:]+ti[:-1])/2., Iti))
return fpulses
def calc_usable(pulses, threashold=0.01):
Is=array([pi[3] for pi in pulses])
sI=0.
for i in range(1,Is.shape[0]):
sI+=Is[i][(Is[:i].sum(axis=0)+Is[i+1:].sum(axis=0))<(threashold*Is[i])].sum()
return sI
def calc_sum(pulses):
Is=array([pi[3] for pi in pulses])
return Is.sum(axis=0)
def tof_plot(det, pulses):
figure(figsize=(18,10))
for i, pi in enumerate(pulses):
semilogy(tp, pi[3], label='%i'%i)
ylim(100,1e5)
def tof_th_plot(det, pulses, threashold=0.01):
figure(figsize=(18,10))
Is=array([pi[3] for pi in pulses])
for i, pi in enumerate(pulses):
fltr=(Is[:i].sum(axis=0)+Is[i+1:].sum(axis=0))<(threashold*Is[i])
semilogy(tp[fltr], pi[3][fltr], label='%i'%i)
ylim(100,1e5)
def lambda_plot(det, pulses):
figure(figsize=(18,10))
for i, pi in enumerate(pulses):
semilogy(lp, pi[1])
semilogy(lp, det.project1d('L', bins=l)[1])
semilogy(lp, Ir)
ylim(100,1e6)
xlim(3,11)
def export(fn, pulses):
It=calc_sum(pulses)
savetxt(fn+'_lpulses.dat', array([pules[0][0], It, Ir]+[pi[1] for pi in pulses]).T)
savetxt(fn+'_tpulses.dat', array([pules[0][2]]+[pi[3] for pi in pulses]).T)
ref=mr.McSim('../results/ersa_ref/')['tof_detector']
Ir=ref.project1d('L', bins=l)[1]
freq_offsets={245: 0.00007, 280: 0.002, 308: 0.0015, 350: 0.003}
rfreq_offsets={245: 0.0000, 280: 0.0002, 308: 0.0015, 350: 0.0007}
results={}
for key,val in freq_offsets.items():
results['%i'%key]=get_pulses('../results/ersa_%iHz/'%key, 2.*key, val)
try:
results['r300-%i'%key]=get_pulses('../results/ersa_r300_%iHz/'%key, 2.*key, val)
except:
continue
for key,val in rfreq_offsets.items():
results['rev-%i'%key]=get_pulses('../results/ersa_reverse_%iHz/'%key, 2.*key, val)
+2 -2
View File
@@ -294,9 +294,9 @@ class Dataset(object):
limits=map(float, self.info['xylimits'].split())
if log:
ax.imshow(self.data, extent=limits, aspect='auto', norm=LogNorm())
ax.imshow(self.data, origin='lower', extent=limits, aspect='auto', norm=LogNorm())
else:
ax.imshow(self.data, extent=limits, aspect='auto')
ax.imshow(self.data, origin='lower', extent=limits, aspect='auto')
ax.set_xlabel(self.info['xlabel'])
ax.set_ylabel(self.info['ylabel'])
ax.set_title(self.info['component'])
+38 -164
View File
@@ -52,9 +52,11 @@ DEFINE INSTRUMENT ESS_reflectometer_Estia
(double omegaa = 1.2, int sample = 0, double sample_length = 0.01, double sample_height = 0.01,
int operationmode = 0, double over_illumination = 0.0, double theta_resolution = 0.04,
double lambda_min = 3.75, double lambda_start = 3.0, double lambda_end = 12.0,
int enable_chopper = 0, int enable_gravity=0, int enable_windows=1,
int enable_polarizer = 0, int enable_analyzer = 0,
double source_power = 2, double frame_usage = 0.98
int enable_chopper = 0, int enable_gravity=0, int enable_windows=1,
double wfm_freq=0.0, double wfm_shift=0.0, double wfm_radius=0.15,
int enable_polarizer = 0, int enable_analyzer = 0,
double source_power = 2,
double selene1_foot1y =0.0, double selene1_foot2y = 0.0
)
DECLARE
@@ -70,8 +72,12 @@ double iscs_z=0.0547095; // ISCS in McStas coordinates using mod-view-opt.instr
double iscs_rot_y=0.4; // ISCS is at TCS-35.6 degree, McStas at -36 degree
double iscs_rot_x=0.7; // downward tilt of Estia axis
// Selene 1 (Will follow later with detailed parameters)
// Selene 1 geometry parameters (optics parameters in Estia_selene1.instr)
double selene1_foot1 = 4.20; // distance of first foot to VS focus
double selene1_foot2 = 7.00; // distance of second foot to VS focus
double selene1_center;
double selene1_shift;
double selene1_rot;
// Selene 2
@@ -97,9 +103,6 @@ double slit_distance = 1.775; //m, distance slit to sample
double selene_theta = 1.25;
double velocity_max ; // m/s neutron velocity of lambda_min
/* flags */
int PP_small = 1 ; // plot all PSD with small area
int PP_large = 1 ; // plot all PSD with large area
double analyzer_max_height = 0.01; // Maximum sample height to be covered by te polarization analyzers
double analyzer1_start = 0.65; // Distance to sample to start the first analyzer
@@ -115,7 +118,7 @@ INITIALIZE
// chopper coupling together
velocity_max= 3.956034E3 / lambda_min; // h/m_n over lambda - ((3.9..e-7m^2/s))/(1e-10m)
chopper_phase=360.0*(chopper_pos*chopper_freq/velocity_max+(pulse_zero-opening_time)*chopper_freq);
chopper_open =360.0*(chopper_pos/(total_length+detector_arm))*frame_usage;
chopper_open = 98.0;
/* print out some calculated parameter for checking purposes */
@@ -130,7 +133,11 @@ Theta2_analyzer1=atan((dist_ana_vfocus+analyzer1_start+analyzer1_length)/dist_an
Theta1_analyzer2=atan((dist_ana_vfocus+analyzer2_start)/dist_ana_vfocus*analyzer_max_height/2.0/analyzer2_start)*180.0/PI;
Theta2_analyzer2=atan((dist_ana_vfocus+analyzer2_start+analyzer2_length)/dist_ana_vfocus*analyzer_max_height/2.0/(analyzer2_start+analyzer2_length))*180.0/PI;
selene1_center=(selene1_foot1+selene1_foot2)/2.;
selene1_shift=(selene1_foot1y+selene1_foot2y)/1000.0;
selene1_rot=atan((selene1_foot2y-selene1_foot1y)/(selene1_foot2-selene1_foot1)/1000.0)*180.0/PI;
printf(" Selene 1 rotation = %.4f deg\n", selene1_rot);
printf(" Selene 1 shift = %.1f mm\n", selene1_shift*1e3);
%}
TRACE
@@ -157,10 +164,12 @@ COMPONENT ISCS = Arm() // rotate around y-axis (slight downward tilt)
COMPONENT arm_feeder = Arm()
AT (0, 0, 0) RELATIVE ISCS
ROTATED (0, 0, 0) RELATIVE ISCS
COMPONENT arm_selene1_center = Arm()
AT (selene1_shift, 0, 2*NBOA_c+selene1_center) RELATIVE ISCS
ROTATED (0, selene1_rot, 0) RELATIVE ISCS
// Axes starting at focus of feeder parallel to c-axis of Selene guides
COMPONENT arm_selene1 = Arm()
AT (0, 0, 2*NBOA_c) RELATIVE ISCS
ROTATED (0, 0, 0) RELATIVE ISCS
AT (0, 0, -selene1_center) RELATIVE arm_selene1_center
COMPONENT arm_selene2 = Arm()
AT (0, 0, 2*NBOA_c+2*selene_c) RELATIVE ISCS
ROTATED (0, 0, 0) RELATIVE ISCS
@@ -199,15 +208,6 @@ COMPONENT moderator = ESS_butterfly(
AT (0, 0, 0) RELATIVE origin
// monitor used for brilliance transfer calculations
COMPONENT tof_source = Monitor_nD(
filename = "tof_source",
options = "xdiv limits=[-0.75 0.75] bins=15 ydiv limits=[-2.0 2.0] bins=40 time limits=[0 0.6] bins=600 lambda limits=[0 35] bins=350",
xwidth=0.01, yheight = 0.01)
WHEN sample==4
AT (0, 0, 0.03) RELATIVE ISCS
ROTATED (0,-selene_theta,0) RELATIVE ISCS
/***************************************
* Geometry of neutron feeder separate *
***************************************/
@@ -216,11 +216,14 @@ COMPONENT tof_source = Monitor_nD(
/****************************************************
* Beam manipulation area around the virtual source *
****************************************************/
/* Absorber to reduce beam to needed size an for shielding purposes (FeNi in CAD model) */
COMPONENT FeNi_in = Slit(xwidth=0.028, yheight=0.078)
AT (0, 0, -0.800) RELATIVE arm_virtual_source_beam
/* Absorber to reduce beam to needed size an for shielding purposes (CPC1 in CAD model) */
COMPONENT CPC1_in = Slit(xwidth=0.0303, yheight=0.0792)
AT (0, 0, -0.890) RELATIVE arm_virtual_source_beam
COMPONENT FeNi_out = Slit(xwidth=0.013, yheight=0.038)
COMPONENT CPC1_monitor = Slit(xwidth=0.016344, yheight=0.05547)
AT (0, 0, -0.3573) RELATIVE arm_virtual_source_beam
COMPONENT CPC1_out = Slit(xwidth=0.013, yheight=0.038)
AT (0, 0, -0.220) RELATIVE arm_virtual_source_beam
@@ -231,30 +234,6 @@ COMPONENT chopper = DiskChopper(radius=chopper_diameter/2.0, yheight=0.02,
WHEN enable_chopper==1
AT (0, 0, chopper_pos-2*NBOA_c) RELATIVE arm_virtual_source_beam
// Beam monitoring before VS
COMPONENT DP_before_virtual_source = DivPos_monitor(
filename = "DP_before_virtual_source" ,
nh = 200, ndiv = 200,
xwidth = 0.2, yheight = 0.2,
maxdiv_h = 6,
restore_neutron = 1)
AT (0, 0, -0.022) RELATIVE arm_virtual_source_beam
COMPONENT PP_small_before_virtual_source = PSD_monitor(
filename = "PP_small_before_virtual_source",
nx = 200, ny = 200,
xwidth = 0.1, yheight = 0.1,
restore_neutron = 1)
WHEN PP_small ==1
AT (0, 0, -0.021) RELATIVE arm_virtual_source_beam
COMPONENT PP_large_before_virtual_source = PSD_monitor(
filename = "PP_large_before_virtual_source",
nx = 200, ny = 200,
xmin = -1.8, xmax = 1.8, ymin = -1.8, ymax = 1.8,
restore_neutron = 1)
WHEN PP_large ==1
AT (0, 0, 2*NBOA_c-0.02) RELATIVE ISCS
/* The actual virtual source mask, two L-shaped absorbers (first top-right) */
@@ -263,66 +242,26 @@ COMPONENT virtual_source_TR = Slit(
WHEN sample!=4
AT (-over_illumination, 0, -0.5*sample_length) RELATIVE arm_virtual_source
// window to absorb neutron not passing through heavy collimation
COMPONENT virtual_source_HC = Slit(
xmin = -0.015, xmax = 0.015,
ymin = -0.015, ymax = 0.015)
// window to cut down to defined size for test setting
COMPONENT virtual_source_HC = Slit(xwidth=sample_length, yheight=sample_height)
WHEN sample==4
AT (0, 0, 0) RELATIVE arm_virtual_source
//
/* The actual virtual source mask, two L-shaped absorbers (second bottom-left) */
COMPONENT virtual_source_BL = Slit(
xmin = -1.0, xmax = 0.0, ymin = -sample_height/2-over_illumination*5, ymax = 1.0)
WHEN sample!=4
AT (over_illumination, 0, 0.5*sample_length) RELATIVE arm_virtual_source
// Beam monitoring after VS
COMPONENT tof_virtual_source = Monitor_nD(
filename = "tof_virtual_source",
options = "xdiv limits=[-0.75 0.75] bins=15 ydiv limits=[-2.0 2.0] bins=40 time limits=[0 0.6] bins=600 lambda limits=[0 35] bins=350",
xwidth=sample_length, yheight = sample_height)
WHEN sample==4
AT (0, 0, 0.019) RELATIVE arm_virtual_source_beam
COMPONENT DL_behind_virtual_source = DivLambda_monitor(
filename = "DL_behind_virtual_source",
nL = 200, nh = 200,
xwidth = 0.02, yheight = 0.06,
maxdiv_h = 3,
Lmin = lambda_start, Lmax = lambda_end,
restore_neutron = 1)
AT (0, 0, 0.02) RELATIVE arm_virtual_source_beam
COMPONENT DP_behind_virtual_source = DivPos_monitor(
filename = "DP_behind_virtual_source" ,
nh = 200, ndiv = 200,
xwidth = 0.2, yheight = 0.2,
maxdiv_h = 3,
restore_neutron = 1)
AT (0, 0, 0.021) RELATIVE arm_virtual_source_beam
COMPONENT PP_small_behind_virtual_source = PSD_monitor(
filename = "PP_small_behind_virtual_source",
nx = 200, ny = 200,
xwidth = 0.02, yheight = 0.05,
restore_neutron = 1)
WHEN PP_small ==1
AT (0, 0, 0.022) RELATIVE arm_virtual_source_beam
COMPONENT PP_large_behind_virtual_source = PSD_monitor(
filename = "PP_large_behind_virtual_source",
nx = 200, ny = 200,
xmin = -1.8, xmax = 1.8, ymin = -1.8, ymax = 1.8,
restore_neutron = 1)
WHEN PP_large ==1
AT (0, 0, 2*NBOA_c+0.023) RELATIVE ISCS
/*************************************
* Geometry of Selene guide separate *
*************************************/
%include "Estia_selene.instr"
%include "Estia_selene1.instr"
%include "Estia_mf.instr"
%include "Estia_selene2.instr"
/**********************************
* Optical components within cave *
@@ -343,60 +282,14 @@ COMPONENT ac_slit = Slit(
/***************
* Sample area *
***************/
// Beam monitoring before sample
COMPONENT PP_small_sample_focus = PSD_monitor(
filename = "PP_small_sample_focus",
nx = 200, ny = 200,
xwidth = 0.02, yheight = 0.05,
restore_neutron = 1)
WHEN PP_small ==1
AT (0, 0, -0.001) RELATIVE arm_sample
COMPONENT DL_sample_focus = DivLambda_monitor(
filename = "DL_sample_focus",
nL = 200, nh = 200,
xwidth = 0.02, yheight = 0.06,
maxdiv_h = 3,
Lmin = lambda_start, Lmax = lambda_end,
restore_neutron = 1)
AT (0, 0, -0.001) RELATIVE arm_sample
/* monitor at sample position for footprint */
COMPONENT PL_sample = Monitor_nD(
filename = "PL_sample",
options = "lambda limits=[0 35] bins=350 x limits=[-0.02 0.02] bins 160",
xwidth = 0.04, yheight = 0.02, restore_neutron=1)
WHEN sample!=4
AT (0, 0, -0.001) RELATIVE arm_sample
ROTATED (0, -90, 0) RELATIVE arm_sample
COMPONENT PSD_sample = PSD_monitor(
filename = "PSD_sample",
nx = 100, ny = 100,
xwidth = 4*sample_length, yheight = 4*sample_height,
restore_neutron = 1)
WHEN sample!=4
AT (0, 0, -0.001) RELATIVE arm_sample
ROTATED (0, -90, 0) RELATIVE arm_sample
COMPONENT tof_sample = Monitor_nD(
filename = "tof_sample",
options = "xdiv limits=[-0.75 0.75] bins=15 ydiv limits=[-2.0 2.0] bins=40 time limits=[0 0.6] bins=600 lambda limits=[0 35] bins=350",
xwidth=sample_length, yheight = sample_height)
options = "x limits=[-0.025 0.025] bins=1000 y limits=[-0.025 0.025] bins=1000 xdiv limits=[-0.75 0.75] bins=150 ydiv limits=[-2.0 2.0] bins=400 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 list all",
xwidth=0.05, yheight = 0.05)
WHEN sample==4
AT (0, 0, -0.0005) RELATIVE arm_sample_beam
ROTATED (0, 0, 0) RELATIVE arm_sample_beam
/* monitor for full beam at sample position ( n/cm²/s ) */
COMPONENT PP_small_sample = PSD_monitor(
filename = "PSD_sample_perp",
nx = 200, ny = 200,
xwidth = 0.05, yheight = 0.05,
restore_neutron = 1)
WHEN sample==4
AT (0, 0, 0) RELATIVE arm_sample_beam
ROTATED (0, 0, 0) RELATIVE arm_sample_beam
/* NiTi multilayer sample */
COMPONENT sample = Mirror(
@@ -463,31 +356,12 @@ COMPONENT analyzer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, ref
ROTATED (0,0,0.0) RELATIVE arm_analyzer
/* detector */
// Beam monitoring at detector position
COMPONENT PP_large_detector = PSD_monitor(
filename = "PP_large_detector",
nx = 200, ny = 200,
xmin = -1.8, xmax = 1.8, ymin = -1.8, ymax = 1.8,
restore_neutron = 1)
AT (0, 0, total_length+detector_arm-0.001) RELATIVE ISCS
// Detector with 0.5x1mm² resolution
COMPONENT PP_detector = PSD_monitor(
filename = "PP_detector",
nx = 1000, ny = 500,
xwidth = 0.5, yheight = 0.5,
restore_neutron = 1)
AT (0, 0, detector_arm) RELATIVE arm_detector
COMPONENT tof_detector = Monitor_nD(
filename = "tof_detector",
options = "x limits=[-0.25 0.25] bins=1000 y limits=[-0.5 0.5] bins=1000 time limits=[0 0.6] bins=6000 lambda limits=[0 35] bins=3500 sx limits=[-1 1] bins=1000 sy limits=[-1 1] bins=1000, list all",
xwidth = 0.5, yheight = 1.0)
WHEN sample!=4
AT (0, 0, detector_arm+0.00001) RELATIVE arm_detector
/*COMPONENT detector = PSD_TOF_monitor(xwidth=0.5, yheight=0.5,
nx=1000, ny=250, tmin=35000, tmax=115000, nt=160, filename="detector")
AT (0, 0, detector_arm) RELATIVE arm_detector
/***********************************************************************/
+172
View File
@@ -0,0 +1,172 @@
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Estia_mf
*
* %Identification
* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
* Date: 01. 03. 2018
* Origin: PSI
* Release: McStas 2.4.1
* Version: 1.0
* %INSTRUMENT_SITE: ESS (E02)
*
* Estia is a vertical sample, focusing reflectometer for small sample
*
* %Description
* These are the components near the middle focus between the Selene guide 1 and 2.
* This file is not intended for stand alone use but is included in the Estia instrument model.
* The removable components, however, allow for a stand alone use for debugging purpose.
*
* %Parameters
* enable_gravity: [0/1] Use gravity in elliptical guide model.
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT Estia_selene(int enable_gravity=0)
DECLARE
%{
double polarizer_max_width = 0.0015; // Maximum sample height to be covered by te polarization analyzers
double polarizer_start = 0.3; // Distance to sample to start the first analyzer
double polarizer_length = 0.6; // length of first analyzer
double Theta1_polarizer, Theta2_polarizer, dist_pol_vfocus; // quantities calculated out of values above and 1.5 degree covered divergence
%}
INITIALIZE
%{
dist_pol_vfocus=polarizer_max_width/2.0/tan(1.5*PI/360.0); // the virtual focus point infront of the actual sample focus where the beams furthest out meet
Theta1_polarizer=atan((dist_pol_vfocus+polarizer_start)/dist_pol_vfocus*polarizer_max_width/2.0/polarizer_start)*180.0/PI;
Theta2_polarizer=atan((dist_pol_vfocus+polarizer_start+polarizer_length)/dist_pol_vfocus*polarizer_max_width/2.0/(polarizer_start+polarizer_length))*180.0/PI;
%}
TRACE
REMOVABLE COMPONENT origin = Progress_bar()
AT (0,0,0) ABSOLUTE
REMOVABLE COMPONENT ISCS = Arm()
AT (0, 0, 0) RELATIVE origin
REMOVABLE COMPONENT arm_selene1 = Arm()
AT (0, 0, 0) RELATIVE ISCS
REMOVABLE COMPONENT arm_selene2 = Arm()
AT (0, 0, 2*selene_c) RELATIVE ISCS
REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002,
target_index=3, Emin=1, Emax=15)
AT (0,0,0) RELATIVE ISCS
ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
COMPONENT arm_ersa = Arm()
AT (0, 0, 0) RELATIVE arm_selene2
ROTATED (0, -selene_theta, 0) RELATIVE arm_selene2
/**************************************
* Middle focus between Selene guides *
**************************************/
COMPONENT arm_polarizer = Arm()
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
ROTATED (selene_theta, -selene_theta, 0) RELATIVE ISCS
COMPONENT polarizer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=-(polarizer_start+polarizer_length), length=polarizer_length,
delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
WHEN enable_polarizer
AT (0, 0, -(polarizer_start+polarizer_length)) RELATIVE arm_polarizer
ROTATED (0,0,90.0) RELATIVE arm_polarizer
COMPONENT ersa_01 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=60+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, -0.08) RELATIVE arm_ersa
COMPONENT ersa_02 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=75+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, -0.07) RELATIVE arm_ersa
COMPONENT ersa_03 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=90+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, -0.06) RELATIVE arm_ersa
COMPONENT ersa_04 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=105+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, -0.05) RELATIVE arm_ersa
COMPONENT ersa_05 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=120+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, -0.04) RELATIVE arm_ersa
COMPONENT ersa_05b = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=45, phase=-22.5+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq<0
AT (0, 0, -0.04) RELATIVE arm_ersa
COMPONENT ersa_11b = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=45, phase=22.5+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq<0
AT (0, 0, 0.04) RELATIVE arm_ersa
COMPONENT ersa_11 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=120+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, 0.04) RELATIVE arm_ersa
COMPONENT ersa_12 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=135+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, 0.05) RELATIVE arm_ersa
COMPONENT ersa_13 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=150+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, 0.06) RELATIVE arm_ersa
COMPONENT ersa_14 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=165+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, 0.07) RELATIVE arm_ersa
COMPONENT ersa_15 = DiskChopper(radius=wfm_radius, yheight=0.05,
theta_0=120, phase=180+wfm_shift*180,
nu=wfm_freq, nslit=2)
WHEN wfm_freq>0
AT (0, 0, 0.08) RELATIVE arm_ersa
/* The proximal polariser comes next */
COMPONENT polarizer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=polarizer_start, length=polarizer_length,
delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
WHEN (enable_polarizer>1)
AT (0, 0, polarizer_start) RELATIVE arm_polarizer
ROTATED (0.0,0,90.0) RELATIVE arm_polarizer
FINALLY
%{
%}
END
+204
View File
@@ -0,0 +1,204 @@
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: ESS_reflectometer_Estia
*
* %Identification
* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
* Date: 01. 03. 2018
* Origin: PSI
* Release: McStas 2.4.1
* Version: 1.0
* %INSTRUMENT_SITE: ESS (E02)
*
* Estia is a vertical sample, focusing reflectometer for small sample
*
* %Description
* The instrument consists of a two part elliptical feeder that focuses onto
* a slit mask that defines the sample footprint, called virtual source.
* The virtual source is imaged on the sample position with a Selene
* type neutron guide (two ellipses).
* This version of the instrument, used in the ESS butterfly moderator,
* has two vertical Selene guide systems that share the same focal points.
* The two guides are only implemented in the feeder as the selene mirrors
* are part of the upgrade program.
*
* %Parameters
* omegaa: [deg] sample rotation omega
* sample: flag to switch between (0) sample Ni/Ti-multilayer
* (1) reference ( R(q_z) = 1 for all q_z )
* (2) Ni-film on silicon
* (3) Natural SiO2 on silicon
* (4) monitor instead of sample (normal to beam)
* sample_length: [m] Size of sample in beam direction, also controls virtual source opening
* sample_height: [m] Size of sample in vertical direction, also controls virtual source opening
*
* operationmode: operation mode (0) high-intensity specular reflectivity
* (1) almost conventional TOF
* over_illumination: [m] Extra opening of virtual source compared to sample size
* theta_resolution: Delta theta / theta adjusted with the slit (operationmode=1)
*
* lambda_start: [A] Beginning of simulated wavelength range
* lambda_end: [A] End of simulated wavelength range
* enable_gravity: [0/1] Use gravity in elliptical guide model.
* enable_chopper: [0/1/2/3] Activate chopper component if !=0. Numbers larger than
* 1 define the number of pulses per chopper opening.
* (pulses skipped=enable_chopper-1)
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT ESS_reflectometer_Estia
(double lambda_start = 3.0, double lambda_end = 12.0,
int enable_gravity=1, int enable_windows=1, direct_beam=0,
double source_power = 5, foil_thickness=0.00001)
DECLARE
%{
/* Geometrical parameters from CAD model of Estia (ESS-0050413)
* TCS coordinate and directional rotation first focus point
* refered to as focus_moderator_y_rot
*/
//TCS position of ISCS: (110,-105,137)
double iscs_x=0.0199717;
double iscs_y=0.0;
double iscs_z=0.0547095; // ISCS in McStas coordinates using mod-view-opt.instr
double iscs_rot_y=0.4; // ISCS is at TCS-35.6 degree, McStas at -36 degree
double iscs_rot_x=0.7; // downward tilt of Estia axis
// Selene 1 geometry parameters (optics parameters in Estia_selene1.instr)
double selene1_foot1 = 4.20; // distance of first foot to VS focus
double selene1_foot2 = 7.00; // distance of second foot to VS focus
double selene1_center;
double selene1_shift;
double selene1_rot;
// Selene 2
/* general instrument geometry parameters */
double total_length = 35.0 ; // m distance moderator-sample (2*c_feeder+4*c_Selene)
double detector_arm = 4.0 ; // m distance sample-detector
/* chopper parameters */
double chopper_diameter = 0.7;
double chopper_pos = 10.895; // m distance source-chopper
double chopper_phase ; // deg phase between pulse and chopper opening
double chopper_open ; // deg of opening angle in the chopper
double pulse_zero = 0.00175; // ms intensity weighted average time of emitted neutron pulse
double opening_time = 0.0015; // ms time to reach full intensity on detector, used to adjust phase to get full intensity at beginning of selected band
double chopper_freq = 14.0 ; // Hz chopper frequency
double slit_distance = 1.775; //m, distance slit to sample
/* derived quantities */
double selene_theta = 1.25;
%}
INITIALIZE
%{
%}
TRACE
/********************************************************************************
* Initial geometry coordinate axes for important components of the simulation. *
********************************************************************************/
COMPONENT origin = Progress_bar()
AT (0,0,0) ABSOLUTE
/* The ISCS is the instrument coordinate system with directions where y points upwards and z lies on the instrument axes. */
COMPONENT ISCS_rot1 = Arm() // position correctly and rotate around z-axis (x points in beam direction)
AT (iscs_x,iscs_y,iscs_z) RELATIVE origin
ROTATED (0,iscs_rot_y,0) RELATIVE origin
COMPONENT ISCS = Arm() // rotate around y-axis (slight downward tilt)
AT (0,0,0) RELATIVE ISCS_rot1
ROTATED (iscs_rot_x,0,0) RELATIVE ISCS_rot1
/***********************
* Instrument Skeleton *
***********************/
// Axes system parallel to the c-axis of the feeder
COMPONENT arm_feeder = Arm()
AT (0, 0, 0) RELATIVE ISCS
ROTATED (0, 0, 0) RELATIVE ISCS
COMPONENT arm_selene1_center = Arm()
AT (selene1_shift, 0, 2*NBOA_c+selene1_center) RELATIVE ISCS
ROTATED (0, selene1_rot, 0) RELATIVE ISCS
// Axes starting at focus of feeder parallel to c-axis of Selene guides
COMPONENT arm_selene1 = Arm()
AT (0, 0, -selene1_center) RELATIVE arm_selene1_center
// Axes of the center of usable beam passing through the virtual source
COMPONENT arm_virtual_source_beam = Arm()
AT (0, 0, 0) RELATIVE arm_selene1
ROTATED (0, selene_theta, 0) RELATIVE ISCS
// 90 degree to main beam is the monitors detector position
COMPONENT arm_monitor = Arm()
AT (0, 0, -0.320) RELATIVE arm_virtual_source_beam
ROTATED (-90, 0, 0) RELATIVE arm_virtual_source_beam
/**********
* Source *
**********/
COMPONENT moderator = ESS_butterfly(
sector = "E", beamline = 2, yheight = 0.03, cold_frac = 0.9,
focus_xw = E02_01_01_Cu_in_xmax-E02_01_01_Cu_in_xmin+0.002,
focus_yh = E02_01_01_Cu_in_yheight+0.002,
target_index=4, Lmin = lambda_start, Lmax = lambda_end,
n_pulses = 1, acc_power=source_power)
AT (0, 0, 0) RELATIVE origin
/***************************************
* Geometry of neutron feeder separate *
***************************************/
%include "Estia_feeder.instr"
/****************************************************
* Beam manipulation area around the virtual source *
****************************************************/
/* Absorber to reduce beam to needed size an for shielding purposes (CPC1 in CAD model) */
COMPONENT CPC1_in = Slit(xwidth=0.0303, yheight=0.0792)
AT (0, 0, -0.890) RELATIVE arm_virtual_source_beam
COMPONENT CPC1_monitor = Slit(xwidth=0.016344, yheight=0.05547)
AT (0, 0, -0.3573) RELATIVE arm_virtual_source_beam
// As TOF detector is rectangular, use focus_r size to limit to actual
// detector size of 0.5'' diameter cylinder
COMPONENT Vanadium_Foil = Incoherent(focus_r=0.00635, p_interact=0.9,
xwidth=0.018, yheight=0.055, zdepth=foil_thickness,
sigma_abs=5.08, sigma_inc=5.08, Vc=13.827,
target_index=1)
WHEN direct_beam<2
AT (0, 0, -0.320) RELATIVE arm_virtual_source_beam
ROTATED (34.3, 0, 0) RELATIVE arm_virtual_source_beam
COMPONENT tof_monitor = TOFLambda_monitor(
filename = "monitor",
tmin=0, tmax=120000, nt=1200,
Lmin=0,Lmax=35,nL=350,
xwidth = 0.025, yheight = 0.025)
WHEN direct_beam==0
AT (0, 0, 0.1) RELATIVE arm_monitor
COMPONENT tof_direct = TOFLambda_monitor(
filename = "VS",
tmin=0, tmax=120000, nt=1200,
Lmin=0,Lmax=35,nL=350,
xwidth = 0.05, yheight = 0.05)
WHEN direct_beam>0
AT (0, 0, 0.08) RELATIVE arm_virtual_source_beam
FINALLY
%{
%}
END
-240
View File
@@ -1,240 +0,0 @@
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Estia_selene
*
* %Identification
* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
* Date: 01. 03. 2018
* Origin: PSI
* Release: McStas 2.4.1
* Version: 1.0
* %INSTRUMENT_SITE: ESS (E02)
*
* Estia is a vertical sample, focusing reflectometer for small sample
*
* %Description
* This is the beam extraction part of the instrument, it is not intended for use but is
* included in the Estia instrument model. The removable components, however,
* allow for a stand alone use for debugging purpose.
*
* %Parameters
* enable_gravity: [0/1] Use gravity in elliptical guide model.
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT Estia_selene(int enable_gravity=0)
DECLARE
%{
// Selene 1 (Will follow later with detailed parameters)
// Selene 2
/* selene parameters */
double selene_xi = 0.6 ; // length of neutron guide over 2*c of ellipses, xi
double selene_b = 0.1047; // m selene short half axis c
double selene_c = 6.0000; // m selene focus half distance c
double selene_coating = 3.8; // m-value of selene coating
/* other variables */
double selene_entry ;
double selene_length ;
double selene_distance ;
double polarizer_max_width = 0.0015; // Maximum sample height to be covered by te polarization analyzers
double polarizer_start = 0.3; // Distance to sample to start the first analyzer
double polarizer_length = 0.6; // length of first analyzer
double Theta1_polarizer, Theta2_polarizer, dist_pol_vfocus; // quantities calculated out of values above and 1.5 degree covered divergence
%}
INITIALIZE
%{
/* derived parameters for the selene guides */
selene_length = 2.0 * selene_c * selene_xi;
selene_distance = ( 1.0 - selene_xi ) * selene_c;
selene_entry = selene_b * sqrt( 1.0 - pow(selene_xi,2.0) );
dist_pol_vfocus=polarizer_max_width/2.0/tan(1.5*PI/360.0); // the virtual focus point infront of the actual sample focus where the beams furthest out meet
Theta1_polarizer=atan((dist_pol_vfocus+polarizer_start)/dist_pol_vfocus*polarizer_max_width/2.0/polarizer_start)*180.0/PI;
Theta2_polarizer=atan((dist_pol_vfocus+polarizer_start+polarizer_length)/dist_pol_vfocus*polarizer_max_width/2.0/(polarizer_start+polarizer_length))*180.0/PI;
%}
TRACE
REMOVABLE COMPONENT origin = Progress_bar()
AT (0,0,0) ABSOLUTE
REMOVABLE COMPONENT ISCS = Arm()
AT (0, 0, 0) RELATIVE origin
REMOVABLE COMPONENT arm_selene1 = Arm()
AT (0, 0, 0) RELATIVE ISCS
REMOVABLE COMPONENT arm_selene2 = Arm()
AT (0, 0, 2*selene_c) RELATIVE ISCS
REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002,
target_index=3, Emin=1, Emax=15)
AT (0,0,0) RELATIVE ISCS
ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
/**************************
* Selene 1 neutron guide *
**************************/
/* Absorber to cut direct view beam (Copper in CAD model) */
COMPONENT slit_before_selene_guide_1 = Slit(
xmin = 0, xmax = selene_entry+0.001,
ymin = 0, ymax = selene_entry+0.005)
AT (0, 0, selene_distance-0.01) RELATIVE arm_selene1
COMPONENT block_before_selene_guide_1 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.005,
ymax=selene_entry/4.0-0.005,
zmin=0.0, zmax=0.001)
AT (0, 0, selene_distance-0.0095) RELATIVE arm_selene1
/* Selene 1 elliptic guide */
COMPONENT selene_guide_1 = Elliptic_guide_gravity(
l=selene_length, dimensionsAt = "mid",
linyh = selene_distance, loutyh= selene_distance,
linxw = selene_distance, loutxw= selene_distance,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=selene_coating, mtop=selene_coating, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance) RELATIVE arm_selene1
/* Absorber to cut direct view beam (Copper in CAD model) */
COMPONENT slit_after_selene_guide_1 = Slit(
xmin = 0, xmax = selene_entry+0.005,
ymin = 0, ymax = selene_entry+0.01)
AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene1
COMPONENT block_after_selene_guide_1 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.015,
ymax=selene_entry/4.0-0.015,
zmin=0.0, zmax=0.001)
AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene1
/**************************************
* Middle focus between Selene guides *
**************************************/
COMPONENT arm_polarizer = Arm()
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
ROTATED (selene_theta, -selene_theta, 0) RELATIVE ISCS
COMPONENT polarizer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=-(polarizer_start+polarizer_length), length=polarizer_length,
delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
WHEN enable_polarizer
AT (0, 0, -(polarizer_start+polarizer_length)) RELATIVE arm_polarizer
ROTATED (0,0,90.0) RELATIVE arm_polarizer
COMPONENT polmon_mfocus=PolLambda_monitor(Lmin=lambda_start, Lmax=lambda_end, nL=30, npol=101, my=1, filename="polmon_mfocus")
WHEN enable_polarizer
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
ROTATED (0.0,0.0,0.0) RELATIVE arm_selene1
COMPONENT PP_small_before_middle_focus = PSD_monitor(
filename = "PP_small_before_middle_focus",
nx = 200, ny = 200,
xwidth = 0.1, yheight = 0.1,
restore_neutron = 1)
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
/* The proximal polariser comes next */
COMPONENT polarizer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=polarizer_start, length=polarizer_length,
delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
WHEN (enable_polarizer>1)
AT (0, 0, polarizer_start) RELATIVE arm_polarizer
ROTATED (0.0,0,90.0) RELATIVE arm_polarizer
/**************************
* Selene 2 neutron guide *
**************************/
/* Absorber to cut direct view beam (Copper in CAD model) */
COMPONENT slit_before_selene_guide_2 = Slit(
xmin = -selene_entry-0.005, xmax=0,
ymin = -selene_entry-0.01, ymax = 0)
AT (0, 0, selene_distance-0.002) RELATIVE arm_selene2
COMPONENT block_before_selene_guide_2 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.015,
ymax=selene_entry/4.0-0.015,
zmin=0.0, zmax=0.001)
AT (0, 0, selene_distance-0.0015) RELATIVE arm_selene2
/* Selene 2 elliptic guide first half */
COMPONENT selene_guide_21 = Elliptic_guide_gravity(
l=0.5*selene_length-0.001, dimensionsAt = "mid",
linyh = selene_distance, loutyh= selene_c+0.001,
linxw = selene_distance, loutxw= selene_c+0.001,
xwidth=selene_b*2, yheight=selene_b*2,
mright=selene_coating, mleft=0, mtop=0, mbottom=selene_coating,
enableGravity=enable_gravity)
AT (0, 0, selene_distance) RELATIVE arm_selene2
/* Absorber to cut direct view beam (Copper in CAD model) */
COMPONENT slit_within_selene_guide_2 = Slit(
xmin = -selene_b, xmax=-selene_b*0.45+0.005,
ymin = -selene_b, ymax=-selene_b*0.45+0.005)
AT (0, 0, 1.0*selene_c) RELATIVE arm_selene2
/*COMPONENT block_within_selene_guide_2 = Absorber(
xmin =-1, xmax=1,
ymin =-selene_b*0.25+0.001,
ymax = selene_b*0.25-0.001,
zmin=0.0, zmax=0.0001)
AT (0, 0, 3.0*selene_c+0.0001) RELATIVE arm_selene2*/
/* Selene 2 elliptic guide first half */
COMPONENT selene_guide_22 = Elliptic_guide_gravity(
l=0.5*selene_length-0.001, dimensionsAt = "mid",
linyh = selene_c+0.001, loutyh= selene_distance,
linxw = selene_c+0.001, loutxw= selene_distance,
xwidth=selene_b*2, yheight=selene_b*2,
mright=selene_coating, mleft=0, mtop=0, mbottom=selene_coating,
enableGravity=enable_gravity)
AT (0, 0, 1.0*selene_c+0.001) RELATIVE arm_selene2
/* Absorber to cut direct view beam (Copper in CAD model) */
COMPONENT slit_after_selene_guide_2 = Slit(
xmin = -selene_entry-0.001, xmax=0, xmax = 0.0,
ymin = -selene_entry-0.005, ymax = 0)
AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene2
COMPONENT block_after_selene_guide_2 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.015,
ymax=selene_entry/4.0-0.015,
zmin=0.0, zmax=0.001)
AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene2
REMOVABLE COMPONENT FocusMonitor = PSD_monitor(
filename = "Focus",
nx = 100, ny = 100,
xwidth = 0.005, yheight = 0.005,
restore_neutron = 1)
AT (0, 0, 4*selene_c) RELATIVE ISCS
ROTATED (0, 1.25, 0) RELATIVE ISCS
FINALLY
%{
%}
END
+270
View File
@@ -0,0 +1,270 @@
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Estia_selene1
*
* %Identification
* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
* Date: 01. 03. 2018
* Origin: PSI
* Release: McStas 2.4.1
* Version: 1.0
* %INSTRUMENT_SITE: ESS (E02)
*
* Estia is a vertical sample, focusing reflectometer for small sample
*
* %Description
* This is the Selene guide 1 part of the instrument, it is not intended for use but is
* included in the Estia instrument model. The removable components, however,
* allow for a stand alone use for debugging purpose.
*
* %Parameters
* enable_gravity: [0/1] Use gravity in elliptical guide model.
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT Estia_selene(int enable_gravity=0)
DECLARE
%{
// Selene 1
// Selene 2
/* selene parameters */
double selene_xi = 0.6 ; // length of neutron guide over 2*c of ellipses, xi
double selene_b = 0.1047; // m selene short half axis c
double selene_c = 6.0000; // m selene focus half distance c
/* Coating along the elliptic guide to allow sharp wavelenth cut-off */
double scoating_01 = 3.58; // coating at Selene 1 start/ Selene 2 end
double scoating_02 = 3.36;
double scoating_03 = 3.19;
double scoating_04 = 3.07;
double scoating_05 = 2.98;
double scoating_06 = 2.92;
double scoating_07 = 2.88;
double scoating_08 = 2.86; // coating in center
double scoating_09 = 2.86;
double scoating_10 = 2.87;
double scoating_11 = 2.91;
double scoating_12 = 2.98;
double scoating_13 = 3.07;
double scoating_14 = 3.20;
double scoating_15 = 3.38; // coating at Selene 1 end/ Selene 2 start
/* other variables */
double selene_entry ;
double selene_length ;
double selene_distance ;
double selene_segment ;
%}
INITIALIZE
%{
/* derived parameters for the selene guides */
selene_length = 2.0 * selene_c * selene_xi;
selene_distance = ( 1.0 - selene_xi ) * selene_c;
selene_entry = selene_b * sqrt( 1.0 - pow(selene_xi,2.0) );
selene_segment = selene_length / 15.;
%}
TRACE
REMOVABLE COMPONENT origin = Progress_bar()
AT (0,0,0) ABSOLUTE
REMOVABLE COMPONENT ISCS = Arm()
AT (0, 0, 0) RELATIVE origin
REMOVABLE COMPONENT arm_selene1 = Arm()
AT (0, 0, 0) RELATIVE ISCS
ROTATED (0, 0, 0) RELATIVE ISCS
REMOVABLE COMPONENT arm_selene2 = Arm()
AT (0, 0, 2*selene_c) RELATIVE ISCS
REMOVABLE COMPONENT source = Moderator(radius = 0.001, focus_xw = selene_entry+0.002, focus_yh = selene_entry+0.002,
target_index=3, Emin=1, Emax=15)
AT (0,0,0) RELATIVE ISCS
ROTATED (-1.25, 1.25, 0) RELATIVE ISCS
/**************************
* Selene 1 neutron guide *
**************************/
/* Absorber to cut direct view beam (Copper in CAD model) */
COMPONENT slit_before_selene_guide_1 = Slit(
xmin = 0, xmax = selene_entry+0.001,
ymin = 0, ymax = selene_entry+0.005)
AT (0, 0, selene_distance-0.01) RELATIVE arm_selene1
COMPONENT block_before_selene_guide_1 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.005,
ymax=selene_entry/4.0-0.005,
zmin=0.0, zmax=0.001)
AT (0, 0, selene_distance-0.0095) RELATIVE arm_selene1
/* Selene 1 elliptic guide */
COMPONENT E02_03_01 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance, loutyh= selene_distance+14*selene_segment,
linxw = selene_distance, loutxw= selene_distance+14*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_01, mtop=scoating_01, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance) RELATIVE arm_selene1
COMPONENT E02_03_02 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+1*selene_segment, loutyh= selene_distance+13*selene_segment,
linxw = selene_distance+1*selene_segment, loutxw= selene_distance+13*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_02, mtop=scoating_02, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_03 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+2*selene_segment, loutyh= selene_distance+12*selene_segment,
linxw = selene_distance+2*selene_segment, loutxw= selene_distance+12*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_03, mtop=scoating_03, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_04 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+3*selene_segment, loutyh= selene_distance+11*selene_segment,
linxw = selene_distance+3*selene_segment, loutxw= selene_distance+11*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_04, mtop=scoating_04, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_05 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+4*selene_segment, loutyh= selene_distance+10*selene_segment,
linxw = selene_distance+4*selene_segment, loutxw= selene_distance+10*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_05, mtop=scoating_05, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_06 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+5*selene_segment, loutyh= selene_distance+9*selene_segment,
linxw = selene_distance+5*selene_segment, loutxw= selene_distance+9*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_06, mtop=scoating_06, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_07 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+6*selene_segment, loutyh= selene_distance+8*selene_segment,
linxw = selene_distance+6*selene_segment, loutxw= selene_distance+8*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_07, mtop=scoating_07, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_08 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+7*selene_segment, loutyh= selene_distance+7*selene_segment,
linxw = selene_distance+7*selene_segment, loutxw= selene_distance+7*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_08, mtop=scoating_08, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_09 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+8*selene_segment, loutyh= selene_distance+6*selene_segment,
linxw = selene_distance+8*selene_segment, loutxw= selene_distance+6*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_09, mtop=scoating_09, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_10 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+9*selene_segment, loutyh= selene_distance+5*selene_segment,
linxw = selene_distance+9*selene_segment, loutxw= selene_distance+5*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_10, mtop=scoating_10, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_11 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+10*selene_segment, loutyh= selene_distance+4*selene_segment,
linxw = selene_distance+10*selene_segment, loutxw= selene_distance+4*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_11, mtop=scoating_11, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_12 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+11*selene_segment, loutyh= selene_distance+3*selene_segment,
linxw = selene_distance+11*selene_segment, loutxw= selene_distance+3*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_12, mtop=scoating_12, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_13 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+12*selene_segment, loutyh= selene_distance+2*selene_segment,
linxw = selene_distance+12*selene_segment, loutxw= selene_distance+2*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_13, mtop=scoating_13, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_14 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+13*selene_segment, loutyh= selene_distance+1*selene_segment,
linxw = selene_distance+13*selene_segment, loutxw= selene_distance+1*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_14, mtop=scoating_14, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene1
COMPONENT E02_03_15 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+14*selene_segment, loutyh= selene_distance,
linxw = selene_distance+14*selene_segment, loutxw= selene_distance,
xwidth=selene_b*2, yheight=selene_b*2,
mright=0, mleft=scoating_15, mtop=scoating_15, mbottom=0,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+14*selene_segment) RELATIVE arm_selene1
/* Absorber to cut direct view beam (Copper in CAD model) */
COMPONENT slit_after_selene_guide_1 = Slit(
xmin = 0, xmax = selene_entry+0.005,
ymin = 0, ymax = selene_entry+0.01)
AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene1
COMPONENT block_after_selene_guide_1 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.015,
ymax=selene_entry/4.0-0.015,
zmin=0.0, zmax=0.001)
AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene1
FINALLY
%{
%}
END
+210
View File
@@ -0,0 +1,210 @@
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Estia_selene2
*
* %Identification
* Written by: Artur Glavic (artur.glavic@psi.ch); Jochen Stahn (jochen.stahn@psi.ch); Christine Klauser (christine.klauser@psi.ch)
* Date: 01. 03. 2018
* Origin: PSI
* Release: McStas 2.4.1
* Version: 1.0
* %INSTRUMENT_SITE: ESS (E02)
*
* Estia is a vertical sample, focusing reflectometer for small sample
*
* %Description
* This is the second Selene guide part of the instrument, it is not intended for use but is
* included in the Estia instrument model. This file does not work no its own, as it
* requires the definitions specified in the Estia_selene1.instr DECLAR and INITIALIZE.
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT Estia_selene2()
DECLARE
%{
%}
INITIALIZE
%{
%}
TRACE
/**************************
* Selene 2 neutron guide *
**************************/
/* Absorber to cut direct view beam (Bor-Al in CAD model) */
COMPONENT slit_before_selene_guide_2 = Slit(
xmin = -selene_entry-0.005, xmax=0,
ymin = -selene_entry-0.01, ymax = 0)
AT (0, 0, selene_distance-0.002) RELATIVE arm_selene2
COMPONENT block_before_selene_guide_2 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.015,
ymax=selene_entry/4.0-0.015,
zmin=0.0, zmax=0.001)
AT (0, 0, selene_distance-0.0015) RELATIVE arm_selene2
/* Selene 2 elliptic guide first half */
COMPONENT E02_04_01 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance, loutyh= selene_distance+14*selene_segment,
linxw = selene_distance, loutxw= selene_distance+14*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_15, mleft=0, mtop=0, mbottom=scoating_15,
enableGravity=enable_gravity)
AT (0, 0, selene_distance) RELATIVE arm_selene2
COMPONENT E02_04_02 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+1*selene_segment, loutyh= selene_distance+13*selene_segment,
linxw = selene_distance+1*selene_segment, loutxw= selene_distance+13*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_14, mleft=0, mtop=0, mbottom=scoating_14,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_03 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+2*selene_segment, loutyh= selene_distance+12*selene_segment,
linxw = selene_distance+2*selene_segment, loutxw= selene_distance+12*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_13, mleft=0, mtop=0, mbottom=scoating_13,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_04 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+3*selene_segment, loutyh= selene_distance+11*selene_segment,
linxw = selene_distance+3*selene_segment, loutxw= selene_distance+11*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_12, mleft=0, mtop=0, mbottom=scoating_12,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_05 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+4*selene_segment, loutyh= selene_distance+10*selene_segment,
linxw = selene_distance+4*selene_segment, loutxw= selene_distance+10*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_11, mleft=0, mtop=0, mbottom=scoating_11,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_06 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+5*selene_segment, loutyh= selene_distance+9*selene_segment,
linxw = selene_distance+5*selene_segment, loutxw= selene_distance+9*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_10, mleft=0, mtop=0, mbottom=scoating_10,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_07 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+6*selene_segment, loutyh= selene_distance+8*selene_segment,
linxw = selene_distance+6*selene_segment, loutxw= selene_distance+8*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_09, mleft=0, mtop=0, mbottom=scoating_09,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_08 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+7*selene_segment, loutyh= selene_distance+7*selene_segment,
linxw = selene_distance+7*selene_segment, loutxw= selene_distance+7*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_08, mleft=0, mtop=0, mbottom=scoating_08,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_09 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+8*selene_segment, loutyh= selene_distance+6*selene_segment,
linxw = selene_distance+8*selene_segment, loutxw= selene_distance+6*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_07, mleft=0, mtop=0, mbottom=scoating_07,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_10 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+9*selene_segment, loutyh= selene_distance+5*selene_segment,
linxw = selene_distance+9*selene_segment, loutxw= selene_distance+5*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_06, mleft=0, mtop=0, mbottom=scoating_06,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_11 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+10*selene_segment, loutyh= selene_distance+4*selene_segment,
linxw = selene_distance+10*selene_segment, loutxw= selene_distance+4*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_05, mleft=0, mtop=0, mbottom=scoating_05,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_12 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+11*selene_segment, loutyh= selene_distance+3*selene_segment,
linxw = selene_distance+11*selene_segment, loutxw= selene_distance+3*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_04, mleft=0, mtop=0, mbottom=scoating_04,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_13 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+12*selene_segment, loutyh= selene_distance+2*selene_segment,
linxw = selene_distance+12*selene_segment, loutxw= selene_distance+2*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_03, mleft=0, mtop=0, mbottom=scoating_03,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_14 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+13*selene_segment, loutyh= selene_distance+1*selene_segment,
linxw = selene_distance+13*selene_segment, loutxw= selene_distance+1*selene_segment,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_02, mleft=0, mtop=0, mbottom=scoating_02,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene2
COMPONENT E02_04_15 = Elliptic_guide_gravity(
l=selene_segment, dimensionsAt = "mid",
linyh = selene_distance+14*selene_segment, loutyh= selene_distance,
linxw = selene_distance+14*selene_segment, loutxw= selene_distance,
xwidth=selene_b*2, yheight=selene_b*2,
mright=scoating_01, mleft=0, mtop=0, mbottom=scoating_01,
enableGravity=enable_gravity)
AT (0, 0, selene_distance+14*selene_segment) RELATIVE arm_selene2
/* Absorber to cut direct view beam (Bor-Al in CAD model) */
COMPONENT slit_after_selene_guide_2 = Slit(
xmin = -selene_entry-0.001, xmax=0, xmax = 0.0,
ymin = -selene_entry-0.005, ymax = 0)
AT (0, 0, 2*selene_c-selene_distance+0.001) RELATIVE arm_selene2
COMPONENT block_after_selene_guide_2 = Absorber(
xmin=-1, xmax=1,
ymin=-selene_entry/4.0+0.015,
ymax=selene_entry/4.0-0.015,
zmin=0.0, zmax=0.001)
AT (0, 0, 2*selene_c-selene_distance+0.0015) RELATIVE arm_selene2
FINALLY
%{
%}
END
+4 -2
View File
@@ -2,8 +2,10 @@
if [ Estia_baseline.instr -nt Estia_baseline.out ] || [ ! -f Estia_baseline.out ] \
|| [ Estia_feeder.instr -nt Estia_baseline.out ] \
|| [ Estia_selene.instr -nt Estia_baseline.out ]; then
|| [ Estia_selene1.instr -nt Estia_baseline.out ] \
|| [ Estia_mf.instr -nt Estia_baseline.out ] \
|| [ Estia_selene2.instr -nt Estia_baseline.out ]; then
rm Estia_baseline.c Estia_baseline.out
mcstas -o Estia_baseline.c Estia_baseline.instr
mpicc -O3 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
mpicc -O3 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI
fi
Binary file not shown.
+12
View File
@@ -0,0 +1,12 @@
#!/bin/bash
mcstas -o Estia_monitor.c Estia_monitor.instr
mpicc -O3 -o Estia_monitor.out Estia_monitor.c -lm -DUSE_MPI
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_ref --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=2
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_trans --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=1
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor100 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.0001
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor001 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.000001
+41 -35
View File
@@ -2,45 +2,51 @@
#-*- coding: utf8 -*-
import sys, os
import numpy as np
from subprocess import call
from numpy import *
from scipy.optimize import leastsq
FOLDER='../results2/'
PREFIX='selene_geo_'
SEED=3798
ITEMS=1000
CALL='mpirun -np 6 Selene_geometry.out --dir=%s --ncount=4e8 enable_gravity=0 sample_length=0.003 '
CALL+='tx_1=%f tz_1=%f ry_1=%f rz_1=%f tx_2=%f tz_2=%f ry_2=%f rz_2=%f '
CALL='/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out '
CALL=+'--dir=../results/selene_geo --ncount=1e9 '
CALL=+'omegaa=1.0 sample=4 sample_length=0.001 sample_height=0.01 '
CALL=+'lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 '
CALL=+'lambda_min=3.75 selene1_foot1y=%.4f selene1_foot2y=%.4f '
TX=[-0.001, 0.001]
TZ=[-0.005, 0.005]
RY=[-0.001, 0.001]
RZ=[-2.5, 2.5]
def B(x,w):
# box function with full width w
return float32(abs(x)<=(w/2.))
def G(x,I0,x0,sigma):
# Gaussian with intensity I0, center x0 and standard deviation sigma
return I0*exp(-0.5*(x-x0)**2/sigma**2)
def Intensity(x, p):
I0, x0, sigma, w=p
return convolve(B(xc,w), G(xc, I0, x0, sigma), mode='same')
def Beam(x,I0,x0,sigma,w):
I=I0*where((x-x0)<(-w/2.), exp(-0.5*(x-x0+w/2.)**2/sigma**2),
where((x-x0)>(w/2.), exp(-0.5*(x-x0-w/2.)**2/sigma**2), 1.))
return I
def residuals(p, x, y):
I0, x0, sigma, w=p
return y-Beam(x, I0, x0, sigma, w)
def FWHM(pi):
rng=xc[where(Beam(xc, *pi)>=(pi[0]/2.))[0]]
return rng[-1]-rng[0]
x=linspace(-10.0005, 10.0005, 20001)
xc=(x[1:]+x[:-1])/2.
def analyze(fi):
sim=mr.McSim(fi)
det=sim['tof_sample']
ignore,y=det.project1d('x', bins=x/1000.)
p,res=leastsq(residuals, (y.max(), (xc*y).sum()/y.sum(), 0.01, 0.1), (xc, y))
return p,det
if __name__=='__main__':
call('mcstas -t -o Selene_geometry.c Selene_geometry.instr'.split())
call('mpicc -O2 -o Selene_geometry.out Selene_geometry.c -lm -DUSE_MPI'.split())
np.random.seed(SEED)
fh=open(os.path.join(FOLDER, PREFIX+'sims.dat'), 'w')
fh.write('# index tx_1 tz_1 ry_1 rz_1 tx_2 tz_2 ry_2 rz_2\n')
fh.write('# seed = %i\n'%SEED)
for i in range(ITEMS):
tx_1=np.random.random()*(TX[1]-TX[0])+TX[0]
tx_2=np.random.random()*(TX[1]-TX[0])+TX[0]
tz_1=np.random.random()*(TZ[1]-TZ[0])+TZ[0]
tz_2=np.random.random()*(TZ[1]-TZ[0])+TZ[0]
ry_1=np.random.random()*(RY[1]-RY[0])+RY[0]
ry_2=np.random.random()*(RY[1]-RY[0])+RY[0]
rz_1=np.random.random()*(RZ[1]-RZ[0])+RZ[0]
rz_2=np.random.random()*(RZ[1]-RZ[0])+RZ[0]
ln=' '.join(['%f'%val for val in [tx_1, tz_1, ry_1, rz_1, tx_2, tz_2, ry_2, rz_2]])
fh.write(ln+'\n')
print i, ':', ln
call((CALL%(os.path.join(FOLDER, PREFIX+'%05i/'%i), tx_1, tz_1, ry_1, rz_1, tx_2, tz_2, ry_2, rz_2)).split())
fh.close()
pass
+177
View File
@@ -0,0 +1,177 @@
#!/bin/bash
DEST=../results
ncount=2e10
use_cores=$SLURM_NPROCS
sample=1
omega=5.0
sample_length=0.01
sample_height=0.01
lambda_start=3.00
lambda_end=15.0
######################## 1% WFM for 30cm and 60cm device #######################
DESTi=$DEST/ersa_245Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=245 wfm_shift=0.0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
DESTi=$DEST/ersa_280Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=280 wfm_shift=0.0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
DESTi=$DEST/ersa_308Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=308 wfm_shift=0.0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
DESTi=$DEST/ersa_350Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=350 wfm_shift=0.0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
DESTi=$DEST/ersa_r300_245Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=245 wfm_shift=0.0 wfm_radius=0.3 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
DESTi=$DEST/ersa_r300_280Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=280 wfm_shift=0.0 wfm_radius=0.3 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
DESTi=$DEST/ersa_r300_308Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=308 wfm_shift=0.0 wfm_radius=0.3 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
DESTi=$DEST/ersa_r300_350Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=350 wfm_shift=0.0 wfm_radius=0.3 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
###################### 0.5% WFM ########################
ncount=4e10
DESTi=$DEST/ersa_reverse_280Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=-280 wfm_shift=0.0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1
DESTi=$DEST/ersa_reverse_308Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=-308 wfm_shift=0.0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1
DESTi=$DEST/ersa_reverse_350Hz
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=-350 wfm_shift=0.0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1
##################### Reference run ###############################
ncount=1e10
DESTi=$DEST/ersa_ref
if [ -e "$DESTi" ]; then
rm -r "$DESTi"
fi
mpirun -np $use_cores Estia_baseline.out \
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.001 \
sample=$sample sample_length=$sample_length sample_height=$sample_height \
enable_windows=0 wfm_freq=0 \
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
+24
View File
@@ -0,0 +1,24 @@
#!/bin/tcsh
#SBATCH -J wfm_McEstia
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout.log
#SBATCH -e stderr.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
module load mcstas
bash compile_if_needed.sh
bash run_wfm.sh
+40
View File
@@ -0,0 +1,40 @@
#!/bin/bash
#SBATCH -J mcEstia_1
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout.log
#SBATCH -e stderr.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq -0.009 0.001 0.009)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
+40
View File
@@ -0,0 +1,40 @@
#!/bin/bash
#SBATCH -J mcEstia_2
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout2.log
#SBATCH -e stderr2.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq -0.09 0.01 -0.01)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
+40
View File
@@ -0,0 +1,40 @@
#!/bin/bash
#SBATCH -J mcEstia_3
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout3.log
#SBATCH -e stderr3.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq 0.01 0.01 0.09)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
+40
View File
@@ -0,0 +1,40 @@
#!/bin/bash
#SBATCH -J mcEstia_4
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout4.log
#SBATCH -e stderr4.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq -0.5 0.1 -0.1)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"
+40
View File
@@ -0,0 +1,40 @@
#!/bin/bash
#SBATCH -J mcEstia_5
#SBATCH -N 2
#SBATCH --ntasks-per-node=24
#SBATCH --time=1-00:00:00
#SBATCH --mail-type=fail
#SBATCH --mail-user=artur.glavic@psi.ch
#SBATCH -o stdout5.log
#SBATCH -e stderr5.log
#SBATCH --partition=ll_long
echo "Starting at `date`"
echo "Running on hosts: $SLURM_NODELIST"
echo "Running on $SLURM_NNODES nodes."
echo "Running on $SLURM_NPROCS processors."
echo "Current working directory is `pwd`"
/usr/bin/modulecmd tcsh load mcstas
bash compile_if_needed.sh
for y1 in $(seq 0.1 0.1 0.5)
do
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
do
printf -v Y1 %.3f $y1
printf -v Y2 %.3f $y2
echo $Y1 $Y2
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
done
done
echo "Program finished with exit code $? at: `date`"