Compare commits
10 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 glavic_a
|
f7bc0d7e56 | ||
|
glavic_a
|
da90617aa4 | ||
|
glavic_a
|
95344a175f | ||
|
glavic_a
|
e85cf16d42 | ||
|
glavic_a
|
83d051ff2f | ||
|
glavic_a
|
80ff1e6785 | ||
|
glavic_a
|
8434fecdde | ||
|
glavic_a
|
9697c9a584 | ||
|
glavic_a
|
7d32fb2313 | ||
|
glavic_a
|
d3964c2883 |
Executable
+10
@@ -0,0 +1,10 @@
|
||||
#!/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/tof_VS_list_p_x_y_hd_vd_t_L/events:GZIP=5 -l /entry1/data/tof_VS_list_p_x_y_hd_vd_t_L/events:CHUNK=3072x7
|
||||
h5repack -i $fi -o $fi.c -f /entry1/data/tof_sample_list_p_x_y_hd_vd_t_L/events:GZIP=5 -l /entry1/data/tof_sample_list_p_x_y_hd_vd_t_L/events:CHUNK=3072x7
|
||||
mv $fi.c $fi
|
||||
done
|
||||
@@ -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
+218
-113
@@ -48,18 +48,27 @@
|
||||
* %End
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE INSTRUMENT ESS_reflectometer_Estia
|
||||
(double omegaa = 1.2, int sample = 0, double sample_length = 0.01, double sample_height = 0.01,
|
||||
int operationmode = 0, double over_illumination = 0.0, double theta_resolution = 0.04,
|
||||
double lambda_min = 3.75, double lambda_start = 3.0, double lambda_end = 12.0,
|
||||
int enable_chopper = 0, int enable_gravity=0, int enable_windows=1,
|
||||
int enable_polarizer = 0, int enable_analyzer = 0,
|
||||
double source_power = 2,
|
||||
double selene1_foot1y =0.0, double selene1_foot2y = 0.0
|
||||
)
|
||||
DEFINE INSTRUMENT ESS_reflectometer_Estia(
|
||||
int enable_gravity=0, int enable_polarizer = 0, int enable_analyzer = 0)
|
||||
|
||||
DECLARE
|
||||
%{
|
||||
// Fix parameters that are variables in normal simulation
|
||||
double omegaa = 0.0;
|
||||
int sample = 4;
|
||||
double sample_length = 0.020;
|
||||
double sample_height = 0.015;
|
||||
int operationmode = 0;
|
||||
double over_illumination = 0.0;
|
||||
double theta_resolution = 0.04;
|
||||
int enable_windows=1;
|
||||
double source_power = 5;
|
||||
int enable_chopper = 0;
|
||||
|
||||
double lambda_start = 0.1;
|
||||
double lambda_end = 30.1;
|
||||
double lambda_min = 3.75;
|
||||
|
||||
/* Geometrical parameters from CAD model of Estia (ESS-0050413)
|
||||
* TCS coordinate and directional rotation first focus point
|
||||
* refered to as focus_moderator_y_rot
|
||||
@@ -77,6 +86,8 @@ double selene1_foot2 = 7.00; // distance of second foot to VS focus
|
||||
double selene1_center;
|
||||
double selene1_shift;
|
||||
double selene1_rot;
|
||||
double selene1_foot1y =0.0;
|
||||
double selene1_foot2y = 0.0;
|
||||
|
||||
// Selene 2
|
||||
|
||||
@@ -216,41 +227,51 @@ COMPONENT moderator = ESS_butterfly(
|
||||
* 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_in = Slit(xwidth=0.027, yheight=0.067)
|
||||
AT (0, 0, -0.810) 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
|
||||
//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
|
||||
COMPONENT CPC1_out = Slit(xwidth=0.0137, yheight=0.03373)
|
||||
AT (0, 0, -0.270) RELATIVE arm_virtual_source_beam
|
||||
|
||||
COMPONENT VS_PSD = PSD_monitor(
|
||||
filename = "VS_PSD", nx=100, ny=100,
|
||||
xwidth=0.05, yheight = 0.05, restore_neutron=1)
|
||||
AT (0, 0, 0) RELATIVE arm_virtual_source_beam
|
||||
|
||||
/* bandwidth definition chopper */
|
||||
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==1
|
||||
AT (0, 0, chopper_pos-2*NBOA_c) RELATIVE arm_virtual_source_beam
|
||||
COMPONENT VS_L = L_monitor(
|
||||
filename = "VS_L", nL=300,
|
||||
Lmin=lambda_start, Lmax=lambda_end,
|
||||
xwidth=0.05, yheight = 0.05, restore_neutron=1)
|
||||
AT (0, 0, 0) RELATIVE arm_virtual_source_beam
|
||||
|
||||
// COMPONENT VS_MCPL = MCPL_output(filename="VS")
|
||||
// AT (0, 0, 0) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT CPC2_in = Slit(xwidth=0.00692, yheight=0.01506)
|
||||
AT (0, 0, 0.077) RELATIVE arm_virtual_source_beam
|
||||
|
||||
COMPONENT CPC2_out = Slit(xwidth=0.0198, yheight=0.04914)
|
||||
AT (0, 0, 0.595) RELATIVE arm_virtual_source_beam
|
||||
|
||||
COMPONENT VS_MCPL = MCPL_output(filename="selene1")
|
||||
AT (0, 0, 2.2) RELATIVE arm_selene1
|
||||
|
||||
|
||||
/* 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
|
||||
/******* Start of scatter logger. Logging scatterings and absorption in the guide system. *********/
|
||||
COMPONENT log_P_start=Shielding_logger()
|
||||
AT (0,0,0) RELATIVE arm_virtual_source_beam
|
||||
EXTEND
|
||||
%{
|
||||
#ifdef scatter_logger_stop
|
||||
#undef scatter_logger_stop
|
||||
#endif
|
||||
#define scatter_logger_stop log_P_start
|
||||
%}
|
||||
|
||||
|
||||
// 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
|
||||
|
||||
|
||||
/*************************************
|
||||
@@ -262,6 +283,156 @@ COMPONENT virtual_source_BL = Slit(
|
||||
|
||||
%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*/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**********************************
|
||||
* Optical components within cave *
|
||||
**********************************/
|
||||
@@ -281,86 +452,20 @@ 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
|
||||
COMPONENT SMPL_PSD = PSD_monitor(
|
||||
filename = "SMPL_PSD", nx=100, ny=100,
|
||||
xwidth=0.05, yheight = 0.05, restore_neutron=1)
|
||||
AT (0, 0, 0) RELATIVE arm_sample_beam
|
||||
ROTATED (0, 0, 0) RELATIVE arm_sample_beam
|
||||
|
||||
COMPONENT SMPL_L = L_monitor(
|
||||
filename = "SMPL_L", nL=300,
|
||||
Lmin=lambda_start, Lmax=lambda_end,
|
||||
xwidth=0.05, yheight = 0.05, restore_neutron=1)
|
||||
AT (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(
|
||||
xwidth = sample_length, yheight = sample_height,
|
||||
center = 1, transmit = 0,
|
||||
R0 = 0.999, alpha = 0.001, m = 50, center = 1, transmit = 0
|
||||
)
|
||||
WHEN sample==1
|
||||
AT (0, 0, 0) RELATIVE arm_sample
|
||||
ROTATED (0, 90, 0) RELATIVE arm_sample
|
||||
|
||||
/* Nickel film on silicon */
|
||||
COMPONENT ni_sample = Mirror(
|
||||
xwidth = sample_length, yheight = sample_height,
|
||||
center = 1, transmit = 0,
|
||||
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
|
||||
|
||||
/***********************************************************************/
|
||||
|
||||
|
||||
@@ -34,7 +34,7 @@ DECLARE
|
||||
// NBOA
|
||||
// Entrance Window
|
||||
double NBOA_Al_entrance_start= 1.850;
|
||||
double NBOA_Al_entrance_length= 0.001;
|
||||
double NBOA_Al_entrance_length= 0.0015;
|
||||
// First segment with Cu-shielding in beginning
|
||||
double E02_01_01_Cu_start= 1.854;
|
||||
double E02_01_01_Cu_length= 0.345;
|
||||
@@ -75,7 +75,7 @@ DECLARE
|
||||
|
||||
// Exit Window
|
||||
double NBOA_Al_exit_start= 5.3396;
|
||||
double NBOA_Al_exit_length= 0.0049;
|
||||
double NBOA_Al_exit_length= 0.0060;
|
||||
|
||||
|
||||
// in-bunker feeder
|
||||
@@ -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,
|
||||
@@ -281,6 +281,19 @@ COMPONENT NBOA_Al_window_exit=Al_window(
|
||||
WHEN enable_windows
|
||||
AT (0,0,NBOA_Al_exit_start) RELATIVE ISCS
|
||||
|
||||
COMPONENT BBG_PSD = PSD_monitor(
|
||||
filename = "BBG_PSD", nx=500, ny=500,
|
||||
xwidth=0.25, yheight = 0.25, restore_neutron=1)
|
||||
AT (0,0,NBOA_Al_exit_start+0.1) RELATIVE ISCS
|
||||
|
||||
COMPONENT BBG_L = L_monitor(
|
||||
filename = "BBG_L", nL=300,
|
||||
Lmin=lambda_start, Lmax=lambda_end,
|
||||
xwidth=0.25, yheight = 0.25, restore_neutron=1)
|
||||
AT (0,0,NBOA_Al_exit_start+0.1) RELATIVE ISCS
|
||||
|
||||
|
||||
|
||||
// Aluminum Window at light shutter flight tube entrance
|
||||
COMPONENT LS_Al_window_in=Al_window(
|
||||
thickness=LS_Al_in_length)
|
||||
@@ -340,7 +353,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 +381,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 +418,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,
|
||||
|
||||
@@ -78,6 +78,19 @@ COMPONENT polarizer1 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, re
|
||||
ROTATED (0,0,90.0) RELATIVE arm_polarizer
|
||||
|
||||
|
||||
COMPONENT MF_PSD = PSD_monitor(
|
||||
filename = "MF_PSD", nx=100, ny=100,
|
||||
xwidth=0.05, yheight = 0.05, restore_neutron=1)
|
||||
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT MF_L = L_monitor(
|
||||
filename = "MF_L", nL=300,
|
||||
Lmin=lambda_start, Lmax=lambda_end,
|
||||
xwidth=0.05, yheight = 0.05, restore_neutron=1)
|
||||
AT (0, 0, 2*selene_c) RELATIVE arm_selene1
|
||||
|
||||
|
||||
|
||||
/* The proximal polariser comes next */
|
||||
COMPONENT polarizer2 = Polariser(nIncRefr=1, d_substrate = 5e-4, reflect_d=0, reflect_u=0, lin=polarizer_start, length=polarizer_length,
|
||||
delta_theta=(Theta1_polarizer+Theta2_polarizer)*PI/180.0, h2=0.1, h1=0.05, abs_ref=1, m_u=4.0, m_d=0.65, both_coated=1, alpha=2.3, W = 0.0014)
|
||||
|
||||
@@ -49,7 +49,7 @@
|
||||
*******************************************************************************/
|
||||
|
||||
DEFINE INSTRUMENT ESS_reflectometer_Estia
|
||||
(double lambda_start = 3.0, double lambda_end = 12.0,
|
||||
(double lambda_start = 0.1, double lambda_end = 30.1,
|
||||
int enable_gravity=1, int enable_windows=1, direct_beam=0,
|
||||
double source_power = 5, foil_thickness=0.00001)
|
||||
|
||||
@@ -164,34 +164,66 @@ COMPONENT moderator = ESS_butterfly(
|
||||
* 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_in = Slit(xwidth=0.027, yheight=0.067)
|
||||
AT (0, 0, -0.810) 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
|
||||
COMPONENT CPC1_monitor = Slit(xwidth=0.01574, yheight=0.03683)
|
||||
AT (0, 0, -0.3536) 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,
|
||||
COMPONENT Vanadium_Foil = Incoherent(focus_r=0.05, 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
|
||||
ROTATED (45.0, 0, 0) RELATIVE arm_virtual_source_beam
|
||||
|
||||
COMPONENT monitor_shielding = Slit(radius=0.01)
|
||||
WHEN direct_beam==0
|
||||
AT (0, 0, 0.030) RELATIVE arm_monitor
|
||||
|
||||
COMPONENT monitor_cossection = Slit(radius=0.025/2.0)
|
||||
WHEN direct_beam==0
|
||||
AT (0, 0, 0.030+0.01449) RELATIVE arm_monitor
|
||||
|
||||
COMPONENT tof_monitor = TOFLambda_monitor(
|
||||
filename = "monitor",
|
||||
tmin=0, tmax=120000, nt=1200,
|
||||
Lmin=0,Lmax=35,nL=350,
|
||||
filename = "monitor", restore_neutron=1,
|
||||
tmin=0, tmax=120000, nt=1201,
|
||||
Lmin=0.1,Lmax=30.1,nL=301,
|
||||
xwidth = 0.025, yheight = 0.025)
|
||||
WHEN direct_beam==0
|
||||
AT (0, 0, 0.1) RELATIVE arm_monitor
|
||||
AT (0, 0, 0.030+0.0145) RELATIVE arm_monitor
|
||||
|
||||
// two additional monitors at center and end of He tube, can be used to model
|
||||
// tof resolution for wavelength dependent absorption length
|
||||
COMPONENT monitor_cossection_2 = Slit(radius=0.025/2.0)
|
||||
WHEN direct_beam==0
|
||||
AT (0, 0, 0.030+0.01449+0.1778/2.0) RELATIVE arm_monitor
|
||||
|
||||
COMPONENT tof_monitor_2 = TOFLambda_monitor(
|
||||
filename = "monitor_2", restore_neutron=1,
|
||||
tmin=0, tmax=120000, nt=1201,
|
||||
Lmin=0.1,Lmax=30.1,nL=301,
|
||||
xwidth = 0.025, yheight = 0.025)
|
||||
WHEN direct_beam==0
|
||||
AT (0, 0, 0.030+0.0145+0.1778/2.0) RELATIVE arm_monitor
|
||||
|
||||
COMPONENT monitor_cossection_3 = Slit(radius=0.025/2.0)
|
||||
WHEN direct_beam==0
|
||||
AT (0, 0, 0.030+0.01449+0.1778) RELATIVE arm_monitor
|
||||
|
||||
COMPONENT tof_monitor_3 = TOFLambda_monitor(
|
||||
filename = "monitor_3",
|
||||
tmin=0, tmax=120000, nt=1201,
|
||||
Lmin=0.1,Lmax=30.1,nL=301,
|
||||
xwidth = 0.025, yheight = 0.025)
|
||||
WHEN direct_beam==0
|
||||
AT (0, 0, 0.030+0.0145+0.1778) RELATIVE arm_monitor
|
||||
|
||||
COMPONENT tof_direct = TOFLambda_monitor(
|
||||
filename = "VS",
|
||||
tmin=0, tmax=120000, nt=1200,
|
||||
Lmin=0,Lmax=35,nL=350,
|
||||
tmin=0, tmax=120000, nt=1201,
|
||||
Lmin=0.1,Lmax=30.1,nL=301,
|
||||
xwidth = 0.05, yheight = 0.05)
|
||||
WHEN direct_beam>0
|
||||
AT (0, 0, 0.08) RELATIVE arm_virtual_source_beam
|
||||
|
||||
@@ -113,7 +113,7 @@ COMPONENT block_before_selene_guide_1 = Absorber(
|
||||
AT (0, 0, selene_distance-0.0095) RELATIVE arm_selene1
|
||||
|
||||
/* Selene 1 elliptic guide */
|
||||
COMPONENT E02_03_01 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -122,7 +122,7 @@ COMPONENT E02_03_01 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_02 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -131,7 +131,7 @@ COMPONENT E02_03_02 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_03 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -140,7 +140,7 @@ COMPONENT E02_03_03 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_04 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -149,7 +149,7 @@ COMPONENT E02_03_04 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_05 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -158,7 +158,7 @@ COMPONENT E02_03_05 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_06 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -167,7 +167,7 @@ COMPONENT E02_03_06 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_07 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -176,7 +176,7 @@ COMPONENT E02_03_07 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_08 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -185,7 +185,7 @@ COMPONENT E02_03_08 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_09 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -194,7 +194,7 @@ COMPONENT E02_03_09 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_10 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -203,7 +203,7 @@ COMPONENT E02_03_10 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_11 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -212,7 +212,7 @@ COMPONENT E02_03_11 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_12 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -221,7 +221,7 @@ COMPONENT E02_03_12 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_13 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -230,7 +230,7 @@ COMPONENT E02_03_13 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_14 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -239,7 +239,7 @@ COMPONENT E02_03_14 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene1
|
||||
|
||||
COMPONENT E02_03_15 = Elliptic_guide_gravity(
|
||||
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,
|
||||
|
||||
@@ -54,7 +54,7 @@ COMPONENT block_before_selene_guide_2 = Absorber(
|
||||
|
||||
|
||||
/* Selene 2 elliptic guide first half */
|
||||
COMPONENT E02_04_01 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -63,7 +63,7 @@ COMPONENT E02_04_01 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_02 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -72,7 +72,7 @@ COMPONENT E02_04_02 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+1*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_03 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -81,7 +81,7 @@ COMPONENT E02_04_03 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+2*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_04 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -90,7 +90,7 @@ COMPONENT E02_04_04 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+3*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_05 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -99,7 +99,7 @@ COMPONENT E02_04_05 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+4*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_06 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -108,7 +108,7 @@ COMPONENT E02_04_06 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+5*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_07 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -117,7 +117,7 @@ COMPONENT E02_04_07 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+6*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_08 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -126,7 +126,7 @@ COMPONENT E02_04_08 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+7*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_09 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -135,7 +135,7 @@ COMPONENT E02_04_09 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+8*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_10 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -144,7 +144,7 @@ COMPONENT E02_04_10 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+9*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_11 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -153,7 +153,7 @@ COMPONENT E02_04_11 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+10*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_12 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -162,7 +162,7 @@ COMPONENT E02_04_12 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+11*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_13 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -171,7 +171,7 @@ COMPONENT E02_04_13 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+12*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_14 = Elliptic_guide_gravity(
|
||||
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,
|
||||
@@ -180,7 +180,7 @@ COMPONENT E02_04_14 = Elliptic_guide_gravity(
|
||||
enableGravity=enable_gravity)
|
||||
AT (0, 0, selene_distance+13*selene_segment) RELATIVE arm_selene2
|
||||
|
||||
COMPONENT E02_04_15 = Elliptic_guide_gravity(
|
||||
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,
|
||||
|
||||
@@ -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
|
||||
@@ -7,5 +7,14 @@ if [ Estia_baseline.instr -nt Estia_baseline.out ] || [ ! -f 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
|
||||
mpicc -O3 -o Estia_baseline.out Estia_baseline.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus \
|
||||
-I/afs/psi.ch/project/sinq/sl6-64/mcstas2.4/mcstas/2.4/libs/mcpl \
|
||||
-L/afs/psi.ch/project/sinq/sl6-64/mcstas2.4/mcstas/2.4/libs/mcpl -lmcpl
|
||||
fi
|
||||
|
||||
if [ Estia_monitor.instr -nt Estia_monitor.out ] || [ ! -f Estia_monitor.out ] \
|
||||
|| [ Estia_feeder.instr -nt Estia_monitor.out ]; then
|
||||
rm Estia_monitor.c Estia_monitor.out
|
||||
mcstas -o Estia_monitor.c Estia_monitor.instr
|
||||
mpicc -O3 -o Estia_monitor.out Estia_monitor.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
|
||||
fi
|
||||
|
||||
@@ -1,4 +0,0 @@
|
||||
foreach z ( `seq -0.5 0.125 0.5` )
|
||||
echo $z
|
||||
./run_simu.sh $z
|
||||
end
|
||||
Executable
+16
@@ -0,0 +1,16 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
ncount=1e9
|
||||
use_cores=$1
|
||||
|
||||
###################### Brilliance Transfer 10x10mm² and 1x1mm² VS ####################
|
||||
|
||||
DESTi=$DEST/mcnp_reference
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
enable_gravity=0 enable_polarizer=0 enable_analyzer=0
|
||||
@@ -0,0 +1,24 @@
|
||||
#!/bin/tcsh
|
||||
#SBATCH -J McEstia
|
||||
#SBATCH -N 3
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#SBATCH --time=12:00:00
|
||||
#SBATCH --mail-type=fail
|
||||
#SBATCH --mail-user=artur.glavic@psi.ch
|
||||
|
||||
#SBATCH -o stdout.log
|
||||
#SBATCH -e stderr.log
|
||||
|
||||
#SBATCH --partition=ll_medium
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
|
||||
module load mcstas
|
||||
|
||||
bash compile_if_needed.sh
|
||||
bash run_mcnpref.sh $SLURM_NPROCS
|
||||
@@ -1,12 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
mcstas -o Estia_monitor.c Estia_monitor.instr
|
||||
mpicc -O3 -o Estia_monitor.out Estia_monitor.c -lm -DUSE_MPI
|
||||
|
||||
|
||||
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0
|
||||
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_ref --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=2
|
||||
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor_trans --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=1
|
||||
|
||||
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor100 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.0001
|
||||
mpirun -np 6 Estia_monitor.out --ncount=1e8 --dir=../results/monitor001 --gravitation enable_windows=1 lambda_start=1.0 lambda_end=35.0 direct_beam=0 foil_thickness=0.000001
|
||||
@@ -1,213 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
ncount=2e7
|
||||
use_cores=4
|
||||
sample=4
|
||||
omega=0.8
|
||||
sample_length=0.01
|
||||
sample_height=0.01
|
||||
lambda_start=3.0
|
||||
lambda_end=32.0
|
||||
|
||||
###################### Brilliance Transfer 10x10mm² and 1x1mm² VS ####################
|
||||
if [ Estia_baseline.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
|
||||
fi
|
||||
#
|
||||
# if [ Estia_baseline_ana1.instr -nt Estia_baseline_ana1.out ]; then
|
||||
# rm Estia_baseline_ana1.c Estia_baseline_ana1.out
|
||||
# mcstas -o Estia_baseline_ana1.c Estia_baseline_ana1.instr
|
||||
# mpicc -O3 -o Estia_baseline_ana1.out Estia_baseline_ana1.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
|
||||
# fi
|
||||
#
|
||||
# if [ Estia_baseline_ana2.instr -nt Estia_baseline_ana2.out ]; then
|
||||
# rm Estia_baseline_ana2.c Estia_baseline_ana2.out
|
||||
# mcstas -o Estia_baseline_ana2.c Estia_baseline_ana2.instr
|
||||
# mpicc -O3 -o Estia_baseline_ana2.out Estia_baseline_ana2.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
|
||||
# fi
|
||||
|
||||
|
||||
###################### Reference and Ni-layer measurement 10x10mm² sample ####################
|
||||
# ncount=1e10
|
||||
sample_length=0.01
|
||||
sample_height=0.01
|
||||
lambda_start=3.5
|
||||
lambda_end=30.0
|
||||
sample=1
|
||||
|
||||
# omega=1.0
|
||||
# DESTi=$DEST/pol_ref_10x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=0 enable_analyzer=0
|
||||
#
|
||||
# DESTi=$DEST/pol1_10x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=1 enable_analyzer=0
|
||||
#
|
||||
# DESTi=$DEST/pol2_10x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=2 enable_analyzer=0
|
||||
# #
|
||||
# omega=7.0
|
||||
# DESTi=$DEST/pol1_10x10_70
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=1 enable_analyzer=0
|
||||
#
|
||||
# DESTi=$DEST/pol2_10x10_70
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=2 enable_analyzer=0
|
||||
#
|
||||
|
||||
# ncount=1e8
|
||||
omega=0.0
|
||||
# DESTi=$DEST/pol1_10x50_70
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --ncount=$ncount --gravitation \
|
||||
# omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0025 \
|
||||
# sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
# lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=1 enable_analyzer=0
|
||||
|
||||
DESTi=$DEST/pol2_10x50_70
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --ncount=$ncount --gravitation \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0025 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=0 enable_polarizer=2 enable_analyzer=0
|
||||
#
|
||||
# ncount=5e9
|
||||
# sample_length=0.01
|
||||
# omega=1.0
|
||||
# sample_length=0.003
|
||||
# DESTi=$DEST/pol_ref_3x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=0 enable_analyzer=0
|
||||
#
|
||||
# DESTi=$DEST/pol1_3x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=1 enable_analyzer=0
|
||||
#
|
||||
# DESTi=$DEST/pol2_3x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=2 enable_analyzer=0
|
||||
#
|
||||
#
|
||||
# sample_height=0.01
|
||||
# sample_length=0.01
|
||||
# DESTi=$DEST/ana1_10x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=0 enable_analyzer=1
|
||||
|
||||
# sample_height=0.001
|
||||
# DESTi=$DEST/ana1_10x1_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline_ana1.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=0 enable_analyzer=1
|
||||
#
|
||||
# sample_height=0.01
|
||||
# DESTi=$DEST/ana2_10x10_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=0 enable_analyzer=2
|
||||
|
||||
# sample_height=0.001
|
||||
# DESTi=$DEST/ana2_10x1_10
|
||||
# if [ -e "$DESTi" ]; then
|
||||
# rm -r "$DESTi"
|
||||
# fi
|
||||
#
|
||||
# mpirun -np $use_cores Estia_baseline_ana2.out \
|
||||
# --dir="$DESTi" --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 enable_polarizer=0 enable_analyzer=2
|
||||
|
||||
|
||||
###################### Reference and Ni-layer measurement 1x1mm² sample ####################
|
||||
|
||||
|
||||
@@ -1,52 +0,0 @@
|
||||
#!/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
|
||||
@@ -1,121 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
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=12.75
|
||||
sample=1
|
||||
|
||||
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" --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
|
||||
|
||||
# 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" --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 ####################
|
||||
#
|
||||
#
|
||||
@@ -1,86 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
use_cores=6
|
||||
|
||||
###################### Compile model if necessary ####################
|
||||
bash compile_if_neede.sh
|
||||
|
||||
######## Reference and Ni-layer conventional measurement 10x10mm² sample #############
|
||||
ncount=1e10
|
||||
sample_length=0.01
|
||||
sample_height=0.01
|
||||
lambda_start=3.25
|
||||
lambda_min=3.75
|
||||
lambda_end=12.75
|
||||
sample=1
|
||||
|
||||
omega=1.0
|
||||
DESTi=$DEST/tof_reference_10x10_10
|
||||
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=1 theta_resolution=0.01 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
# Ni Sample
|
||||
ncount=2e10
|
||||
sample=2
|
||||
omega=0.5
|
||||
DESTi=$DEST/tof_nickle_10x10_05
|
||||
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=1 theta_resolution=0.03 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
ncount=1e10
|
||||
omega=1.2
|
||||
DESTi=$DEST/tof_nickle_10x10_12
|
||||
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=1 theta_resolution=0.03 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
ncount=5e9
|
||||
omega=3.0
|
||||
DESTi=$DEST/tof_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=1 theta_resolution=0.03 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
ncount=2e9
|
||||
omega=7.5
|
||||
DESTi=$DEST/tof_nickle_10x10_75
|
||||
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=1 theta_resolution=0.03 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_min=$lambda_min lambda_end=$lambda_end enable_gravity=1 enable_chopper=1
|
||||
|
||||
###################### Reference and Ni-layer measurement 1x1mm² sample ####################
|
||||
|
||||
|
||||
@@ -1,46 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
ncount=1e9
|
||||
use_cores=6
|
||||
sample=4
|
||||
omega=0.8
|
||||
sample_length=0.01
|
||||
sample_height=0.01
|
||||
lambda_start=3.0
|
||||
lambda_end=32.0
|
||||
|
||||
bash compile_if_neede.sh
|
||||
|
||||
|
||||
################# Reference and Ni-layer measurement 2 pulse skipping ####################
|
||||
ncount=4e9
|
||||
sample_length=0.01
|
||||
sample_height=0.01
|
||||
lambda_start=3.9
|
||||
lambda_end=24.0
|
||||
sample=1
|
||||
omega=2.0
|
||||
|
||||
DESTi=$DEST/single_skip_reference
|
||||
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.0002 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3
|
||||
|
||||
sample=2
|
||||
DESTi=$DEST/single_skip_nickle
|
||||
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.0002 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=1 enable_chopper=3
|
||||
@@ -1,119 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
ncount=1e9
|
||||
use_cores=6
|
||||
sample=1
|
||||
omega=15.0
|
||||
sample_length=0.02
|
||||
sample_height=0.01
|
||||
|
||||
lambda_start=3.5
|
||||
lambda_min=5.0
|
||||
lambda_end=15.5
|
||||
frame_usage=0.92
|
||||
|
||||
###################### Brilliance Transfer 10x10mm² and 1x1mm² VS ####################
|
||||
bash compile_if_neede.sh
|
||||
|
||||
###################### Reference and Ni-layer measurement 10x10mm² sample ####################
|
||||
|
||||
frame_usage=0.92
|
||||
DESTi=$DEST/chopper_f1_092
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1 \
|
||||
frame_usage=$frame_usage lambda_min=$lambda_min
|
||||
|
||||
frame_usage=0.94
|
||||
DESTi=$DEST/chopper_f1_094
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1 \
|
||||
frame_usage=$frame_usage lambda_min=$lambda_min
|
||||
|
||||
frame_usage=0.96
|
||||
DESTi=$DEST/chopper_f1_096
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1 \
|
||||
frame_usage=$frame_usage lambda_min=$lambda_min
|
||||
|
||||
frame_usage=0.974
|
||||
DESTi=$DEST/chopper_f1_097s
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1 \
|
||||
frame_usage=$frame_usage lambda_min=$lambda_min
|
||||
|
||||
frame_usage=0.98
|
||||
DESTi=$DEST/chopper_f1_098
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1 \
|
||||
frame_usage=$frame_usage lambda_min=$lambda_min
|
||||
|
||||
frame_usage=0.99
|
||||
DESTi=$DEST/chopper_f1_099
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1 \
|
||||
frame_usage=$frame_usage lambda_min=$lambda_min
|
||||
|
||||
|
||||
frame_usage=0.98
|
||||
lambda_start=2.75
|
||||
lambda_min=3.75
|
||||
lambda_end=12.75
|
||||
DESTi=$DEST/chopper_f1_375A
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_baseline.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0001 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
lambda_start=$lambda_start lambda_end=$lambda_end enable_gravity=0 enable_chopper=1 \
|
||||
frame_usage=$frame_usage lambda_min=$lambda_min
|
||||
|
||||
|
||||
###################### Reference and Ni-layer measurement 1x1mm² sample ####################
|
||||
|
||||
|
||||
@@ -1,32 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
DEST=../results
|
||||
ncount=1e8
|
||||
use_cores=6
|
||||
sample=1
|
||||
omega=0.8
|
||||
sample_length=0.01
|
||||
sample_height=0.015
|
||||
|
||||
if [ Estia_vs.instr -nt Estia_vs.out ]; then
|
||||
rm Estia_vs.c Estia_vs.out
|
||||
mcstas -t -o Estia_vs.c Estia_vs.instr
|
||||
mpicc -O3 -o Estia_vs.out Estia_vs.c -lm -DUSE_MPI -DUSE_NEXUS -lNeXus
|
||||
fi
|
||||
|
||||
|
||||
omega=2.0
|
||||
DESTi=$DEST/vs_data
|
||||
if [ -e "$DESTi" ]; then
|
||||
rm -r "$DESTi"
|
||||
fi
|
||||
|
||||
mpirun -np $use_cores Estia_vs.out \
|
||||
--dir="$DESTi" --format=NeXuS --ncount=$ncount --gravitation \
|
||||
omegaa=$omega operationmode=0 theta_resolution=0.04 over_illumination=0.0025 \
|
||||
sample=$sample sample_length=$sample_length sample_height=$sample_height \
|
||||
enable_gravity=1 enable_chopper=0 source_power=5
|
||||
|
||||
###################### Reference and Ni-layer measurement 1x1mm² sample ####################
|
||||
|
||||
|
||||
@@ -1,40 +0,0 @@
|
||||
#!/bin/bash
|
||||
#SBATCH -J mcEstia_1
|
||||
#SBATCH -N 2
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#SBATCH --time=1-00:00:00
|
||||
#SBATCH --mail-type=fail
|
||||
#SBATCH --mail-user=artur.glavic@psi.ch
|
||||
|
||||
#SBATCH -o stdout.log
|
||||
#SBATCH -e stderr.log
|
||||
|
||||
#SBATCH --partition=ll_long
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
/usr/bin/modulecmd tcsh load mcstas
|
||||
bash compile_if_needed.sh
|
||||
|
||||
for y1 in $(seq -0.009 0.001 0.009)
|
||||
do
|
||||
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
|
||||
do
|
||||
printf -v Y1 %.3f $y1
|
||||
printf -v Y2 %.3f $y2
|
||||
echo $Y1 $Y2
|
||||
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
|
||||
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
|
||||
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
|
||||
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
|
||||
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
|
||||
done
|
||||
done
|
||||
|
||||
|
||||
echo "Program finished with exit code $? at: `date`"
|
||||
|
||||
@@ -0,0 +1,31 @@
|
||||
#!/bin/tcsh
|
||||
#SBATCH -J mcEstia_S
|
||||
#SBATCH -N 2
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#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_short
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
#module unload intel
|
||||
#module load intel
|
||||
module load mcstas
|
||||
bash compile_if_needed.sh
|
||||
|
||||
mpirun -np $SLURM_NPROCS Estia_baseline.out \
|
||||
--dir=../results/shielding_2 --ncount=5e9 \
|
||||
enable_gravity=0
|
||||
|
||||
|
||||
echo "Program finished with exit code $? at: `date`"
|
||||
|
||||
@@ -1,40 +0,0 @@
|
||||
#!/bin/bash
|
||||
#SBATCH -J mcEstia_2
|
||||
#SBATCH -N 2
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#SBATCH --time=1-00:00:00
|
||||
#SBATCH --mail-type=fail
|
||||
#SBATCH --mail-user=artur.glavic@psi.ch
|
||||
|
||||
#SBATCH -o stdout2.log
|
||||
#SBATCH -e stderr2.log
|
||||
|
||||
#SBATCH --partition=ll_long
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
/usr/bin/modulecmd tcsh load mcstas
|
||||
bash compile_if_needed.sh
|
||||
|
||||
for y1 in $(seq -0.09 0.01 -0.01)
|
||||
do
|
||||
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
|
||||
do
|
||||
printf -v Y1 %.3f $y1
|
||||
printf -v Y2 %.3f $y2
|
||||
echo $Y1 $Y2
|
||||
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
|
||||
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
|
||||
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
|
||||
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
|
||||
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
|
||||
done
|
||||
done
|
||||
|
||||
|
||||
echo "Program finished with exit code $? at: `date`"
|
||||
|
||||
@@ -1,40 +0,0 @@
|
||||
#!/bin/bash
|
||||
#SBATCH -J mcEstia_3
|
||||
#SBATCH -N 2
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#SBATCH --time=1-00:00:00
|
||||
#SBATCH --mail-type=fail
|
||||
#SBATCH --mail-user=artur.glavic@psi.ch
|
||||
|
||||
#SBATCH -o stdout3.log
|
||||
#SBATCH -e stderr3.log
|
||||
|
||||
#SBATCH --partition=ll_long
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
/usr/bin/modulecmd tcsh load mcstas
|
||||
bash compile_if_needed.sh
|
||||
|
||||
for y1 in $(seq 0.01 0.01 0.09)
|
||||
do
|
||||
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
|
||||
do
|
||||
printf -v Y1 %.3f $y1
|
||||
printf -v Y2 %.3f $y2
|
||||
echo $Y1 $Y2
|
||||
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
|
||||
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
|
||||
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
|
||||
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
|
||||
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
|
||||
done
|
||||
done
|
||||
|
||||
|
||||
echo "Program finished with exit code $? at: `date`"
|
||||
|
||||
@@ -1,40 +0,0 @@
|
||||
#!/bin/bash
|
||||
#SBATCH -J mcEstia_4
|
||||
#SBATCH -N 2
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#SBATCH --time=1-00:00:00
|
||||
#SBATCH --mail-type=fail
|
||||
#SBATCH --mail-user=artur.glavic@psi.ch
|
||||
|
||||
#SBATCH -o stdout4.log
|
||||
#SBATCH -e stderr4.log
|
||||
|
||||
#SBATCH --partition=ll_long
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
/usr/bin/modulecmd tcsh load mcstas
|
||||
bash compile_if_needed.sh
|
||||
|
||||
for y1 in $(seq -0.5 0.1 -0.1)
|
||||
do
|
||||
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
|
||||
do
|
||||
printf -v Y1 %.3f $y1
|
||||
printf -v Y2 %.3f $y2
|
||||
echo $Y1 $Y2
|
||||
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
|
||||
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
|
||||
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
|
||||
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
|
||||
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
|
||||
done
|
||||
done
|
||||
|
||||
|
||||
echo "Program finished with exit code $? at: `date`"
|
||||
|
||||
@@ -1,40 +0,0 @@
|
||||
#!/bin/bash
|
||||
#SBATCH -J mcEstia_5
|
||||
#SBATCH -N 2
|
||||
#SBATCH --ntasks-per-node=24
|
||||
#SBATCH --time=1-00:00:00
|
||||
#SBATCH --mail-type=fail
|
||||
#SBATCH --mail-user=artur.glavic@psi.ch
|
||||
|
||||
#SBATCH -o stdout5.log
|
||||
#SBATCH -e stderr5.log
|
||||
|
||||
#SBATCH --partition=ll_long
|
||||
|
||||
echo "Starting at `date`"
|
||||
echo "Running on hosts: $SLURM_NODELIST"
|
||||
echo "Running on $SLURM_NNODES nodes."
|
||||
echo "Running on $SLURM_NPROCS processors."
|
||||
echo "Current working directory is `pwd`"
|
||||
|
||||
/usr/bin/modulecmd tcsh load mcstas
|
||||
bash compile_if_needed.sh
|
||||
|
||||
for y1 in $(seq 0.1 0.1 0.5)
|
||||
do
|
||||
for y2 in $(seq -0.009 0.001 0.009) $(seq -0.09 0.01 -0.01) $(seq 0.01 0.01 0.09) $(seq -0.5 0.1 -0.1) $(seq 0.1 0.1 0.5)
|
||||
do
|
||||
printf -v Y1 %.3f $y1
|
||||
printf -v Y2 %.3f $y2
|
||||
echo $Y1 $Y2
|
||||
/afs/psi.ch/project/sinq/sl6-64/bin/mpirun -np $SLURM_NPROCS Estia_baseline.out \
|
||||
--dir=../results/selene1_geo_$Y1\_$Y2 --ncount=1e9 \
|
||||
omegaa=1.0 sample=4 sample_length=0.00005 sample_height=0.01 \
|
||||
lambda_start=2.5 lambda_end=15.0 enable_gravity=1 enable_chopper=1 \
|
||||
lambda_min=3.75 selene1_foot1y=$Y1 selene1_foot2y=$Y2
|
||||
done
|
||||
done
|
||||
|
||||
|
||||
echo "Program finished with exit code $? at: `date`"
|
||||
|
||||
Reference in New Issue
Block a user