Compare commits
8 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 glavic_a
|
8434fecdde | ||
|
glavic_a
|
9697c9a584 | ||
|
glavic_a
|
0cc9f0ae5d | ||
|
glavic_a
|
00bd8b18c0 | ||
 ll_glavic
|
7a49a0f22c | ||
|
glavic_a
|
a58b2d3030 | ||
|
glavic_a
|
6dee3c98c1 | ||
|
glavic_a
|
8c10d63987 |
@@ -1,9 +0,0 @@
|
||||
#!/bin/bash
|
||||
# compress McStas hdf5 files to increase performance and save disk space
|
||||
|
||||
for fi in $*
|
||||
do
|
||||
echo $fi
|
||||
h5repack -i $fi -o $fi.c -f /entry1/data/of_detector_list_p_x_y_t_L_sx_sy/events:GZIP=5 -l /entry1/data/of_detector_list_p_x_y_t_L_sx_sy/events:CHUNK=3072x7
|
||||
mv $fi.c $fi
|
||||
done
|
||||
+137
-26
@@ -10,36 +10,147 @@ seterr(all='ignore')
|
||||
import estia_help as eh
|
||||
import mcstas_reader as mr
|
||||
|
||||
# work around for wrong wavelength calculation in pulse skipping mode
|
||||
eh.l_code='((t-0.0335)%%(%f/14.)+0.0335-0.0015)*101.436605'
|
||||
|
||||
# Scaling factor of source monitor after normalizing by area to get to brilliance.
|
||||
# Converts from 1.5x4.0 degree² acceptance range to 1/steradian
|
||||
# Factor 10 due to 0.1 Angstrom binning.
|
||||
BT_SOURCE_SCALE=10.*(180.**2/pi**2)/1.5/4.0
|
||||
BT_SAMPLE_SCALE=10.*(180.**2/pi**2)/1.5/1.5
|
||||
|
||||
|
||||
if __name__=='__main__':
|
||||
if not os.path.exists('../results/mooc_analyzed'):
|
||||
os.mkdir('../results/mooc_analyzed')
|
||||
if not os.path.exists('../results/analyzed'):
|
||||
os.mkdir('../results/analyzed')
|
||||
|
||||
if len(sys.argv)>1:
|
||||
only_items=map(int, sys.argv[1:])
|
||||
else:
|
||||
only_items=range(100)
|
||||
|
||||
print "pulse skipping mode"
|
||||
opts=dict(qres=0.02, qmin=0.007, qmax=0.2, mindq=2e-5, crop_overlap=False)
|
||||
names=['q_z', 'Intensity(ToF)', 'Reflectivity(ToF)',
|
||||
'1s-Stat(ToF)', '1pulse-Stat(ToF)', '10s-Statistics(ToF)',
|
||||
'Intensity(λ)', 'Reflectivity(λ)']
|
||||
units=['A^{-1}', 'cts/s', '1', 'cts/s', '1']
|
||||
info='Reflectivity simulation for 50x10mm² sample, 2-pulse skipping'
|
||||
if 1 in only_items:
|
||||
print "Brilliance Transfer 5x10"
|
||||
bt_10=mr.McSim('../results/brilliance_5x10/mccode.h5')
|
||||
d=bt_10['tof_sample.L']
|
||||
vs=bt_10['tof_virtual_source.L']
|
||||
r=bt_10['tof_source.L']
|
||||
x, y, e=eh.calc_brilliance_transfer(d, BT_SAMPLE_SCALE*2., r, BT_SOURCE_SCALE)
|
||||
savetxt('../results/analyzed/brilliance_transfer_5x10.dat', array([x, y, e]).T)
|
||||
savetxt('../results/analyzed/brilliance_5x10.dat', array([x, d.data*BT_SAMPLE_SCALE*2.,
|
||||
d.errors*BT_SAMPLE_SCALE*2.]).T)
|
||||
x, y, e=eh.calc_brilliance_transfer(vs, BT_SAMPLE_SCALE, r, BT_SOURCE_SCALE)
|
||||
savetxt('../results/analyzed/brilliance_transfer_vs_5x10.dat', array([x, y, e]).T)
|
||||
|
||||
btn_10=mr.McSim('../results/brilliance_nowindow_5x10/mccode.h5')
|
||||
d=btn_10['tof_sample.L']
|
||||
vs=btn_10['tof_virtual_source.L']
|
||||
r=btn_10['tof_source.L']
|
||||
x, y, e=eh.calc_brilliance_transfer(d, BT_SAMPLE_SCALE*2., r, BT_SOURCE_SCALE)
|
||||
savetxt('../results/analyzed/brilliance_transfer_nowindow_5x10.dat', array([x, y, e]).T)
|
||||
x, y, e=eh.calc_brilliance_transfer(vs, BT_SAMPLE_SCALE, r, BT_SOURCE_SCALE)
|
||||
savetxt('../results/analyzed/brilliance_transfer_nowindow_vs_5x10.dat', array([x, y, e]).T)
|
||||
|
||||
ref=mr.McSim('../results/mooc_reference_ps')['tof_detector']
|
||||
|
||||
for i in range(6):
|
||||
print "dataset",i
|
||||
sample=mr.McSim('../results/mooc_model_ps_%i'%i)['tof_detector']
|
||||
q, I, R=eh.calcR(sample, ref, 1.5, 1.5, skip_pulses=2, **opts)
|
||||
q, I_real, R_real=eh.calcR(sample, ref, 1.5, 1.5, skip_pulses=2, use_tof=False, **opts)
|
||||
idx_start, idx_end=where((I>0.))[0][[0,-1]]
|
||||
eh.save_w_header('../results/mooc_analyzed/model%i_skip_reflectivity.dat'%i,
|
||||
[q[idx_start:idx_end+1],
|
||||
I[idx_start:idx_end+1], R[idx_start:idx_end+1],
|
||||
random.poisson(I[idx_start:idx_end+1])/I[idx_start:idx_end+1]*R[idx_start:idx_end+1],
|
||||
random.poisson(I[idx_start:idx_end+1]*3./14.)/I[idx_start:idx_end+1]*14./3.*R[idx_start:idx_end+1],
|
||||
random.poisson(I[idx_start:idx_end+1]*10.)/I[idx_start:idx_end+1]/10.*R[idx_start:idx_end+1],
|
||||
|
||||
if 2 in only_items:
|
||||
print "Brilliance Transfer 1x1"
|
||||
bt_1=mr.McSim('../results/brilliance_1x1/mccode.h5')
|
||||
d=bt_1['tof_sample.L']
|
||||
r=bt_1['tof_source.L']
|
||||
x, y, e=eh.calc_brilliance_transfer(d, BT_SAMPLE_SCALE*100, r, BT_SOURCE_SCALE)
|
||||
savetxt('../results/analyzed/brilliance_transfer_1x1.dat', array([x, y, e]).T)
|
||||
savetxt('../results/analyzed/brilliance_1x1.dat', array([x, d.data*BT_SAMPLE_SCALE*100,
|
||||
d.errors*BT_SAMPLE_SCALE*100]).T)
|
||||
|
||||
if 3 in only_items:
|
||||
print "Refelctiviy 10x10 sample"
|
||||
opts=dict(qres=0.01, qmin=0.003, qmax=0.5, mindq=5e-4)
|
||||
names=['q_z', 'Intensity(ToF)', 'Reflectivity(ToF)', 'dR (1s ToF)', 'Intensity(λ)', 'Reflectivity(λ)']
|
||||
units=['A^{-1}', 'cts/s', '1', 'cts/s', '1']
|
||||
info='Reflectivity simulation for 10x10mm² sample, omega=%.1f deg'
|
||||
|
||||
print "0.8 degree"
|
||||
ref=mr.McSim('../results/reference_10x10_10/mccode.h5')['tof_detector']
|
||||
ni_08=mr.McSim('../results/nickle_10x10_08/mccode.h5')['tof_detector']
|
||||
q, I, R=eh.calcR(ni_08, ref, 0.8, 1.0, detcorr=True, **opts)
|
||||
q, I_real, R_real=eh.calcR(ni_08, ref, 0.8, 1.0, use_tof=False, **opts)
|
||||
idx_start, idx_end=where((I>0.)&(I_real>0.))[0][[0,-1]]
|
||||
eh.save_w_header('../results/analyzed/reflectivity_10x10_08.dat',
|
||||
[q[idx_start:idx_end+1],
|
||||
I[idx_start:idx_end+1], R[idx_start:idx_end+1],
|
||||
(R/sqrt(I))[idx_start:idx_end+1],
|
||||
I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
|
||||
info, names, units)
|
||||
info, names, units)
|
||||
del(ni_08)
|
||||
|
||||
print "3.0 degree"
|
||||
ni_30=mr.McSim('../results/nickle_10x10_30/mccode.h5')['tof_detector']
|
||||
q, I, R=eh.calcR(ni_30, ref, 3.0, 1.0, detcorr=True, **opts)
|
||||
q, I_real, R_real=eh.calcR(ni_30, ref, 3.0, 1.0, use_tof=False, **opts)
|
||||
idx_start, idx_end=where((I>0.)&(I_real>0.))[0][[0,-1]]
|
||||
eh.save_w_header('../results/analyzed/reflectivity_10x10_30.dat',
|
||||
[q[idx_start:idx_end+1],
|
||||
I[idx_start:idx_end+1], R[idx_start:idx_end+1],
|
||||
(R/sqrt(I))[idx_start:idx_end+1],
|
||||
I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
|
||||
info, names, units)
|
||||
del(ni_30)
|
||||
|
||||
print "8.0 degree"
|
||||
ni_80=mr.McSim('../results/nickle_10x10_80/mccode.h5')['tof_detector']
|
||||
q, I, R=eh.calcR(ni_80, ref, 8.0, 1.0, detcorr=True, **opts)
|
||||
q, I_real, R_real=eh.calcR(ni_80, ref, 8.0, 1.0, use_tof=False, **opts)
|
||||
idx_start, idx_end=where((I>0.)&(I_real>0.))[0][[0,-1]]
|
||||
eh.save_w_header('../results/analyzed/reflectivity_10x10_80.dat',
|
||||
[q[idx_start:idx_end+1],
|
||||
I[idx_start:idx_end+1], R[idx_start:idx_end+1],
|
||||
(R/sqrt(I))[idx_start:idx_end+1],
|
||||
I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
|
||||
info, names, units)
|
||||
|
||||
if 4 in only_items:
|
||||
print "10x10 sample pulse skipping mode."
|
||||
opts=dict(qres=0.04, qmin=0.003, qmax=0.2, mindq=2e-4)
|
||||
names=['q_z', 'Intensity(ToF)', 'Reflectivity(ToF)', 'Intensity(λ)', 'Reflectivity(λ)']
|
||||
units=['A^{-1}', 'cts/s', '1', 'cts/s', '1']
|
||||
info='Reflectivity simulation for 10x10mm² sample, 2-pulse skipping'
|
||||
|
||||
ref=mr.McSim('../results/single_skip_reference/mccode.h5')['tof_detector']
|
||||
ni=mr.McSim('../results/single_skip_nickle/mccode.h5')['tof_detector']
|
||||
q, I, R=eh.calcR(ni, ref, 2.0, skip_pulses=2, **opts)
|
||||
q, I_real, R_real=eh.calcR(ni, ref, 2.0, skip_pulses=2, use_tof=False, **opts)
|
||||
idx_start, idx_end=where((I>0.))[0][[0,-1]]
|
||||
eh.save_w_header('../results/analyzed/single_skip_reflectivity.dat',
|
||||
[q[idx_start:idx_end+1],
|
||||
I[idx_start:idx_end+1], R[idx_start:idx_end+1],
|
||||
I_real[idx_start:idx_end+1], R_real[idx_start:idx_end+1]],
|
||||
info, names, units)
|
||||
|
||||
if 5 in only_items:
|
||||
print "Extracting event statistics."
|
||||
fh=open('../results/analyzed/event_stats.dat', 'w')
|
||||
fh.write('# event statistics for Estia\n')
|
||||
fh.write('# item avg. total peak total avg. ROI peak ROI\n')
|
||||
fh.write('# [cps] [cps] [cps/mm²] [cps/mm²]\n')
|
||||
line='%(name)s %(total) 12e %(peak) 12e %(roi) 12e %(roi_peak) 12e\n'
|
||||
|
||||
ds=mr.McSim('../results/reference_10x10_10/mccode.h5')['tof_detector']
|
||||
res=eh.calcStat(ds)
|
||||
res['name']='Reference'
|
||||
fh.write(line%res)
|
||||
|
||||
ds=mr.McSim('../results/nickle_10x10_08/mccode.h5')['tof_detector']
|
||||
res=eh.calcStat(ds)
|
||||
res['name']='Ni-0.8deg'
|
||||
fh.write(line%res)
|
||||
|
||||
ds=mr.McSim('../results/nickle_10x10_30/mccode.h5')['tof_detector']
|
||||
res=eh.calcStat(ds)
|
||||
res['name']='Ni-3.0deg'
|
||||
fh.write(line%res)
|
||||
|
||||
ds=mr.McSim('../results/nickle_10x10_80/mccode.h5')['tof_detector']
|
||||
res=eh.calcStat(ds)
|
||||
res['name']='Ni-8.0deg'
|
||||
fh.write(line%res)
|
||||
|
||||
fh.close()
|
||||
|
||||
|
||||
@@ -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'])
|
||||
|
||||
@@ -0,0 +1,753 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Dose_calculator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Rodion Kolevatov
|
||||
* Date: August 2018
|
||||
* Version: $Revision: 1.0 $
|
||||
* Release: ---
|
||||
* Origin: IFE
|
||||
*
|
||||
* Calculating dose rate at the outer surface of lateral shielding of uniform thickness
|
||||
* made of material specified in the component input (concrete, iron or lead).
|
||||
* Shielding starts at a distance 0.5*Innerspace from the center of the guide.
|
||||
*
|
||||
* %D
|
||||
* A Shielding_logger together with a set of Shielding_iterators is required.
|
||||
* Iterators should contain Monitor_nD components for recording capture rates.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
* Shielding type: Concrete, Iron or Lead
|
||||
* Names of text files where the coating capture rates are recorded
|
||||
* Inner space in the shielding given by parameter Innerspace (in meters).
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Dose_calculator
|
||||
DEFINITION PARAMETERS ()
|
||||
SETTING PARAMETERS (string Material = "Fe", double Innerspace=0.3, double Thickness=0.2, string SteelTubing="", double TubingThickness=0.0, string NiCaptureFile="NiCapture.dat",
|
||||
string TiCaptureFile = "TiCapture.dat", string TotalCaptureFile = "TotalCapture.dat", string OutputFile="Dose.dat")
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#ifndef LIN_INT_ROUTINE
|
||||
#define LIN_INT_ROUTINE 1
|
||||
//Multi-d linear interpolation routine. Array of args, number of args, sizes of args in interpolation grid, argument grid in single line, data in single line
|
||||
double lint (double * args, int Narg, int * sizearg, double * argtable, double * datatable)
|
||||
{
|
||||
if (Narg==1) {
|
||||
// printf ("Narg=1, arg = %g\n", *args);
|
||||
int point=0;
|
||||
if ((*args)<(*argtable) ) return *datatable;
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return * (datatable+ *sizearg-1); // if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) point++;
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup = (*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*(*(datatable+point))+weightlow*(*(datatable+point-1));
|
||||
}//arg within range of values
|
||||
}// if Narg==1
|
||||
else if (Narg >1){//if more than one argument
|
||||
//lookup how large is the data (Narg-1)D slice for fixed value of the first argument
|
||||
int slicesize=1;
|
||||
int i,point=0;
|
||||
for ( i=1;i<Narg;i++) slicesize*=sizearg[i];
|
||||
// printf("Narg > 1, slicesize = %d, arggridstart = %g, argument = %g, arggridend = %g\n",slicesize,*argtable,*args, *(argtable+ *sizearg-1) );
|
||||
//search weights for the first argument and get results with one argument less
|
||||
if ((*args)<(*argtable) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable);
|
||||
// if the value is too low -- return what is on the lower bound
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), (datatable+ slicesize*(*sizearg)));
|
||||
// if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) {
|
||||
// printf ("*(argtable+point) = %g\n", *(argtable+point));
|
||||
point++;}
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup =(*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+point*slicesize )
|
||||
+weightlow*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+(point-1)*slicesize);
|
||||
} // argument within the ragne
|
||||
} //Narg>1
|
||||
};
|
||||
#endif
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
#include <dirent.h>
|
||||
#include <unistd.h>
|
||||
|
||||
//*** NEUTRON CONVERSION DATA *****//
|
||||
#ifndef NEUTRON_GAMMA_DATA
|
||||
#define NEUTRON_GAMMA_DATA 1
|
||||
//Photons per capture in NiTi coating
|
||||
double fraction[]={0.003565789, 0.009694611, 0.028838269, 0.149117282, 0.117348519, 0.042269932, 0.386434136, 0.092896131, 0.046767054, 0.032512648, 0.043547581, 0.040537585, 0.510133557, 0.320364465};
|
||||
double energy [] = {0.150,
|
||||
0.200,
|
||||
0.300,
|
||||
0.400,
|
||||
0.500,
|
||||
0.600,
|
||||
0.800,
|
||||
1.,
|
||||
1.50,
|
||||
2.000,
|
||||
3.0,
|
||||
4.0,
|
||||
5.000,
|
||||
6.000,
|
||||
7.0,
|
||||
8.000,
|
||||
9.0,
|
||||
10.000,
|
||||
11.000};
|
||||
int nEgroup= 19;
|
||||
//Spectrum of capture photons from Ni and Ti
|
||||
//Ni gamma capture specrum
|
||||
double fractionNi[]={0.003722,//150
|
||||
0.015513,//200
|
||||
0.048432,//300
|
||||
0.047485,//400
|
||||
0.198059,//500
|
||||
0.002809,//600
|
||||
0.003142,//800
|
||||
0.06579,//1000
|
||||
0.054534,//1500
|
||||
0.028872,//2000
|
||||
0.063882,//3000
|
||||
0.050286,//4000
|
||||
0.035567,//5000
|
||||
0.080326,//6000
|
||||
0.156571,//7000
|
||||
0.130355,//8000
|
||||
0.534153,//9000
|
||||
0.000558,//10000
|
||||
0.000260};//11000
|
||||
|
||||
|
||||
//Ti capture gamma spectrum
|
||||
double fractionTi[]={0.009326, //0.15
|
||||
0.0, //0.2
|
||||
0.001938,//0.3
|
||||
0.302632,//0.4
|
||||
0.000882,//0.5
|
||||
0.000668, //0.6
|
||||
0.000914, //0.8
|
||||
0.023046,//1.0
|
||||
0.947104,//1.5
|
||||
0.232257,//2.0
|
||||
0.086418,//3.0
|
||||
0.162878,//4.0
|
||||
0.11209,//5.0
|
||||
0.013419,//6.0
|
||||
0.0875051, //7.0
|
||||
0.009083,//8.0
|
||||
0.001923, //9.0
|
||||
0.002189,//10.0
|
||||
0.000332};//11.0
|
||||
//double fractionB[]={0.0, 0.0, 0.0, 0.0, 0.93, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0,0.0};
|
||||
double fractionB[]={0,//0.15
|
||||
0,//0.2
|
||||
0,//0.3
|
||||
2.32983E-06,//0.4
|
||||
0.935769511,//0.5
|
||||
1.94926E-06,//0.6
|
||||
2.22356E-05,//0.8
|
||||
3.29279E-05,//1.0
|
||||
0.000248357,//1.5
|
||||
0.000398671,//2.0
|
||||
0.000410038,//3.0
|
||||
0.001042575,//4.0
|
||||
0.000961964,//5.0
|
||||
8.714E-05,//6.0
|
||||
0.000248803,//7.0
|
||||
0.00017111,//8.0
|
||||
3.9074E-05,//9.0
|
||||
7.28E-07,//10.0
|
||||
6.27e-6};//11.0
|
||||
|
||||
/*
|
||||
//conversion test input
|
||||
double fraction [] = {0.5,0.5};
|
||||
double energy [] = {2.0,2.0};
|
||||
int nEgroup= 2;
|
||||
*/
|
||||
//*** SHIELDING DATA TABLES ****//
|
||||
// Mashkovich set + NIST data for 15 MeV
|
||||
// Linear attenuation for concrete
|
||||
double AttenuationArgsConc[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataConc[]={31.7, 28.5, 24.6, 21.9, 20.0, 18.5, 16.3, 14.6, 13.1, 11.9, 10.3, 8.74, 8.37, 7.34, 6.65, 6.19, 5.61, 5.29, 4.8208}; // in m^-1
|
||||
int AttenuationSizeConc[]={18};
|
||||
// Iron
|
||||
double AttenuationArgsFe[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataFe[]={139., 106., 83.3, 71.7, 64.6, 59.5, 52.0, 46.7, 42.2, 38.1, 33.3, 29.1, 28.4, 26.0, 24.8, 24.0, 23.4, 23.4, 24.3}; //in m^-1
|
||||
int AttenuationSizeFe[]={19};
|
||||
double AttenuationArgsPb[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataPb[]={2180., 1070., 425., 244., 170., 133., 95.2, 77.1, 65.8, 57.7, 50.8, 47.6, 46.8, 47.2, 48.1, 49.4, 52.0, 55.0, 64.2}; //in m^-1
|
||||
int AttenuationSizePb[]={19};
|
||||
/*
|
||||
//NIST set
|
||||
//Concrete
|
||||
double AttenuationArgsConc[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataConc[]={
|
||||
33.028,
|
||||
29.486,
|
||||
25.231,
|
||||
22.5009,
|
||||
20.5045,
|
||||
18.9428,
|
||||
16.6221,
|
||||
14.9385,
|
||||
13.3561,
|
||||
12.1624,
|
||||
10.4811,
|
||||
8.5123,
|
||||
7.3991,
|
||||
6.6884,
|
||||
6.2031,
|
||||
5.5936,
|
||||
5.2394,
|
||||
4.8208};
|
||||
int AttenuationSizeConc[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsFe[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataFe[]={
|
||||
1.55E+02,
|
||||
1.15E+02,
|
||||
8.65E+01,
|
||||
7.40E+01,
|
||||
6.63E+01,
|
||||
6.07E+01,
|
||||
5.27E+01,
|
||||
4.72E+01,
|
||||
4.21E+01,
|
||||
3.84E+01,
|
||||
3.36E+01,
|
||||
2.85E+01,
|
||||
2.61E+01,
|
||||
2.48E+01,
|
||||
2.41E+01,
|
||||
2.36E+01,
|
||||
2.36E+01,
|
||||
2.43E+01};
|
||||
int AttenuationSizeFe[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsPb[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in meV
|
||||
double AttenuationDataPb[]={
|
||||
2.28E+03,
|
||||
1.13E+03,
|
||||
4.57E+02,
|
||||
2.63E+02,
|
||||
1.83E+02,
|
||||
1.42E+02,
|
||||
1.01E+02,
|
||||
8.05E+01,
|
||||
6.66E+01,
|
||||
5.92E+01,
|
||||
5.22E+01,
|
||||
4.80E+01,
|
||||
4.76E+01,
|
||||
4.84E+01,
|
||||
4.98E+01,
|
||||
5.30E+01,
|
||||
5.64E+01,
|
||||
6.42E+01
|
||||
}; //in m^-1
|
||||
int AttenuationSizePb[]={18};
|
||||
*/
|
||||
|
||||
/*
|
||||
//Linear attenuation test input
|
||||
double AttenuationArgs[]={0.15, 10.0}; // in meV
|
||||
double AttenuationData[]={7.0, 7.0 }; // in m^-1
|
||||
int AttenuationSize[]={2};
|
||||
*/
|
||||
|
||||
// Dose buildup factors concrete
|
||||
double BuildupDataConc[]={1., 1.74, 2.26, 2.95, 3.79, 4.51, 5.57, 6.51, 3.18,
|
||||
1., 2.82, 5.13, 11.2, 24.2, 42.7, 87.6, 153., 353.,
|
||||
1., 2.52, 4.66, 10.8, 25.6, 48.2, 107., 198., 497.,
|
||||
1., 2.27, 4.03, 8.97, 20.2, 30.4, 75.6, 131, 292,
|
||||
1., 1.98, 3.24, 6.42, 12.7, 20.7, 37.2, 57.1, 106,
|
||||
1., 1.77, 2.65, 4.61, 7.97, 11.7, 18.6, 26.0, 42.2,
|
||||
1., 1.67, 2.38, 3.84, 6.20, 8.71, 13.1, 17.7, 27.4,
|
||||
1., 1.61, 2.18, 3.37, 5.23, 7.15, 10.5, 13.9, 20.9,
|
||||
1., 1.49, 1.93, 2.80, 4.14, 5.52, 7.86, 10.2, 15.5,
|
||||
1., 1.41, 1.76, 2.45, 3.51, 4.59, 6.43, 8.31, 12.2,
|
||||
1., 1.35, 1.64, 2.22, 3.10, 4.01, 5.57, 7.19, 10.6,
|
||||
1., 1.26, 1.46, 1.86, 2.50, 3.16, 4.34, 5.59, 8.27};
|
||||
|
||||
int BuildupSizeConc[] ={12, //energy
|
||||
9}; // mu d
|
||||
double BuildupArgsConc[]={
|
||||
0.05, 0.15, 0.3, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0 // energy in MeV
|
||||
,
|
||||
0., 1., 2., 4., 7., 10., 15., 20., 30. //mu d
|
||||
};
|
||||
//Dose buildup factors Fe
|
||||
double BuildupDataFe[]={
|
||||
1., 1.5, 2.2, 3.1, 4.1, 4.6, 5.4, 5.9, //0.1
|
||||
1., 2.0, 3.1, 5.3, 8.9, 14., 22., 31., //0.2
|
||||
1., 2.1, 3.3, 6.0, 12., 23., 49., 84., //0.4
|
||||
1., 1.98, 3.09, 5.98, 11.7, 19.2, 35.4, 55.6, //0.5
|
||||
1., 1.87, 2.89, 5.39, 10.2, 16.2, 28.3, 42.7, //1.0
|
||||
1., 1.76, 2.43, 4.13, 7.25, 10.9, 17.6, 25.1, //2.0
|
||||
1., 1.55, 2.15, 3.51, 5.85, 8.51, 13.5, 19.1, //3.0
|
||||
1., 1.45, 1.94, 3.03, 4.91, 7.11, 11.2, 16.0, //4.0
|
||||
1., 1.34, 1.72, 2.58, 4.14, 6.02, 9.89, 14.7, //6.0
|
||||
1., 1.27, 1.56, 2.23, 3.49, 5.07, 8.50, 13.0, //8.0
|
||||
1., 1.20, 1.42, 1.95, 2.99, 4.35, 7.54, 12.4, //9.0
|
||||
1., 1.48, 1.86, 2.72, 4.30 , 6.37, 11.4, 19.1};//15.0
|
||||
int BuildupSizeFe[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsFe[]={0.1, 0.2, 0.4, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
//Dose buildup factors Pb
|
||||
double BuildupDataPb[]={1., 1.01, 1.03, 1.06, 1.15, 1.16, 1.18, 1.19,
|
||||
1., 1.11, 1.17, 1.25, 1.34, 1.41, 1.5, 1.56,
|
||||
1., 1.17, 1.29, 1.46, 1.58, 1.72, 1.89, 2.02,
|
||||
1., 1.24, 1.42, 1.69, 2.00, 2.27, 2.65, 2.73,
|
||||
1., 1.37, 1.69, 2.26, 3.02, 3.74, 4.81, 5.86,
|
||||
1., 1.39, 1.76, 2.51, 3.66, 4.84, 6.87, 9.00,
|
||||
1., 1.34, 1.68, 2.43, 3.75, 5.30, 8.44, 12.3,
|
||||
1., 1.27, 1.56, 2.25, 3.61, 5.44, 9.80, 16.3,
|
||||
1., 1.21, 1.46, 2.08, 3.44, 5.55, 11.7, 23.6,
|
||||
1., 1.18, 1.40, 1.97, 3.34, 5.69, 13.8, 32.7,
|
||||
1., 1.14, 1.30, 1.74, 2.89, 5.07, 14.1, 44.6,
|
||||
1., 1.11, 1.23, 1.58, 2.52, 4.34, 12.5, 39.2};
|
||||
|
||||
int BuildupSizePb[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsPb[]={0.15, 0.30, 0.40, 0.5, 1.0, 2.0, 3.0, 4.0, 5.1, 6.0, 8.0, 10.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
|
||||
|
||||
// Flux-to-dose conversion
|
||||
/* //NRB-99?
|
||||
double FtoDArgs[]= {0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.80, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*0.752E-10, 3600*1.00E-10, 3600*1.51E-10, 3600*2.00E-10, 3600*2.47E-10, 3600*2.91E-10, 3600*3.73E-10, 3600*4.48E-10, 3600*7.49E-10, 3600*12.0E-10, 3600*16.0E-10, 3600*19.9E-10, 3600*23.8E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 13;
|
||||
int FtoDSize[] = {13};
|
||||
*/
|
||||
//ESS
|
||||
double FtoDArgs[]= {
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.511,
|
||||
0.6,
|
||||
0.662,
|
||||
0.8,
|
||||
1.,
|
||||
1.12,
|
||||
1.33,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
6.13,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in MeV
|
||||
double FtoDData[] = {
|
||||
2.69E-07,
|
||||
3.60E-07,
|
||||
5.44E-07,
|
||||
7.20E-07,
|
||||
8.89E-07,
|
||||
9.07E-07,
|
||||
1.05E-06,
|
||||
1.14E-06,
|
||||
1.34E-06,
|
||||
1.62E-06,
|
||||
1.76E-06,
|
||||
2.01E-06,
|
||||
2.20E-06,
|
||||
2.69E-06,
|
||||
3.51E-06,
|
||||
4.21E-06,
|
||||
4.82E-06,
|
||||
5.40E-06,
|
||||
5.47E-06,
|
||||
6.70E-06,
|
||||
7.92E-06,
|
||||
1.09E-05
|
||||
}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 22;
|
||||
int FtoDSize[] = {22};
|
||||
|
||||
|
||||
|
||||
// test input flux to dose
|
||||
/*
|
||||
double FtoDArgs[]= {0.15, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*7.0E-10,3600*7.0E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 2;
|
||||
int FtoDSize[] = {2};
|
||||
*/
|
||||
#endif //end of insertion of datatables.
|
||||
/*** end of insertion to DECLARE section **/
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
|
||||
/** Define material for dose calculation **/
|
||||
|
||||
int imat=0, itubing=0;
|
||||
if(strcmp(SteelTubing,"yes")==0||strcmp(SteelTubing,"Yes")==0||strcmp(SteelTubing,"YES")==0||strcmp(SteelTubing,"y")==0||strcmp(SteelTubing,"Y")==0||strcmp(SteelTubing,"1")==0){
|
||||
itubing=1;
|
||||
imat=0;
|
||||
printf("Calculating dose for STEEL TUBING\n Outer material will be CONCRETE\n");}
|
||||
else if((strcmp(Material,"Fe")==0)||(strcmp(Material,"fe")==0)||(strcmp(Material,"Iron")==0)||(strcmp(Material,"iron")==0)||(strcmp(Material,"IRON")==0)){
|
||||
itubing=0;
|
||||
imat=1;
|
||||
printf("Calculating dose for IRON\n");}
|
||||
else if ((strcmp(Material,"Pb")==0)||(strcmp(Material,"pb")==0)||(strcmp(Material,"Lead")==0)||(strcmp(Material,"lead")==0)||(strcmp(Material,"LEAD")==0))
|
||||
{itubing=0;
|
||||
imat=2;
|
||||
printf("Calculating dose for LEAD\n");}
|
||||
else if((strcmp(Material,"Conc")==0)||(strcmp(Material,"Concrete")==0)||(strcmp(Material,"conc")==0)||(strcmp(Material,"concrete")==0))
|
||||
{itubing=0;
|
||||
imat=0;
|
||||
printf("Calculating dose for CONCRETE\n");}
|
||||
else
|
||||
{
|
||||
printf ("No/wrong material specified for dose calculation, possible options: Fe, Pb, concrete\n");
|
||||
printf ("Using default: concrete\n");
|
||||
}
|
||||
/**Reading datafiles with capture per bin**/
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
double d1, d2, d3, d4;
|
||||
int ibin=0,NBINS;
|
||||
char line[1000],line1[1000];
|
||||
char dirname[1000];
|
||||
memset(dirname,0,sizeof(dirname));
|
||||
strcat(strcat(dirname,mcdirname),separator);
|
||||
|
||||
char filename[1000];
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FILE *dataIn;
|
||||
|
||||
memset(filename,0,sizeof(filename));//resetting filename to NULL string
|
||||
strcat(strcat(filename,dirname),NiCaptureFile);
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",NiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4) ibin++;
|
||||
};
|
||||
|
||||
|
||||
NBINS=ibin;
|
||||
|
||||
//printf("Shielding calculator:\n Input file contains %d bins.\n", NBINS);
|
||||
|
||||
double* zhat = malloc(NBINS*sizeof(double));
|
||||
double* Ihat= malloc(NBINS*sizeof(double));
|
||||
double* IhatNi= malloc(NBINS*sizeof(double));
|
||||
double* IhatTi= malloc(NBINS*sizeof(double));
|
||||
rewind(dataIn);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatNi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
memset(filename,0,sizeof(filename));// setting filename to NULL string
|
||||
strcat(strcat(filename,dirname),TiCaptureFile); //setting filename to TiCapture file
|
||||
//printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatTi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
fclose(dataIn);
|
||||
|
||||
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),TotalCaptureFile);
|
||||
//printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TotalCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
Ihat[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
int i,iz,iEgroup; // Auxiliary variables
|
||||
|
||||
double zbinlength = zhat[2]-zhat[1];
|
||||
|
||||
|
||||
|
||||
double RRR = (0.5*Innerspace)+Thickness+itubing*TubingThickness;
|
||||
double* doseUnshielded= malloc(NBINS*sizeof(double));
|
||||
double* doseShielded= malloc(NBINS*sizeof(double));
|
||||
|
||||
|
||||
|
||||
|
||||
/*Moving along the bins, calculating the dose */
|
||||
|
||||
|
||||
|
||||
for ( iz=0; iz< NBINS; iz++)
|
||||
{double zc=zhat[iz]; // the z-value for the calculation
|
||||
double z;
|
||||
doseShielded[iz]=0.0;
|
||||
doseUnshielded[iz]=0.0;
|
||||
|
||||
/*computing the integral*/
|
||||
double dmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
doseUnshielded[iz]+=(fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926/(RRR*RRR+z*z)*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData);
|
||||
if (imat==2){//LEAD
|
||||
dmfp=lint(&En,1,AttenuationSizePb, AttenuationArgsPb, AttenuationDataPb)*Thickness*sqrt(RRR*RRR+z*z)/RRR;
|
||||
double args[]={En,dmfp};
|
||||
doseShielded[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926
|
||||
*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)
|
||||
*exp(-dmfp)/(RRR*RRR+z*z)*lint(args, 2, BuildupSizePb,BuildupArgsPb,BuildupDataPb);}
|
||||
else if (imat==1){//IRON
|
||||
dmfp=lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*Thickness*sqrt(RRR*RRR+z*z)/RRR;
|
||||
double args[]={En,dmfp};
|
||||
doseShielded[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926
|
||||
*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)
|
||||
*exp(-dmfp)/(RRR*RRR+z*z)*lint(args, 2, BuildupSizeFe,BuildupArgsFe,BuildupDataFe);}
|
||||
else {//CONCRETE, imat=0, default
|
||||
dmfp=lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*Thickness*sqrt(RRR*RRR+z*z)/RRR + itubing*lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*TubingThickness*sqrt(RRR*RRR+z*z)/RRR;
|
||||
double args[]={En,dmfp};
|
||||
doseShielded[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926
|
||||
*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)
|
||||
*exp(-dmfp)/(RRR*RRR+z*z)*lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);}
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
} /* Position scan iz */
|
||||
|
||||
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
FILE *dose_out;
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),OutputFile);
|
||||
printf("Dose calculator:\nWriting to file %s\n",filename);
|
||||
dose_out=fopen(filename,"w");
|
||||
|
||||
|
||||
fprintf(dose_out,"#Dose at the outer surface of lateral %s shielding of thickness %g m\n#Guide is in the center of the innner space of width %g m\n#Presense of Steel Tubing = %s; of %g m\n#Position z,m\tOverall loss, n/s/m\tNi capture, n/s/m,\tTi capture, n/s/m\tDose unshielded, uSv/h\tDose shielded, uSv/h\n",Material, Thickness, Innerspace, SteelTubing,TubingThickness);
|
||||
|
||||
|
||||
for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
fprintf( dose_out,"%g\t%.3g\t%.3g\t%.3g\t%.3g\t%.3g\n",zhat[i], Ihat[i]/zbinlength, IhatNi[i]/zbinlength, IhatTi[i]/zbinlength,doseUnshielded[i],doseShielded[i]);
|
||||
}
|
||||
fclose(dose_out);
|
||||
free(zhat); free(Ihat); free(IhatNi); free(IhatTi);
|
||||
free(doseShielded); free(doseUnshielded);
|
||||
|
||||
printf("Done with writing dose rates\n");
|
||||
%}
|
||||
|
||||
END
|
||||
File diff suppressed because it is too large
Load Diff
+262
-39
@@ -49,13 +49,13 @@
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE INSTRUMENT ESS_reflectometer_Estia
|
||||
(double omegaa = 1.2, int sample = 1, string sample_file="alessandra_model_2_nores_000.dat",
|
||||
double sample_length = 0.01, double sample_height = 0.01,
|
||||
(double omegaa = 25, int sample = 0, double sample_length = 0.02, double sample_height = 0.015,
|
||||
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,
|
||||
double lambda_min = 3.75, double lambda_start = 0.1, 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
|
||||
double source_power = 5,
|
||||
double selene1_foot1y =0.0, double selene1_foot2y = 0.0
|
||||
)
|
||||
|
||||
DECLARE
|
||||
@@ -71,8 +71,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
|
||||
|
||||
@@ -90,7 +94,7 @@ double chopper_phase ; // deg phase between pulse and chopper ope
|
||||
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 chopper_freq = 14.0 ; // Hz chopper frequency
|
||||
|
||||
double slit_distance = 1.775; //m, distance slit to sample
|
||||
|
||||
@@ -98,9 +102,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
|
||||
@@ -113,12 +114,10 @@ double Theta1_analyzer1, Theta1_analyzer2, Theta2_analyzer1, Theta2_analyzer2, d
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
|
||||
chopper_freq/=fmax(1.0, enable_chopper);
|
||||
// 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 */
|
||||
@@ -133,7 +132,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
|
||||
@@ -160,10 +163,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
|
||||
@@ -201,6 +206,7 @@ COMPONENT moderator = ESS_butterfly(
|
||||
n_pulses = 1+enable_chopper, acc_power=source_power)
|
||||
AT (0, 0, 0) RELATIVE origin
|
||||
|
||||
|
||||
/***************************************
|
||||
* Geometry of neutron feeder separate *
|
||||
***************************************/
|
||||
@@ -209,11 +215,14 @@ COMPONENT moderator = ESS_butterfly(
|
||||
/****************************************************
|
||||
* 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
|
||||
|
||||
|
||||
@@ -221,33 +230,200 @@ COMPONENT FeNi_out = Slit(xwidth=0.013, yheight=0.038)
|
||||
COMPONENT chopper = DiskChopper(radius=chopper_diameter/2.0, yheight=0.02,
|
||||
theta_0=chopper_open, phase=chopper_phase+chopper_open/2.0,
|
||||
nu=chopper_freq, nslit=1)
|
||||
WHEN enable_chopper!=0
|
||||
WHEN enable_chopper==1
|
||||
AT (0, 0, chopper_pos-2*NBOA_c) RELATIVE arm_virtual_source_beam
|
||||
|
||||
|
||||
|
||||
/* The actual virtual source mask, two L-shaped absorbers (first top-right) */
|
||||
COMPONENT virtual_source_TR = Slit(
|
||||
xmin = 0.0, xmax = 1.0, ymin = -1.0, ymax = sample_height/2+over_illumination*5)
|
||||
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
|
||||
|
||||
|
||||
/******* Start of scatter logger. Logging scatterings and absorption in the guide system. *********/
|
||||
COMPONENT log_P_start=Shielding_logger()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_logger_stop
|
||||
#undef scatter_logger_stop
|
||||
#endif
|
||||
#define scatter_logger_stop log_P_start
|
||||
%}
|
||||
|
||||
|
||||
|
||||
|
||||
/*************************************
|
||||
* Geometry of Selene guide separate *
|
||||
*************************************/
|
||||
%include "Estia_selene.instr"
|
||||
%include "Estia_selene1.instr"
|
||||
|
||||
%include "Estia_mf.instr"
|
||||
|
||||
%include "Estia_selene2.instr"
|
||||
|
||||
|
||||
|
||||
COMPONENT Lambda_EndGuide = L_monitor(
|
||||
nL = 1000, filename = "Lambda_EndGuide", Lmin=0.5, Lmax=100,
|
||||
restore_neutron = 1)
|
||||
AT (0, 0, -slit_distance-0.01) RELATIVE arm_sample_beam
|
||||
|
||||
|
||||
|
||||
/*Stop of scatter logger. Record scatterings in what is between start and stop, that is, parabolic feeder. */
|
||||
COMPONENT log_P_stop=Shielding_logger_stop(logger=log_P_start)
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
|
||||
/**IIIIIIIIIIIIIIIIIII ITERATOR SECTION. PROCESSING STORED EVENTS IIIIIIIIIIIIIIII***/
|
||||
/* Insert this into McStas instrument file to do the costs evaluation */
|
||||
COMPONENT arm_iter_P1_start=Arm()
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
COMPONENT iter_P1_start = Shielding_log_iterator_Ni_new()
|
||||
AT (0,0,0) RELATIVE moderator //ABSOLUTE
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iter_P1_start
|
||||
%}
|
||||
JUMP arm_iter_P1_stop WHEN(optics_not_hit)
|
||||
|
||||
/*Monitoring the tracks stored by the scatter logger*/
|
||||
/*Putting dummy arm to register all neutrons to ensure that monitors_nD with shape "previous" will process them */
|
||||
COMPONENT arm_iter_P1_dummy=Arm ()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT mndP01=Monitor_nD (
|
||||
restore_neutron=1, zmin=0.0,
|
||||
zmax=35.0,
|
||||
bins=350, options="previous no slit z ", filename="NiCapture.dat")
|
||||
AT(0,0,0) RELATIVE moderator //RELATIVE source
|
||||
|
||||
COMPONENT arm_iter_P1_stop=Arm()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT iter_P1_stop = Shielding_log_iterator_stop(iterator=iter_P1_start)
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
/*Moving again to the reference point of iterator start
|
||||
when there are still some tracks stored to perform iterations with,
|
||||
checked by the function MC_GETPAR */
|
||||
COMPONENT a11i = Arm()
|
||||
AT (0,0,0) RELATIVE Lambda_EndGuide
|
||||
JUMP arm_iter_P1_start WHEN(MC_GETPAR(iter_P1_stop,loop))
|
||||
|
||||
/*IIIIIIIIIIIIIIIIIIII END OF PROCESSING. END OF ITERATOR SECTION IIIIIIIIIIIIIIIIIIIIIIIIIIIII*/
|
||||
|
||||
|
||||
|
||||
|
||||
/**IIIIIIIIIIIIIIIIIII ITERATOR SECTION2. PROCESSING STORED EVENTS IIIIIIIIIIIIIIII***/
|
||||
/* Insert this into McStas instrument file to do the costs evaluation */
|
||||
COMPONENT arm_iter_P2_start=Arm()
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
COMPONENT iter_P2_start = Shielding_log_iterator_Ti_new()
|
||||
AT (0,0,0) RELATIVE moderator //ABSOLUTE
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iter_P2_start
|
||||
%}
|
||||
JUMP arm_iter_P2_stop WHEN(optics_not_hit)
|
||||
|
||||
/*Monitoring the tracks stored by the scatter logger*/
|
||||
/*Putting dummy arm to register all neutrons to ensure that monitors_nD with shape "previous" will process them */
|
||||
COMPONENT arm_iter_P2_dummy=Arm ()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT mndP02=Monitor_nD (
|
||||
restore_neutron=1, zmin=0.0,
|
||||
zmax=35.0,
|
||||
bins=350, options="previous no slit z ", filename="TiCapture.dat")
|
||||
AT(0,0,0) RELATIVE moderator //RELATIVE source
|
||||
|
||||
COMPONENT arm_iter_P2_stop=Arm()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT iter_P2_stop = Shielding_log_iterator_stop(iterator=iter_P2_start)
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
/*Moving again to the reference point of iterator start
|
||||
when there are still some tracks stored to perform iterations with,
|
||||
checked by the function MC_GETPAR */
|
||||
COMPONENT a12i = Arm()
|
||||
AT (0,0,0) RELATIVE Lambda_EndGuide
|
||||
JUMP arm_iter_P2_start WHEN(MC_GETPAR(iter_P2_stop,loop))
|
||||
|
||||
/*IIIIIIIIIIIIIIIIIIII END OF PROCESSING. END OF ITERATOR2 SECTION IIIIIIIIIIIIIIIIIIIIIIIIIIIII*/
|
||||
|
||||
/**IIIIIIIIIIIIIIIIIII ITERATOR SECTION3. PROCESSING STORED EVENTS IIIIIIIIIIIIIIII***/
|
||||
/* Insert this into McStas instrument file to do the costs evaluation */
|
||||
COMPONENT arm_iter_P3_start=Arm()
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
COMPONENT iter_P3_start = Shielding_log_iterator_total()
|
||||
AT (0,0,0) RELATIVE moderator //ABSOLUTE
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iter_P3_start
|
||||
%}
|
||||
JUMP arm_iter_P3_stop WHEN(optics_not_hit)
|
||||
|
||||
/*Monitoring the tracks stored by the scatter logger*/
|
||||
/*Putting dummy arm to register all neutrons to ensure that monitors_nD with shape "previous" will process them */
|
||||
COMPONENT arm_iter_P3_dummy=Arm ()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT mndP03=Monitor_nD (
|
||||
restore_neutron=1, zmin=0.0,
|
||||
zmax=35.0,
|
||||
bins=350, options="previous no slit z ", filename="TotalCapture.dat")
|
||||
AT(0,0,0) RELATIVE moderator //RELATIVE source
|
||||
|
||||
COMPONENT arm_iter_P3_stop=Arm()
|
||||
AT (0,0,0) RELATIVE PREVIOUS
|
||||
|
||||
COMPONENT iter_P3_stop = Shielding_log_iterator_stop(iterator=iter_P3_start, last=1)
|
||||
AT(0,0,0) RELATIVE moderator
|
||||
|
||||
/*Moving again to the reference point of iterator start
|
||||
when there are still some tracks stored to perform iterations with,
|
||||
checked by the function MC_GETPAR */
|
||||
COMPONENT a13i = Arm()
|
||||
AT (0,0,0) RELATIVE Lambda_EndGuide
|
||||
JUMP arm_iter_P3_start WHEN(MC_GETPAR(iter_P3_stop,loop))
|
||||
|
||||
/*IIIIIIIIIIIIIIIIIIII END OF PROCESSING. END OF ITERATOR3 SECTION IIIIIIIIIIIIIIIIIIIIIIIIIIIII*/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**********************************
|
||||
@@ -269,6 +445,24 @@ COMPONENT ac_slit = Slit(
|
||||
/***************
|
||||
* Sample area *
|
||||
***************/
|
||||
COMPONENT tof_sample = Monitor_nD(
|
||||
filename = "tof_sample",
|
||||
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) RELATIVE arm_sample_beam
|
||||
ROTATED (0, 0, 0) RELATIVE arm_sample_beam
|
||||
|
||||
|
||||
/* NiTi multilayer sample */
|
||||
COMPONENT sample = Mirror(
|
||||
xwidth = sample_length, yheight = sample_height,
|
||||
center = 1, transmit = 0,
|
||||
reflect = "NiTiML.ref"
|
||||
)
|
||||
WHEN sample==0
|
||||
AT (0, 0, 0) RELATIVE arm_sample
|
||||
ROTATED (0, 90, 0) RELATIVE arm_sample
|
||||
|
||||
/* ideal reflector as reference */
|
||||
COMPONENT reference_sample = Mirror(
|
||||
@@ -280,28 +474,57 @@ COMPONENT reference_sample = Mirror(
|
||||
AT (0, 0, 0) RELATIVE arm_sample
|
||||
ROTATED (0, 90, 0) RELATIVE arm_sample
|
||||
|
||||
/* User defined sample file */
|
||||
COMPONENT mooc_sample = Mirror(
|
||||
/* Nickel film on silicon */
|
||||
COMPONENT ni_sample = Mirror(
|
||||
xwidth = sample_length, yheight = sample_height,
|
||||
center = 1, transmit = 0,
|
||||
reflect = sample_file
|
||||
reflect = "Si-Ni.ref"
|
||||
)
|
||||
WHEN sample==2
|
||||
AT (0, 0, 0) RELATIVE arm_sample
|
||||
ROTATED (0, 90, 0) RELATIVE arm_sample
|
||||
|
||||
/* Silicon with natural oxide */
|
||||
COMPONENT si_sample = Mirror(
|
||||
xwidth = sample_length, yheight = sample_height,
|
||||
center = 1, transmit = 0,
|
||||
reflect = "Si-SiO2.ref"
|
||||
)
|
||||
WHEN sample==3
|
||||
AT (0, 0, 0) RELATIVE arm_sample
|
||||
ROTATED (0, 90, 0) RELATIVE arm_sample
|
||||
|
||||
|
||||
COMPONENT arm_analyzer = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_detector
|
||||
ROTATED (-selene_theta+(Theta1_analyzer1-Theta2_analyzer1)/2.0, 0, 0) RELATIVE arm_detector
|
||||
|
||||
COMPONENT arm_analyzer2 = Arm()
|
||||
AT (0, 0, 0) RELATIVE arm_detector
|
||||
ROTATED (-selene_theta+(Theta1_analyzer2-Theta2_analyzer2)/2.0, 0, 0) RELATIVE arm_detector
|
||||
|
||||
/* polarization analyser */
|
||||
COMPONENT analyzer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer1_start, length=analyzer1_length,
|
||||
delta_theta=(Theta1_analyzer1+Theta2_analyzer1+0.05)*PI/180.0, h2=0.14, h1=0.05, abs_ref=0,
|
||||
m_u=5.5, m_d=0.45, both_coated=0, alpha=2.3, W = 0.0014)
|
||||
WHEN enable_analyzer
|
||||
AT (0, 0.0, analyzer1_start) RELATIVE arm_analyzer
|
||||
ROTATED (0,0,0.0) RELATIVE arm_analyzer
|
||||
|
||||
COMPONENT analyzer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=analyzer2_start, length=analyzer2_length,
|
||||
delta_theta=(Theta1_analyzer2+Theta2_analyzer2+0.05)*PI/180.0, h2=0.2, h1=0.05, abs_out=0,
|
||||
m_u=5.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
|
||||
WHEN enable_analyzer==2
|
||||
AT (0, 0.0, analyzer2_start) RELATIVE arm_analyzer
|
||||
ROTATED (0,0,0.0) RELATIVE arm_analyzer
|
||||
|
||||
/* 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
|
||||
// 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
|
||||
|
||||
/***********************************************************************/
|
||||
|
||||
|
||||
@@ -201,7 +201,7 @@ COMPONENT E02_01_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
%}
|
||||
|
||||
/* Insert part of the feeder before the monolith exit window*/
|
||||
COMPONENT E02_01_01 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_01_01 = Elliptic_guide_gravity_custom(
|
||||
l=E02_01_01_length-NBOA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_01_01_start, loutyh= NBOA_c*2-E02_01_01_start-E02_01_01_length,
|
||||
linxw = E02_01_01_start, loutxw= NBOA_c*2-E02_01_01_start-E02_01_01_length,
|
||||
@@ -229,7 +229,7 @@ COMPONENT E02_01_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_01_02 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_01_02 = Elliptic_guide_gravity_custom(
|
||||
l=E02_01_02_length-NBOA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_01_02_start, loutyh= NBOA_c*2-E02_01_02_start-E02_01_02_length,
|
||||
linxw = E02_01_02_start, loutxw= NBOA_c*2-E02_01_02_start-E02_01_02_length,
|
||||
@@ -266,7 +266,7 @@ COMPONENT E02_01_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_01_03 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_01_03 = Elliptic_guide_gravity_custom(
|
||||
l=E02_01_03_length, dimensionsAt = "mid",
|
||||
linyh = E02_01_03_start, loutyh= NBOA_c*2-E02_01_03_start-E02_01_03_length,
|
||||
linxw = E02_01_03_start, loutxw= NBOA_c*2-E02_01_03_start-E02_01_03_length,
|
||||
@@ -340,7 +340,7 @@ COMPONENT E02_02_012_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_02_01 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_02_01 = Elliptic_guide_gravity_custom(
|
||||
l=E02_02_01_length-NFGA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_02_01_start, loutyh= NFGA_c*2-E02_02_01_start-E02_02_01_length,
|
||||
linxw = E02_02_01_start, loutxw= NFGA_c*2-E02_02_01_start-E02_02_01_length,
|
||||
@@ -368,7 +368,7 @@ COMPONENT E02_02_022_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_02_02 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_02_02 = Elliptic_guide_gravity_custom(
|
||||
l=E02_02_02_length-NFGA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_02_02_start, loutyh= NFGA_c*2-E02_02_02_start-E02_02_02_length,
|
||||
linxw = E02_02_02_start, loutxw= NFGA_c*2-E02_02_02_start-E02_02_02_length,
|
||||
@@ -405,7 +405,7 @@ COMPONENT E02_02_033_wedge = AbsorberFixed(xmin=-0.1, xmax=0.0,
|
||||
PROP_Z0;
|
||||
%}
|
||||
|
||||
COMPONENT E02_02_03 = Elliptic_guide_gravity(
|
||||
COMPONENT E02_02_03 = Elliptic_guide_gravity_custom(
|
||||
l=E02_02_03_length-NFGA_gap, dimensionsAt = "mid",
|
||||
linyh = E02_02_03_start, loutyh= NFGA_c*2-E02_02_03_start-E02_02_03_length,
|
||||
linxw = E02_02_03_start, loutxw= NFGA_c*2-E02_02_03_start-E02_02_03_length,
|
||||
|
||||
@@ -0,0 +1,96 @@
|
||||
/*******************************************************************************
|
||||
* 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
|
||||
|
||||
|
||||
/**************************************
|
||||
* 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
|
||||
|
||||
|
||||
/* 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
|
||||
@@ -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
|
||||
@@ -1,204 +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 *
|
||||
**************************************/
|
||||
|
||||
|
||||
|
||||
/**************************
|
||||
* 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
|
||||
|
||||
/* 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
|
||||
@@ -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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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
|
||||
@@ -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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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_custom(
|
||||
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
|
||||
+18001
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,949 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Shielding_calculator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Rodion Kolevatov
|
||||
* Date: August 2018
|
||||
* Version: $Revision: 1.0 $
|
||||
* Release: ---
|
||||
* Origin: IFE
|
||||
*
|
||||
* Calculating lateral shielding thickness required to attenuate dose rate from gammas
|
||||
* arising from coating capture to a level specified by MaxRate.
|
||||
*
|
||||
* %D
|
||||
* Start of the trace-region in which SCATTER events should be logged.
|
||||
* Whenever a SCATTER occurs in components between this one and its counterpart Scatter_logger_stop the neutron state is logged.
|
||||
* The log is kept in memory - but only for one single
|
||||
* numerical neutron at a time, so there should be no real danger of memory overrun.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
* Allowed dose rate beyond shielding, MaxRate
|
||||
* Names of text files where the coating capture rates are recorded
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_calculator
|
||||
DEFINITION PARAMETERS ()
|
||||
SETTING PARAMETERS (double MaxRate = 3.0, double Innerspace = 0.3, string NiCaptureFile = "NiCapture.dat",
|
||||
string TiCaptureFile = "TiCapture.dat", string TotalCaptureFile = "TotalCapture.dat", string OutputFile="Shielding.dat")
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#ifndef LIN_INT_ROUTINE
|
||||
#define LIN_INT_ROUTINE 1
|
||||
//Multi-d linear interpolation routine. Array of args, number of args, sizes of args in interpolation grid, argument grid in single line, data in single line
|
||||
double lint (double * args, int Narg, int * sizearg, double * argtable, double * datatable)
|
||||
{
|
||||
if (Narg==1) {
|
||||
// printf ("Narg=1, arg = %g\n", *args);
|
||||
int point=0;
|
||||
if ((*args)<(*argtable) ) return *datatable;
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return * (datatable+ *sizearg-1); // if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) point++;
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup = (*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*(*(datatable+point))+weightlow*(*(datatable+point-1));
|
||||
}//arg within range of values
|
||||
}// if Narg==1
|
||||
else if (Narg >1){//if more than one argument
|
||||
//lookup how large is the data (Narg-1)D slice for fixed value of the first argument
|
||||
int slicesize=1;
|
||||
int i,point=0;
|
||||
for ( i=1;i<Narg;i++) slicesize*=sizearg[i];
|
||||
// printf("Narg > 1, slicesize = %d, arggridstart = %g, argument = %g, arggridend = %g\n",slicesize,*argtable,*args, *(argtable+ *sizearg-1) );
|
||||
//search weights for the first argument and get results with one argument less
|
||||
if ((*args)<(*argtable) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable);
|
||||
// if the value is too low -- return what is on the lower bound
|
||||
else if ((*args)>*(argtable+ *sizearg-1) ) return lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), (datatable+ slicesize*(*sizearg)));
|
||||
// if the value is too large return what corresponds to largest point on the grid available
|
||||
else { // now argument is withing the range of values
|
||||
//interval lookup
|
||||
while (*args>*(argtable+point)) {
|
||||
// printf ("*(argtable+point) = %g\n", *(argtable+point));
|
||||
point++;}
|
||||
//weights
|
||||
double interval = *(argtable + point) - *(argtable+point-1);
|
||||
double weightup =(*args - *(argtable+point-1))/interval;
|
||||
double weightlow =(*(argtable + point) - *args)/interval;
|
||||
return weightup*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+point*slicesize )
|
||||
+weightlow*lint(args+1,Narg-1, sizearg+1, argtable+(*sizearg), datatable+(point-1)*slicesize);
|
||||
} // argument within the ragne
|
||||
} //Narg>1
|
||||
};
|
||||
#endif
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
#include <dirent.h>
|
||||
#include <unistd.h>
|
||||
/******Parameters relevant for calculations of the shielding*******/
|
||||
#define stepPb 0.001 // Step for calculations with steel shielding (2mm)
|
||||
#define stepFe 0.001 // Step for calculations with steel shielding (5mm)
|
||||
#define stepbt 0.005 // Step for calculations with concrete shielding. (1cm)
|
||||
#define Rtubing 0.05 //in meters. Vacuum tubing inner radius, used for calculation of combined shielding in Febt case.
|
||||
#define TFe 0.1 // in meters. Thickness of steel vacuum tubing. (Assuming guide model as guide in evacuated steel pipe with wall thickness TFe).
|
||||
#define stepcomb 0.005 // Step for calculation of combined shielding, TFe layer of steel plus concrete. (2.5 cm)
|
||||
|
||||
//*** NEUTRON CONVERSION DATA *****//
|
||||
#ifndef NEUTRON_GAMMA_DATA
|
||||
#define NEUTRON_GAMMA_DATA 1
|
||||
//Photons per capture in NiTi coating
|
||||
double fraction[]={0.003565789, 0.009694611, 0.028838269, 0.149117282, 0.117348519, 0.042269932, 0.386434136, 0.092896131, 0.046767054, 0.032512648, 0.043547581, 0.040537585, 0.510133557, 0.320364465};
|
||||
double energy [] = {0.150,
|
||||
0.200,
|
||||
0.300,
|
||||
0.400,
|
||||
0.500,
|
||||
0.600,
|
||||
0.800,
|
||||
1.,
|
||||
1.50,
|
||||
2.000,
|
||||
3.0,
|
||||
4.0,
|
||||
5.000,
|
||||
6.000,
|
||||
7.0,
|
||||
8.000,
|
||||
9.0,
|
||||
10.000,
|
||||
11.000};
|
||||
int nEgroup= 19;
|
||||
//Spectrum of capture photons from Ni and Ti
|
||||
//Ni gamma capture specrum
|
||||
double fractionNi[]={0.003722,//150
|
||||
0.015513,//200
|
||||
0.048432,//300
|
||||
0.047485,//400
|
||||
0.198059,//500
|
||||
0.002809,//600
|
||||
0.003142,//800
|
||||
0.06579,//1000
|
||||
0.054534,//1500
|
||||
0.028872,//2000
|
||||
0.063882,//3000
|
||||
0.050286,//4000
|
||||
0.035567,//5000
|
||||
0.080326,//6000
|
||||
0.156571,//7000
|
||||
0.130355,//8000
|
||||
0.534153,//9000
|
||||
0.000558,//10000
|
||||
0.000260};//11000
|
||||
|
||||
|
||||
//Ti capture gamma spectrum
|
||||
double fractionTi[]={0.009326, //0.15
|
||||
0.0, //0.2
|
||||
0.001938,//0.3
|
||||
0.302632,//0.4
|
||||
0.000882,//0.5
|
||||
0.000668, //0.6
|
||||
0.000914, //0.8
|
||||
0.023046,//1.0
|
||||
0.947104,//1.5
|
||||
0.232257,//2.0
|
||||
0.086418,//3.0
|
||||
0.162878,//4.0
|
||||
0.11209,//5.0
|
||||
0.013419,//6.0
|
||||
0.0875051, //7.0
|
||||
0.009083,//8.0
|
||||
0.001923, //9.0
|
||||
0.002189,//10.0
|
||||
0.000332};//11.0
|
||||
//double fractionB[]={0.0, 0.0, 0.0, 0.0, 0.93, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0,0.0};
|
||||
double fractionB[]={0,//0.15
|
||||
0,//0.2
|
||||
0,//0.3
|
||||
2.32983E-06,//0.4
|
||||
0.935769511,//0.5
|
||||
1.94926E-06,//0.6
|
||||
2.22356E-05,//0.8
|
||||
3.29279E-05,//1.0
|
||||
0.000248357,//1.5
|
||||
0.000398671,//2.0
|
||||
0.000410038,//3.0
|
||||
0.001042575,//4.0
|
||||
0.000961964,//5.0
|
||||
8.714E-05,//6.0
|
||||
0.000248803,//7.0
|
||||
0.00017111,//8.0
|
||||
3.9074E-05,//9.0
|
||||
7.28E-07,//10.0
|
||||
6.27e-6};//11.0
|
||||
|
||||
/*
|
||||
//conversion test input
|
||||
double fraction [] = {0.5,0.5};
|
||||
double energy [] = {2.0,2.0};
|
||||
int nEgroup= 2;
|
||||
*/
|
||||
//*** SHIELDING DATA TABLES ****//
|
||||
// Mashkovich set + NIST data for 15 MeV
|
||||
// Linear attenuation for concrete
|
||||
double AttenuationArgsConc[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataConc[]={31.7, 28.5, 24.6, 21.9, 20.0, 18.5, 16.3, 14.6, 13.1, 11.9, 10.3, 8.74, 8.37, 7.34, 6.65, 6.19, 5.61, 5.29, 4.8208}; // in m^-1
|
||||
int AttenuationSizeConc[]={18};
|
||||
// Iron
|
||||
double AttenuationArgsFe[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataFe[]={139., 106., 83.3, 71.7, 64.6, 59.5, 52.0, 46.7, 42.2, 38.1, 33.3, 29.1, 28.4, 26.0, 24.8, 24.0, 23.4, 23.4, 24.3}; //in m^-1
|
||||
int AttenuationSizeFe[]={19};
|
||||
double AttenuationArgsPb[]={0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0, 2.75, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15.0}; // in meV
|
||||
double AttenuationDataPb[]={2180., 1070., 425., 244., 170., 133., 95.2, 77.1, 65.8, 57.7, 50.8, 47.6, 46.8, 47.2, 48.1, 49.4, 52.0, 55.0, 64.2}; //in m^-1
|
||||
int AttenuationSizePb[]={19};
|
||||
/*
|
||||
//NIST set
|
||||
//Concrete
|
||||
double AttenuationArgsConc[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataConc[]={
|
||||
33.028,
|
||||
29.486,
|
||||
25.231,
|
||||
22.5009,
|
||||
20.5045,
|
||||
18.9428,
|
||||
16.6221,
|
||||
14.9385,
|
||||
13.3561,
|
||||
12.1624,
|
||||
10.4811,
|
||||
8.5123,
|
||||
7.3991,
|
||||
6.6884,
|
||||
6.2031,
|
||||
5.5936,
|
||||
5.2394,
|
||||
4.8208};
|
||||
int AttenuationSizeConc[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsFe[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.};
|
||||
|
||||
double AttenuationDataFe[]={
|
||||
1.55E+02,
|
||||
1.15E+02,
|
||||
8.65E+01,
|
||||
7.40E+01,
|
||||
6.63E+01,
|
||||
6.07E+01,
|
||||
5.27E+01,
|
||||
4.72E+01,
|
||||
4.21E+01,
|
||||
3.84E+01,
|
||||
3.36E+01,
|
||||
2.85E+01,
|
||||
2.61E+01,
|
||||
2.48E+01,
|
||||
2.41E+01,
|
||||
2.36E+01,
|
||||
2.36E+01,
|
||||
2.43E+01};
|
||||
int AttenuationSizeFe[]={18};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
double AttenuationArgsPb[]={
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.6,
|
||||
0.8,
|
||||
1.,
|
||||
1.25,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in meV
|
||||
double AttenuationDataPb[]={
|
||||
2.28E+03,
|
||||
1.13E+03,
|
||||
4.57E+02,
|
||||
2.63E+02,
|
||||
1.83E+02,
|
||||
1.42E+02,
|
||||
1.01E+02,
|
||||
8.05E+01,
|
||||
6.66E+01,
|
||||
5.92E+01,
|
||||
5.22E+01,
|
||||
4.80E+01,
|
||||
4.76E+01,
|
||||
4.84E+01,
|
||||
4.98E+01,
|
||||
5.30E+01,
|
||||
5.64E+01,
|
||||
6.42E+01
|
||||
}; //in m^-1
|
||||
int AttenuationSizePb[]={18};
|
||||
*/
|
||||
|
||||
/*
|
||||
//Linear attenuation test input
|
||||
double AttenuationArgs[]={0.15, 10.0}; // in meV
|
||||
double AttenuationData[]={7.0, 7.0 }; // in m^-1
|
||||
int AttenuationSize[]={2};
|
||||
*/
|
||||
|
||||
// Dose buildup factors concrete
|
||||
double BuildupDataConc[]={1., 1.74, 2.26, 2.95, 3.79, 4.51, 5.57, 6.51, 3.18,
|
||||
1., 2.82, 5.13, 11.2, 24.2, 42.7, 87.6, 153., 353.,
|
||||
1., 2.52, 4.66, 10.8, 25.6, 48.2, 107., 198., 497.,
|
||||
1., 2.27, 4.03, 8.97, 20.2, 30.4, 75.6, 131, 292,
|
||||
1., 1.98, 3.24, 6.42, 12.7, 20.7, 37.2, 57.1, 106,
|
||||
1., 1.77, 2.65, 4.61, 7.97, 11.7, 18.6, 26.0, 42.2,
|
||||
1., 1.67, 2.38, 3.84, 6.20, 8.71, 13.1, 17.7, 27.4,
|
||||
1., 1.61, 2.18, 3.37, 5.23, 7.15, 10.5, 13.9, 20.9,
|
||||
1., 1.49, 1.93, 2.80, 4.14, 5.52, 7.86, 10.2, 15.5,
|
||||
1., 1.41, 1.76, 2.45, 3.51, 4.59, 6.43, 8.31, 12.2,
|
||||
1., 1.35, 1.64, 2.22, 3.10, 4.01, 5.57, 7.19, 10.6,
|
||||
1., 1.26, 1.46, 1.86, 2.50, 3.16, 4.34, 5.59, 8.27};
|
||||
|
||||
int BuildupSizeConc[] ={12, //energy
|
||||
9}; // mu d
|
||||
double BuildupArgsConc[]={
|
||||
0.05, 0.15, 0.3, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0 // energy in MeV
|
||||
,
|
||||
0., 1., 2., 4., 7., 10., 15., 20., 30. //mu d
|
||||
};
|
||||
//Dose buildup factors Fe
|
||||
double BuildupDataFe[]={
|
||||
1., 1.5, 2.2, 3.1, 4.1, 4.6, 5.4, 5.9, //0.1
|
||||
1., 2.0, 3.1, 5.3, 8.9, 14., 22., 31., //0.2
|
||||
1., 2.1, 3.3, 6.0, 12., 23., 49., 84., //0.4
|
||||
1., 1.98, 3.09, 5.98, 11.7, 19.2, 35.4, 55.6, //0.5
|
||||
1., 1.87, 2.89, 5.39, 10.2, 16.2, 28.3, 42.7, //1.0
|
||||
1., 1.76, 2.43, 4.13, 7.25, 10.9, 17.6, 25.1, //2.0
|
||||
1., 1.55, 2.15, 3.51, 5.85, 8.51, 13.5, 19.1, //3.0
|
||||
1., 1.45, 1.94, 3.03, 4.91, 7.11, 11.2, 16.0, //4.0
|
||||
1., 1.34, 1.72, 2.58, 4.14, 6.02, 9.89, 14.7, //6.0
|
||||
1., 1.27, 1.56, 2.23, 3.49, 5.07, 8.50, 13.0, //8.0
|
||||
1., 1.20, 1.42, 1.95, 2.99, 4.35, 7.54, 12.4, //9.0
|
||||
1., 1.48, 1.86, 2.72, 4.30 , 6.37, 11.4, 19.1};//15.0
|
||||
int BuildupSizeFe[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsFe[]={0.1, 0.2, 0.4, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
//Dose buildup factors Pb
|
||||
double BuildupDataPb[]={1., 1.01, 1.03, 1.06, 1.15, 1.16, 1.18, 1.19,
|
||||
1., 1.11, 1.17, 1.25, 1.34, 1.41, 1.5, 1.56,
|
||||
1., 1.17, 1.29, 1.46, 1.58, 1.72, 1.89, 2.02,
|
||||
1., 1.24, 1.42, 1.69, 2.00, 2.27, 2.65, 2.73,
|
||||
1., 1.37, 1.69, 2.26, 3.02, 3.74, 4.81, 5.86,
|
||||
1., 1.39, 1.76, 2.51, 3.66, 4.84, 6.87, 9.00,
|
||||
1., 1.34, 1.68, 2.43, 3.75, 5.30, 8.44, 12.3,
|
||||
1., 1.27, 1.56, 2.25, 3.61, 5.44, 9.80, 16.3,
|
||||
1., 1.21, 1.46, 2.08, 3.44, 5.55, 11.7, 23.6,
|
||||
1., 1.18, 1.40, 1.97, 3.34, 5.69, 13.8, 32.7,
|
||||
1., 1.14, 1.30, 1.74, 2.89, 5.07, 14.1, 44.6,
|
||||
1., 1.11, 1.23, 1.58, 2.52, 4.34, 12.5, 39.2};
|
||||
|
||||
int BuildupSizePb[]={12, // energy
|
||||
8}; // mu d
|
||||
double BuildupArgsPb[]={0.15, 0.30, 0.40, 0.5, 1.0, 2.0, 3.0, 4.0, 5.1, 6.0, 8.0, 10.0, // energy in MeV
|
||||
0., 1., 2., 4., 7., 10., 15., 20.0};
|
||||
|
||||
|
||||
|
||||
|
||||
// Flux-to-dose conversion
|
||||
/* //NRB-99?
|
||||
double FtoDArgs[]= {0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.80, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*0.752E-10, 3600*1.00E-10, 3600*1.51E-10, 3600*2.00E-10, 3600*2.47E-10, 3600*2.91E-10, 3600*3.73E-10, 3600*4.48E-10, 3600*7.49E-10, 3600*12.0E-10, 3600*16.0E-10, 3600*19.9E-10, 3600*23.8E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 13;
|
||||
int FtoDSize[] = {13};
|
||||
*/
|
||||
//ESS
|
||||
double FtoDArgs[]= {
|
||||
0.15,
|
||||
0.2,
|
||||
0.3,
|
||||
0.4,
|
||||
0.5,
|
||||
0.511,
|
||||
0.6,
|
||||
0.662,
|
||||
0.8,
|
||||
1.,
|
||||
1.12,
|
||||
1.33,
|
||||
1.5,
|
||||
2.,
|
||||
3.,
|
||||
4.,
|
||||
5.,
|
||||
6.,
|
||||
6.13,
|
||||
8.,
|
||||
10.,
|
||||
15.}; // in MeV
|
||||
double FtoDData[] = {
|
||||
2.69E-07,
|
||||
3.60E-07,
|
||||
5.44E-07,
|
||||
7.20E-07,
|
||||
8.89E-07,
|
||||
9.07E-07,
|
||||
1.05E-06,
|
||||
1.14E-06,
|
||||
1.34E-06,
|
||||
1.62E-06,
|
||||
1.76E-06,
|
||||
2.01E-06,
|
||||
2.20E-06,
|
||||
2.69E-06,
|
||||
3.51E-06,
|
||||
4.21E-06,
|
||||
4.82E-06,
|
||||
5.40E-06,
|
||||
5.47E-06,
|
||||
6.70E-06,
|
||||
7.92E-06,
|
||||
1.09E-05
|
||||
}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 22;
|
||||
int FtoDSize[] = {22};
|
||||
|
||||
|
||||
|
||||
// test input flux to dose
|
||||
/*
|
||||
double FtoDArgs[]= {0.15, 10.0}; // in meV
|
||||
double FtoDData[] = {3600*7.0E-10,3600*7.0E-10}; // from phot/s/m2 to mkSv/hr conversion
|
||||
int nFtoDgroups= 2;
|
||||
int FtoDSize[] = {2};
|
||||
*/
|
||||
#endif //end of insertion of datatables.
|
||||
/*** end of insertion to DECLARE section **/
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
/**Reading datafiles with capture per bin**/
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
double d1, d2, d3, d4;
|
||||
int ibin=0,NBINS;
|
||||
char line[1000],line1[1000];
|
||||
char dirname[1000];
|
||||
memset(dirname,0,sizeof(dirname));
|
||||
strcat(strcat(dirname,mcdirname),separator);
|
||||
|
||||
char filename[1000];
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
FILE *dataIn;
|
||||
|
||||
memset(filename,0,sizeof(filename));//resetting filename to NULL string
|
||||
strcat(strcat(filename,dirname),NiCaptureFile);
|
||||
|
||||
printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",NiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4) ibin++;
|
||||
};
|
||||
|
||||
|
||||
NBINS=ibin;
|
||||
|
||||
//printf("Shielding calculator:\n Input file contains %d bins.\n", NBINS);
|
||||
|
||||
double* zhat = malloc(NBINS*sizeof(double));
|
||||
double* Ihat= malloc(NBINS*sizeof(double));
|
||||
double* IhatNi= malloc(NBINS*sizeof(double));
|
||||
double* IhatTi= malloc(NBINS*sizeof(double));
|
||||
rewind(dataIn);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatNi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
memset(filename,0,sizeof(filename));// setting filename to NULL string
|
||||
strcat(strcat(filename,dirname),TiCaptureFile); //setting filename to TiCapture file
|
||||
printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TiCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
IhatTi[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
fclose(dataIn);
|
||||
|
||||
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),TotalCaptureFile);
|
||||
printf("Shielding calculator:\n Next input: %s\n",filename);
|
||||
|
||||
if( access( filename, R_OK ) != -1 ) {
|
||||
// file exists
|
||||
dataIn=fopen(filename,"r");
|
||||
} else {
|
||||
// file doesn't exist
|
||||
printf("Shielding calculator could not find file\n%s\nexiting.",TotalCaptureFile);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
|
||||
ibin=0;
|
||||
if (dataIn==NULL) {printf("Can't open file\n"); exit(0);};
|
||||
while (fgets(line, 1000,dataIn)!=NULL)
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
if (ibin!=NBINS){printf("Shielding calculator:\nMismatch in number of bins between input files.\nNo shielding calculation performed, exiting\n");
|
||||
exit(0);}
|
||||
|
||||
rewind(dataIn);
|
||||
|
||||
//printf("Shielding calculator:\n Reading file %s\n",filename);
|
||||
ibin=0;
|
||||
while (fgets(line, 1000,dataIn))
|
||||
{
|
||||
if (*line=='#') continue; // ignore comment line
|
||||
if(sscanf(line,"%le %le %le %le",&d1, &d2, &d3, &d4)==4)
|
||||
{
|
||||
zhat[ibin]=d1;
|
||||
Ihat[ibin]=d2;
|
||||
ibin++;
|
||||
};
|
||||
};
|
||||
//printf("Shielding calculator:\n Closing file %s...",filename);
|
||||
fclose(dataIn);
|
||||
//printf("done\n");
|
||||
|
||||
int i,ilayerConc=1,ilayerFe=1,ilayerFeConc=1,ilayerPb=1,iz,iEgroup; // Auxiliary variables
|
||||
|
||||
double zbinlength = zhat[2]-zhat[1];
|
||||
|
||||
|
||||
|
||||
double* RvaluesPb = malloc(NBINS*sizeof(double));
|
||||
double* RvaluesFe= malloc(NBINS*sizeof(double));
|
||||
double* Rvaluesbt= malloc(NBINS*sizeof(double));
|
||||
double* thicknessPb= malloc(NBINS*sizeof(double));
|
||||
double* thicknessFe= malloc(NBINS*sizeof(double));
|
||||
double* thicknessbt= malloc(NBINS*sizeof(double));
|
||||
double* RvaluesFebt= malloc(NBINS*sizeof(double));
|
||||
double* thicknessFebt= malloc(NBINS*sizeof(double));
|
||||
|
||||
double* doseFe= malloc(NBINS*sizeof(double));
|
||||
double* dosebt= malloc(NBINS*sizeof(double));
|
||||
double* doseFebt= malloc(NBINS*sizeof(double));
|
||||
double* dosePb= malloc(NBINS*sizeof(double));
|
||||
|
||||
|
||||
|
||||
|
||||
/* for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
doseFe[i]=0.0; dosebt[i]=0.0;doseFebt[i]=0.0;dosePb[i]=0.0;
|
||||
}
|
||||
*/
|
||||
|
||||
for ( iz=0; iz< NBINS; iz++)
|
||||
{
|
||||
double zc=zhat[iz]; // the z-value for the calculation
|
||||
double z;
|
||||
|
||||
//CONCRETE
|
||||
// ilayerConc may have already some value calculated for an upstream piece.
|
||||
// We start from it and if shielding is too thick, reduce this number.
|
||||
do {
|
||||
double thickbt = (ilayerConc)*stepbt;
|
||||
double Rbt = (0.5*Innerspace) + thickbt;
|
||||
Rvaluesbt[iz]=Rbt;
|
||||
dosebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
dBtmfp=lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*(Rbt-(0.5*Innerspace))*sqrt(Rbt*Rbt+z*z)/Rbt;
|
||||
double args[]={En,dBtmfp};
|
||||
dosebt[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dBtmfp)/(Rbt*Rbt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerConc--;
|
||||
} while ((dosebt[iz] < MaxRate)&&(ilayerConc>0));
|
||||
|
||||
//increase ilayerConc until the dose rate requirements are satisfied.
|
||||
do {
|
||||
double thickbt = (ilayerConc)*stepbt;
|
||||
double Rbt = (0.5*Innerspace) + thickbt;
|
||||
Rvaluesbt[iz]=Rbt;
|
||||
dosebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
dBtmfp=lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*(Rbt-(0.5*Innerspace))*sqrt(Rbt*Rbt+z*z)/Rbt;
|
||||
double args[]={En,dBtmfp};
|
||||
dosebt[iz] += (fracNi* IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dBtmfp)/(Rbt*Rbt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerConc++;
|
||||
} while ((dosebt[iz] >=MaxRate));
|
||||
|
||||
|
||||
|
||||
//CONCRETE+TUBING
|
||||
do{
|
||||
double thickFebt= (ilayerFeConc)*stepcomb;
|
||||
double RFebt=(0.5*Innerspace)+TFe+thickFebt;
|
||||
RvaluesFebt[iz]=RFebt;
|
||||
doseFebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
|
||||
dFebtmfp=(lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*thickFebt+
|
||||
lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*TFe)
|
||||
*sqrt(RFebt*RFebt+z*z)/RFebt;
|
||||
double args[]={En,dFebtmfp};
|
||||
|
||||
doseFebt[iz]+= (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFebtmfp)/(RFebt*RFebt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFeConc--;
|
||||
} while ((doseFebt[iz] < MaxRate)&&(ilayerFeConc>0));
|
||||
|
||||
|
||||
do {
|
||||
double thickFebt= (ilayerFeConc)*stepcomb;
|
||||
double RFebt=(0.5*Innerspace)+TFe+thickFebt;
|
||||
RvaluesFebt[iz]=RFebt;
|
||||
doseFebt[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp, dPbmfp, dBtmfp,dFebtmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup], fracB=fractionB[iEgroup];
|
||||
|
||||
dFebtmfp=(lint(&En,1,AttenuationSizeConc, AttenuationArgsConc, AttenuationDataConc)*thickFebt+
|
||||
lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*TFe)
|
||||
*sqrt(RFebt*RFebt+z*z)/RFebt;
|
||||
double args[]={En,dFebtmfp};
|
||||
|
||||
doseFebt[iz]+= (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatTi[i]-IhatNi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFebtmfp)/(RFebt*RFebt+z*z)*
|
||||
lint(args, 2, BuildupSizeConc,BuildupArgsConc,BuildupDataConc);
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFeConc++;
|
||||
} while ((doseFebt[iz] >=MaxRate));
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
//IRON
|
||||
do {
|
||||
double thickFe = (ilayerFe)*stepFe;
|
||||
double RFe = (0.5*Innerspace) + thickFe;
|
||||
RvaluesFe[iz]=RFe;
|
||||
doseFe[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
|
||||
|
||||
dFemfp=lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*thickFe*sqrt(RFe*RFe+z*z)/RFe;
|
||||
double args[]={En,dFemfp};
|
||||
doseFe[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFemfp)/(RFe*RFe+z*z)*
|
||||
lint(args, 2, BuildupSizeFe,BuildupArgsFe,BuildupDataFe);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFe--;
|
||||
} while ((doseFe[iz] <MaxRate)&&(ilayerFe>0));
|
||||
|
||||
|
||||
|
||||
do {
|
||||
double thickFe = (ilayerFe)*stepFe;
|
||||
double RFe = (0.5*Innerspace) + thickFe;
|
||||
RvaluesFe[iz]=RFe;
|
||||
doseFe[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dFemfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
dFemfp=lint(&En,1,AttenuationSizeFe, AttenuationArgsFe, AttenuationDataFe)*(RFe-(0.5*Innerspace))*sqrt(RFe*RFe+z*z)/RFe;
|
||||
double args[]={En,dFemfp};
|
||||
doseFe[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dFemfp)/(RFe*RFe+z*z)*
|
||||
lint(args, 2, BuildupSizeFe,BuildupArgsFe,BuildupDataFe);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerFe++;
|
||||
} while ((doseFe[iz] >=MaxRate));
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
//LEAD
|
||||
do {
|
||||
double thickPb = (ilayerPb)*stepPb;
|
||||
double RPb = (0.5*Innerspace) + thickPb;
|
||||
RvaluesPb[iz]=RPb;
|
||||
dosePb[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dPbmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
|
||||
dPbmfp=lint(&En,1,AttenuationSizePb, AttenuationArgsPb, AttenuationDataPb)*(RPb-(0.5*Innerspace))*sqrt(RPb*RPb+z*z)/RPb;
|
||||
double args[]={En,dPbmfp};
|
||||
dosePb[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dPbmfp)/(RPb*RPb+z*z)*
|
||||
lint(args, 2, BuildupSizePb,BuildupArgsPb,BuildupDataPb);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerPb--;
|
||||
} while ((dosePb[iz] <MaxRate)&&(ilayerPb>0));
|
||||
|
||||
|
||||
do {
|
||||
double thickPb = (ilayerPb)*stepPb;
|
||||
double RPb = (0.5*Innerspace) + thickPb;
|
||||
RvaluesPb[iz]=RPb;
|
||||
dosePb[iz]=0.0;
|
||||
/*computing the integral*/
|
||||
double dPbmfp;
|
||||
for (i=0;i<NBINS;i++){
|
||||
z=zc-zhat[i];
|
||||
|
||||
for (iEgroup=0; iEgroup<nEgroup;iEgroup++){
|
||||
double En = energy[iEgroup], fracNi=fractionNi[iEgroup],fracTi=fractionTi[iEgroup],fracB=fractionB[iEgroup];
|
||||
|
||||
dPbmfp=lint(&En,1,AttenuationSizePb, AttenuationArgsPb, AttenuationDataPb)*(RPb-(0.5*Innerspace))*sqrt(RPb*RPb+z*z)/RPb;
|
||||
double args[]={En,dPbmfp};
|
||||
dosePb[iz] += (fracNi*IhatNi[i]+fracTi*IhatTi[i]+fracB*(Ihat[i]-IhatNi[i]-IhatTi[i])) *0.25/3.1415926*
|
||||
lint(&En, 1, FtoDSize, FtoDArgs,FtoDData)* exp(-dPbmfp)/(RPb*RPb+z*z)*
|
||||
lint(args, 2, BuildupSizePb,BuildupArgsPb,BuildupDataPb);
|
||||
|
||||
} /* summing energy groups */
|
||||
} /* calculation of integral i=0...NBINS */
|
||||
ilayerPb++;
|
||||
} while ((dosePb[iz] >=MaxRate));
|
||||
|
||||
|
||||
} /* Position scan iz */
|
||||
|
||||
/* Selection of appropriate radius for the gamma shielding */
|
||||
for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
int index=0;
|
||||
thicknessFe[i]=RvaluesFe[i]-(0.5*Innerspace);
|
||||
thicknessbt[i]=Rvaluesbt[i]-(0.5*Innerspace);
|
||||
thicknessPb[i]=RvaluesPb[i]-(0.5*Innerspace);
|
||||
thicknessFebt[i]=RvaluesFebt[i]-(0.5*Innerspace)-TFe;
|
||||
}
|
||||
|
||||
#ifdef _WIN32
|
||||
#define separator "\\"
|
||||
#else
|
||||
#define separator "/"
|
||||
#endif
|
||||
|
||||
FILE *cost_out;
|
||||
memset(filename,0,sizeof(filename));
|
||||
strcat(strcat(filename,dirname),OutputFile);
|
||||
printf("Shielding calculator:\nWriting to file %s\n",filename);
|
||||
cost_out=fopen(filename,"w");
|
||||
|
||||
|
||||
fprintf(cost_out,"#Required shielding dimensions for the dose rate not to exceed %g uSv/hr at the surface with shielding inner dimensions %g meters\n#Position z,m\tNeutrons Ihat(z),n/s/bin\tNeutrons I(z),n/s/m\tCaptured by Ni /s/bin \tCapturedby Ni /s/m\tCaptured by Ti /s/bin\tCaptured by Ti /s/m\tThickness Pb, m\tThickness Fe, m\tThickness concrete, m (no tubing)\tThickness concrete, m (steel tubing %g cm)\tOuter radius Pb,m\tOuter radius Fe,m\tOuter radius concrete, m\tOuter radius concrete,m (with tubing)\n",MaxRate,Innerspace,100*TFe);
|
||||
|
||||
|
||||
for ( i=0; i< NBINS; i++)
|
||||
{
|
||||
fprintf( cost_out,"%g\t%.4g\t%.4g\t%.4g\t%.4g\t%.4g\t%.4g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",zhat[i],Ihat[i], Ihat[i]/zbinlength, IhatNi[i], IhatNi[i]/zbinlength, IhatTi[i],IhatTi[i]/zbinlength,thicknessPb[i], thicknessFe[i],thicknessbt[i],thicknessFebt[i], RvaluesPb[i],RvaluesFe[i],Rvaluesbt[i],RvaluesFebt[i]);
|
||||
}
|
||||
fclose(cost_out);
|
||||
free(zhat); free(Ihat); free(IhatNi); free(IhatTi);
|
||||
free(RvaluesPb);
|
||||
free(RvaluesFe);
|
||||
free(Rvaluesbt);
|
||||
free(thicknessPb);
|
||||
free(thicknessFe);
|
||||
free(thicknessbt);
|
||||
free(RvaluesFebt);
|
||||
free(thicknessFebt);
|
||||
|
||||
free(doseFe);
|
||||
free(dosebt);
|
||||
free(doseFebt);
|
||||
free(dosePb);
|
||||
|
||||
printf("Done with writing shielding\n");
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,210 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Updated by: Rodon Kolevatov
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the beginning of the region in trace in which pseudo-neutrons
|
||||
* are to be propagated. Pseudo-neutrons are neutrons which are generated by the function
|
||||
* <i>compute_func</i> from before and after SCATTER neutron states which have been logged by
|
||||
* a set of Scatter_logger/Scatter_logger_stop components.
|
||||
* N.B. This component should be immediately preceeded by an Arm. Any components between this one
|
||||
* and a subsequent Scatter_log_iterator_stop-component will be visited by the set of pseudo-neutrons
|
||||
* as if they were regular neutrons in the classical McStas-manner.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* compute_func [ ] : Address of the function that computes a psuedo neutron from before and after scatter states.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_Ni_new
|
||||
DEFINITION PARAMETERS (compute_func=NULL)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS (nstate_initial,s0,s1)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#ifndef M1THICKNESS
|
||||
#define M1THICKNESS 1500.0 //thickness of the m=1 coating
|
||||
#endif
|
||||
/*This is the specialized pseudo-neutron function that computes
|
||||
an escaping neutron from logged before and after SCATTER neutron states
|
||||
with weight corresponding to capture in the Nickel layers*/
|
||||
int exit_neutron_Ni(double *ns_tilde, struct Generalized_State_t *S0, struct Generalized_State_t *S1){
|
||||
/*!!Note that the transformation into global coordinate system must be done while logging
|
||||
as we do not have access to neither the component name nor can get the component rotation by index.*/
|
||||
Coords c1,c2;
|
||||
Rotation R1,R2;
|
||||
|
||||
/*so now compute the pseudo neutron state and possibly user variables*/
|
||||
/*momentum transfer at reflection. ASSUMES NO GRAVITY???*/
|
||||
double velocity=sqrt((S1->_vx)*(S1->_vx)+(S1->_vy)*(S1->_vy)+(S1->_vz)*(S1->_vz));
|
||||
double qtransf = V2Q*sqrt((S1->_vx-S0->_vx)*(S1->_vx-S0->_vx)+(S1->_vy-S0->_vy)*(S1->_vy-S0->_vy)+(S1->_vz-S0->_vz)*(S1->_vz-S0->_vz));
|
||||
double qinm = qtransf/0.0218;
|
||||
/*position comes from "new" state*/
|
||||
ns_tilde[0]=S1->_x;ns_tilde[1]=S1->_y;ns_tilde[2]=S1->_z;
|
||||
/*velocity is the "old" state*/
|
||||
ns_tilde[3]=S0->_vx;ns_tilde[4]=S0->_vy;ns_tilde[5]=S0->_vz;
|
||||
/*time from new*/
|
||||
ns_tilde[6]=S1->_t;
|
||||
/*spin comes from "new" state*/
|
||||
ns_tilde[7]=S1->_sx;ns_tilde[8]=S1->_sy;ns_tilde[9]=S1->_sz;
|
||||
/*weight of capture in Ni is determined analytically*/
|
||||
double captureprob;
|
||||
if ((S1->_mvalue<0)||(qinm<0.0001)) {ns_tilde[10]=0.0; captureprob=0.0;}
|
||||
else if (qinm<=1.01){
|
||||
/*q<qc(Ni), loss due to diffusion beyond the coating cutoff*/
|
||||
/* captureprob = 4.0582E-09*2200.0/velocity/((0.085787L*18.5 + 0.008321L*5.7)*1.0E-08+4.0582E-09*2200.0/velocity)*(S0->_p-S1->_p)/S0->_p; */
|
||||
double sigmadif ;//diffuse scattering cross section
|
||||
if (velocity > 1950.0){ sigmadif = 18.5; //lambda < 2A, totally incoherent scattering
|
||||
} else if (velocity < 1300.0){sigmadif = 5.2; // lambda >3A, totally coherent regime, only incoherent part contributes
|
||||
} else sigmadif=5.2+13.3*(1950-velocity)/(1950.0-1300.0);
|
||||
/*Now checking if the coating is m=1 or has higher m*/
|
||||
|
||||
if(S1->_mvalue<1.05) /*m=1 coating=> "triple thickness approximation" with Im k determined for reflectivity dip*/
|
||||
{
|
||||
double imk; /*imaginary part of momentum in the layer*/
|
||||
imk=1.e-8*(sigmadif*0.09121+0.41*2200.0/velocity)*velocity/3956.0*M1THICKNESS; /*v(m/s)x\lambda(AA)=3956*/
|
||||
captureprob = 4.49*2200.0/velocity/(sigmadif+4.49*2200.0/velocity)*(S0->_p-S1->_p)*(1-exp(-2.0*imk*M1THICKNESS*3.0))/S0->_p;
|
||||
} else { // reflection loss below qcNi in case of multilayer.
|
||||
/*Taking the minimum of what is for m=1 coating and what is when assume that neutron is physically lost beyond the cutoff.*/
|
||||
double imk; /*imaginary part of momentum in the layer*/
|
||||
imk=1.e-8*(sigmadif*0.09121+0.41*2200.0/velocity)*velocity/3956.0*M1THICKNESS; /*v(m/s)x\lambda(AA)=3956*/
|
||||
double c1, c2;
|
||||
c1 = 4.49*2200.0/velocity/(sigmadif+4.49*2200.0/velocity)*(S0->_p-S1->_p)*(1-exp(-2.0*imk*M1THICKNESS*3.0))/S0->_p;
|
||||
c2 = 0.005*(S1->_mvalue+0.1)*(S0->_p-S1->_p)/S0->_p; // for the matters of conservative estimate might be divided by 1-R at mvalue+0.1.
|
||||
captureprob = (c1>c2) ? c1 : c2 ;
|
||||
}
|
||||
} //end of "if (qinm<=1.02)"
|
||||
else if (qinm<=S1->_mvalue+0.1) { /*q>q_c(Ni) and reflection, absorption per hit */
|
||||
captureprob=0.005*qinm;
|
||||
} else { /* transmission beyond smirrorcutoff, some reserve for m=1 coating*/
|
||||
if(S1->_mvalue<1.05) captureprob=0.0025*1.3*1.3/qinm;
|
||||
else captureprob=0.0025*(S1->_mvalue+0.1)*(S1->_mvalue+0.1)/qinm;}
|
||||
/*check change in weight*/
|
||||
if((S0->_p-S1->_p)>1.e-5*S0->_p) ns_tilde[10]=S0->_p*captureprob; else ns_tilde[10]=0.0;
|
||||
/* if "mvalue" is -1, then absorption happened not on the coating but somewhere else. The weight of capture in nickel should be set to zero.*/
|
||||
|
||||
if(ns_tilde[10]>(S0->_p-S1->_p)) { printf("%f\t%f\t%f\t%E\t%E\t%f\t%E\n",velocity, qinm, S1->_mvalue, S0->_p, S1->_p, captureprob, ns_tilde[10]); exit(0);}
|
||||
return 0;
|
||||
}
|
||||
|
||||
#define NOABS \
|
||||
do {/* Nothing*/} while(0)
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int (*pseudo_neutron_state_function_Ni) (double *, struct Generalized_State_t *, struct Generalized_State_t *);
|
||||
|
||||
struct Generalized_State_t *s1,*s0;
|
||||
|
||||
double *nstate_initial;
|
||||
|
||||
int optics_not_hit;
|
||||
|
||||
/*need a pointer to the structure set up by the logger*/
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
if (compute_func) {
|
||||
pseudo_neutron_state_function_Ni=compute_func;
|
||||
}else{
|
||||
pseudo_neutron_state_function_Ni=exit_neutron_Ni;
|
||||
}
|
||||
nstate_initial=NULL;
|
||||
|
||||
optics_not_hit=0;
|
||||
|
||||
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
|
||||
//printf("Entering Iterator_Ni\n");
|
||||
/*I am the start of the pseudo neutron iterator*/
|
||||
if (nstate_initial==NULL){
|
||||
optics_not_hit=0; /* Fresh start, resetting variable */
|
||||
double *ns=nstate_initial=calloc(11,sizeof(double));
|
||||
ns[0]=x;ns[1]=y; ns[2]=z;
|
||||
ns[3]=vx;ns[4]=vy;ns[5]=vz;
|
||||
ns[6]=t;
|
||||
ns[7]=sx;ns[8]=sy;ns[9]=sz;
|
||||
ns[10]=p;
|
||||
|
||||
s0=Bounce_store;
|
||||
s1=Bounce_store+1;
|
||||
/* Remove std. ABSORB to avoid breaking analysis loop */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_iterator_stop
|
||||
|
||||
|
||||
}
|
||||
|
||||
//if (s1->_p!=-1){
|
||||
if (s1->_p!=-1 && s1->_p!=-2){
|
||||
/*a neutron weight of -1 is nonsensical. I.e. if s1->p>=0, it means there are two states in the log from which we can compute a pseudo-neutron*/
|
||||
double nstate[11];
|
||||
if ( pseudo_neutron_state_function_Ni(nstate,s0,s1) ){
|
||||
printf("Warning: (%s): error reported when computing pseudo neutron\n",NAME_CURRENT_COMP);
|
||||
}
|
||||
/*set neutron state for subsequent components*/
|
||||
x=nstate[0];y=nstate[1];z=nstate[2];
|
||||
vx=nstate[3];vy=nstate[4];vz=nstate[5];
|
||||
t=nstate[6];
|
||||
sx=nstate[7];sy=nstate[8];sz=nstate[9];
|
||||
p=nstate[10];
|
||||
s0++;
|
||||
s1++;
|
||||
// printf("Ni: z=%g\tp=%g\n",z,p);
|
||||
}else if (Bounce_store[1]._p==-1){
|
||||
/*This can happen now only in the case optics was not hit, i.e. no reflection, no absorption*/
|
||||
x=s0->_x;y=s0->_y;z=s0->_z;
|
||||
vx=s0->_vx;vy=s0->_vy;vz=s0->_vz;
|
||||
t=s0->_t;
|
||||
sx=s0->_sx;sy=s0->_sy;sz=s0->_sz;
|
||||
p=s0->_p;
|
||||
optics_not_hit = 1;
|
||||
// fprintf(stdout,"No interactions with optics. Use \"optics_not_hit\" variable to JUMP at scatter_logger_stop and avoid sending to ScatterLogIterator those neutrons which were actually not scattered.\n");
|
||||
} else {
|
||||
printf("Here happens what should not happen: s1->_p=%g, Bounce_store[1]._p=%g\n",s1->_p,Bounce_store[1]._p);
|
||||
fprintf(stderr,"This should not happen. Period.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,186 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the beginning of the region in trace in which pseudo-neutrons
|
||||
* are to be propagated. Pseudo-neutrons are neutrons which are generated by the function
|
||||
* <i>compute_func</i> from before and after SCATTER neutron states which have been logged by
|
||||
* a set of Scatter_logger/Scatter_logger_stop components.
|
||||
* N.B. This component should be immediately preceeded by an Arm. Any components between this one
|
||||
* and a subsequent Scatter_log_iterator_stop-component will be visited by the set of pseudo-neutrons
|
||||
* as if they were regular neutrons in the classical McStas-manner.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* compute_func [ ] : Address of the function that computes a psuedo neutron from before and after scatter states.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_Ti_new
|
||||
DEFINITION PARAMETERS (compute_func=NULL)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS (nstate_initial,s0,s1)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
/*This is the specialized pseudo-neutron function that computes
|
||||
an escaping neutron from logged before and after SCATTER neutron states
|
||||
with weights corresponding to capture in Ti layers*/
|
||||
int exit_neutron_Ti(double *ns_tilde, struct Generalized_State_t *S0, struct Generalized_State_t *S1){
|
||||
/*!!Note that the transformation into global coordinate system must be done while logging
|
||||
as we do not have access to neither the component name nor can get the component rotation by index.*/
|
||||
Coords c1,c2;
|
||||
Rotation R1,R2;
|
||||
|
||||
/*so now compute the pseudo neutron state and possibly user variables*/
|
||||
/*momentum transfer at reflection*/
|
||||
double velocity=sqrt((S1->_vx)*(S1->_vx)+(S1->_vy)*(S1->_vy)+(S1->_vz)*(S1->_vz));
|
||||
double qtransf = V2Q*sqrt((S1->_vx-S0->_vx)*(S1->_vx-S0->_vx)+(S1->_vy-S0->_vy)*(S1->_vy-S0->_vy)+(S1->_vz-S0->_vz)*(S1->_vz-S0->_vz));
|
||||
double qinm = qtransf/0.0218;
|
||||
/*position comes from "new" state*/
|
||||
ns_tilde[0]=S1->_x;ns_tilde[1]=S1->_y;ns_tilde[2]=S1->_z;
|
||||
/*velocity is the "old" state*/
|
||||
ns_tilde[3]=S0->_vx;ns_tilde[4]=S0->_vy;ns_tilde[5]=S0->_vz;
|
||||
/*time from new*/
|
||||
ns_tilde[6]=S1->_t;
|
||||
/*spin comes from "new" state*/
|
||||
ns_tilde[7]=S1->_sx;ns_tilde[8]=S1->_sy;ns_tilde[9]=S1->_sz;
|
||||
/*weight of capture in Ni is determined analytically*/
|
||||
double captureprob;
|
||||
/* if "mvalue" is -1, than absorption happened not on the coating but somewhere else. The weight of capture in nickel should be set to zero.*/
|
||||
if ((S1->_mvalue<0)||(qinm<0.0001)) {ns_tilde[10]=0.0; captureprob=0.0;}
|
||||
else if (S1->_mvalue<1.02){captureprob=0.0; //negligible probability for capture in Ti for m=1 coatings
|
||||
}
|
||||
else // now for coatings with m>1
|
||||
{if (qinm<=1.02){
|
||||
/*q<qc(Ni) means diffusion to higher divergences and interaction with coating there*/
|
||||
captureprob=0.0045*(S1->_mvalue-0.9)*(S0->_p-S1->_p)/S0->_p;
|
||||
} else if (qinm<=S1->_mvalue+0.1) { /*q>q_c(Ni) and reflection, absorption per hit */
|
||||
captureprob=0.0045*(qinm-1.0);
|
||||
} else { /* transmission beyond smirrorcutoff*/
|
||||
captureprob=0.00225*(S1->_mvalue-0.9)*(S1->_mvalue+0.1)/qinm;}}
|
||||
/*check change in weight*/
|
||||
if((S0->_p-S1->_p)>1.e-5*S0->_p) ns_tilde[10]=S0->_p*captureprob; else ns_tilde[10]=0.0;
|
||||
//printf("%f\t%f\t%f\t%E\t%E\t%f\t%E\n",velocity, qinm, S1->_mvalue, S0->_p, S1->_p, captureprob, ns_tilde[10]);
|
||||
return 0;
|
||||
}
|
||||
|
||||
#define NOABS \
|
||||
do {/* Nothing*/} while(0)
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int (*pseudo_neutron_state_function_Ti) (double *, struct Generalized_State_t *, struct Generalized_State_t *);
|
||||
|
||||
struct Generalized_State_t *s1,*s0;
|
||||
|
||||
double *nstate_initial;
|
||||
|
||||
int optics_not_hit;
|
||||
|
||||
/*need a pointer to the structure set up by the logger*/
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
if (compute_func) {
|
||||
pseudo_neutron_state_function_Ti=compute_func;
|
||||
}else{
|
||||
pseudo_neutron_state_function_Ti=exit_neutron_Ti;
|
||||
}
|
||||
nstate_initial=NULL;
|
||||
|
||||
optics_not_hit=0;
|
||||
|
||||
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
|
||||
//printf("Entering Iterator_Ti\n");
|
||||
/*I am the start of the pseudo neutron iterator*/
|
||||
if (nstate_initial==NULL){
|
||||
optics_not_hit=0; /* Fresh start, resetting variable */
|
||||
double *ns=nstate_initial=calloc(11,sizeof(double));
|
||||
ns[0]=x;ns[1]=y; ns[2]=z;
|
||||
ns[3]=vx;ns[4]=vy;ns[5]=vz;
|
||||
ns[6]=t;
|
||||
ns[7]=sx;ns[8]=sy;ns[9]=sz;
|
||||
ns[10]=p;
|
||||
|
||||
s0=Bounce_store;
|
||||
s1=Bounce_store+1;
|
||||
/* Remove std. ABSORB to avoid breaking analysis loop */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_iterator_stop
|
||||
|
||||
|
||||
}
|
||||
|
||||
//if (s1->_p!=-1){
|
||||
if (s1->_p!=-1 && s1->_p!=-2){
|
||||
/*a neutron weight of -1 is nonsensical. I.e. if s1->p>=0, it means there are two states in the log from which we can compute a pseudo-neutron*/
|
||||
double nstate[11];
|
||||
if ( pseudo_neutron_state_function_Ti(nstate,s0,s1) ){
|
||||
printf("Warning: (%s): error reported when computing pseudo neutron\n",NAME_CURRENT_COMP);
|
||||
}
|
||||
/*set neutron state for subsequent components*/
|
||||
x=nstate[0];y=nstate[1];z=nstate[2];
|
||||
vx=nstate[3];vy=nstate[4];vz=nstate[5];
|
||||
t=nstate[6];
|
||||
sx=nstate[7];sy=nstate[8];sz=nstate[9];
|
||||
p=nstate[10];
|
||||
s0++;
|
||||
s1++;
|
||||
// printf("Ti: z=%g\tp=%g\n",z,p);
|
||||
}else if (Bounce_store[1]._p==-1){
|
||||
/*This can happen now only in the case optics was not hit, i.e. no reflection, no absorption*/
|
||||
x=s0->_x;y=s0->_y;z=s0->_z;
|
||||
vx=s0->_vx;vy=s0->_vy;vz=s0->_vz;
|
||||
t=s0->_t;
|
||||
sx=s0->_sx;sy=s0->_sy;sz=s0->_sz;
|
||||
p=s0->_p;
|
||||
optics_not_hit = 1;
|
||||
// fprintf(stdout,"No interactions with optics. Use \"optics_not_hit\" variable to JUMP at scatter_logger_stop and avoid sending to ScatterLogIterator those neutrons which were actually not scattered.\n");
|
||||
} else {
|
||||
printf("Here happens what should not happen: s1->_p=%g, Bounce_store[1]._p=%g\n",s1->_p,Bounce_store[1]._p);
|
||||
fprintf(stderr,"This should not happen. Period.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,116 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator_stop.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration stop element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the end of the trace-region in which pseudo-neutrons are handled. Please see the Scatter_log_iterator-component for more details.
|
||||
* N.B. This component should be immediately followed by a construction like:
|
||||
* COMPONENT a1 = Arm()
|
||||
* AT (0,0,0) ABSOLUTE
|
||||
* JUMP a0 WHEN(MC_GETPAR(iterator_stop_comp_name,loop))
|
||||
*
|
||||
* This is to extract the value of the loop variable from the innards of this component
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* iterator [ ] Instance name of the Scatter_log_iterator log component preceeding this one.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_stop
|
||||
DEFINITION PARAMETERS (iterator)
|
||||
SETTING PARAMETERS (int last=0)
|
||||
OUTPUT PARAMETERS (loop)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int loop;
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
loop=1;
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
#ifdef scatter_iterator_stop
|
||||
#undef scatter_iterator_stop
|
||||
#endif
|
||||
#define scatter_iterator_stop iterator
|
||||
scatter_iterator_stop:
|
||||
loop=1;
|
||||
struct Generalized_State_t *s1=MC_GETPAR(iterator,s1);
|
||||
|
||||
if (s1->_p==-1){
|
||||
/*we have reached the end - unset loop and reset neutron state to whatever it was before we entered the pseudo neutron iterator*/
|
||||
loop=0;
|
||||
double *ns=MC_GETPAR(iterator,nstate_initial);
|
||||
x=ns[0];y=ns[1];z=ns[2];
|
||||
vx=ns[3];vy=ns[4];vz=ns[5];
|
||||
t=ns[6];
|
||||
sx=ns[7];sy=ns[8];sz=ns[9];
|
||||
p=ns[10];
|
||||
|
||||
free(ns);
|
||||
MC_GETPAR(iterator,nstate_initial)=NULL;
|
||||
/* Restore std ABSORB */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb mcabsorbAll
|
||||
// printf("No ABSORB after the iterator\n");
|
||||
} else if (s1->_p==-2) {
|
||||
/*if the weight equals -2 it means that the neutron was absorbed in components within the scatter_logger (see scatter logger definition)
|
||||
and should not be propagated further after this iterator*/
|
||||
// printf("Performing ABSORB after the iterator\n");
|
||||
/*we have reached the end - unset loop and reset neutron state to whatever it was before we entered the pseudo neutron iterator*/
|
||||
loop=0;
|
||||
double *ns=MC_GETPAR(iterator,nstate_initial);
|
||||
x=ns[0];y=ns[1];z=ns[2];
|
||||
vx=ns[3];vy=ns[4];vz=ns[5];
|
||||
t=ns[6];
|
||||
sx=ns[7];sy=ns[8];sz=ns[9];
|
||||
p=ns[10];
|
||||
|
||||
free(ns);
|
||||
MC_GETPAR(iterator,nstate_initial)=NULL;
|
||||
/* Restore std ABSORB */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb mcabsorbAll
|
||||
if (last>0) ABSORB; //Will sample new neutron.
|
||||
}
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,165 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_log_iterator.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen
|
||||
* Date: November 2012
|
||||
* Version: $Revision: 1.21 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics
|
||||
*
|
||||
* Iteration element for a Scatter_log
|
||||
*
|
||||
* %D
|
||||
*
|
||||
* This component marks the beginning of the region in trace in which pseudo-neutrons
|
||||
* are to be propagated. Pseudo-neutrons are neutrons which are generated by the function
|
||||
* <i>compute_func</i> from before and after SCATTER neutron states which have been logged by
|
||||
* a set of Scatter_logger/Scatter_logger_stop components.
|
||||
* N.B. This component should be immediately preceeded by an Arm. Any components between this one
|
||||
* and a subsequent Scatter_log_iterator_stop-component will be visited by the set of pseudo-neutrons
|
||||
* as if they were regular neutrons in the classical McStas-manner.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* compute_func [ ] : Address of the function that computes a psuedo neutron from before and after scatter states.
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_log_iterator_total
|
||||
DEFINITION PARAMETERS (compute_func=NULL)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS (nstate_initial,s0,s1)
|
||||
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
|
||||
|
||||
SHARE
|
||||
%{
|
||||
/*This is the specialized pseudo-neutron function that computes
|
||||
an escaping neutron from logged before and after SCATTER neutron states*/
|
||||
int exit_neutron(double *ns_tilde, struct Generalized_State_t *S0, struct Generalized_State_t *S1){
|
||||
/*!!Note that the transformation into global coordinate system must be done while logging
|
||||
as we do not have access to neither the component name nor can get the component rotation by index.*/
|
||||
Coords c1,c2;
|
||||
Rotation R1,R2;
|
||||
|
||||
/*so now compute the pseudo neutron state and possibly user variables*/
|
||||
/*position comes from "new" state*/
|
||||
ns_tilde[0]=S1->_x;ns_tilde[1]=S1->_y;ns_tilde[2]=S1->_z;
|
||||
/*velocity is the "old" state*/
|
||||
ns_tilde[3]=S0->_vx;ns_tilde[4]=S0->_vy;ns_tilde[5]=S0->_vz;
|
||||
/*time from new*/
|
||||
ns_tilde[6]=S1->_t;
|
||||
/*spin comes from "new" state*/
|
||||
ns_tilde[7]=S1->_sx;ns_tilde[8]=S1->_sy;ns_tilde[9]=S1->_sz;
|
||||
/*weight is difference old-new to mean the neutrons "deposited" in the guide wall*/
|
||||
ns_tilde[10]=S0->_p-S1->_p;
|
||||
return 0;
|
||||
}
|
||||
|
||||
#define NOABS \
|
||||
do {/* Nothing*/} while(0)
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int (*pseudo_neutron_state_function) (double *, struct Generalized_State_t *, struct Generalized_State_t *);
|
||||
|
||||
struct Generalized_State_t *s1,*s0;
|
||||
|
||||
double *nstate_initial;
|
||||
|
||||
int optics_not_hit;
|
||||
|
||||
/*need a pointer to the structure set up by the logger*/
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
if (compute_func) {
|
||||
pseudo_neutron_state_function=compute_func;
|
||||
}else{
|
||||
pseudo_neutron_state_function=exit_neutron;
|
||||
}
|
||||
nstate_initial=NULL;
|
||||
|
||||
optics_not_hit=0;
|
||||
|
||||
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
//printf("Entering Iterator_total\n");
|
||||
/*I am the start of the pseudo neutron iterator*/
|
||||
if (nstate_initial==NULL){
|
||||
optics_not_hit=0; /* Fresh start, resetting variable */
|
||||
double *ns=nstate_initial=calloc(11,sizeof(double));
|
||||
ns[0]=x;ns[1]=y; ns[2]=z;
|
||||
ns[3]=vx;ns[4]=vy;ns[5]=vz;
|
||||
ns[6]=t;
|
||||
ns[7]=sx;ns[8]=sy;ns[9]=sz;
|
||||
ns[10]=p;
|
||||
|
||||
s0=Bounce_store;
|
||||
s1=Bounce_store+1;
|
||||
/* Remove std. ABSORB to avoid breaking analysis loop */
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_iterator_stop
|
||||
|
||||
|
||||
}
|
||||
|
||||
//if (s1->_p!=-1){
|
||||
if (s1->_p!=-1 && s1->_p!=-2){
|
||||
/*a neutron weight of -1 is nonsensical. I.e. if s1->p>=0, it means there are two states in the log from which we can compute a pseudo-neutron*/
|
||||
double nstate[11];
|
||||
if ( pseudo_neutron_state_function(nstate,s0,s1) ){
|
||||
printf("Warning: (%s): error reported when computing pseudo neutron\n",NAME_CURRENT_COMP);
|
||||
}
|
||||
/*set neutron state for subsequent components*/
|
||||
x=nstate[0];y=nstate[1];z=nstate[2];
|
||||
vx=nstate[3];vy=nstate[4];vz=nstate[5];
|
||||
t=nstate[6];
|
||||
sx=nstate[7];sy=nstate[8];sz=nstate[9];
|
||||
p=nstate[10];
|
||||
s0++;
|
||||
s1++;
|
||||
// printf("Tot: z=%g\tp=%g\n",z,p);
|
||||
}else if (Bounce_store[1]._p==-1){
|
||||
/*This can happen now only in the case optics was not hit, i.e. no reflection, no absorption*/
|
||||
x=s0->_x;y=s0->_y;z=s0->_z;
|
||||
vx=s0->_vx;vy=s0->_vy;vz=s0->_vz;
|
||||
t=s0->_t;
|
||||
sx=s0->_sx;sy=s0->_sy;sz=s0->_sz;
|
||||
p=s0->_p;
|
||||
optics_not_hit = 1;
|
||||
// fprintf(stdout,"No interactions with optics. Use \"optics_not_hit\" variable to JUMP at scatter_logger_stop and avoid sending to ScatterLogIterator those neutrons which were actually not scattered.\n");
|
||||
} else {
|
||||
printf("Here happens what should not happen: s1->_p=%g, Bounce_store[1]._p=%g\n",s1->_p,Bounce_store[1]._p);
|
||||
fprintf(stderr,"This should not happen. Period.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
%}
|
||||
|
||||
MCDISPLAY
|
||||
%{
|
||||
/* A bit ugly; hard-coded dimensions. */
|
||||
magnify("");
|
||||
line(0,0,0,0.2,0,0);
|
||||
line(0,0,0,0,0.2,0);
|
||||
line(0,0,0,0,0,0.2);
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,217 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_logger.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen, Peter K Willendrup & Esben Klinkby
|
||||
* Date: January 2013
|
||||
* Version: $Revision: 1.12 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics / DTU Nutech
|
||||
*
|
||||
* Logging iteractions of neutrons with components
|
||||
*
|
||||
* %D
|
||||
* Start of the trace-region in which SCATTER events should be logged.
|
||||
* Whenever a SCATTER occurs in components between this one and its counterpart Scatter_logger_stop the neutron state is logged.
|
||||
* The log is kept in memory - but only for one single
|
||||
* numerical neutron at a time, so there should be no real danger of memory overrun.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* (none)
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_logger
|
||||
DEFINITION PARAMETERS ()
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
|
||||
// m_local_refl is to be used in shielding applications and accessed by shielding logger
|
||||
#ifndef MVALUE_LOCAL_IS_DEF
|
||||
#define MVALUE_LOCAL_IS_DEF 1
|
||||
double m_local_refl=-1;
|
||||
#endif
|
||||
|
||||
struct Generalized_State_t {
|
||||
double _x,_y,_z,_vx,_vy,_vz;
|
||||
double _p,_t,_sx,_sy,_sz;
|
||||
double _mvalue; // coating m-value in the place of the hit
|
||||
long long int _nid;
|
||||
int _comp, _idx;
|
||||
};
|
||||
|
||||
#define SCATTER_LOG do { \
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE-1){\
|
||||
struct Generalized_State_t *bp=&(Bounce_store[bounce_store_index]);\
|
||||
if( (bp-1)->_p!=p || (bp-1)->_vx!=vx || (bp-1)->_vy!=vy || (bp-1)->_vz!=vz ){\
|
||||
Coords ctmp=POS_A_CURRENT_COMP;\
|
||||
Coords _r = coords_set(x,y,z);\
|
||||
Coords _v = coords_set(vx,vy,vz);\
|
||||
Coords _s = coords_set(sx,sy,sz);\
|
||||
Coords _rg,_vg,_sg;\
|
||||
Rotation _Rt;\
|
||||
rot_transpose(ROT_A_CURRENT_COMP,_Rt);\
|
||||
_rg=coords_add(rot_apply(_Rt,_r),ctmp);\
|
||||
_vg=rot_apply(_Rt,_v);\
|
||||
_sg=rot_apply(_Rt,_s);\
|
||||
coords_get(_rg,&(bp->_x),&(bp->_y),&(bp->_z));\
|
||||
coords_get(_vg,&(bp->_vx),&(bp->_vy),&(bp->_vz));\
|
||||
coords_get(_sg,&(bp->_sx),&(bp->_sy),&(bp->_sz));\
|
||||
bp->_t=t;\
|
||||
bp->_p=p;\
|
||||
bp->_nid=mcget_run_num();\
|
||||
bp->_comp=INDEX_CURRENT_COMP;\
|
||||
bp->_idx=bounce_store_index;\
|
||||
/* New: registering m-value at reflection. it is set to -1 by component the if we missed coating */\
|
||||
bp->_mvalue=m_local_refl;\
|
||||
/* printf("Recording scattering, writing state (%d) to the buffer, r: %g %g %g v: %g %g %g p:%g\n",\
|
||||
bounce_store_index, bp->_x, bp->_y, bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_p);*/\
|
||||
bounce_store_index++;\
|
||||
}\
|
||||
}else if(bounce_store_index==(BOUNCE_LOG_SIZE-1) && !bounce_store_overrun){\
|
||||
printf("Warning (%s): Scatter_log overrun at %llu - not logging any more events\n",NAME_CURRENT_COMP,mcget_run_num());\
|
||||
bounce_store_overrun=1;\
|
||||
}\
|
||||
do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcScattered++;} while(0);\
|
||||
} while(0)
|
||||
|
||||
#define ABSORB_LOG do { /*printf("DOING ABSORB_LOG\n");*/ \
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE-1){\
|
||||
struct Generalized_State_t *bp=&(Bounce_store[bounce_store_index]);\
|
||||
/* if( (bp-1)->_p!=p || (bp-1)->_vx!=vx || (bp-1)->_vy!=vy || (bp-1)->_vz!=vz )*//*<-- Check removed, works wrong if neutron is absorbed at the entrance*/{\
|
||||
Coords ctmp=POS_A_CURRENT_COMP;\
|
||||
Coords _r = coords_set(x,y,z);\
|
||||
Coords _v = coords_set(vx,vy,vz);\
|
||||
Coords _s = coords_set(sx,sy,sz);\
|
||||
Coords _rg,_vg,_sg;\
|
||||
Rotation _Rt;\
|
||||
rot_transpose(ROT_A_CURRENT_COMP,_Rt);\
|
||||
_rg=coords_add(rot_apply(_Rt,_r),ctmp);\
|
||||
_vg=rot_apply(_Rt,_v);\
|
||||
_sg=rot_apply(_Rt,_s);\
|
||||
coords_get(_rg,&(bp->_x),&(bp->_y),&(bp->_z));\
|
||||
coords_get(_vg,&(bp->_vx),&(bp->_vy),&(bp->_vz));\
|
||||
/*bp->_vx=0.; bp->_vy=0.; bp->_vz=0.;*/\
|
||||
/*vx=0; vy=0;vz=0;*/\
|
||||
coords_get(_sg,&(bp->_sx),&(bp->_sy),&(bp->_sz));\
|
||||
bp->_t=t;\
|
||||
bp->_p=0.;\
|
||||
bp->_nid=mcget_run_num();\
|
||||
bp->_comp=INDEX_CURRENT_COMP;\
|
||||
bp->_idx=bounce_store_index;\
|
||||
bp->_mvalue=m_local_refl;\
|
||||
/* printf("Recording absorption event, writing state (%d) to the buffer, r: %g %g %g v: %g %g %g p:%g\n",\
|
||||
bounce_store_index, bp->_x, bp->_y, bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_p);*/\
|
||||
bounce_store_index++;\
|
||||
}\
|
||||
}else if(bounce_store_index==(BOUNCE_LOG_SIZE-1) && !bounce_store_overrun){\
|
||||
printf("Warning (%s): Scatter_log overrun at %llu - not logging any more events\n",NAME_CURRENT_COMP,mcget_run_num());\
|
||||
bounce_store_overrun=1;\
|
||||
}\
|
||||
absorbed_in_optics=1;\
|
||||
do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcDEBUG_ABSORB(); MAGNET_OFF; goto mcabsorb;} while(0);\
|
||||
} while(0)
|
||||
|
||||
|
||||
#define SCATTER0\
|
||||
do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcScattered++;} while(0)
|
||||
|
||||
#define ABSORB0\
|
||||
do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcDEBUG_ABSORB(); MAGNET_OFF; goto mcabsorb;} while(0)
|
||||
|
||||
|
||||
|
||||
const int BOUNCE_LOG_SIZE=10000000;
|
||||
struct Generalized_State_t *Bounce_store;
|
||||
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int bounce_store_index;
|
||||
int bounce_store_overrun;
|
||||
int absorbed_in_optics;
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
bounce_store_index=0;
|
||||
absorbed_in_optics=0;
|
||||
Bounce_store=malloc(sizeof(*Bounce_store)*BOUNCE_LOG_SIZE);
|
||||
Bounce_store[BOUNCE_LOG_SIZE-1]._p=-1;
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
#undef SCATTER
|
||||
#define SCATTER SCATTER_LOG
|
||||
|
||||
#undef ABSORB
|
||||
#define ABSORB ABSORB_LOG
|
||||
|
||||
/*
|
||||
#ifdef scatter_logger_stop
|
||||
#undef scatter_logger_stop
|
||||
#endif
|
||||
|
||||
#define scatter_logger_stop logger
|
||||
*/
|
||||
#undef mcabsorb
|
||||
#define mcabsorb scatter_logger_stop
|
||||
|
||||
bounce_store_index=0;/*we are now starting logging so we should start afresh*/
|
||||
absorbed_in_optics=0;/*we are now starting logging so we should start afresh*/
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE){
|
||||
struct Generalized_State_t *bp=&(Bounce_store[bounce_store_index]);
|
||||
Coords ctmp=POS_A_CURRENT_COMP;
|
||||
Coords r = coords_set(x,y,z);
|
||||
Coords v = coords_set(vx,vy,vz);
|
||||
Coords s = coords_set(sx,sy,sz);
|
||||
|
||||
Coords rg,vg,sg;
|
||||
|
||||
Rotation Rt;
|
||||
rot_transpose(ROT_A_CURRENT_COMP,Rt);
|
||||
rg=coords_add(rot_apply(Rt,r),ctmp);
|
||||
vg=rot_apply(Rt,v);
|
||||
sg=rot_apply(Rt,s);
|
||||
coords_get(rg,&(bp->_x),&(bp->_y),&(bp->_z));
|
||||
coords_get(vg,&(bp->_vx),&(bp->_vy),&(bp->_vz));
|
||||
coords_get(sg,&(bp->_sx),&(bp->_sy),&(bp->_sz));
|
||||
bp->_t=t;
|
||||
bp->_p=p;
|
||||
bp->_nid=mcget_run_num();
|
||||
bp->_comp=INDEX_CURRENT_COMP;
|
||||
bp->_idx=bounce_store_index;
|
||||
bp->_mvalue=0.0; // This is the incoming particle. No collision, setting mvalue=0.
|
||||
// printf("Starting scatter logger, writing state (%d) to the buffer, r: %g %g %g v: %g %g %g p: %g\n", bounce_store_index, bp->_x, bp->_y, bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_p);
|
||||
bounce_store_index++;
|
||||
}else if(bounce_store_index==BOUNCE_LOG_SIZE && !bounce_store_overrun){
|
||||
printf("Warning (%s): Scatter_log overrun - not logging any more SCATTER events\n",NAME_CURRENT_COMP);
|
||||
bounce_store_overrun=1;
|
||||
}
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
%}
|
||||
|
||||
END
|
||||
@@ -0,0 +1,121 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* McStas, neutron ray-tracing package
|
||||
* Copyright 1997-2002, All rights reserved
|
||||
* Risoe National Laboratory, Roskilde, Denmark
|
||||
* Institut Laue Langevin, Grenoble, France
|
||||
*
|
||||
* Component: Scatter_logger_stop.comp
|
||||
*
|
||||
* %I
|
||||
*
|
||||
* Written by: Erik B Knudsen, Peter K Willendrup & Esben Klinkby
|
||||
* Date: January 2013
|
||||
* Version: $Revision: 1.12 $
|
||||
* Release: McStas 2.1
|
||||
* Origin: DTU Physics / DTU Nutech
|
||||
*
|
||||
* Stop logging iteractions of neutrons with components
|
||||
*
|
||||
* %D
|
||||
* Component which marks the end of the region where SCATTER events should be logged.
|
||||
*
|
||||
* %P
|
||||
* Input parameters:
|
||||
*
|
||||
* logger: The Scatter_logger.comp which began the logging region. This is necessary to allow communication between the components. ( )
|
||||
*
|
||||
* %E
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE COMPONENT Shielding_logger_stop
|
||||
DEFINITION PARAMETERS (logger)
|
||||
SETTING PARAMETERS ()
|
||||
OUTPUT PARAMETERS ()
|
||||
|
||||
SHARE
|
||||
%{
|
||||
#define SCATTER0\
|
||||
do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcScattered++;} while(0)
|
||||
#define ABSORB0\
|
||||
do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
|
||||
mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcDEBUG_ABSORB(); MAGNET_OFF; goto mcabsorb;} while(0)
|
||||
%}
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
int bounce_store_index;
|
||||
int bounce_store_overrun;
|
||||
int logger_buffer_cleared;
|
||||
%}
|
||||
|
||||
INITIALIZE
|
||||
%{
|
||||
#ifndef logger
|
||||
fprintf(stderr,"Error(%s): Logger undefined - can't stop noexisting logger\n", NAME_CURRENT_COMP);
|
||||
#endif
|
||||
logger_buffer_cleared=0;
|
||||
%}
|
||||
|
||||
|
||||
TRACE
|
||||
%{
|
||||
|
||||
#undef SCATTER
|
||||
#define SCATTER SCATTER0
|
||||
|
||||
#undef ABSORB
|
||||
#define ABSORB ABSORB0
|
||||
|
||||
#undef mcabsorb
|
||||
#define mcabsorb mcabsorbAll
|
||||
|
||||
#ifdef scatter_logger_stop
|
||||
#undef scatter_logger_stop
|
||||
#endif
|
||||
|
||||
#define scatter_logger_stop logger
|
||||
scatter_logger_stop:
|
||||
|
||||
if (bounce_store_index<BOUNCE_LOG_SIZE){
|
||||
struct Generalized_State_t *bp;
|
||||
bp=&(Bounce_store[bounce_store_index]);
|
||||
bp->_x=0;
|
||||
bp->_y=0;
|
||||
bp->_z=0;
|
||||
bp->_vx=0;
|
||||
bp->_vy=0;
|
||||
bp->_vz=0;
|
||||
bp->_sx=0;
|
||||
bp->_sy=0;
|
||||
bp->_sz=0;
|
||||
bp->_t=0;
|
||||
if(absorbed_in_optics==1) bp->_p=-2; else bp->_p=-1;
|
||||
bp->_nid=0;
|
||||
bp->_comp=0;
|
||||
bounce_store_index++;
|
||||
}else if(bounce_store_index==BOUNCE_LOG_SIZE && !bounce_store_overrun){
|
||||
printf("Warning (%s): Scatter_log overrun - cannot set stop bit. Aborting\n",NAME_CURRENT_COMP);
|
||||
exit(1);
|
||||
}
|
||||
%}
|
||||
|
||||
FINALLY
|
||||
%{
|
||||
/*If we are a logger stop we should also free the memory*/
|
||||
int i;
|
||||
struct Generalized_State_t *bp;
|
||||
for (i=0;i<BOUNCE_LOG_SIZE;i++){
|
||||
bp=&(Bounce_store[i]);
|
||||
/* maybe this should be in the share block - must be able to discern between different loggers*/
|
||||
/* printf("SCATTERLOG: %d %lld %g %g %g %g %g %g %g %g %g %g %g %d\n", \*/
|
||||
/* i,bp->_nid,bp->_x,bp->_y,bp->_z, bp->_vx, bp->_vy, bp->_vz, bp->_t, \*/
|
||||
/* bp->_sx, bp->_sy, bp->_sz, bp->_p, bp->_comp);*/
|
||||
}
|
||||
if (!logger_buffer_cleared) { free(Bounce_store); logger_buffer_cleared=1;}
|
||||
|
||||
/*do nothing if the memory has already been freed*/
|
||||
%}
|
||||
|
||||
END
|
||||
+18000
File diff suppressed because it is too large
Load Diff
+18000
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@@ -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.
Executable
+52
@@ -0,0 +1,52 @@
|
||||
#!/usr/bin/env python
|
||||
#-*- coding: utf8 -*-
|
||||
|
||||
import sys, os
|
||||
from subprocess import call
|
||||
from numpy import *
|
||||
from scipy.optimize import leastsq
|
||||
|
||||
|
||||
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 '
|
||||
|
||||
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__':
|
||||
pass
|
||||
+108
-104
@@ -1,117 +1,121 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
ncount=4e9
|
||||
use_cores=${SLURM_NPROCS:-6}
|
||||
sample_length=0.05
|
||||
ncount=1e9
|
||||
use_cores=6
|
||||
sample=4
|
||||
omega=0.8
|
||||
sample_length=0.005
|
||||
sample_height=0.01
|
||||
lambda_start=3.0
|
||||
lambda_end=32.0
|
||||
|
||||
###################### Brilliance Transfer 10x10mm² and 1x1mm² VS ####################
|
||||
bash compile_if_needed.sh
|
||||
|
||||
DESTi=$DEST/brilliance_nowindow_5x10
|
||||
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.000 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
enable_windows=0 \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0
|
||||
|
||||
DESTi=$DEST/brilliance_5x10
|
||||
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.000 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0
|
||||
|
||||
|
||||
ncount=1e10
|
||||
sample_length=0.001
|
||||
sample_height=0.001
|
||||
|
||||
DESTi=$DEST/brilliance_1x1
|
||||
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.000 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0
|
||||
|
||||
|
||||
# ###################### Reference and Ni-layer measurement 10x10mm² sample ####################
|
||||
ncount=1e10
|
||||
sample_length=0.01
|
||||
sample_height=0.01
|
||||
lambda_start=3.25
|
||||
lambda_end=25.5
|
||||
lambda_end=12.75
|
||||
sample=1
|
||||
|
||||
bash compile_if_needed.sh
|
||||
echo "Running with $use_cores threads"
|
||||
|
||||
# ###################### Reference and sample measurement pulse skipping ####################
|
||||
|
||||
#
|
||||
# Reference in pulse skipping mode
|
||||
omega=1.5
|
||||
DESTi=$DEST/mooc_reference_ps
|
||||
echo "Reference, pulse skipping"
|
||||
echo "Results saved to $DESTi"
|
||||
omega=1.0
|
||||
DESTi=$DEST/reference_10x10_10
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0002 \
|
||||
sample=1 sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3
|
||||
#
|
||||
# Sample models in pulse skipping mode
|
||||
for MODEL in 0 1 2 3 4 5
|
||||
do
|
||||
omega=1.5
|
||||
DESTi="$DEST/mooc_model_ps_$MODEL"
|
||||
echo "Model $MODEL, pulse skipping"
|
||||
echo "Results saved to $DESTi"
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0002 \
|
||||
sample=2 sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3 \
|
||||
sample_file="alessandra_model_2_nores_00$MODEL.dat"
|
||||
|
||||
done
|
||||
|
||||
#
|
||||
# Sample models in pulse skipping mode, larger angle
|
||||
for MODEL in 0 1 2 3 4 5
|
||||
do
|
||||
omega=4.0
|
||||
DESTi="$DEST/mooc_model_ps2_$MODEL"
|
||||
echo "Model $MODEL, pulse skipping"
|
||||
echo "Results saved to $DESTi"
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0002 \
|
||||
sample=2 sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3 \
|
||||
sample_file="alessandra_model_2_nores_00$MODEL.dat"
|
||||
|
||||
done
|
||||
|
||||
# ###################### Reference and samples in conventional mode ####################
|
||||
#
|
||||
#
|
||||
#
|
||||
# Reference in conventional mode
|
||||
ncount=4e9
|
||||
lambda_end=13.5
|
||||
theta_resolution=0.1
|
||||
|
||||
omega=1.5
|
||||
DESTi=$DEST/mooc_reference_tof
|
||||
echo "Reference, omega=$omega, Delta omega=$theta_resolution"
|
||||
echo "Results saved to $DESTi"
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus \
|
||||
omegaa=$omega operationmode=1 theta_resolution=$theta_resolution over_illumination=0.0002 \
|
||||
sample=1 sample_length=$sample_length sample_height=$sample_height \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.000 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
#
|
||||
# Sample models in pulse skipping mode
|
||||
for MODEL in 0 1 2 3 4 5
|
||||
do
|
||||
for omega in 0.3 0.75 1.875 4.6875 11.71875
|
||||
do
|
||||
ncount=$(bc <<< "5000000000/$omega/$omega")
|
||||
theta_resolution=$(bc <<< "scale=10;0.05*$omega")
|
||||
DESTi="$DEST/mooc_model_tof_$MODEL-$omega"
|
||||
echo "Model $MODEL, omega=$omega, Delta omega=$theta_resolution, n=$ncount"
|
||||
echo "Results saved to $DESTi"
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
# Ni Sample
|
||||
|
||||
ncount=6e9
|
||||
sample=2
|
||||
omega=0.8
|
||||
DESTi=$DEST/nickle_10x10_08
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --ncount=$ncount --gravitation --format=NeXus \
|
||||
omegaa=$omega operationmode=1 theta_resolution=$theta_resolution over_illumination=0.0002 \
|
||||
sample=2 sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1 \
|
||||
sample_file="alessandra_model_2_nores_00$MODEL.dat"
|
||||
done
|
||||
done
|
||||
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.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
ncount=2e9
|
||||
omega=3.0
|
||||
DESTi=$DEST/nickle_10x10_30
|
||||
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.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
|
||||
ncount=1e9
|
||||
omega=8.0
|
||||
DESTi=$DEST/nickle_10x10_80
|
||||
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.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
# ###################### Reference and Ni-layer measurement 1x1mm² sample ####################
|
||||
#
|
||||
#
|
||||
|
||||
Executable → Regular
+13
-5
@@ -1,15 +1,15 @@
|
||||
#!/bin/tcsh
|
||||
#SBATCH -J mooc_McEstia
|
||||
#SBATCH -J mcEstia_S
|
||||
#SBATCH -N 2
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#SBATCH --time=1-00:00:00
|
||||
#SBATCH --time=2: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
|
||||
#SBATCH --partition=ll_short
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
@@ -17,7 +17,15 @@ echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
|
||||
#module unload intel
|
||||
#module load intel
|
||||
module load mcstas
|
||||
bash compile_if_needed.sh
|
||||
|
||||
mpirun -np $SLURM_NPROCS Estia_baseline.out \
|
||||
--dir=../results/shielding --ncount=1e9 \
|
||||
sample=0
|
||||
|
||||
|
||||
echo "Program finished with exit code $? at: `date`"
|
||||
|
||||
bash run_simu.sh
|
||||
Reference in New Issue
Block a user