McStas-3.3: No neutron can be detected after the Guide component with an MCPL_input source

[zTENGz]

I use SSW of MCNP6.2 to create a mono-directional (parallel to +z direction), single energy (14.87meV) neutron plane source which is uniformly distributed with x ranging from -5cm to 5cm and y ranging from -10 cm to 10cm. The number of neutron is 1e4. Then the wssa file is transformed to mcpl file via ssw2mcpl. A simple instrument file in created. It firstly reads mcpl source, watches source's wavelength, energy, and xy-distribution via the Monitor_nD component, secondly a 5m neutron guide is placed after the source via the Guide component, finally another group of Monitor_nD components are placed after the neutron guide. There are signals in the Monitor_nD after the mcpl source. But there are no signals in the Monitor_nD after the Guide component. I don't know what is wrong.

All related files can be downloaded from the link below https://drive.google.com/drive/folders/1lbrspix2gbRsKrr5m9G3smGKJBsw9KY-?usp=sharing

[willend]

Hello @zTENGz,

Thank you for the files, I have now had a closer look:

1) The main challenge when inputting an MCPL input file is "knowing where the particles are and where they are going". Looking at the output of mcpltool (similar shown when using the MCPL_input component) reveals that the neutrons are in fact at 500.1 cm (5.001 m) from the coordinate system of theMCPL_inputcomponent:

mcpltool 14.87meV_noTR_SSW.mcpl 
Opened MCPL file 14.87meV_noTR_SSW.mcpl:

  Basic info
    Format             : MCPL-3
    No. of particles   : 10000
(removed some lines for clarity...)

index     pdgcode   ekin[MeV]       x[cm]       y[cm]       z[cm]          ux          uy          uz    time[ms]      weight  userflags
    0        2112   1.487e-08      2.3867      7.2552       500.1           0           0           1  0.00059289           1 0x00000010
    1        2112   1.487e-08     -4.8566    -0.64316       500.1           0           0           1  0.00059289           1 0x00000010
    2        2112   1.487e-08       2.117     -9.6387       500.1           0           0           1  0.00059289           1 0x00000010
    3        2112   1.487e-08     -1.1762     -9.6445       500.1           0           0           1  0.00059289           1 0x00000010

This means that from the point of view of the Guide component which is displaced 0.1 m from your source_arm/source components, which in turn means that the neutrons "already passed" the Guide entry plane - and will hence not be considered. (In fact they are ABSORB'ed)

In the attached instrument file (rename .instr.txt -> .instr) I have displaced MCPL_input 5.001 m "back" wrt. the source_arm so that the neutrons rather become positioned at/near the source_arm z=0 plane. I have further added a couple of Monitor_nD's that measure the z-location, mostly for illustration purposes.

2) To know where your neutrons are heading, I have also added a Monitor_nD component measuring xdiv vs ydiv. This shows that all of your neutrons are effectively of 0 divergence, i.e. travelling perfectly parallel to z - which means that they will not reflect on the Guide walls, but simply pass unaltered through. This is seen in the attached dataset (plot and zip file attached) where the signal on PSD monitors m3 and m7 is identical.

Hope this initial help allows you tp progress further. :-)

Best Peter Willendrup

mcpl_source_NG.instr.txt mcpl_source_NG_20230710_141149.zip

[zTENGz]

Hi, Willen,

Thanks for your help.

This time I change the +z direction of mono-directional neutron plane source. It is modified by rotating 12 degrees counterclockwise from the +z direction via TR card in MCNP input (*TR1 0 0 0 12 90 -78 90 0 90 102 90 12). The settings of single energy (14.87meV) and uniformly distributed with x ranging from -5 cm to 5 cm and y ranging from -10 cm to 10 cm are not changed. A new SSW source is created. I add ROTATED ( 0, 0, -12) ABSOLUTE in the block of COMPONENT source_arm = Arm() and follow your approach to use Monitor_nD to measure x,y,z-location of source neutron. Then the offset values of x and z for MCPL_input are determined. I have two questions. 1. I'm not sure ROTATED ( 0, 0, -12) ABSOLUTE is right or not. 2. There are no signals in the Monitor_nD after the Guide component again.

Sincerely Yung-Hung Teng

All related files are updated and can be downloaded from the link below https://drive.google.com/drive/folders/1lbrspix2gbRsKrr5m9G3smGKJBsw9KY-?usp=sharing

[willend]

