[mcstas-users] ALLOW_BACKPROP, PROP_Z0 vs --gravitation

Richard Wagner wagnerrichard at ill.fr
Mon Aug 9 18:24:58 CEST 2021


Hi Peter,

find enclosed my example of an instrument where the back-propagation fails.

The source is an MCPL_input component and I added a small mcpl-input 
file, too.
The particles in the MCPL file all start at z=1.84m

Without gravitation back-propagation works, with gravitation it fails.

Any ideas, what could cause such a behavior?

Best,

Richard

PS: The PSI_source instrument you enclosed works for me (== 
back-propagation with gravity turned on)


On 26/07/2021 16:58, Peter Kjær Willendrup wrote:
> Hi Richard,
>
>
> I just came back from vacation and had a look at the issue you report. 
> ALLOW_BACKPROP ought to work irrespective of the gravity setting in fact.
>
> I tried putting your snippet of code into the PSI_source instrument 
> “behind” the source (see attached) - and I don’t seem to reproduce 
> your finding?
>
> Could the issue potentially be specific to the instrument you are 
> using it in? If you send me a copy (plus need inputfiles etc.) I will 
> have a closer look.
>
> (Feel free to do this outside the mailinglist if you prefer. :) )
>
>
> Best,
> Peter
>
>
>> On 18 Jul 2021, at 00.14, Richard Wagner <wagnerrichard at ill.fr 
>> <mailto:wagnerrichard at ill.fr>> wrote:
>>
>> Dear All,
>>
>> I have the following issue with a Shape Component in a simulation. It 
>> is used to filter neutrons that have not are not originating from the 
>> moderator window.
>>
>> COMPONENT Mod_Face = Shape(xwidth=.30,yheight=.24)
>> AT (0,0,0.16) RELATIVE ABSOLUTE
>> EXTEND %{
>>   /* Propagate back to a small rectangle in front of moderator */
>>   printf("(before) x=%g y=%g z=%g\n", x, y, z);
>>
>>   ALLOW_BACKPROP;
>>   PROP_Z0;
>>   SCATTER;
>>
>>   printf("(after) x=%g y=%g z=%g\n", x, y, z);
>>   printf("(if) %d\n", (fabs(x)> xwidth/2 || fabs(y)> yheight/2));
>>   /* Remove neutrons that are not from around the moderators */
>>   if (fabs(x)> xwidth/2 || fabs(y)> yheight/2)
>>   {
>>     printf("(absorbed) x=%g y=%g z=%g\n", x, y, z);
>>     ABSORB;
>>   }
>> %}
>>
>>
>> If I run the McStas simulation with the --gravitation flag the back 
>> propagation is not working.
>>
>> Log entries will exemplaryly read:
>> (if) 0
>> (before) x=0.131648 y=0.00568935 z=1.84
>> (after) x=0.131648 y=0.00568935 z=1.84
>>
>> And no neutrons will be filtered out
>>
>> With gravity set to 'off' I get:
>>
>> (before) x=0.110096 y=0.118489 z=1.84
>> (after) x=-0.118996 y=-0.00294855 z=0
>> (if) 0
>> (before) x=-0.0518426 y=-0.0952835 z=1.84
>> (after) x=0.519784 y=-0.514964 z=0
>> (if) 1
>> (absorbed) x=0.519784 y=-0.514964 z=0
>>
>> So it works as intended.
>>
>> Is this behavior as expected? Or differently is backprop not working 
>> with gravity turned on?
>> I use McStas 2.7 and need to make the simulation with gravity.
>>
>> Any ideas?
>>
>> Best,
>>
>> Richard
>>
>> --
>> *Richard Wagner*
>> Research Engineer
>> Nuclear and Particle Physics Group
>> Institut Laue-Langevin - ILL
>> 71, avenue des Martyrs
>> CS 20156
>> 38042 Grenoble Cedex 9
>> France
>>
>> www.ill.eu
>> _______________________________________________
>> mcstas-users mailing list
>> mcstas-users at mcstas.org <mailto:mcstas-users at mcstas.org>
>> https://mailman2.mcstas.org/mailman/listinfo/mcstas-users 
>> <https://mailman2.mcstas.org/mailman/listinfo/mcstas-users>
>
-- 
*Richard Wagner*
Research Engineer
Nuclear and Particle Physics Group
Institut Laue-Langevin - ILL
71, avenue des Martyrs
CS 20156
38042 Grenoble Cedex 9
France

