Skip to content

Commit

Permalink
Merge pull request #1764 from McStasMcXtrace/single_crystal_order_fix
Browse files Browse the repository at this point in the history
Fix of order system in Single_crystal
  • Loading branch information
willend authored Nov 19, 2024
2 parents d2da55d + 49b90c5 commit 68aba44
Show file tree
Hide file tree
Showing 2 changed files with 130 additions and 45 deletions.
100 changes: 100 additions & 0 deletions mcstas-comps/examples/Tests_samples/Test_Powders/Test_Powders.instr
Original file line number Diff line number Diff line change
@@ -0,0 +1,100 @@
/*******************************************************************************
* McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Test PowderN output
*
* %Identification
* Written by: M. Bertelsen
* Date: November 19 2024
* Origin: ESS DMSC
* %INSTRUMENT_SITE: Tests_samples
*
* Test output of PowderN and Single_crystal on a spherical monitor / PSD.
*
* %Description
* A test instrument for Powder output from different sample components.
* Currently does not validate, but Single_crystal should be able to
* replicate PowderN results.
*
* %Example: lambda=2.5 directbeam=0 comp=1 SPLITS=1 Detector: Sph_mon_I=5.67059e+08
* %Example: lambda=2.5 directbeam=0 comp=1 SPLITS=1 Detector: psd_mon_I=1.46997e+06
* %Example: lambda=2.5 directbeam=0 comp=2 SPLITS=1 Detector: Sph_mon_I=5.67059e+08
* %Example: lambda=2.5 directbeam=0 comp=2 SPLITS=1 Detector: psd_mon_I=1.46997e+06
*
* %Parameters
* lambda: [Angs] Wavelength emitted from source, 1% wl-interval around.
* L1: [m] Source-sample distance
* directbeam: [1] Suppress direct beam or not
* reflections_powderN: [str] List of powder reflections for powderN
* reflections_Single_crystal: [str] List of powder reflections for Single_crystal
* SPLITS: [1] Number of SPLIT's before sample
* frac_t: [1] Fraction of stats assigned to unscattered, "direct beam"
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT Test_PowderN(int comp=1, lambda=2.5, L1=10, int directbeam=0, string reflections_powderN="Ge.laz", string reflections_Single_crystal="Ge.lau", int SPLITS=1, frac_t=0.1)

DECLARE %{
%}

/* The INITIALIZE section is executed when the simulation starts */
/* (C code). You may use them as component parameter values. */
INITIALIZE
%{
%}

/* Here comes the TRACE section, where the actual */
/* instrument is defined as a sequence of components. */
TRACE

REMOVABLE COMPONENT Origin = Progress_bar()
AT (0,0,0) ABSOLUTE

/* source with constant flux */
REMOVABLE COMPONENT Source = Source_gen(
radius = 0.11, dist = L1, focus_xw = 0.01, focus_yh = 0.1,
lambda0 = lambda, dlambda = lambda*0.01,
T1=229.6,I1=5.32e13,T2=1102, I2=4.37e12, T3=437.1,I3=3.84e13)
AT (0, 0, 0) RELATIVE Origin

/* TIP: monochromator cradle */
SPLIT SPLITS COMPONENT sample_cradle = Arm()
AT (0, 0, L1) RELATIVE PREVIOUS

COMPONENT Pow_PowderN = PowderN(
radius=0.005, yheight=0.1, reflections=reflections_powderN, p_transmit=frac_t)
WHEN (comp==1)
AT (0, 0, 0) RELATIVE sample_cradle
EXTEND %{
if(INSTRUMENT_GETPAR(directbeam) == 0) {
if (!SCATTERED) {
ABSORB;
}
}
%}

SPLIT SPLITS COMPONENT Pow_Single_crystal = Single_crystal(
radius=0.005, yheight=0.1, reflections=reflections_Single_crystal, p_transmit=frac_t,
powder=1, mosaic=5, order=1, barns=0)
WHEN (comp==2)
AT (0, 0, 0) RELATIVE sample_cradle
EXTEND %{
if(INSTRUMENT_GETPAR(directbeam) == 0) {
if (!SCATTERED) {
ABSORB;
}
}
%}

COMPONENT Sph_mon = PSD_monitor_4PI(nx=360,ny=180, radius=1, restore_neutron=1, filename="Sphere")
AT (0, 0, 0) RELATIVE PREVIOUS

