Work on ESS source description / comparison with MCPL inputs

mcstas-comps/examples/ESS/ESS_MCPL_generate_filtered/ESS_MCPL_generate_filtered.instr

@@ -0,0 +1,241 @@
+/*******************************************************************************
+* Instrument: ESS_MCPL_generate_filtered
+*
+* %I
+* Written by: Peter Willendrup <pkwi@fysik.dtu.dk>
+* Date: 2016-09-26
+* Origin: ESS
+* %INSTRUMENT_SITE: ESS
+*
+* Instrument meant to 'filter' legacy target division MCPL files.
+*
+* %D
+* Instrument meant to 'filter' legacy target division MCPL files (i.e. S2.mcpl.gz -> S2_filtered.mcpl.gz) for use with 
+* e.g. the ESS_butterfly_MCPL_test instrument and filter=1.
+*
+* To function, the instrument needs a local copy of the MCPL description you want to work with, e.g. from
+* https://public.esss.dk/users/willend/MCPL/
+*
+* The instrument assumes that binary MCPL datasets are available in . named [sector][beamline].mcpl.gz, i.e. W8.mcpl.gz.
+*
+* %P
+* sector:        [str]  Defines the 'sector' of your instrument position. Valid values are "N","S","E" and "W"
+* beamline:      [1]    Defines the 'beamline number' of your instrument position. Valid values are 1..10 or 1..11 depending on sector
+* Lmin:          [AA]   Minimum wavelength simulated
+* Lmax:          [AA]   Maximum wavelength simulated
+* c_performance: [1]    Cold brilliance scalar performance multiplicator c_performance > 0
+* t_performance: [1]    Thermal brilliance scalar performance multiplicator t_performance > 0
+* index:         [1]    Target index for source focusing. Defaults to illuminate the "cold collimated" brilliance monitor, thereby suppressing "dist"
+* dist:          [m]    Distance from origin to focusing rectangle; at (0,0,dist) - alternatively use target_index
+* cold:          [1]    Defines the statistical fraction of events emitted from the cold part of the moderator
+* Yheight:       [m]    Defines the moderator height. Valid values are 0.03 m and 0.06 m
+* delta:         [m]    Parameter that allows to scan "collimator" position
+* allmons:       [1]    Flag to include also debugging-monitor output
+*
+* %L
+* <reference/HTML link>
+* Benchmarking website available at <a href="http://ess_butterfly.mcstas.org">http://ess_butterfly.mcstas.org</a>
+* %E
+*******************************************************************************/
+DEFINE INSTRUMENT ESS_MCPL_generate_filtered(string sector="S",
+  int beamline=2,Lmin=0.2,Lmax=20,c_performance=1,t_performance=1,
+  int index=0,dist=2,cold=0.5,Yheight=0.03,delta=0,int allmons=0)
+
+
+DECLARE %{
+    double calcAlpha(double length, double radius) {
+    // calculate angle of arm after curved guide
+    return RAD2DEG * length/radius;
+  }
+
+  double calcX(double length, double radius) {
+    // calculate position and angle of arm after curved guide
+    double alpha = DEG2RAD * calcAlpha(length, radius);
+    return radius*(1.0-cos(alpha));
+  }
+
+  double calcZ(double length, double radius) {
+    // calculate position and angle of arm after curved guide
+    double alpha = DEG2RAD * calcAlpha(length, radius);
+    return radius*sin(alpha);
+  }
+
+  double XW, YH;
+  char options1[256],options2[256],options3[256],options4[256];
+  char srcdef[128];
+  double WidthC=0.072,WidthT=0.108;
+  double lambdamin, lambdamax;
+  double TCollmin;
+  double TCollmax;
+  #pragma acc declare create(TCollmin,TCollmax)
+  double EminTh=20, EmaxTh=100, EminC=0, EmaxC=20;
+  double weightmultiplier;
+  #pragma acc declare create(EminTh,EmaxTh,EminC,EmaxC,weightmultiplier)
+  /* 10 beamlines in sector N and E  - plus one location added for drawing */
+  double iBeamlinesN[] = { 30.0,  36.0,  42.0,  48.0,  54.0,  60.0,  66.0,  72.0,  78.0,  84.0,  90.0};
+  double iBeamlinesE[] = {-30.0, -36.0, -42.0, -48.0, -54.0, -60.0, -66.0, -72.0, -78.0, -84.0, -90.0};
+  /* 11 beamlines in sector S and W - plus one location added for drawing */
+  double iBeamlinesW[] = { 150.0,  144.7,  138.0,  132.7,  126.0,  120.7,  114.0,  108.7,  102.0,  96.7,  90.0,  84.0};
+  double iBeamlinesS[] = {-150.0, -144.7, -138.0, -132.7, -126.0, -120.7, -114.0, -108.7, -102.0, -96.7, -90.0, -84.0};
+  double* iBeamlines;
+  double ANGLE;
+  double DeltaX,DeltaZ;
+  char MCPLfile[128];
+  char FILTERfile[128];
+  int Divisions;
+%}
+
+USERVARS %{
+  int IsCold;
+  double SrcX;
+  double SrcY;
+  double SrcZ;
+  double E_min;
+  double E_max;
+  double Eneutron;
+  double T0;
+  double L0;
+  /* particle data in original MCNP / TCS units */
+  double xorig;
+  double yorig;
+  double zorig;
+  double vxorig;
+  double vyorig;
+  double vzorig;
+  double torig;
+%}
+
+INITIALIZE
+%{
+  lambdamin=Lmin;
+  lambdamax=Lmax;
+  XW=1.05*(WidthC+2*WidthT);
+  YH=1.05*Yheight;
+  sprintf(options1,"user1 bins=201 limits=[-%g,%g]",XW/2,XW/2);
+  sprintf(options4,"user1 bins=201 limits=[-%g,%g]",YH/2,YH/2);
+  sprintf(options2,"user1 bins=201 limits=[-%g,%g], user2 bins=201 limits=[-%g,%g]",XW/2,XW/2,YH/2,YH/2);
+  sprintf(options3,"user1 bins=201 limits=[-%g,%g], user2 bins=201 limits=[-%g,%g]",1.05*(WidthC/2),1.05*(WidthC/2),1.05*Yheight/2,1.05*Yheight/2);
+  sprintf(srcdef,"2015");
+  if (beamline==1) {
+    TCollmin=0;
+    TCollmax=0.058;
+  } else if (beamline==2) {
+    TCollmin=0;
+    TCollmax=0.06;
+  }
+  else {
+    TCollmin=0.011;
+    TCollmax=0.071;
+  }
+  #pragma acc update device(TCollmin,TCollmax)
+  if (strcasestr(sector,"N")) {
+    iBeamlines=iBeamlinesN;
+    DeltaX=-0.0585; DeltaZ=0.0925;
+  } else if (strcasestr(sector,"W")) {
+    iBeamlines=iBeamlinesW;
+    DeltaX=0.0585; DeltaZ=0.0925;
+  } else if (strcasestr(sector,"S")) {
+    iBeamlines=iBeamlinesS;
+    DeltaX=0.0585; DeltaZ=-0.0925;
+  } else if (strcasestr(sector,"E")) {
+    iBeamlines=iBeamlinesE;
+    DeltaX=-0.0585; DeltaZ=-0.0925;
+  }
+  ANGLE=iBeamlines[beamline-1]-90;
+  sprintf(MCPLfile,"%s%i.mcpl.gz",sector,beamline);
+  sprintf(FILTERfile,"%s%i_filtered.mcpl.gz",sector,beamline);
+  printf("MCPLfile is %s\n",MCPLfile);
+  printf("Filtered MCPLfile is %s\n",FILTERfile);
+
+  weightmultiplier=1.56e16/1e5; /* (ESS protons/s) / MCNP nps */
+  #pragma acc update device(EminTh,EmaxTh,EminC,EmaxC,weightmultiplier)
+%}
+
+TRACE
+
+COMPONENT Origin = Progress_bar()
+AT (0, 0, 0) ABSOLUTE
+
+/* 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_once(filename=MCPLfile,verbose=1)
+AT(0,0,-0.137) RELATIVE PREVIOUS
+EXTEND %{
+  SCATTER;
+  xorig=x; yorig=y; zorig=z;
+  vxorig=vx; vyorig=vy; vzorig=vz;
+  torig=t;
+%}
+
+COMPONENT Sphere1 = PSD_monitor_4PI(filename="nonrotated", radius=2.2,restore_neutron=1, nowritefile=!allmons)
+WHEN (allmons) AT (0,0,0) RELATIVE PREVIOUS
+
+COMPONENT Source = ESS_butterfly(sector=sector,beamline=beamline,Lmin=Lmin,Lmax=Lmax,c_performance=c_performance,t_performance=t_performance,dist=dist,target_index=index,cold_frac=cold, yheight=Yheight,
+				   focus_xw=0.12, focus_yh=0.12)
+  WHEN (0==1) AT (DeltaX,0,DeltaZ) ABSOLUTE
+  ROTATED (0, ANGLE, 0) ABSOLUTE
+
+COMPONENT Sphere0 = PSD_monitor_4PI(filename="rotated", radius=2.2,restore_neutron=1, nowritefile=!allmons)
+WHEN (allmons) AT (0,0,0) RELATIVE Source
+
+COMPONENT Focus_cut=Shape(xwidth=0.01,yheight=0.01)
+AT(0,0,2) RELATIVE Source
+EXTEND %{
+  ALLOW_BACKPROP;
+  PROP_Z0;
+  SCATTER;
+  if (fabs(x)>0.06 || fabs(y)>0.06) {
+    ABSORB;
+  } 
+%}
+
+COMPONENT ATSource=PSD_monitor(xwidth=0.31,yheight=0.31, nowritefile=!allmons)
+  WHEN (allmons) AT (0,0,0.08) RELATIVE Source
+
+COMPONENT BackTrace = Shape(xwidth=0.3,yheight=0.3)
+AT (0,0,0.08) RELATIVE Source
+EXTEND %{
+  /* Propagate back to a small rectangle in front of moderators */
+  SCATTER;
+  double myL = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz);
+  ALLOW_BACKPROP;
+  PROP_Z0;
+  SCATTER;
+  /* Remove neutrons that are not from around the moderators */
+  if (fabs(x)>0.12 || fabs(y)>0.03) {
+   ABSORB;
+  }
+  if ( myL < INSTRUMENT_GETPAR(Lmin) || myL > INSTRUMENT_GETPAR(Lmax) ) {
+    ABSORB;
+  }
+%}
+
+COMPONENT POSTATSource=PSD_monitor(xwidth=0.31,yheight=0.31)
+AT (0,0,0) RELATIVE BackTrace
+
+COMPONENT Restore = Arm()
+  AT (0,0,0) RELATIVE vin
+EXTEND %{
+  x=xorig; y=yorig; z=zorig;
+  vx=vxorig; vy=vyorig; vz=vzorig;
+  t=torig;
+%}
+
+COMPONENT vout = MCPL_output(filename=FILTERfile)
+  AT (0,0,0) RELATIVE vin
+
+
+FINALLY
+%{
+%}
+
+END