Hello again Yung-Hung,

There are a few important details about placing components in McStas that are good to understand well. I suggest looking at some of the basic McStas slides in e.g. the attached PDF taken from our "School" repo https://github.com/McStasMcXtrace/Schools/ (Intro slides from one of the last schools held at ESS).

ROTATED (alpha, beta, gamma) RELATIVE somewhere means rotated first alpha around x, then beta around y, finally gamma around z.

(We are luckily in the position of only one rotation, which makes things relatively simple - but if you need multiple and in another order, simply use multiple arms referring to each other. :) )

There is an illustration of placement of a Guide on page 26-30 which is good to look at, also have the remarks on page 48-51 in mind.

Since you are rotating -12 degrees "away" from original z, a first guess is that we need a rotation around y (vertical), I added an input variable Rguide=-12 deg as default anddGuide=0.1 m as default for rotating and placing the Guide.

I have implemented this in the attached, but as a rotation not of the first ABSOLUTE component (which will have very little effect) but rather on the optics. The result is not perfect, but at least there is transport through your Guide.

What typically helps me a lot in the debugging of such geometries is to use the mcdisplay tool "aka TRACE mode in mcgui" that can be used to visualise the instrument geometry, also PSD_monitor_4PI is useful, a perfect, spherical monitor.

(NB: I noticed that there are issues with mcdisplay if the instr file has multiple '.' in the filename, so the attached instr only has one suffix of '.txt' -> rename to '.instr' after downloading)

Best of luck simulating, I may return later with a further improved version of the instrument. Peter

mcpl_rotated_source_NG_PW.txt

1_McStas_intro_and_general_concepts.pdf

[willend]

Here is another version that I think is close to what you want to achieve.

I have added a rotation of Rsrc=-12 deg by default, dGuide=5.1 m by default and an Rguide=0 deg by default.

I have then applied the rotation -Rsrc around y to source_arm so that the beam becomes rectified along "typical" beam-direction z. I have further implemented an Arm that "propagates neutrons back" to the plane of the moderator so that they are no longer at the surface where they were recorded in MCNP but rather "close to the moderator".

The neutrons of course still have a perfect direction (0 divergence), so a natural next step would be to add a "per-particle" direction-cosine in MCNP (random number) - but if you want to experiment with the Guide reflecting them already now, we can e.g. do a scan of Rguide +/- 1 degrees

mcrun mcpl_rotated_source_NG_PW2.instr Rsrc=-12 Rguide=-1,1 -N21 --autoplot on the command-line. In mcgui use a run dialogue like this picture: (I have marked scan range, number of steps and how to get a plot after the scan series)

Once such a simulation is done --autoplot will bring up this plot:

If you then ctrl+click on the image as indicated ^ , you will see how the transported beam changes with rotation Rguide: (illustrating the PSD "m7" after the Guide as function of the scan)

Best and good luck simulating on. Peter

mcpl_rotated_source_NG_PW2.txt

[zTENGz]

Hello again Yung-Hung,

There are a few important details about placing components in McStas that are good to understand well. I suggest looking at some of the basic McStas slides in e.g. the attached PDF taken from our "School" repo https://github.com/McStasMcXtrace/Schools/ (Intro slides from one of the last schools held at ESS).

ROTATED (alpha, beta, gamma) RELATIVE somewhere means rotated first alpha around x, then beta around y, finally gamma around z.

(We are luckily in the position of only one rotation, which makes things relatively simple - but if you need multiple and in another order, simply use multiple arms referring to each other. :) )

There is an illustration of placement of a Guide on page 26-30 which is good to look at, also have the remarks on page 48-51 in mind.

Since you are rotating -12 degrees "away" from original z, a first guess is that we need a rotation around y (vertical), I added an input variable Rguide=-12 deg as default anddGuide=0.1 m as default for rotating and placing the Guide.

I have implemented this in the attached, but as a rotation not of the first ABSOLUTE component (which will have very little effect) but rather on the optics. The result is not perfect, but at least there is transport through your Guide.

What typically helps me a lot in the debugging of such geometries is to use the mcdisplay tool "aka TRACE mode in mcgui" that can be used to visualise the instrument geometry, also PSD_monitor_4PI is useful, a perfect, spherical monitor.

(NB: I noticed that there are issues with mcdisplay if the instr file has multiple '.' in the filename, so the attached instr only has one suffix of '.txt' -> rename to '.instr' after downloading)