COMPONENT peak_dir = Arm()
AT (0,0,0) RELATIVE sample_cradle
ROTATED (0, 100, 0) RELATIVE sample_cradle

COMPONENT psd_mon = PSD_monitor(xwidth=0.08, yheight=0.05, nx=100, ny=100, restore_neutron=1, filename="PSD")
AT (0,0,1) RELATIVE peak_dir

/* The END token marks the instrument definition end */
END
75 changes: 30 additions & 45 deletions mcstas-comps/samples/Single_crystal.comp
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
* Modified by: PW, June 2017: Doc updates
* Modified by: PW, Feb 2018: GPU edits
* Modified by: EF, June 2023: can now read CIF files
* Modified by: MB, November 2024: Order 0.5 optimization with large split numbers
* Modified by: MB, November 2024: Fixed issue with leaving crystal due to order
*
* Mosaic single crystal with multiple scattering vectors, optimised for speed
* with large crystals and many reflections.
Expand All @@ -37,14 +37,12 @@
* In order to dramatically improve the simulation efficiency, we recommend to
* use a SPLIT keyword on this component (or prior to it), as well as to disable
* the multiple scattering handling by setting order=1. This is especially powerful
* for large reflection lists such as with macromolecular proteins. With order=1
* the component still calculates extinction as the ray leaves the crystal, setting
* order=0.5 removes this step, further increasing the gains with large reflection
* lists and split.When an incoming particle is identical to the preceeding,
* reciprocal space initialisation is skipped, and a Monte Carlo choice is done on
* available reflections from the last repciprocal space calculation! To assist the
* user in choosing a "relevant" value of the SPLIT, a rolling average of the number
* of available reflections is calculated and presented in the component output.
* for large reflection lists such as with macromolecular proteins. When an incoming
* particle is identical to the preceeding, reciprocal space initialisation is
* skipped, and a Monte Carlo choice is done on available reflections from the last
* repciprocal space calculation! To assist the user in choosing a "relevant" value
* of the SPLIT, a rolling average of the number of available reflections is
* calculated and presented in the component output.
*
* <b>Mosacitiy modes:</b>
* The component features three independent ways of parametrising mosaicity:
Expand Down Expand Up @@ -183,7 +181,7 @@
* cy: [] c on y axis
* cz: [] c on z axis
* reflections: [string] File name containing structure factors of reflections (LAZ LAU CIF, FullProf, ShelX). Use empty ("") or NULL for incoherent scattering only
* order: [1] Limit multiple scattering up to given order (0: all, 1: first, 2: second, ...) Order 0.5 special case that skips coherent scattering when leaving crystal
* order: [1] Limit multiple scattering up to given order (0: all, 1: first, 2: second, ...)
*
* Optional input parameters
*
Expand Down Expand Up @@ -1084,7 +1082,6 @@ DECLARE
off_struct offdata;
struct hkl_data *hkl_list;
#ifndef OPENACC
struct tau_data tau_list_initial[MCSX_REFL_SLIST_SIZE];
struct tau_data tau_list[MCSX_REFL_SLIST_SIZE];
#endif
%}
Expand Down Expand Up @@ -1172,13 +1169,13 @@ INITIALIZE
exit(fprintf(stderr,"Single_crystal: %s: powder and PG modes can not be used together!\n"
"ERROR Please use EITHER powder or PG mode.\n", NAME_CURRENT_COMP));