www.ill.eu <www.ill.eu>


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20210809/ee90e190/attachment-0001.html>
-------------- next part --------------
/*******************************************************************************
*
* Instrument: Test_Sources
*
* %Identification
* Written by: Erik B Knudsen <erkn at fysik.dtu.dk>
* Date: Mar 2016
* Origin: DTU
* %INSTRUMENT_SITE: Tests_MCPL_etc
*
* A test instrument for MCPL_input
*
* %Description
*
* This is a unit test for the MCPL_input component.
*
* %Example: -n1e3 repeat=1 Detector: m1_I=2.42284e+11
*
* %Parameters
*
* repeat: [1]       Repeat the contents of the inputfile this many times. NB: When using MPI you implicitly repeat by #mpi processes
* E_smear: [1]      When repeating, make Gaussian MC choice on particle energy with spread E_smear * particle energy
* pos_smear: [m]    When repeating, make uniform MC choice on sphere of radius pos_spear around particle position
* dir_smear: [deg]  When repeating, make Gaussian MC choice in cone of opening dir_smear around particle direction
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT Test_MCPL_input(string fn="out_first_1000.mcpl.gz", repeat=1, E_smear=0.1, pos_smear=0.001,dir_smear=0.01)

DEPENDENCY "-DMCPLPATH=@MCCODE_LIB@/libs/mcpl/voutput.mcpl" 

DECLARE 
%{
  long long ncount_i;
  /**/
%}

INITIALIZE
%{
//#define QUOTE(name) #name
//#define STR(macro) QUOTE(macro)

//#ifndef MCPLPATH
//#define MCPLPATH=/usr/share/mcstas/2.3/MCPL/voutput.mcpl
//#endif

//#define MCPLFILE STR(MCPLPATH)

//printf("Using the input file: %s\n", MCPLFILE);

 ncount_i=0;
%}

TRACE

COMPONENT Origin = Progress_bar()
  AT (0, 0, 0) ABSOLUTE/* read neutrons from an mcpl file*/

/* read neutrons from an mcpl file*/
/*
COMPONENT vinROT2 = Arm()
AT(0,0,0) RELATIVE PREVIOUS
  ROTATED (0,90,0) RELATIVE PREVIOUS

COMPONENT vinROT1 = Arm()
AT(0,0,0) RELATIVE PREVIOUS
  ROTATED (-90,0,0) RELATIVE PREVIOUS
*/
COMPONENT vin = MCPL_input(filename=fn,verbose=1,repeat_count=repeat,E_smear=E_smear,pos_smear=pos_smear,dir_smear=dir_smear)
// AT(0,0,0) RELATIVE PREVIOUS 
AT(0,0+0.242,0) RELATIVE PREVIOUS 
/*AT(0,0+0.147,0) RELATIVE PREVIOUS */
ROTATED (0,90,0) RELATIVE PREVIOUS
EXTEND %{
  //y = y + 0.242; 
  // renormalization of intensity
  p*=1.56e16;
  p/=1e7;
  //printf("(0) x=%g y=%g z=%g\n", x, y, z);
%}

COMPONENT Mod_Face = Shape(xwidth=.30,yheight=.24)
  AT (0,0,0.16) RELATIVE ABSOLUTE
EXTEND %{
  /* Propagate back to a small rectangle in front of moderators */
  printf("(before) x=%g y=%g z=%g\n", x, y, z);
  ALLOW_BACKPROP;
  PROP_Z0;
  SCATTER;
  printf("(after) x=%g y=%g z=%g\n", x, y, z);  
  printf("(if) %d\n", (fabs(x)> xwidth/2 || fabs(y)> yheight/2));
  /* Remove neutrons that are not from around the moderators */
  if (fabs(x)> xwidth/2 || fabs(y)> yheight/2) 
  //if (fabs(x)> 0.15 || fabs(y)> 0.12) 
  {
    printf("(absorbed) x=%g y=%g z=%g\n", x, y, z);   
    ABSORB;   
  }
  /* Optional filters on lambda and time pulse
  /*double myL = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz);
  if (myL<mcipLmin || myL>mcipLmax) {
    ABSORB;
  }
  t+=rand01()*ESS_SOURCE_DURATION;
  if (t>3*ESS_SOURCE_DURATION) {
    ABSORB;
  }*/
%}