Best of luck simulating, I may return later with a further improved version of the instrument. Peter

mcpl_rotated_source_NG_PW.txt

1_McStas_intro_and_general_concepts.pdf

Hi, Peter

Thanks for your help.

About the PSD_monitor_4PI, the webpage of The PSD_monitor_4PI Component said it is "An (n times m) pixel spherical PSD monitor using a cylindrical projection". Is the idea of cylindrical projection similar to Mercator projection method? I found a figure about Mercator projection in wiki as below. I can't understand the relationship of the radius R of the sphere and the n times m pixel matrix. Does it mean that 2πR divided by n in longitude and 2πR divided by m in latitude?

About the "aka TRACE mode in mcgui", does it mean the Trace function? The screenshot of html file is below. This figure help me to realize that the definition of +x direction and -x direction in McStas are different from the definition in MCNP.

Regards Yung-Hung

[willend]

Hi, Peter

Thanks for your help.

About the PSD_monitor_4PI, the webpage of The PSD_monitor_4PI Component said it is "An (n times m) pixel spherical PSD monitor using a cylindrical projection". Is the idea of cylindrical projection similar to Mercator projection method? I found a figure about Mercator projection in wiki as below.

Yes, I believe that is right.

I can't understand the relationship of the radius R of the sphere and the n times m pixel matrix. Does it mean that 2πR divided by n in longitude and 2πR divided by m in latitude?

As I remember this is 2πR divided by n in longitude ∈ [-π π] and πR divided by m in latitude ∈ [-π/2 π/2]. Of course a sphere of any R will have the same angular sub-division in n and m, R is only there since the neutrons are propagated to the surface before being assigned to the monitor bins.

About the "aka TRACE mode in mcgui", does it mean the Trace function?

Exactly. There are a multitude of variants of "mcdisplay" available that stronger / weaker for different types of visualisation, see https://github.com/McStasMcXtrace/McCode/wiki/mcdisplay-variants---table-overview . They way they work is that instrument geometry and particle ray geometry is communicated via stdout from the simulation binary during the course of simulation (ala what you get from ./my_instrument.out --trace --no-output-files)

This figure help me to realize that the definition of +x direction and -x direction in McStas are different from the definition in MCNP.

Great, this is much better than "flying blind". (And when two frameworks are in use that can in principle both make "any needed orientation" this easily becomes a little complex. :-) ) For a McStas-provided example of a non-trivial instrument that uses MCPL input from MCNP(X), see your $MCSTAS/examples/ESS_butterfly_MCPL_test.instr - to simulate the beam out of an ESS beamport it will need the corresponding input file from https://public.esss.dk/users/willend/MCPL/

[zTENGz]

Here is another version that I think is close to what you want to achieve.

I have added a rotation of Rsrc=-12 deg by default, dGuide=5.1 m by default and an Rguide=0 deg by default.

I have then applied the rotation -Rsrc around y to source_arm so that the beam becomes rectified along "typical" beam-direction z. I have further implemented an Arm that "propagates neutrons back" to the plane of the moderator so that they are no longer at the surface where they were recorded in MCNP but rather "close to the moderator".

The neutrons of course still have a perfect direction (0 divergence), so a natural next step would be to add a "per-particle" direction-cosine in MCNP (random number) - but if you want to experiment with the Guide reflecting them already now, we can e.g. do a scan of Rguide +/- 1 degrees

mcrun mcpl_rotated_source_NG_PW2.instr Rsrc=-12 Rguide=-1,1 -N21 --autoplot on the command-line. In mcgui use a run dialogue like this picture: (I have marked scan range, number of steps and how to get a plot after the scan series)

Once such a simulation is done --autoplot will bring up this plot:

If you then ctrl+click on the image as indicated ^ , you will see how the transported beam changes with rotation Rguide: (illustrating the PSD "m7" after the Guide as function of the scan)

Best and good luck simulating on. Peter

mcpl_rotated_source_NG_PW2.txt

Hi, Peter

Thanks for your help.

