/******************************************************************************* * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: TOFRes_monitor * * %I * Written by: Kristian Nielsen * Date: 1999 * Origin: Risoe * Modified by: EF, 16th Apr 2003: imported from Monitor_nD to enable many shapes * Modified by: T. Weber, Nov 2020: a) added live calculations; b) updated for McStas 3 / OpenACC * * Monitor for resolution calculations, TOFRes version, for use with TOFRes_sample in 3.x * * %D * A single detector/monitor, used together with the TOFRes_sample component to * compute instrument resolution functions. Outputs a list of neutron * scattering events in the sample along with their intensities in the * detector. The output file may be analyzed with the mcresplot front-end. * * Example: TOFRes_monitor(filename="Output.res", res_sample_comp="RSample", xmin=-0.1, xmax=0.1, ymin=-0.1, ymax=0.1) * * Setting the monitor geometry. * The optional parameter 'options' may be set as a string with the * following keywords. Default is rectangular ('square'): * box Box of size xwidth, yheight, zdepth * cylinder To get a cylindrical monitor (diameter is xwidth, height is yheight). * banana Same as cylinder, without top/bottom, on restricted angular area * disk Disk flat xy monitor. diameter is xwidth. * sphrere To get a spherical monitor (e.g. a 4PI) (diameter is xwidth). * square Square flat xy monitor (xwidth, yheight) * * %P * INPUT PARAMETERS: * * xmin: [m] Lower x bound of detector opening * xmax: [m] Upper x bound of detector opening * ymin: [m] Lower y bound of detector opening * ymax: [m] Upper y bound of detector opening * zmin: [m] Lower z bound of detector opening * zmax: [m] Upper z bound of detector opening * filename: [string] Name of output file. If unset, use automatic name * res_sample_comp: [WITH quotes] Name of TOFRes_sample component in the instrument definition * bufsize: [1] Number of events to store. Use 0 to store all * * OPTIONAL PARAMETERS (derived from Monitor_nD) * * xwidth: [m] Width/diameter of detector * yheight: [m] Height of detector * zdepth: [m] Thichness of detector * radius: [m] Radius of sphere/cylinder monitor * options: [str] String that specifies the geometry of the monitor * restore_neutron: [1] If set, the monitor does not influence the neutron state * live_calc: [1] If set, the monitor directly outputs the resolution matrix * * CALCULATED PARAMETERS: * * DEFS: [struct] structure containing Monitor_nD Defines * Vars: [struct] structure containing Monitor_nD variables * num_events: [long] number of recorded events * * %E *******************************************************************************/ DEFINE COMPONENT TOFRes_monitor SETTING PARAMETERS (string res_sample_comp, string filename=0, string options=0, xwidth=.1, yheight=.1, zdepth=0, radius=0, xmin=0, xmax=0, ymin=0, ymax=0, zmin=0, zmax=0, bufsize=0, int restore_neutron=0, int live_calc=1) /* these are protected C variables */ /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "monitor_nd-lib" %include "cov-lib" %} DECLARE%{ MonitornD_Defines_type DEFS; MonitornD_Variables_type Vars; unsigned long num_events; double* reso_events; double* reso_probabilities; char res_pi_var[20]; char res_ki_x_var[20]; char res_ki_y_var[20]; char res_ki_z_var[20]; char res_kf_x_var[20]; char res_kf_y_var[20]; char res_kf_z_var[20]; char res_rx_var[20]; char res_ry_var[20]; char res_rz_var[20]; /* Arrays for storing resolution matrix */ DArray2d Covar_p; DArray2d Covar_p2; DArray2d Covar_N; DArray2d Res_p; DArray2d Res_p2; DArray2d Res_N; int num_cores; %} INITIALIZE %{ /* Use instance name for monitor output if no input was given */ if (!strcmp(filename, "\0")) sprintf(filename, "%s", NAME_CURRENT_COMP); num_events = 0; reso_events = 0; reso_probabilities = 0; if(live_calc) { reso_events = malloc(mcget_ncount()*4*sizeof(double)); reso_probabilities = malloc(mcget_ncount()*sizeof(double)); } int i; char tmp[1024]; strcpy(Vars.compcurname, NAME_CURRENT_COMP); if (options != NULL) strncpy(tmp, options, 1024); else strcpy(tmp, ""); if (strstr(tmp, "list")) exit(fprintf(stderr, "TOFRes_monitor: %s: Error: Only use geometry keywords (remove list from 'option').\n", NAME_CURRENT_COMP)); if (!bufsize) sprintf(Vars.option, "%s borders list all, ud1 ud2 ud3 ud4 ud5 ud6 ud7 ud8 ud9 ud10 n", tmp); else sprintf(Vars.option, "%s borders list=%f, ud1 ud2 ud3 ud4 ud5 ud6 ud7 ud8 ud9 ud10 n", tmp, bufsize); if (radius) xwidth = 2*radius; Monitor_nD_Init(&DEFS, &Vars, xwidth, yheight, zdepth, xmin, xmax, ymin, ymax, zmin, zmax, 0, 0); Vars.Coord_Type[0] = DEFS.COORD_USERDOUBLE0; /* otherwise p is always the first variable */ if (Vars.Coord_Number != 11) exit(fprintf(stderr,"TOFRes_monitor: %s: Error: Invalid number of variables to monitor (%li).\n", NAME_CURRENT_COMP, Vars.Coord_Number+1)); /* set the labels */ /* we have to record ki_x ki_y ki_z kf_x kf_y kf_z x y z p_i p_f */ int idx = 0; strcpy(tmp,"ki_x"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"ki_y"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"ki_z"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"kf_x"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"kf_y"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"kf_z"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"x"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"y"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"z"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"p_i"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); strcpy(tmp,"p_f"); strcpy(Vars.Coord_Label[idx], tmp); strcpy(Vars.Coord_Var[idx++], tmp); if (filename != NULL) strncpy(Vars.Mon_File, filename, 128); /* Initialize uservar strings */ int *index_ptr=COMP_GETPAR3(TOFRes_sample, res_sample_comp, compindex); int index = *index_ptr; sprintf(res_pi_var, "res_pi_%i", index); sprintf(res_ki_x_var, "res_ki_x_%i", index); sprintf(res_ki_y_var, "res_ki_y_%i", index); sprintf(res_ki_z_var, "res_ki_z_%i", index); sprintf(res_kf_x_var, "res_kf_x_%i", index); sprintf(res_kf_y_var, "res_kf_y_%i", index); sprintf(res_kf_z_var, "res_kf_z_%i", index); sprintf(res_rx_var, "res_rx_%i", index); sprintf(res_ry_var, "res_ry_%i", index); sprintf(res_rz_var, "res_rz_%i", index); if(live_calc) { /* Allocate arrays to save cov and res matrix */ Covar_p = create_darr2d(4, 4); Covar_p2 = create_darr2d(4, 4); Covar_N = create_darr2d(4, 4); Res_p = create_darr2d(4, 4); Res_p2 = create_darr2d(4, 4); Res_N = create_darr2d(4, 4); #ifdef USE_MPI num_cores=mpi_node_count; #else num_cores=1; #endif } %} TRACE %{ double t0 = 0; double t1 = 0; int intersect = 0; unsigned long event_idx = 0; double event_ki[3], event_kf[3]; double event_pi, event_pf; double event_pos[3]; if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) /* square xy */ { PROP_Z0; intersect = (x>=Vars.mxmin && x<=Vars.mxmax && y>=Vars.mymin && y<=Vars.mymax); } else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK) /* disk xy */ { PROP_Z0; intersect = ((x*x + y*y) <= Vars.Sphere_Radius*Vars.Sphere_Radius); } 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); if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA) && (intersect != 1)) intersect = 0; /* remove top/bottom for banana */ } 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)); } 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)) { 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 sphere. */ if ((Vars.Flag_Shape < 0) && (t1 > 0)) PROP_DT(t1); /* t1 outgoing beam */ else PROP_DT(t0); /* t0 incoming beam */ } /* Now fetch data from the TOFRes_sample. */ if(p > 0. && (!bufsize || num_events < bufsize)) { /* old behaviour not supported by openacc: struct Res_sample_struct *s = (struct Res_sample_struct *)(COMP_GETPAR3(Res_sample, res_sample_comp, res_struct));*/ /* pi */ event_pi = particle_getvar(_particle, res_pi_var, NULL); if(event_pi > 0.) { /* pf */ event_pf = p / event_pi; /* ki */ event_ki[0] = particle_getvar(_particle, res_ki_x_var, NULL); event_ki[1] = particle_getvar(_particle, res_ki_y_var, NULL); event_ki[2] = particle_getvar(_particle, res_ki_z_var, NULL); /* kf */ event_kf[0] = particle_getvar(_particle, res_kf_x_var, NULL); event_kf[1] = particle_getvar(_particle, res_kf_y_var, NULL); event_kf[2] = particle_getvar(_particle, res_kf_z_var, NULL); /* pos */ event_pos[0] = particle_getvar(_particle, res_rx_var, NULL); event_pos[1] = particle_getvar(_particle, res_ry_var, NULL); event_pos[2] = particle_getvar(_particle, res_rz_var, NULL); #pragma acc atomic capture { event_idx = num_events++; } /* variables for Monitor_nD */ Vars.UserDoubles[0] = event_ki[0]; Vars.UserDoubles[1] = event_ki[1]; Vars.UserDoubles[2] = event_ki[2]; Vars.UserDoubles[3] = event_kf[0]; Vars.UserDoubles[4] = event_kf[1]; Vars.UserDoubles[5] = event_kf[2]; Vars.UserDoubles[6] = event_pos[0]; Vars.UserDoubles[7] = event_pos[1]; Vars.UserDoubles[8] = event_pos[2]; Vars.UserDoubles[9] = event_pi; Vars.UserDoubles[10] = event_pf; Monitor_nD_Trace(&DEFS, &Vars, _particle); /* live calculation */ if(live_calc) { reso_events[event_idx*4 + 0] = event_ki[0] - event_kf[0]; reso_events[event_idx*4 + 1] = event_ki[1] - event_kf[1]; reso_events[event_idx*4 + 2] = event_ki[2] - event_kf[2]; reso_events[event_idx*4 + 3] = tl2_k_to_E( event_ki[0], event_ki[1], event_ki[2], event_kf[0], event_kf[1], event_kf[2]); reso_probabilities[event_idx] = event_pi * event_pf; } } SCATTER; } /* if p */ } /* end if intersection */ if (restore_neutron) { RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } %} SAVE %{ int dbg_save_events = 0; /* save results, but do not free pointers */ Monitor_nD_Save(&DEFS, &Vars); /* live calculation */ if(live_calc && num_events > 0) { if(dbg_save_events) { /* save individual neutron events */ #ifndef USE_MPI const char* event_filename = "reso_events.dat"; #else char event_filename[256]; sprintf(event_filename, "reso_events_%d.dat", mpi_node_rank); #endif tl2_save_events(reso_events, reso_probabilities, event_filename, num_events); } double cov[4*4], res[4*4]; if(tl2_reso(reso_events, reso_probabilities, cov, res, num_events)) { printf("Resolution calculation used %d neutron events.\n", num_events); tl2_print_mat(cov, "Covariance matrix", 4, 4); tl2_print_mat(res, "Resolution matrix", 4, 4); printf("Please run \"mcresplot %s\" for a full analysis.\n", filename); } else { printf("Error: Resolution matrix could not be calculated."); } int i,j; for(i=0;i<4;i++) { for(j=0;j<4;j++) { /* "Events" */ Covar_N[i][j]=1; Res_N[i][j]=1; /* Covariance/resolution matrixces (potentiall pr. MPI node) */ Covar_p[i][j]=cov[4*i+j]/num_cores; Res_p[i][j]=res[4*i+j]/num_cores; /* "errors" */ Covar_p2[i][j]=(cov[4*i+j]/num_cores)*(cov[4*i+j]/num_cores); Res_p2[i][j]=(res[4*i+j]/num_cores)*(res[4*i+j]/num_cores); } } /* Nasty call to DETECTOR_OUT_2D to store covariance matrix */ char covar_fname[1024]; char res_fname[1024]; sprintf(covar_fname,"%s_%s", _comp->_name, "covar"); sprintf(res_fname,"%s_%s", _comp->_name, "resol"); DETECTOR_OUT_2D("Covariance", "Columns", "Rows", 0.0, 4.0, 0.0, 4.0, 4, 4, &Covar_N[0][0],&Covar_p[0][0],&Covar_p2[0][0], covar_fname); DETECTOR_OUT_2D("Resolution", "Columns", "Rows", 0.0, 4.0, 0.0, 4.0, 4, 4, &Res_N[0][0],&Res_p[0][0],&Res_p2[0][0], res_fname); } %} FINALLY %{ /* free pointers */ Monitor_nD_Finally(&DEFS, &Vars); if(live_calc) { free(reso_events); free(reso_probabilities); } destroy_darr2d(Covar_p); destroy_darr2d(Covar_p2); destroy_darr2d(Covar_N); destroy_darr2d(Res_p); destroy_darr2d(Res_p2); destroy_darr2d(Res_N); %} MCDISPLAY %{ Monitor_nD_McDisplay(&DEFS, &Vars); %} END