if (powder && !(order==1 || order==0.5)) {
if (powder && !(order==1)) {
fprintf(stderr,"Single_crystal: %s: powder mode means implicit choice of no multiple scattering!\n"
"WARNING setting order=1\n", NAME_CURRENT_COMP);
order=1;
}

if (PG && !(order==1 || order==0.5)) {
if (PG && !(order==1)) {
fprintf(stderr,"Single_crystal: %s: PG mode means implicit choice of no multiple scattering!\n"
"WARNING setting order=1\n", NAME_CURRENT_COMP);
order=1;
Expand Down Expand Up @@ -1250,10 +1247,6 @@ TRACE
char type; /* type of last event: t=transmit,c=coherent or i=incoherent */
int itype; /* type of last event: t=1,c=2 or i=3 */

int force_transmit; /* Flag to handle cross-section weighting in case of finite order */

force_transmit=0;

#ifdef OPENACC
#ifdef USE_OFF
off_struct thread_offdata = offdata;
Expand Down Expand Up @@ -1298,12 +1291,17 @@ TRACE
itype = 0;

#ifndef OPENACC
T = tau_list_initial;
T = tau_list;
hkl_info.type = type;
#endif
do { /* Loop over multiple scattering events */
/* Angles for powder randomization */
double Alpha, Beta, Gamma;
double lab_vx, lab_vy, lab_vz;

lab_vx = vx;
lab_vy = vy;
lab_vz = vz;

if (hkl_info.shape == 0)
intersect = cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius, yheight);
Expand Down Expand Up @@ -1331,14 +1329,13 @@ TRACE

l_full = t2*v;

if (order==0.5 && force_transmit) {
// With order 0.5, exit before calculating coherent cross section on exit
// Exit due to truncated order, weight with relevant cross-sections to distance l_full
p*=exp(-abs_xlen*l_full);
intersect=0;
break;
}

if (order && event_counter >= order) {
/* Exit due to truncated order, weight with relevant cross-sections to distance l_full */
p*=exp(-abs_xlen*l_full);
intersect=0;
break;
}

/* (1). Compute incoming wave vector ki */
if (powder) { /* orientation of crystallite is random */
Alpha = randpm1()*PI*powder;
Expand All @@ -1351,14 +1348,12 @@ TRACE
PGrotate(&vx, &vy, &vz, Alpha, hkl_info.csx, hkl_info.csy, hkl_info.csz);
}



/* ------------------------------------------------------------------------- */
/* lattice curvature option: rotate neutron velocity */
/* WARNING: cannot be used together with the PG c-rotation! */
curv_xangle = 0;
curv_yangle = 0;

_vx = vx;
_vy = vy;
_vz = vz;
Expand Down Expand Up @@ -1400,12 +1395,8 @@ TRACE
/* in case we use 'SPLIT' then consecutive neutrons can be identical when entering here
and we may skip the hkl_search call. One tau_list is reserved for data for the initial
ray results so that it potentially can be reused later. */
if (event_counter==0) {
T = tau_list_initial;
} else {
T = tau_list;
}
if ((order==1 || order==0.5) && fabs(kix - hkl_info.kix) < deltak
if (order==1 && fabs(kix - hkl_info.kix) < deltak
&& fabs(kiy - hkl_info.kiy) < deltak
&& fabs(kiz - hkl_info.kiz) < deltak) {
hkl_info.nb_reuses++;
Expand All @@ -1414,7 +1405,6 @@ TRACE
coh_refl = hkl_info.coh_refl;
coh_xsect = hkl_info.coh_xsect;
tau_count = hkl_info.tau_count;
T = tau_list_initial;
} else {
#endif
/* Max possible tau for this ki with 5*sigma delta-d/d cutoff. */
Expand Down Expand Up @@ -1469,13 +1459,6 @@ TRACE
if(tot_xlen <= 0){
ABSORB; // Should we really absorb here? If "nothing" can happen we perhaps ought to "pass" instead?
}

if (force_transmit) {
/* Exit due to truncated order, weight with relevant cross-sections to distance l_full */
p*=exp(-abs_xlen*l_full);
intersect=0;
break;
}

/* (5). Transmission */
p_trans = exp(-tot_xlen*l_full);
Expand Down Expand Up @@ -1535,7 +1518,11 @@ TRACE
else
l = -log(1 - rand0max((1 - exp(-tot_xlen*l_full))))/tot_xlen;

/* Propagate to scattering point */
/* Propagate to scattering point in lab space */
vx = lab_vx;
vy = lab_vy;
vz = lab_vz;

PROP_DT(l/v);
event_counter++;

Expand Down Expand Up @@ -1638,8 +1625,6 @@ TRACE
if (PG) { /* orientation of crystallite is longer random */
PGderotate(&vx, &vy, &vz, Alpha, hkl_info.csx, hkl_info.csy, hkl_info.csz);
}
/* exit if multiple scattering order has been reached */
if (order && event_counter >= order) { force_transmit=1; }
/* Repeat loop for next scattering event. */
} while (intersect); /* end do (intersect) (multiple scattering loop) */
} /* if intersect */
Expand Down

0 comments on commit 68aba44

Please sign in to comment.