/*COMPONENT vin = MCPL_input(filename=MCPLFILE,verbose=1,repeat_count=repeat,E_smear=E_smear,pos_smear=pos_smear,dir_smear=dir_smear)
AT( 0,0,0) RELATIVE PREVIOUS */


COMPONENT Lambda_m1 = Monitor_nD(
  xwidth=2, yheight=2,
  options="lambda limits=[0.1 20], parallel, previous", bins=200,
  restore_neutron=1
) AT (0,0,2.7) ABSOLUTE

COMPONENT Source_divergence = Divergence_monitor(
    xwidth=4, 
    yheight=4, 
    nh=90, 
    nv=90, 
    filename="div_source", 
    maxdiv_h=90, 
    maxdiv_v=90, 
    restore_neutron=1)
AT (0, 0, 2.7) ABSOLUTE


/*
COMPONENT m3 = Monitor_nD(
  xwidth=0.2, yheight=0.2,
  options="auto t parallel, previous", bins=40
) AT (0,0,5) ABSOLUTE

COMPONENT m4 = Monitor_nD(
  xwidth=0.2, yheight=0.2,
  options="auto E, parallel, previous", bins=40
) AT (0,0,0) ABSOLUTE

COMPONENT m5 = Monitor_nD(
  xwidth=0.2, yheight=0.2, user1=2112,
  options="parallel, previous, neutron, user1, energy, x, y, z, vx, vy, vz, time, list all neutrons"
) AT (0,0,0) ABSOLUTE
*/

COMPONENT Guide_psd1 = PSD_monitor(
    filename="psd_1", restore_neutron=1, xwidth=.5, yheight=.5)
AT (0, 0, 2.7) ABSOLUTE
/*
COMPONENT psd_monitor_4pi = PSD_monitor_4PI(radius=0.9, nx=180, ny=180, 
filename="Events_1", restore_neutron=1)
AT (0, 0, 0) ABSOLUTE
*/

COMPONENT Guide_psd_3m = PSD_monitor(
    filename="psd_3m", restore_neutron=1, xwidth=2, yheight=2)
AT (0, 0, 3) ABSOLUTE

COMPONENT Guide_psd_10m = PSD_monitor(
    filename="psd_10m", restore_neutron=1, xwidth=4, yheight=4)
AT (0, 0, 10) ABSOLUTE

/*COMPONENT reflector_ellip = Elliptic_guide_gravity(
 l = 9.5, xwidth = .5,
 yheight = .5, linxw = .25,
 loutxw = .25, linyh = .25,
 loutyh = .25, dimensionsAt = "mid",
 m = 6) AT (0,0,.5) ABSOLUTE*/

COMPONENT m2_100m = Monitor_nD(
  xwidth=5, yheight=5,
  options="auto x, auto y", bins=100,
  restore_neutron=1
) AT (0,0,100) ABSOLUTE

COMPONENT Guide_psd_Det_100m = PSD_monitor(
    filename="psd_det", restore_neutron=1, xwidth=5, yheight=5)
AT (0, 0, 0.1) RELATIVE PREVIOUS

COMPONENT event_detector = Monitor_nD(
 xwidth = 10, yheight = 10,
 restore_neutron = 1, options = "list all auto, x y z vx vy vz t",
 filename = "events.dat")
AT (0,0,0) RELATIVE PREVIOUS

END
-------------- next part --------------
A non-text attachment was scrubbed...
Name: out_first_1000.mcpl.gz
Type: application/x-gzip
Size: 28037 bytes
Desc: not available
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20210809/ee90e190/attachment-0001.bin>


More information about the mcstas-users mailing list