/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Monitor_nD
*
* %Identification
* Written by: Emmanuel Farhi
* Date: 14th Feb 2000.
* Origin: ILL
* Release: McStas 1.6
* Version: $Revision$
* Modified by: EF, 29th Feb 2000 : added more options, monitor shape, theta, phi
* Modified by: EF, 01st Feb 2001 : PreMonitor for correlation studies (0.13.6)
* Modified by: EF, 5th Apr 2001 : use global functions (0.14) compile faster
* Modified by: EF, 23th Jul 2001 : log of signal, init arrays to 0, box (0.15)
* Modified by: EF, 04th Sep 2001 : log/abs of variables (0.16)
* Modified by: EF, 24th Oct 2001 : capture flux [p*lambda/1.7985] (0.16.3)
* Modified by: EF, 27th Aug 2002 : monitor a variable in place of I (0.16.5)
* Modified by: EF, 25th Oct 2002 : banana, and auto for each variable (0.16.5)
*
* This component is a general Monitor that can output 0/1/2D signals
* (Intensity or signal vs. [something] and vs. [something] ...)
*
* %Description
* This component is a general Monitor that can output 0/1/2D signals
* It can produce many 1D signals (one for any variable specified in
* option list), or a single 2D output (two variables correlation).
* Also, an additional 'list' of neutron events can be produced.
* By default, monitor is square (in x/y plane). A disk shape is also possible
* The 'cylinder' and 'banana' option will change that for a banana shape
* The 'sphere' option simulates spherical detector. The 'box' is a box.
* The cylinder, sphere and banana should be centered on the scattering point.
* The monitored flux may be per monitor unit area, and weighted by
* a lambda/lambda(2200m/s) factor to obtain standard integrated capture flux.
* In normal configuration, the Monitor_nD measures the current parameters
* of the neutron that is beeing detected. But a PreMonitor_nD component can
* be used in order to study correlations between a neutron being detected in
* a Monitor_nD place, and given parameters that are monitored elsewhere
* (at PreMonitor_nD).
* The monitor can also act as a 3He gas detector, taking into account the
* detection efficiency.
*
* The 'bins' and 'limits' modifiers are to be used after each variable,
* and 'auto','log' and 'abs' come before it. (eg: auto abs log hdiv bins=10
* limits=[-5 5]) When placed after all variables, these two latter modifiers
* apply to the signal (e.g. intensity). Unknown keywords are ignored.
* If no limits are specified for a given observable, reasonable defaults will be
* applied. Note that these implicit limits are even applied in list mode.
*
* Implicit limits for typical variables:
* (consult monitor_nd-lib.c if you don't find your variable here)
* x, y, z: Derived from detection-object geometry
* k: [0 10] Angs-1
* v: [0 1e6] m/s
* t: [0 1] s
* p: [0 FLT_MAX] in intensity-units
* vx, vy: [-1000 1000] m/s
* vz: [0 10000] m/s
* kx, ky: [-1 1] Angs-1
* kz: [-10 10] Angs-1
* energy, omega: [0 100] meV
* lambda,wavelength: [0 100] Angs
* sx, sy, sz: [-1 1] in polarisation-units
* angle: [-50 50] deg
* divergence, vdiv, hdiv, xdiv, ydiv: [-5 5] deg
* longitude, lattitude: [-180 180] deg
* neutron: [0 simulaton_ncount]
* id, pixel id: [0 FLT_MAX]
* uservars u1,u2,u3: [-1e10 1e10]
*
* In the case of multiple components at the same position, the 'parallel'
* keyword must be used in each instance instead of defining a GROUP.
*
* Possible options are
* Variables to record:
* kx ky kz k wavevector [Angs-1] Wavevector on x,y,z and norm
* vx vy vz v [m/s] Velocity on x,y,z and norm
* x y z radius [m] Distance, Position and norm
* xy, yz, xz [m] Radial position in xy, yz and xz plane
* kxy kyz kxz [Angs-1] Radial wavevector in xy, yz and xz plane
* vxy vyz vxz [m/s] Radial velocity in xy, yz and xz plane
* t time [s] Time of Flight
* energy omega [meV] energy of neutron
* lambda wavelength [Angs] wavelength of neutron
* sx sy sz [1] Spin
* vdiv ydiv dy [deg] vertical divergence (y)
* hdiv divergence xdiv [deg] horizontal divergence (x)
* angle [deg] divergence from direction
* theta longitude [deg] longitude (x/z) for sphere and cylinder
* phi lattitude [deg] lattitude (y/z) for sphere and cylinder
*
* user user1 will monitor the [Mon_Name]_Vars.UserVariable{1|2|3}
* user2 user3 to be assigned in an other component (see below)
*
* p intensity flux [n/s or n/cm^2/s]
* ncounts n neutron [1] neutron ID, i.e current event index
* pixel id [1] pixelID in histogram made of preceeding vars, e.g. 'theta y'. To set an offset PixelID use the 'min=value' keyword. Sets event mode.
*
* Other options keywords are:
* abs Will monitor the abs of the following variable or of the signal (if used after all variables)
* auto Automatically set detector limits for one/all
* all {limits|bins|auto} To set all limits or bins values or auto mode
* binary {float|double} with 'source' option, saves in compact files
* bins=[bins=20] Number of bins in the detector along dimension
* borders To also count off-limits neutrons (X < min or X > max)
* capture weight by lambda/lambda(2200m/s) capture flux
* file=string Detector image file name. default is component name, plus date and variable extension.
* incoming Monitor incoming beam in non flat det
* limits=[min max] Lower/Upper limits for axes (see up for the variable unit)
* list=[counts=1000] or all For a long file of neutron characteristics with [counts] or all events
* log Will monitor the log of the following variable or of the signal (if used after all variables)
* min=[min_value] Same as limits, but only sets the min or max
* max=[max_value]
* multiple Create multiple independant 1D monitors files
* no or not Revert next option
* outgoing Monitor outgoing beam (default)
* parallel Use this option when the next component is at the same position (parallel components)
* per cm2 Intensity will be per cm^2 (detector area). Displays beam section.
* per steradian Displays beam solid angle in steradian
* premonitor Will monitor neutron parameters stored previously with PreMonitor_nD.
* signal=[var] Will monitor [var] instead of usual intensity
* slit or absorb Absorb neutrons that are out detector
* source The monitor will save neutron states
* unactivate To unactivate detector (0D detector)
* verbose To display additional informations
* 3He_pressure=[3 in bars] The 3He gas pressure in detector. 3He_pressure=0 is perfect detector (default)
*
* Detector shape options (specified as xwidth,yheight,zdepth or x/y/z/min/max)
* box Box of size xwidth, yheight, zdepth.
* cylinder To get a cylindrical monitor (diameter is xwidth or set radius, height is yheight).
* banana Same as cylinder, without top/bottom, on restricted angular area; use theta variable with limits to define arc. (diameter is xwidth or set radius, height is yheight).
* disk Disk flat xy monitor. diameter is xwidth.
* sphere To get a spherical monitor (e.g. a 4PI) (diameter is xwidth or set radius).
* square Square flat xy monitor (xwidth, yheight).
* previous The monitor uses PREVIOUS component as detector surface. Or use 'geometry' parameter to specify any PLY/OFF geometry file.
*
* EXAMPLES:
*
* - MyMon = Monitor_nD(xwidth = 0.1, yheight = 0.1, zdepth = 0,
* options = "intensity per cm2 angle,limits=[-5 5] bins=10,with
* borders, file = mon1");
* will monitor neutron angle from [z] axis, between -5
* and 5 degrees, in 10 bins, into "mon1.A" output 1D file
*
*
- options = "sphere theta phi outgoing"
* for a sphere PSD detector (out beam) and saves into file "MyMon_[Date_ID].th_ph"
*
*
- options = "banana, theta limits=[10,130], bins=120, y"
* a theta/height banana detector
*
*
- options = "angle radius all auto"
* is a 2D monitor with automatic limits
*
*
- options = "list=1000 kx ky kz energy"
* records 1000 neutron event in a file
*
*
- options = "multiple kx ky kz, auto abs log t, and list all neutrons"
* makes 4 output 1D files and produces a complete list for all neutrons
* and monitor log(abs(tof)) within automatic limits (for t)
*
*
- options = "theta y, sphere, pixel min=100"
* a 4pi detector which outputs an event list with pixelID from the actual
* detector surface, starting from index 100.
*
*
* To dynamically define a number of bins, or limits:
* Use in DECLARE: char op[256];
* Use in INITIALIZE: sprintf(op, "lambda limits=[%g %g], bins=%i", lmin, lmax, lbin);
* Use in TRACE: Monitor_nD(... options=op ...)
*
* How to monitor any instrument/component variable into a Monitor_nD
* Suppose you want to monitor a variable 'age' which you assign somwhere in
* the instrument:
* COMPONENT MyMonitor = Monitor_nD(
* xwidth = 0.1, yheight = 0.1,
* user1="age", username1="Age of the Captain [years]",
* options="user1, auto")
* AT ...
*
* See also the example in PreMonitor_nD to
* monitor neutron parameters cross-correlations.
*
* %BUGS
* The 'auto' option for guessing optimal variable bounds should NOT be used with MPI
* as each process may use different limits.
*
* %Parameters
* INPUT PARAMETERS:
*
* xwidth: [m] Width of detector.
* yheight: [m] Height of detector.
* zdepth: [m] Thickness of detector (z).
* radius: [m] Radius of sphere/banana shape monitor
* options: [str] String that specifies the configuration of the monitor. The general syntax is "[x] options..." (see Descr.).
*
* Optional input parameters (override xwidth yheight zdepth):
* xmin: [m] Lower x bound of opening
* xmax: [m] Upper x bound of opening
* ymin: [m] Lower y bound of opening
* ymax: [m] Upper y bound of opening
* zmin: [m] Lower z bound of opening
* zmax: [m] Upper z bound of opening
* filename: [str] Output file name (overrides file=XX option).
* bins: [1] Number of bins to force for all variables. Use 'bins' keyword in 'options' for heterogeneous bins
* min: [u] Minimum range value to force for all variables. Use 'min' or 'limits' keyword in 'options' for other limits
* max: [u] Maximum range value to force for all variables. Use 'max' or 'limits' keyword in 'options' for other limits
* user1: [str] Variable name of USERVAR to be monitored by user1.
* user2: [str] Variable name of USERVAR to be monitored by user2.
* user3: [str] Variable name of USERVAR to be monitored by user3.
* username1: [str] Name assigned to User1
* username2: [str] Name assigned to User2
* username3: [str] Name assigned to User3
* restore_neutron: [0|1] If set, the monitor does not influence the neutron state. Equivalent to setting the 'parallel' option.
* geometry: [str] Name of an OFF file to specify a complex geometry detector
* nowritefile: [1] If set, monitor will skip writing to disk
*
* OUTPUT PARAMETERS:
*
* DEFS: [struct] structure containing Monitor_nD Defines
* Vars: [struct] structure containing Monitor_nD variables
*
* %Link
* PreMonitor_nD
*
* %End
******************************************************************************/
DEFINE COMPONENT Monitor_nD
DEFINITION PARAMETERS ()
SETTING PARAMETERS (
string user1="", string user2="", string user3="",
xwidth=0, yheight=0, zdepth=0,
xmin=0, xmax=0, ymin=0, ymax=0, zmin=0, zmax=0,
bins=0, min=-1e40, max=1e40, restore_neutron=0, radius=0,
string options="NULL", string filename="NULL",string geometry="NULL", int nowritefile=0,
string username1="NULL", string username2="NULL", string username3="NULL"
)
/* these are protected C variables */
OUTPUT PARAMETERS ()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "monitor_nd-lib"
%include "read_table-lib"
%include "interoff-lib"
%}
DECLARE
%{
MonitornD_Defines_type DEFS;
MonitornD_Variables_type Vars;
MCDETECTOR detector;
off_struct offdata;
%}
INITIALIZE
%{
char tmp[CHAR_BUF_LENGTH];
strcpy(Vars.compcurname, NAME_CURRENT_COMP);
if (options != NULL)
strncpy(Vars.option, options, CHAR_BUF_LENGTH);
else {
strcpy(Vars.option, "x y");
printf("Monitor_nD: %s has no option specified. Setting to PSD ('x y') monitor.\n", NAME_CURRENT_COMP);
}
Vars.compcurpos = POS_A_CURRENT_COMP;
if (strstr(Vars.option, "source"))
strcat(Vars.option, " list, x y z vx vy vz t sx sy sz ");
if (bins) { sprintf(tmp, " all bins=%ld ", (long)bins); strcat(Vars.option, tmp); }
if (min > -FLT_MAX && max < FLT_MAX) { sprintf(tmp, " all limits=[%g %g]", min, max); strcat(Vars.option, tmp); }
else if (min > -FLT_MAX) { sprintf(tmp, " all min=%g", min); strcat(Vars.option, tmp); }
else if (max < FLT_MAX) { sprintf(tmp, " all max=%g", max); strcat(Vars.option, tmp); }
/* transfer, "zero", and check username- and user variable strings to Vars struct*/
strncpy(Vars.UserName1,
username1 && strlen(username1) && strcmp(username1, "0") && strcmp(username1, "NULL") ?
username1 : "", 128);
strncpy(Vars.UserName2,
username2 && strlen(username2) && strcmp(username2, "0") && strcmp(username2, "NULL") ?
username2 : "", 128);
strncpy(Vars.UserName3,
username3 && strlen(username3) && strcmp(username3, "0") && strcmp(username3, "NULL") ?
username3 : "", 128);
if(user1 && strlen(user1) && strcmp(user1, "0") && strcmp(user1, "NULL")){
strncpy(Vars.UserVariable1,user1,128);
int fail;_class_particle testparticle;
particle_getvar(&testparticle,Vars.UserVariable1,&fail);
if(fail){
fprintf(stderr,"Warning (%s): user1=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user1);
}
}
if(user2 && strlen(user2) && strcmp(user2, "0") && strcmp(user2, "NULL")){
strncpy(Vars.UserVariable2,user2,128);
int fail;_class_particle testparticle;
particle_getvar(&testparticle,Vars.UserVariable2,&fail);
if(fail){
fprintf(stderr,"Warning (%s): user2=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user2);
}
}
if(user3 && strlen(user3) && strcmp(user3, "0") && strcmp(user3, "NULL")){
strncpy(Vars.UserVariable3,user3,128);
int fail;_class_particle testparticle;
particle_getvar(&testparticle,Vars.UserVariable3,&fail);
if(fail){
fprintf(stderr,"Warning (%s): user3=%s is unknown. The signal will not be resolved - this is likely not what you intended.\n",NAME_CURRENT_COMP,user3);
}
}
/*sanitize parameters set for curved shapes*/
if(strstr(Vars.option,"cylinder") || strstr(Vars.option,"banana") || strstr(Vars.option,"sphere")){
/*this _is_ an explicit curved shape. Should have a radius. Inherit from xwidth or zdepth (diameters), x has precedence.*/
if (!radius){
if(xwidth){
radius=xwidth/2.0;
}else{
radius=zdepth/2.0;
}
}else{
xwidth=2*radius;
}
if(!yheight){
/*if not set - use the diameter as height for the curved object. This will likely only happen for spheres*/
yheight=2*radius;
}
}else if (radius) {
/*radius is set - this must be a curved shape. Infer shape from yheight, and set remaining values
(xwidth etc. They are used inside monitor_nd-lib.*/
xwidth = zdepth = 2*radius;
if (yheight){
/*a height is given (and no shape explitly set - assume cylinder*/
strcat(Vars.option, " banana");
}else {
strcat(Vars.option, " sphere");
yheight=2*radius;
}
}
int offflag=0;
if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) {
#ifndef USE_OFF
fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n");
exit(-1);
#else
if (!off_init( geometry, xwidth, yheight, zdepth, 1, &offdata )) {
printf("Monitor_nD: %s could not initiate the OFF geometry %s. \n"
" Defaulting to normal Monitor dimensions.\n",
NAME_CURRENT_COMP, geometry);
strcpy(geometry, "");
} else {
offflag=1;
}
#endif
}
if (!radius && !xwidth && !yheight && !zdepth && !xmin && !xmax && !ymin && !ymax &&
!strstr(Vars.option, "previous") && (!geometry || !strlen(geometry)))
exit(printf("Monitor_nD: %s has no dimension specified. Aborting (radius, xwidth, yheight, zdepth, previous, geometry).\n", NAME_CURRENT_COMP));
Monitor_nD_Init(&DEFS, &Vars, xwidth, yheight, zdepth, xmin,xmax,ymin,ymax,zmin,zmax,offflag);
if (Vars.Flag_OFF) {
offdata.mantidflag=Vars.Flag_mantid;
offdata.mantidoffset=Vars.Coord_Min[Vars.Coord_Number-1];
}
if (filename && strlen(filename) && strcmp(filename,"NULL") && strcmp(filename,"0"))
strncpy(Vars.Mon_File, filename, 128);
/* check if user given filename with ext will be used more than once */
if ( ((Vars.Flag_Multiple && Vars.Coord_Number > 1) || Vars.Flag_List) && strchr(Vars.Mon_File,'.') )
{ char *XY; XY = strrchr(Vars.Mon_File,'.'); *XY='_'; }
if (restore_neutron) Vars.Flag_parallel=1;
detector.m = 0;
#ifdef USE_MPI
MPI_MASTER(
if (strstr(Vars.option, "auto") && mpi_node_count > 1)
printf("Monitor_nD: %s is using automatic limits option 'auto' together with MPI.\n"
"WARNING this may create incorrect distributions (but integrated flux will be right).\n", NAME_CURRENT_COMP);
);
#else
#ifdef OPENACC
if (strstr(Vars.option, "auto"))
printf("Monitor_nD: %s is requesting automatic limits option 'auto' together with OpenACC.\n"
"WARNING this feature is NOT supported using OpenACC and has been disabled!\n", NAME_CURRENT_COMP);
#endif
#endif
%}
TRACE
%{
double transmit_he3=1.0;
double multiplier_capture=1.0;
double t0 = 0;
double t1 = 0;
int pp;
int intersect = 0;
char Flag_Restore = 0;
#ifdef OPENACC
#ifdef USE_OFF
off_struct thread_offdata = offdata;
#endif
#else
#define thread_offdata offdata
#endif
/* this is done automatically
STORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
*/
#ifdef USE_OFF
if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL"))
{
/* determine intersections with object */
intersect = off_intersect_all(&t0, &t1, NULL, NULL,
x,y,z, vx, vy, vz, 0, 0, 0, &thread_offdata );
if (Vars.Flag_mantid) {
if(intersect) {
Vars.OFF_polyidx=thread_offdata.nextintersect;
} else {
Vars.OFF_polyidx=-1;
}
}
}
else
#endif
if ( (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE)
|| (abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK) ) /* square xy or disk xy */
{
// propagate to xy plane and find intersection
// make sure the event is recoverable afterwards
t0 = t;
ALLOW_BACKPROP;
PROP_Z0;
if ( (t>=t0) && (z==0.0) ) // forward propagation to xy plane was successful
{
if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE)
{
// square xy
intersect = (x>=Vars.mxmin && x<=Vars.mxmax && y>=Vars.mymin && y<=Vars.mymax);
}
else
{
// disk xy
intersect = (SQR(x) + SQR(y)) <= SQR(Vars.Sphere_Radius);
}
}
else
{
intersect=0;
}
}
else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) /* sphere */
{
intersect = sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius);
/* intersect = (intersect && t0 > 0); */
}
else if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA)) /* cylinder */
{
intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius, Vars.Cylinder_Height);
}
else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) /* box */
{
intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz,
fabs(Vars.mxmax-Vars.mxmin), fabs(Vars.mymax-Vars.mymin), fabs(Vars.mzmax-Vars.mzmin));
}
else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_PREVIOUS) /* previous comp */
{ intersect = 1; }
if (intersect)
{
if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND)
|| (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA)
|| (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) )
{
/* check if we have to remove the top/bottom with BANANA shape */
if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA) {
if (intersect == 1) { // Entered and left through sides
if (t0 < 0 && t1 > 0) {
t0 = t; /* neutron was already inside ! */
}
if (t1 < 0 && t0 > 0) { /* neutron exit before entering !! */
t1 = t;
}
/* t0 is now time of incoming intersection with the detection area */
if ((Vars.Flag_Shape < 0) && (t1 > 0)) {
PROP_DT(t1); /* t1 outgoing beam */
} else {
PROP_DT(t0); /* t0 incoming beam */
}
} else if (intersect == 3 || intersect == 5) { // Entered from top or bottom, left through side
if ((Vars.Flag_Shape < 0) && (t1 > 0)) {
PROP_DT(t1); /* t1 outgoing beam */
} else {
intersect=0;
Flag_Restore=1;
}
} else if (intersect == 9 || intersect == 17) { // Entered through side, left from top or bottom
if ((Vars.Flag_Shape < 0) && (t1 > 0)) {
intersect=0;
Flag_Restore=1;
} else {
PROP_DT(t0); /* t0 incoming beam */
}
} else if (intersect == 13 || intersect == 19) { // Went through top/bottom on entry and exit
intersect=0;
Flag_Restore=1;
} else {
printf("Cylinder_intersect returned unexpected value %l\n", intersect);
}
} else {
// All other shapes than the BANANA
if (t0 < 0 && t1 > 0)
t0 = t; /* neutron was already inside ! */
if (t1 < 0 && t0 > 0) /* neutron exit before entering !! */
t1 = t;
/* t0 is now time of incoming intersection with the detection area */
if ((Vars.Flag_Shape < 0) && (t1 > 0))
PROP_DT(t1); /* t1 outgoing beam */
else
PROP_DT(t0); /* t0 incoming beam */
}
/* Final test if we are on lid / bottom of banana/sphere */
if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA || abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) {
if (Vars.Cylinder_Height && fabs(y) >= Vars.Cylinder_Height/2 - FLT_EPSILON) {
intersect=0;
Flag_Restore=1;
}
}
}
}
if (intersect)
{
/* Now get the data to monitor: current or keep from PreMonitor */
/* if (Vars.Flag_UsePreMonitor != 1)*/
/* {*/
/* Vars.cp = p;*/
/* Vars.cx = x;*/
/* Vars.cvx = vx;*/
/* Vars.csx = sx;*/
/* Vars.cy = y;*/
/* Vars.cvy = vy;*/
/* Vars.csy = sy;*/
/* Vars.cz = z;*/
/* Vars.cvz = vz;*/
/* Vars.csz = sz;*/
/* Vars.ct = t;*/
/* }*/
if ((Vars.He3_pressure > 0) && (t1 != t0) && ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX)))
{
transmit_he3 = exp(-7.417*Vars.He3_pressure*fabs(t1-t0)*2*PI*K2V);
/* will monitor the absorbed part */
p = p * (1-transmit_he3);
}
if (Vars.Flag_capture)
{
multiplier_capture = V2K*sqrt(vx*vx+vy*vy+vz*vz);
if (multiplier_capture != 0) multiplier_capture = 2*PI/multiplier_capture; /* lambda. lambda(2200 m/2) = 1.7985 Angs */
p = p * multiplier_capture/1.7985;
}
pp = Monitor_nD_Trace(&DEFS, &Vars, _particle);
if (pp==0.0)
{
ABSORB;
}
else if(pp==1)
{
SCATTER;
}
/*set weight to undetected part if capture and/or he3_pressure*/
if (Vars.He3_pressure > 0){
/* after monitor, only remains 1-p_detect */
p = p * transmit_he3/(1.0-transmit_he3);
}
if (Vars.Flag_capture){
p = p / multiplier_capture*1.7985;
}
if (Vars.Flag_parallel) /* back to neutron state before detection */
Flag_Restore = 1;
} /* end if intersection */
else {
if (Vars.Flag_Absorb && !Vars.Flag_parallel)
{
// restore neutron ray before absorbing for correct mcdisplay
RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
ABSORB;
}
else Flag_Restore = 1; /* no intersection, back to previous state */
}
if (Flag_Restore)
{
RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
}
%}
SAVE
%{
if (!nowritefile) {
/* save results, but do not free pointers */
detector = Monitor_nD_Save(&DEFS, &Vars);
}
%}
FINALLY
%{
/* free pointers */
Monitor_nD_Finally(&DEFS, &Vars);
%}
MCDISPLAY
%{
if (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL"))
{
off_display(offdata);
} else {
Monitor_nD_McDisplay(&DEFS, &Vars);
}
%}
END