mcstas-comps/examples/ESS/ESS_MCPL_generate_filtered/README.md

@@ -0,0 +1,49 @@
+# The `ESS_MCPL_generate_filtered` Instrument
+
+*McStas: Instrument meant to 'filter' legacy target division MCPL files.*
+
+## Identification
+
+- **Site:** ESS
+- **Author:** Peter Willendrup <pkwi@fysik.dtu.dk>
+- **Origin:** ESS
+- **Date:** 2016-09-26
+
+## Description
+
+```text
+Instrument meant to 'filter' legacy target division MCPL files (i.e. S2.mcpl.gz -> S2_filtered.mcpl.gz) for use with
+e.g. the ESS_butterfly_MCPL_test instrument and filter=1.
+
+To function, the instrument needs a local copy of the MCPL description you want to work with, e.g. from
+https://public.esss.dk/users/willend/MCPL/
+
+The instrument assumes that binary MCPL datasets are available in . named [sector][beamline].mcpl.gz, i.e. W8.mcpl.gz.
+```
+
+## Input parameters
+
+Parameters in **boldface** are required; the others are optional.
+
+| Name | Unit | Description | Default |
+|------|------|-------------|---------|
+| sector | str | Defines the 'sector' of your instrument position. Valid values are "N","S","E" and "W" | "S" |
+| beamline | 1 | Defines the 'beamline number' of your instrument position. Valid values are 1..10 or 1..11 depending on sector | 2 |
+| Lmin | AA | Minimum wavelength simulated | 0.2 |
+| Lmax | AA | Maximum wavelength simulated | 20 |
+| c_performance | 1 | Cold brilliance scalar performance multiplicator c_performance > 0 | 1 |
+| t_performance | 1 | Thermal brilliance scalar performance multiplicator t_performance > 0 | 1 |
+| index | 1 | Target index for source focusing. Defaults to illuminate the "cold collimated" brilliance monitor, thereby suppressing "dist" | 0 |
+| dist | m | Distance from origin to focusing rectangle; at (0,0,dist) - alternatively use target_index | 2 |
+| cold | 1 | Defines the statistical fraction of events emitted from the cold part of the moderator | 0.5 |
+| Yheight | m | Defines the moderator height. Valid values are 0.03 m and 0.06 m | 0.03 |
+| delta | m | Parameter that allows to scan "collimator" position | 0 |
+| allmons | 1 | Flag to include also debugging-monitor output | 0 |
+
+## Links
+
+- [Source code](ESS_MCPL_generate_filtered.instr) for `ESS_MCPL_generate_filtered.instr`.
+- <reference/HTML link>
+- Benchmarking website available at <a href="http://ess_butterfly.mcstas.org">http://ess_butterfly.mcstas.org</a>
+
+---