This time I use SSW to write a non-mono-directional neutron plane source. The surface plane is set by rotating 12 degrees counterclockwise from the +z direction via TR card in MCNP input (*TR1 0 0 0 12 90 -78 90 0 90 102 90 12) and it is at guide of about 1m from the PE moderator. In MCNP calculation, the number of neutrons above 100 meV and the number of neutrons below 100 meV is about 1:1. Then I use mcpl_rotated_source_NG_PW2.txt to read the mcpl source file, set Rsrc=-12 and do the Rguide sweep between -10 and 0. The search result is shown below. The maximum intensity is occurred at about -6.5 degree. Then I do another calculation that set Rsrc=-12 and Rguide=-6.5. The result is shown below. It seems most neutrons from the plane source do not enter the guide. How can I deal with the problem? In SSW card I write SSW 225(48) that a track will be recorded if it crosses the surface in the correct direction and is either entering a cell in the list whose entry is positive. The surface cards of cell 48 is as 201 -202 203 -204 225 -226.
The surface definition is as below. 201 pz -10.0
202 pz 10.0
203 2 px -5.0
204 2 px 5.0
In MCNP calculation without super-mirror effect the neutron flux below 40 meV at about 1m of guide is about 4.6e5 n/cm2/s, and the value at about 6m of guide is about 2.0e5 n/cm2/s. Both values are tallyed in the range of 0 degree to 1 degree from the direction of guide.

The second question is about the neutron intensity normalization. I can get the number of neutrons below 100 meV from MCNP calculation with FM card 2.4968e15 proton/s corresponding to 0.4mA proton current at SSW surface plane. Can I use this value to adjust Mcstas intensity results at the position of plane source and at the position of guide exit?

Sincerely Yung-Hung

All related files are updated and can be downloaded from the link below https://drive.google.com/drive/folders/1lbrspix2gbRsKrr5m9G3smGKJBsw9KY-?usp=sharing

[willend]

Hi @zTENGz,

Just a brief reply for now as I have just started my summer vacation. If time allows I will be back with a more in-depth analysis once I am back around August 7th. :-)

I have been poking around a bit more in your instrument file and experimenting - as it looks to me an x-translation of the guide wrt. the beam is somehow missing. In the attached instrument I have implemented such a translation and now the Guide seems to transport more reasonably. For guidance I used the 3 PSD's of varying size placed at the guide entry and visually inspected that the "spot" is aligned with the guide dimensions.

That being said, to get everything right I think you need to mentally "project" the MCNP geometry and the McStas geometry "onto" each other: Moderator, guide opening and "neutron cloud" to ensure you get the tranformations right. (Think of the neutron cloud likely being on the intersection between the chosen SSW surface and a dXtran sphere.)

As you know, in general neutron guides only transport cold and thermal neutrons effectively, meaning of wavelength ~1Å and colder, I.e. ~81 meV or lower in energy. Also, the neutrons need to hit the guide entry and further, the critical angle of super-mirrors are of the order of ~1 degree. Having only had a brief look in at the data of the MCPL file, these effects all together will dilute your statistics quite a bit. For that reason I have also implemented a "repeat" mechanism: Via Gaussian filters in direction, energy and position, the neutrons may be "re-used" multiple times. (Total intensity/weight is conserved - not multiplied.) See MCPL_input docs for details.

I attach screenshots of the output of 1 run of the MCPL file on 1 core plus a 16 x repeat (done via mpi on 8 cores, i.e. 2 passes of the file pr core)

Please find attached an instrument file that seems to transport your neutrons reasonably - try experimenting / monitoring further yourself to get it even better.

Best of luck simulating until I get back, Peter

mcpl_rotated_source_NG_PW3.txt

[willend]

Dear Yung-Hung,

Just remembered your second question also:

The second question is about the neutron intensity normalization. I can get the number of neutrons below 100 meV from MCNP calculation with FM card 2.4968e15 proton/s corresponding to 0.4mA proton current at SSW surface plane. Can I use this value to adjust Mcstas intensity results at the position of plane source and at the position of guide exit?

My recommendation here would be to take the same approach as used the earlier mentioned ESS_butterfly_MCPL_test.instr, namely to think in terms of "MCNP source particles":

COMPONENT vin = MCPL_input(.....)
AT(0,0,0) RELATIVE PREVIOUS
EXTEND %{
  SCATTER;
  p*=1.56e16;
  p/=1e5;
%}

1) Normalise by MCNP file statistics, i.e. what NPS was used in underlying MCNP simulation ( in the example this is p/=1e5; ) 2) Multiply with number of protons / s at planned accelerator current (in the example p*=1.56e16;)

Good luck, Peter

[zTENGz]

Dear Peter,

Have a nice summer vacation.

Best wishes, Yung-Hung