/******************************************************************************* * * McCode, neutron/xray ray-tracing package * Copyright (C) 1997-2009, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Runtime: share/mccode-r.c * * %Identification * Written by: KN * Date: Aug 29, 1997 * Release: McStas X.Y/McXtrace X.Y * Version: $Revision$ * * Runtime system for McStas and McXtrace. * Embedded within instrument in runtime mode. * Contains SECTIONS: * MPI handling (sum, send, recv) * format definitions * I/O * mcdisplay support * random numbers * coordinates handling * vectors math (solve 2nd order, normals, randvec...) * parameter handling * signal and main handlers * * Usage: Automatically embbeded in the c code whenever required. * * $Id$ * *******************************************************************************/ /******************************************************************************* * The I/O format definitions and functions *******************************************************************************/ /** Include header files to avoid implicit declarations (not allowed on LLVM) */ #include #include // UNIX specific headers (non-Windows) #if defined(__unix__) || defined(__APPLE__) #include #endif #include #ifdef _WIN32 #include # define mkdir( D, M ) _mkdir( D ) #endif #ifndef DANSE #ifdef MC_ANCIENT_COMPATIBILITY int mctraceenabled = 0; int mcdefaultmain = 0; #endif /* else defined directly in the McCode generated C code */ static long mcseed = 0; /* seed for random generator */ static long mcstartdate = 0; /* start simulation time */ static int mcdisable_output_files = 0; /* --no-output-files */ mcstatic int mcgravitation = 0; /* use gravitation flag, for PROP macros */ int mcMagnet = 0; /* magnet stack flag */ mcstatic int mcdotrace = 0; /* flag for --trace and messages for DISPLAY */ int mcallowbackprop = 0; /* flag to enable negative/backprop */ /* Number of particle histories to simulate. */ #ifdef NEUTRONICS mcstatic unsigned long long int mcncount = 1; mcstatic unsigned long long int mcrun_num = 0; #else mcstatic unsigned long long int mcncount = 1000000; mcstatic unsigned long long int mcrun_num = 0; #endif /* NEUTRONICS */ #else #include "mcstas-globals.h" #endif /* !DANSE */ /* SECTION: NeXus compression */ #ifndef NX_COMPRESSION #define NX_COMPRESSION NX_COMP_NONE #endif /* SECTION: MPI handling ==================================================== */ #ifdef USE_MPI /* MPI rank */ static int mpi_node_rank; static int mpi_node_root = 0; /******************************************************************************* * mc_MPI_Reduce: Gathers arrays from MPI nodes using Reduce function. *******************************************************************************/ int mc_MPI_Sum(double *sbuf, long count) { if (!sbuf || count <= 0) return(MPI_SUCCESS); /* nothing to reduce */ else { /* we must cut the buffer into blocks not exceeding the MPI max buffer size of 32000 */ long offset=0; double *rbuf=NULL; int length=MPI_REDUCE_BLOCKSIZE; /* defined in mccode-r.h */ int i=0; rbuf = calloc(count, sizeof(double)); if (!rbuf) exit(-fprintf(stderr, "Error: Out of memory %li (mc_MPI_Sum)\n", count*sizeof(double))); while (offset < count) { if (!length || offset+length > count-1) length=count-offset; else length=MPI_REDUCE_BLOCKSIZE; if (MPI_Reduce((double*)(sbuf+offset), (double*)(rbuf+offset), length, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) != MPI_SUCCESS) return MPI_ERR_COUNT; offset += length; } for (i=0; i count-1) length=count-offset; else length=MPI_REDUCE_BLOCKSIZE; if (MPI_Send((void*)(sbuf+offset*dsize), length, dtype, dest, tag++, MPI_COMM_WORLD) != MPI_SUCCESS) return MPI_ERR_COUNT; offset += length; } return MPI_SUCCESS; } /* mc_MPI_Send */ /******************************************************************************* * mc_MPI_Recv: Receives arrays from MPI nodes by blocks to avoid buffer limit * the buffer must have been allocated previously. *******************************************************************************/ int mc_MPI_Recv(void *sbuf, long count, MPI_Datatype dtype, int source) { int dsize; long offset=0; int tag=1; int length=MPI_REDUCE_BLOCKSIZE; /* defined in mccode-r.h */ if (!sbuf || count <= 0) return(MPI_SUCCESS); /* nothing to recv */ MPI_Type_size(dtype, &dsize); while (offset < count) { if (offset+length > count-1) length=count-offset; else length=MPI_REDUCE_BLOCKSIZE; if (MPI_Recv((void*)(sbuf+offset*dsize), length, dtype, source, tag++, MPI_COMM_WORLD, MPI_STATUS_IGNORE) != MPI_SUCCESS) return MPI_ERR_COUNT; offset += length; } return MPI_SUCCESS; } /* mc_MPI_Recv */ #endif /* USE_MPI */ /* SECTION: parameters handling ============================================= */ /* Instrument input parameter type handling. */ /******************************************************************************* * mcparm_double: extract double value from 's' into 'vptr' *******************************************************************************/ static int mcparm_double(char *s, void *vptr) { char *p; double *v = (double *)vptr; if (!s) { *v = 0; return(1); } *v = strtod(s, &p); if(*s == '\0' || (p != NULL && *p != '\0') || errno == ERANGE) return 0; /* Failed */ else return 1; /* Success */ } /******************************************************************************* * mcparminfo_double: display parameter type double *******************************************************************************/ static char * mcparminfo_double(char *parmname) { return "double"; } /******************************************************************************* * mcparmerror_double: display error message when failed extract double *******************************************************************************/ static void mcparmerror_double(char *parm, char *val) { fprintf(stderr, "Error: Invalid value '%s' for floating point parameter %s (mcparmerror_double)\n", val, parm); } /******************************************************************************* * mcparmprinter_double: convert double to string *******************************************************************************/ static void mcparmprinter_double(char *f, void *vptr) { double *v = (double *)vptr; sprintf(f, "%g", *v); } /******************************************************************************* * mcparm_int: extract int value from 's' into 'vptr' *******************************************************************************/ static int mcparm_int(char *s, void *vptr) { char *p; int *v = (int *)vptr; long x; if (!s) { *v = 0; return(1); } *v = 0; x = strtol(s, &p, 10); if(x < INT_MIN || x > INT_MAX) return 0; /* Under/overflow */ *v = x; if(*s == '\0' || (p != NULL && *p != '\0') || errno == ERANGE) return 0; /* Failed */ else return 1; /* Success */ } /******************************************************************************* * mcparminfo_int: display parameter type int *******************************************************************************/ static char * mcparminfo_int(char *parmname) { return "int"; } /******************************************************************************* * mcparmerror_int: display error message when failed extract int *******************************************************************************/ static void mcparmerror_int(char *parm, char *val) { fprintf(stderr, "Error: Invalid value '%s' for integer parameter %s (mcparmerror_int)\n", val, parm); } /******************************************************************************* * mcparmprinter_int: convert int to string *******************************************************************************/ static void mcparmprinter_int(char *f, void *vptr) { int *v = (int *)vptr; sprintf(f, "%d", *v); } /******************************************************************************* * mcparm_string: extract char* value from 's' into 'vptr' (copy) *******************************************************************************/ static int mcparm_string(char *s, void *vptr) { char **v = (char **)vptr; if (!s) { *v = NULL; return(1); } *v = (char *)malloc(strlen(s) + 1); if(*v == NULL) { exit(-fprintf(stderr, "Error: Out of memory %li (mcparm_string).\n", (long)strlen(s) + 1)); } strcpy(*v, s); return 1; /* Success */ } /******************************************************************************* * mcparminfo_string: display parameter type string *******************************************************************************/ static char * mcparminfo_string(char *parmname) { return "string"; } /******************************************************************************* * mcparmerror_string: display error message when failed extract string *******************************************************************************/ static void mcparmerror_string(char *parm, char *val) { fprintf(stderr, "Error: Invalid value '%s' for string parameter %s (mcparmerror_string)\n", val, parm); } /******************************************************************************* * mcparmprinter_string: convert string to string (including esc chars) *******************************************************************************/ static void mcparmprinter_string(char *f, void *vptr) { char **v = (char **)vptr; char *p; if (!*v) { *f='\0'; return; } strcpy(f, ""); for(p = *v; *p != '\0'; p++) { switch(*p) { case '\n': strcat(f, "\\n"); break; case '\r': strcat(f, "\\r"); break; case '"': strcat(f, "\\\""); break; case '\\': strcat(f, "\\\\"); break; default: strncat(f, p, 1); } } /* strcat(f, "\""); */ } /* mcparmprinter_string */ /* now we may define the parameter structure, using previous functions */ static struct { int (*getparm)(char *, void *); char * (*parminfo)(char *); void (*error)(char *, char *); void (*printer)(char *, void *); } mcinputtypes[] = { { mcparm_double, mcparminfo_double, mcparmerror_double, mcparmprinter_double }, { mcparm_int, mcparminfo_int, mcparmerror_int, mcparmprinter_int }, { mcparm_string, mcparminfo_string, mcparmerror_string, mcparmprinter_string } }; /******************************************************************************* * mcestimate_error: compute sigma from N,p,p2 in Gaussian large numbers approx *******************************************************************************/ double mcestimate_error(double N, double p1, double p2) { double pmean, n1; if(N <= 1) return p1; pmean = p1 / N; n1 = N - 1; /* Note: underflow may cause p2 to become zero; the fabs() below guards against this. */ return sqrt((N/n1)*fabs(p2 - pmean*pmean)); } double (*mcestimate_error_p) (double V2, double psum, double p2sum)=mcestimate_error; /* ========================================================================== */ /* MCCODE_R_IO_C */ /* ========================================================================== */ #ifndef MCCODE_R_IO_C #define MCCODE_R_IO_C "$Revision$" /* SECTION: file i/o handling ================================================ */ #ifndef HAVE_STRCASESTR // from msysgit: https://code.google.com/p/msysgit/source/browse/compat/strcasestr.c char *strcasestr(const char *haystack, const char *needle) { int nlen = strlen(needle); int hlen = strlen(haystack) - nlen + 1; int i; for (i = 0; i < hlen; i++) { int j; for (j = 0; j < nlen; j++) { unsigned char c1 = haystack[i+j]; unsigned char c2 = needle[j]; if (toupper(c1) != toupper(c2)) goto next; } return (char *) haystack + i; next: ; } return NULL; } #endif #ifndef HAVE_STRCASECMP int strcasecmp( const char *s1, const char *s2 ) { int c1, c2; do { c1 = tolower( (unsigned char) *s1++ ); c2 = tolower( (unsigned char) *s2++ ); } while (c1 == c2 && c1 != 0); return c2 > c1 ? -1 : c1 > c2; } #endif /******************************************************************************* * mcfull_file: allocates a full file name=mcdirname+file. Catenate extension if missing. *******************************************************************************/ char *mcfull_file(char *name, char *ext) { int dirlen=0; char *mem =NULL; dirlen = mcdirname ? strlen(mcdirname) : 0; mem = (char*)malloc(dirlen + strlen(name) + CHAR_BUF_LENGTH); if(!mem) { exit(-fprintf(stderr, "Error: Out of memory %li (mcfull_file)\n", (long)(dirlen + strlen(name) + 256))); } strcpy(mem, ""); /* prepend directory name to path if name does not contain a path */ if (dirlen > 0 && !strchr(name, MC_PATHSEP_C)) { strcat(mem, mcdirname); strcat(mem, MC_PATHSEP_S); } /* dirlen */ strcat(mem, name); if (!strchr(name, '.') && ext && strlen(ext)) { /* add extension if not in file name already */ strcat(mem, "."); strcat(mem, ext); } return(mem); } /* mcfull_file */ /******************************************************************************* * mcnew_file: opens a new file within mcdirname if non NULL * the file is opened in "a" (append, create if does not exist) * the extension 'ext' is added if the file name does not include one. * the last argument is set to 0 if file did not exist, else to 1. *******************************************************************************/ FILE *mcnew_file(char *name, char *ext, int *exists) { char *mem; FILE *file=NULL; if (!name || strlen(name) == 0 || mcdisable_output_files) return(NULL); mem = mcfull_file(name, ext); /* create mcdirname/name.ext */ /* check for existence */ file = fopen(mem, "r"); /* for reading -> fails if does not exist */ if (file) { fclose(file); *exists=1; } else *exists=0; /* open the file for writing/appending */ #ifdef USE_NEXUS if (mcformat && strcasestr(mcformat, "NeXus")) { /* NXhandle nxhandle is defined in the .h with USE_NEXUS */ NXaccess mode = (*exists ? NXACC_CREATE5 | NXACC_RDWR : NXACC_CREATE5); if (NXopen(mem, mode, &nxhandle) != NX_OK) file = NULL; else file = (FILE*)&nxhandle; /* to make it non NULL */ } else #endif file = fopen(mem, "a+"); if(!file) fprintf(stderr, "Warning: could not open output file '%s' for %s (mcnew_file)\n", mem, *exists ? "append" : "create"); free(mem); return file; } /* mcnew_file */ /******************************************************************************* * mcdetector_statistics: compute detector statistics, error bars, [x I I_err N] 1D * RETURN: updated detector structure * Used by: mcdetector_import *******************************************************************************/ MCDETECTOR mcdetector_statistics( MCDETECTOR detector) { if (!detector.p1 || !detector.m) return(detector); /* compute statistics and update MCDETECTOR structure ===================== */ double sum_z = 0, min_z = 0, max_z = 0; double fmon_x =0, smon_x = 0, fmon_y =0, smon_y=0, mean_z=0; double Nsum=0, P2sum=0; double sum_xz = 0, sum_yz = 0, sum_x = 0, sum_y = 0, sum_x2z = 0, sum_y2z = 0; int i,j; char hasnan=0, hasinf=0; char israw = ((char*)strcasestr(detector.format,"raw") != NULL); double *this_p1=NULL; /* new 1D McCode array [x I E N]. Freed after writing data */ /* if McCode/PGPLOT and rank==1 we create a new m*4 data block=[x I E N] */ if (detector.rank == 1 && strcasestr(detector.format,"McCode")) { this_p1 = (double *)calloc(detector.m*detector.n*detector.p*4, sizeof(double)); if (!this_p1) exit(-fprintf(stderr, "Error: Out of memory creating %li 1D " MCCODE_STRING " data set for file '%s' (mcdetector_import)\n", detector.m*detector.n*detector.p*4*sizeof(double*), detector.filename)); } max_z = min_z = detector.p1[0]; /* compute sum and moments (not for lists) */ if (!strcasestr(detector.format,"list") && detector.m) for(j = 0; j < detector.n*detector.p; j++) { for(i = 0; i < detector.m; i++) { double x,y,z; double N, E; long index= !detector.istransposed ? i*detector.n*detector.p + j : i+j*detector.m; char hasnaninf=0; if (detector.m) x = detector.xmin + (i + 0.5)/detector.m*(detector.xmax - detector.xmin); else x = 0; if (detector.n && detector.p) y = detector.ymin + (j + 0.5)/detector.n/detector.p*(detector.ymax - detector.ymin); else y = 0; z = detector.p1[index]; N = detector.p0 ? detector.p0[index] : 1; E = detector.p2 ? detector.p2[index] : 0; if (detector.p2 && !israw) detector.p2[index] = (*mcestimate_error_p)(detector.p0[index],detector.p1[index],detector.p2[index]); /* set sigma */ if (detector.rank == 1 && this_p1 && strcasestr(detector.format,"McCode")) { /* fill-in 1D McCode array [x I E N] */ this_p1[index*4] = x; this_p1[index*4+1] = z; this_p1[index*4+2] = detector.p2 ? detector.p2[index] : 0; this_p1[index*4+3] = N; } if (isnan(z) || isnan(E) || isnan(N)) hasnaninf=hasnan=1; if (isinf(z) || isinf(E) || isinf(N)) hasnaninf=hasinf=1; /* compute stats integrals */ if (!hasnaninf) { sum_xz += x*z; sum_yz += y*z; sum_x += x; sum_y += y; sum_z += z; sum_x2z += x*x*z; sum_y2z += y*y*z; if (z > max_z) max_z = z; if (z < min_z) min_z = z; Nsum += N; P2sum += E; } } } /* for j */ /* compute 1st and 2nd moments. For lists, sum_z=0 so this is skipped. */ if (sum_z && detector.n*detector.m*detector.p) { fmon_x = sum_xz/sum_z; fmon_y = sum_yz/sum_z; smon_x = sum_x2z/sum_z-fmon_x*fmon_x; smon_x = smon_x > 0 ? sqrt(smon_x) : 0; smon_y = sum_y2z/sum_z-fmon_y*fmon_y; smon_y = smon_y > 0 ? sqrt(smon_y) : 0; mean_z = sum_z/detector.n/detector.m/detector.p; } /* store statistics into detector */ detector.intensity = sum_z; detector.error = Nsum ? (*mcestimate_error_p)(Nsum, sum_z, P2sum) : 0; detector.events = Nsum; detector.min = min_z; detector.max = max_z; detector.mean = mean_z; detector.centerX = fmon_x; detector.halfwidthX= smon_x; detector.centerY = fmon_y; detector.halfwidthY= smon_y; /* if McCode/PGPLOT and rank==1 replace p1 with new m*4 1D McCode and clear others */ if (detector.rank == 1 && this_p1 && strcasestr(detector.format,"McCode")) { detector.p1 = this_p1; detector.n = detector.m; detector.m = 4; detector.p0 = detector.p2 = NULL; detector.istransposed = 1; } if (detector.n*detector.m*detector.p > 1) snprintf(detector.signal, CHAR_BUF_LENGTH, "Min=%g; Max=%g; Mean=%g;", detector.min, detector.max, detector.mean); else strcpy(detector.signal, "None"); snprintf(detector.values, CHAR_BUF_LENGTH, "%g %g %g", detector.intensity, detector.error, detector.events); switch (detector.rank) { case 1: snprintf(detector.statistics, CHAR_BUF_LENGTH, "X0=%g; dX=%g;", detector.centerX, detector.halfwidthX); break; case 2: case 3: snprintf(detector.statistics, CHAR_BUF_LENGTH, "X0=%g; dX=%g; Y0=%g; dY=%g;", detector.centerX, detector.halfwidthX, detector.centerY, detector.halfwidthY); break; default: strcpy(detector.statistics, "None"); } if (hasnan) printf("WARNING: Nan detected in component/file %s %s\n", detector.component, strlen(detector.filename) ? detector.filename : ""); if (hasinf) printf("WARNING: Inf detected in component/file %s %s\n", detector.component, strlen(detector.filename) ? detector.filename : ""); return(detector); } /* mcdetector_statistics */ /******************************************************************************* * mcdetector_import: build detector structure, merge non-lists from MPI * compute basic stat, write "Detector:" line * RETURN: detector structure. Invalid data if detector.p1 == NULL * Invalid detector sets m=0 and filename="" * Simulation data sets m=0 and filename=mcsiminfo_name * This function is equivalent to the old 'mcdetector_out', returning a structure *******************************************************************************/ MCDETECTOR mcdetector_import( char *format, char *component, char *title, long m, long n, long p, char *xlabel, char *ylabel, char *zlabel, char *xvar, char *yvar, char *zvar, double x1, double x2, double y1, double y2, double z1, double z2, char *filename, double *p0, double *p1, double *p2, Coords position) { time_t t; /* for detector.date */ long date_l; /* date as a long number */ char istransposed=0; char c[CHAR_BUF_LENGTH]; /* temp var for signal label */ MCDETECTOR detector; /* build MCDETECTOR structure ============================================= */ /* make sure we do not have NULL for char fields */ /* these also apply to simfile */ strncpy (detector.filename, filename ? filename : "", CHAR_BUF_LENGTH); strncpy (detector.format, format ? format : "McCode" , CHAR_BUF_LENGTH); /* add extension if missing */ if (strlen(detector.filename) && !strchr(detector.filename, '.')) { /* add extension if not in file name already */ strcat(detector.filename, ".dat"); } strncpy (detector.component, component ? component : MCCODE_STRING " component", CHAR_BUF_LENGTH); snprintf(detector.instrument, CHAR_BUF_LENGTH, "%s (%s)", mcinstrument_name, mcinstrument_source); snprintf(detector.user, CHAR_BUF_LENGTH, "%s on %s", getenv("USER") ? getenv("USER") : MCCODE_NAME, getenv("HOST") ? getenv("HOST") : "localhost"); time(&t); /* get current write time */ date_l = (long)t; /* same but as a long */ snprintf(detector.date, CHAR_BUF_LENGTH, "%s", ctime(&t)); if (strlen(detector.date)) detector.date[strlen(detector.date)-1] = '\0'; /* remove last \n in date */ detector.date_l = date_l; if (!mcget_run_num() || mcget_run_num() >= mcget_ncount()) snprintf(detector.ncount, CHAR_BUF_LENGTH, "%llu", mcget_ncount() #ifdef USE_MPI *mpi_node_count #endif ); else snprintf(detector.ncount, CHAR_BUF_LENGTH, "%g/%g", (double)mcget_run_num(), (double)mcget_ncount()); detector.p0 = p0; detector.p1 = p1; detector.p2 = p2; /* handle transposition (not for NeXus) */ if (!strcasestr(detector.format, "NeXus")) { if (m<0 || n<0 || p<0) istransposed = !istransposed; if (strcasestr(detector.format, "transpose")) istransposed = !istransposed; if (istransposed) { /* do the swap once for all */ long i=m; m=n; n=i; } } m=labs(m); n=labs(n); p=labs(p); /* make sure dimensions are positive */ detector.istransposed = istransposed; /* determine detector rank (dimensionality) */ if (!m || !n || !p || !p1) detector.rank = 4; /* invalid: exit with m=0 filename="" */ else if (m*n*p == 1) detector.rank = 0; /* 0D */ else if (n == 1 || m == 1) detector.rank = 1; /* 1D */ else if (p == 1) detector.rank = 2; /* 2D */ else detector.rank = 3; /* 3D */ /* from rank, set type */ switch (detector.rank) { case 0: strcpy(detector.type, "array_0d"); m=n=p=1; break; case 1: snprintf(detector.type, CHAR_BUF_LENGTH, "array_1d(%ld)", m*n*p); m *= n*p; n=p=1; break; case 2: snprintf(detector.type, CHAR_BUF_LENGTH, "array_2d(%ld, %ld)", m, n*p); n *= p; p=1; break; case 3: snprintf(detector.type, CHAR_BUF_LENGTH, "array_3d(%ld, %ld, %ld)", m, n, p); break; default: m=0; strcpy(detector.type, ""); strcpy(detector.filename, "");/* invalid */ } detector.m = m; detector.n = n; detector.p = p; /* these only apply to detector files ===================================== */ snprintf(detector.position, CHAR_BUF_LENGTH, "%g %g %g", position.x, position.y, position.z); /* may also store actual detector orientation in the future */ strncpy(detector.title, title && strlen(title) ? title : component, CHAR_BUF_LENGTH); strncpy(detector.xlabel, xlabel && strlen(xlabel) ? xlabel : "X", CHAR_BUF_LENGTH); /* axis labels */ strncpy(detector.ylabel, ylabel && strlen(ylabel) ? ylabel : "Y", CHAR_BUF_LENGTH); strncpy(detector.zlabel, zlabel && strlen(zlabel) ? zlabel : "Z", CHAR_BUF_LENGTH); strncpy(detector.xvar, xvar && strlen(xvar) ? xvar : "x", CHAR_BUF_LENGTH); /* axis variables */ strncpy(detector.yvar, yvar && strlen(yvar) ? yvar : detector.xvar, CHAR_BUF_LENGTH); strncpy(detector.zvar, zvar && strlen(zvar) ? zvar : detector.yvar, CHAR_BUF_LENGTH); /* set "variables" as e.g. "I I_err N" */ strcpy(c, "I "); if (strlen(detector.zvar)) strncpy(c, detector.zvar,32); else if (strlen(detector.yvar)) strncpy(c, detector.yvar,32); else if (strlen(detector.xvar)) strncpy(c, detector.xvar,32); if (detector.rank == 1) snprintf(detector.variables, CHAR_BUF_LENGTH, "%s %s %s_err N", detector.xvar, c, c); else snprintf(detector.variables, CHAR_BUF_LENGTH, "%s %s_err N", c, c); /* limits */ detector.xmin = x1; detector.xmax = x2; detector.ymin = y1; detector.ymax = y2; detector.zmin = z1; detector.zmax = z2; if (abs(detector.rank) == 1) snprintf(detector.limits, CHAR_BUF_LENGTH, "%g %g", x1, x2); else if (detector.rank == 2) snprintf(detector.limits, CHAR_BUF_LENGTH, "%g %g %g %g", x1, x2, y1, y2); else snprintf(detector.limits, CHAR_BUF_LENGTH, "%g %g %g %g %g %g", x1, x2, y1, y2, z1, z2); /* if MPI and nodes_nb > 1: reduce data sets when using MPI =============== */ #ifdef USE_MPI if (!strcasestr(detector.format,"list") && mpi_node_count > 1 && m) { /* we save additive data: reduce everything into mpi_node_root */ if (p0) mc_MPI_Sum(p0, m*n*p); if (p1) mc_MPI_Sum(p1, m*n*p); if (p2) mc_MPI_Sum(p2, m*n*p); if (!p0) { /* additive signal must be then divided by the number of nodes */ int i; for (i=0; i CHAR_BUF_LENGTH) break; snprintf(ThisParam, CHAR_BUF_LENGTH, " %s(%s)", mcinputtable[i].name, (*mcinputtypes[mcinputtable[i].type].parminfo) (mcinputtable[i].name)); strcat(Parameters, ThisParam); if (strlen(Parameters) >= CHAR_BUF_LENGTH-64) break; } /* output data ============================================================ */ if (f != stdout) fprintf(f, "%sFile: %s%c%s\n", pre, mcdirname, MC_PATHSEP_C, mcsiminfo_name); else fprintf(f, "%sCreator: %s\n", pre, MCCODE_STRING); fprintf(f, "%sSource: %s\n", pre, mcinstrument_source); fprintf(f, "%sParameters: %s\n", pre, Parameters); fprintf(f, "%sTrace_enabled: %s\n", pre, mctraceenabled ? "yes" : "no"); fprintf(f, "%sDefault_main: %s\n", pre, mcdefaultmain ? "yes" : "no"); fprintf(f, "%sEmbedded_runtime: %s\n", pre, #ifdef MC_EMBEDDED_RUNTIME "yes" #else "no" #endif ); fflush(f); } /* mcinfo_out */ /******************************************************************************* * mcruninfo_out_backend: output simulation tags/info (both in SIM and data files) * Used in: mcsiminfo_init (ascii case), mcdetector_out_xD_ascii, mcinfo(stdout) *******************************************************************************/ static void mcruninfo_out_backend(char *pre, FILE *f, int info) { int i; char Parameters[CHAR_BUF_LENGTH]; if (!f || mcdisable_output_files) return; fprintf(f, "%sFormat: %s%s\n", pre, mcformat && strlen(mcformat) ? mcformat : MCCODE_NAME, mcformat && strcasestr(mcformat,"McCode") ? " with text headers" : ""); fprintf(f, "%sURL: %s\n", pre, "http://www.mccode.org"); fprintf(f, "%sCreator: %s\n", pre, MCCODE_STRING); fprintf(f, "%sInstrument: %s\n", pre, mcinstrument_source); fprintf(f, "%sNcount: %llu\n", pre, mcget_ncount()); fprintf(f, "%sTrace: %s\n", pre, mcdotrace ? "yes" : "no"); fprintf(f, "%sGravitation: %s\n", pre, mcgravitation ? "yes" : "no"); snprintf(Parameters, CHAR_BUF_LENGTH, "%ld", mcseed); fprintf(f, "%sSeed: %s\n", pre, Parameters); fprintf(f, "%sDirectory: %s\n", pre, mcdirname ? mcdirname : "."); #ifdef USE_MPI if (mpi_node_count > 1) fprintf(f, "%sNodes: %i\n", pre, mpi_node_count); #endif /* output parameter string ================================================ */ for(i = 0; i < mcnumipar; i++) { if (!info){ (*mcinputtypes[mcinputtable[i].type].printer)(Parameters, mcinputtable[i].par); fprintf(f, "%sParam: %s=%s\n", pre, mcinputtable[i].name, Parameters); }else{ /*if an info run, some variables might not have values. Flag these by "NULL"*/ if(mcinputtable[i].val && strlen(mcinputtable[i].val)){ /* ... those with defautl values*/ (*mcinputtypes[mcinputtable[i].type].printer)(Parameters, mcinputtable[i].par); fprintf(f, "%sParam: %s=%s\n", pre, mcinputtable[i].name, Parameters); }else{ /* ... and those without */ fprintf(f, "%sParam: %s=NULL\n", pre, mcinputtable[i].name); } } } } /* mcruninfo_out_backend */ /************************ * wrapper function to mcruninfo_out_backend * Regular runs use this whereas the single call from mcinfo is directly to the backend *************************/ static void mcruninfo_out(char *pre, FILE *f){ mcruninfo_out_backend(pre,f,0); } /******************************************************************************* * mcsiminfo_out: wrapper to fprintf(mcsiminfo_file) *******************************************************************************/ void mcsiminfo_out(char *format, ...) { va_list ap; if(mcsiminfo_file && !mcdisable_output_files) { va_start(ap, format); vfprintf(mcsiminfo_file, format, ap); va_end(ap); } } /* mcsiminfo_out */ /******************************************************************************* * mcdatainfo_out: output detector header * mcdatainfo_out(prefix, file_handle, detector) writes info to data file *******************************************************************************/ static void mcdatainfo_out(char *pre, FILE *f, MCDETECTOR detector) { if (!f || !detector.m || mcdisable_output_files) return; /* output data ============================================================ */ fprintf(f, "%sDate: %s (%li)\n", pre, detector.date, detector.date_l); fprintf(f, "%stype: %s\n", pre, detector.type); fprintf(f, "%sSource: %s\n", pre, detector.instrument); fprintf(f, "%scomponent: %s\n", pre, detector.component); fprintf(f, "%sposition: %s\n", pre, detector.position); fprintf(f, "%stitle: %s\n", pre, detector.title); fprintf(f, !mcget_run_num() || mcget_run_num() >= mcget_ncount() ? "%sNcount: %s\n" : "%sratio: %s\n", pre, detector.ncount); if (strlen(detector.filename)) { fprintf(f, "%sfilename: %s\n", pre, detector.filename); } fprintf(f, "%sstatistics: %s\n", pre, detector.statistics); fprintf(f, "%ssignal: %s\n", pre, detector.signal); fprintf(f, "%svalues: %s\n", pre, detector.values); if (detector.rank >= 1) { fprintf(f, "%sxvar: %s\n", pre, detector.xvar); fprintf(f, "%syvar: %s\n", pre, detector.yvar); fprintf(f, "%sxlabel: %s\n", pre, detector.xlabel); fprintf(f, "%sylabel: %s\n", pre, detector.ylabel); if (detector.rank > 1) { fprintf(f, "%szvar: %s\n", pre, detector.zvar); fprintf(f, "%szlabel: %s\n", pre, detector.zlabel); } } fprintf(f, abs(detector.rank)==1 ? "%sxlimits: %s\n" : "%sxylimits: %s\n", pre, detector.limits); fprintf(f, "%svariables: %s\n", pre, strcasestr(detector.format, "list") ? detector.ylabel : detector.variables); fflush(f); } /* mcdatainfo_out */ /* mcdetector_out_array_ascii: output a single array to a file * m: columns * n: rows * p: array * f: file handle (already opened) */ static void mcdetector_out_array_ascii(long m, long n, double *p, FILE *f, char istransposed) { if(f) { int i,j; for(j = 0; j < n; j++) { for(i = 0; i < m; i++) { fprintf(f, "%.10g ", p[!istransposed ? i*n + j : j*m+i]); } fprintf(f,"\n"); } } } /* mcdetector_out_array_ascii */ /******************************************************************************* * mcdetector_out_0D_ascii: called by mcdetector_out_0D for ascii output *******************************************************************************/ MCDETECTOR mcdetector_out_0D_ascii(MCDETECTOR detector) { int exists=0; FILE *outfile = NULL; /* Write data set information to simulation description file. */ MPI_MASTER( mcsiminfo_out("\nbegin data\n"); // detector.component mcdatainfo_out(" ", mcsiminfo_file, detector); mcsiminfo_out("end data\n"); /* Don't write if filename is NULL: mcnew_file handles this (return NULL) */ outfile = mcnew_file(detector.component, "dat", &exists); if(outfile) { /* write data file header and entry in simulation description file */ mcruninfo_out( "# ", outfile); mcdatainfo_out("# ", outfile, detector); /* write I I_err N */ fprintf(outfile, "%g %g %g\n", detector.intensity, detector.error, detector.events); fclose(outfile); } ); /* MPI_MASTER */ return(detector); } /* mcdetector_out_0D_ascii */ /******************************************************************************* * mcdetector_out_1D_ascii: called by mcdetector_out_1D for ascii output *******************************************************************************/ MCDETECTOR mcdetector_out_1D_ascii(MCDETECTOR detector) { int exists=0; FILE *outfile = NULL; MPI_MASTER( /* Write data set information to simulation description file. */ mcsiminfo_out("\nbegin data\n"); // detector.filename mcdatainfo_out(" ", mcsiminfo_file, detector); mcsiminfo_out("end data\n"); /* Loop over array elements, writing to file. */ /* Don't write if filename is NULL: mcnew_file handles this (return NULL) */ outfile = mcnew_file(detector.filename, "dat", &exists); if(outfile) { /* write data file header and entry in simulation description file */ mcruninfo_out( "# ", outfile); mcdatainfo_out("# ", outfile, detector); /* output the 1D array columns */ mcdetector_out_array_ascii(detector.m, detector.n, detector.p1, outfile, detector.istransposed); fclose(outfile); } ); /* MPI_MASTER */ return(detector); } /* mcdetector_out_1D_ascii */ /******************************************************************************* * mcdetector_out_2D_ascii: called by mcdetector_out_2D for ascii output *******************************************************************************/ MCDETECTOR mcdetector_out_2D_ascii(MCDETECTOR detector) { int exists=0; FILE *outfile = NULL; MPI_MASTER( /* Loop over array elements, writing to file. */ /* Don't write if filename is NULL: mcnew_file handles this (return NULL) */ outfile = mcnew_file(detector.filename, "dat", &exists); if(outfile) { /* write header only if file has just been created (not appending) */ if (!exists) { /* Write data set information to simulation description file. */ mcsiminfo_out("\nbegin data\n"); // detector.filename mcdatainfo_out(" ", mcsiminfo_file, detector); mcsiminfo_out("end data\n"); mcruninfo_out( "# ", outfile); mcdatainfo_out("# ", outfile, detector); } fprintf(outfile, "# Data [%s/%s] %s:\n", detector.component, detector.filename, detector.zvar); mcdetector_out_array_ascii(detector.m, detector.n*detector.p, detector.p1, outfile, detector.istransposed); if (detector.p2) { fprintf(outfile, "# Errors [%s/%s] %s_err:\n", detector.component, detector.filename, detector.zvar); mcdetector_out_array_ascii(detector.m, detector.n*detector.p, detector.p2, outfile, detector.istransposed); } if (detector.p0) { fprintf(outfile, "# Events [%s/%s] N:\n", detector.component, detector.filename); mcdetector_out_array_ascii(detector.m, detector.n*detector.p, detector.p0, outfile, detector.istransposed); } fclose(outfile); if (!exists) { if (strcasestr(detector.format, "list")) printf("Events: \"%s\"\n", strlen(detector.filename) ? detector.filename : detector.component); } } /* if outfile */ ); /* MPI_MASTER */ #ifdef USE_MPI if (strcasestr(detector.format, "list") && mpi_node_count > 1) { int node_i=0; /* loop along MPI nodes to write sequentially */ for(node_i=0; node_i strlen(original)) n = strlen(original); else original += strlen(original)-n; strncpy(valid, original, n); for (i=0; i < n; i++) { if ( (valid[i] > 122) || (valid[i] < 32) || (strchr("!\"#$%&'()*+,-.:;<=>?@[\\]^`/ \n\r\t", valid[i]) != NULL) ) { if (i) valid[i] = '_'; else valid[i] = 'm'; } } valid[i] = '\0'; return(valid); } /* strcpy_valid */ /* end ascii output section ================================================= */ #ifdef USE_NEXUS /* ========================================================================== */ /* NeXus output */ /* ========================================================================== */ #define nxprintf(...) nxstr('d', __VA_ARGS__) #define nxprintattr(...) nxstr('a', __VA_ARGS__) /******************************************************************************* * nxstr: output a tag=value data set (char) in NeXus/current group * when 'format' is larger that 1024 chars it is used as value for the 'tag' * else the value is assembled with format and following arguments. * type='d' -> data set * 'a' -> attribute for current data set *******************************************************************************/ static int nxstr(char type, NXhandle *f, char *tag, char *format, ...) { va_list ap; char value[CHAR_BUF_LENGTH]; int i; int ret=NX_OK; if (!tag || !format || !strlen(tag) || !strlen(format)) return(NX_OK); /* assemble the value string */ if (strlen(format) < CHAR_BUF_LENGTH) { va_start(ap, format); ret = vsnprintf(value, CHAR_BUF_LENGTH, format, ap); va_end(ap); i = strlen(value); } else { i = strlen(format); } if (type == 'd') { /* open/put/close data set */ if (NXmakedata (f, tag, NX_CHAR, 1, &i) != NX_OK) return(NX_ERROR); NXopendata (f, tag); if (strlen(format) < CHAR_BUF_LENGTH) ret = NXputdata (f, value); else ret = NXputdata (f, format); NXclosedata(f); } else { if (strlen(format) < CHAR_BUF_LENGTH) ret = NXputattr (f, tag, value, strlen(value), NX_CHAR); else ret = NXputattr (f, tag, format, strlen(format), NX_CHAR); } return(ret); } /* nxstr */ /******************************************************************************* * mcinfo_readfile: read a full file into a string buffer which is allocated * Think to free the buffer after use. * Used in: mcinfo_out_nexus (nexus) *******************************************************************************/ char *mcinfo_readfile(char *filename) { FILE *f = fopen(filename, "rb"); if (!f) return(NULL); fseek(f, 0, SEEK_END); long fsize = ftell(f); rewind(f); char *string = malloc(fsize + 1); if (string) { int n = fread(string, fsize, 1, f); fclose(f); string[fsize] = 0; } return(string); } /******************************************************************************* * mcinfo_out: output instrument/simulation groups in NeXus file * Used in: mcsiminfo_init (nexus) *******************************************************************************/ static void mcinfo_out_nexus(NXhandle f) { FILE *fid; /* for intrument source code/C/IDF */ char *buffer=NULL; time_t t =time(NULL); /* for date */ char entry0[CHAR_BUF_LENGTH]; int count=0; char name[CHAR_BUF_LENGTH]; char class[CHAR_BUF_LENGTH]; if (!f || mcdisable_output_files) return; /* write NeXus NXroot attributes */ /* automatically added: file_name, HDF5_Version, file_time, NeXus_version */ nxprintattr(f, "creator", "%s generated with " MCCODE_STRING, mcinstrument_name); /* count the number of existing NXentry and create the next one */ NXgetgroupinfo(f, &count, name, class); sprintf(entry0, "entry%i", count+1); /* create the main NXentry (mandatory in NeXus) */ if (NXmakegroup(f, entry0, "NXentry") == NX_OK) if (NXopengroup(f, entry0, "NXentry") == NX_OK) { nxprintf(nxhandle, "program_name", MCCODE_STRING); nxprintf(f, "start_time", ctime(&t)); nxprintf(f, "title", "%s%s%s simulation generated by instrument %s", mcdirname && strlen(mcdirname) ? mcdirname : ".", MC_PATHSEP_S, mcsiminfo_name, mcinstrument_name); nxprintattr(f, "program_name", MCCODE_STRING); nxprintattr(f, "instrument", mcinstrument_name); nxprintattr(f, "simulation", "%s%s%s", mcdirname && strlen(mcdirname) ? mcdirname : ".", MC_PATHSEP_S, mcsiminfo_name); /* write NeXus instrument group */ if (NXmakegroup(f, "instrument", "NXinstrument") == NX_OK) if (NXopengroup(f, "instrument", "NXinstrument") == NX_OK) { int i; char *string=NULL; /* write NeXus parameters(types) data =================================== */ string = (char*)malloc(CHAR_BUF_LENGTH); if (string) { strcpy(string, ""); for(i = 0; i < mcnumipar; i++) { char ThisParam[CHAR_BUF_LENGTH]; snprintf(ThisParam, CHAR_BUF_LENGTH, " %s(%s)", mcinputtable[i].name, (*mcinputtypes[mcinputtable[i].type].parminfo) (mcinputtable[i].name)); if (strlen(string) + strlen(ThisParam) < CHAR_BUF_LENGTH) strcat(string, ThisParam); } nxprintattr(f, "Parameters", string); free(string); } nxprintattr(f, "name", mcinstrument_name); nxprintf (f, "name", mcinstrument_name); nxprintattr(f, "Source", mcinstrument_source); nxprintattr(f, "Trace_enabled", mctraceenabled ? "yes" : "no"); nxprintattr(f, "Default_main", mcdefaultmain ? "yes" : "no"); nxprintattr(f, "Embedded_runtime", #ifdef MC_EMBEDDED_RUNTIME "yes" #else "no" #endif ); /* add instrument source code when available */ buffer = mcinfo_readfile(mcinstrument_source); if (buffer && strlen(buffer)) { long length=strlen(buffer); nxprintf (f, "description", buffer); NXopendata(f,"description"); nxprintattr(f, "file_name", mcinstrument_source); nxprintattr(f, "file_size", "%li", length); nxprintattr(f, "MCCODE_STRING", MCCODE_STRING); NXclosedata(f); nxprintf (f,"instrument_source", "%s " MCCODE_NAME " " MCCODE_PARTICLE " Monte Carlo simulation", mcinstrument_name); free(buffer); } else nxprintf (f, "description", "File %s not found (instrument description %s is missing)", mcinstrument_source, mcinstrument_name); /* add Mantid/IDF.xml when available */ char *IDFfile=NULL; IDFfile = (char*)malloc(CHAR_BUF_LENGTH); sprintf(IDFfile,"%s%s",mcinstrument_source,".xml"); buffer = mcinfo_readfile(IDFfile); if (buffer && strlen(buffer)) { NXmakegroup (nxhandle, "instrument_xml", "NXnote"); NXopengroup (nxhandle, "instrument_xml", "NXnote"); nxprintf(f, "data", buffer); nxprintf(f, "description", "IDF.xml file found with instrument %s", mcinstrument_source); nxprintf(f, "type", "text/xml"); NXclosegroup(f); /* instrument_xml */ free(buffer); } free(IDFfile); NXclosegroup(f); /* instrument */ } /* NXinstrument */ /* write NeXus simulation group */ if (NXmakegroup(f, "simulation", "NXnote") == NX_OK) if (NXopengroup(f, "simulation", "NXnote") == NX_OK) { nxprintattr(f, "name", "%s%s%s", mcdirname && strlen(mcdirname) ? mcdirname : ".", MC_PATHSEP_S, mcsiminfo_name); nxprintf (f, "name", "%s", mcsiminfo_name); nxprintattr(f, "Format", mcformat && strlen(mcformat) ? mcformat : MCCODE_NAME); nxprintattr(f, "URL", "http://www.mccode.org"); nxprintattr(f, "program", MCCODE_STRING); nxprintattr(f, "Instrument",mcinstrument_source); nxprintattr(f, "Trace", mcdotrace ? "yes" : "no"); nxprintattr(f, "Gravitation",mcgravitation ? "yes" : "no"); nxprintattr(f, "Seed", "%li", mcseed); nxprintattr(f, "Directory", mcdirname); #ifdef USE_MPI if (mpi_node_count > 1) nxprintf(f, "Nodes", "%i", mpi_node_count); #endif /* output parameter string ================================================ */ if (NXmakegroup(f, "Param", "NXparameters") == NX_OK) if (NXopengroup(f, "Param", "NXparameters") == NX_OK) { int i; char string[CHAR_BUF_LENGTH]; for(i = 0; i < mcnumipar; i++) { if (mcget_run_num() || (mcinputtable[i].val && strlen(mcinputtable[i].val))) { if (mcinputtable[i].par == NULL) strncpy(string, (mcinputtable[i].val ? mcinputtable[i].val : ""), CHAR_BUF_LENGTH); else (*mcinputtypes[mcinputtable[i].type].printer)(string, mcinputtable[i].par); nxprintf(f, mcinputtable[i].name, "%s", string); nxprintattr(f, mcinputtable[i].name, string); } } NXclosegroup(f); /* Param */ } /* NXparameters */ NXclosegroup(f); /* simulation */ } /* NXsimulation */ /* create a group to hold all monitors */ NXmakegroup(f, "data", "NXdetector"); /* leave the NXentry opened (closed at exit) */ } /* NXentry */ } /* mcinfo_out_nexus */ /******************************************************************************* * mcdatainfo_out_nexus: output detector header * mcdatainfo_out_nexus(detector) create group and write info to NeXus data file * open data:NXdetector then filename:NXdata and write headers/attributes * requires: NXentry to be opened *******************************************************************************/ static void mcdatainfo_out_nexus(NXhandle f, MCDETECTOR detector) { char data_name[CHAR_BUF_LENGTH]; if (!f || !detector.m || mcdisable_output_files) return; strcpy_valid(data_name, detector.filename && strlen(detector.filename) ? detector.filename : detector.component); /* the NXdetector group has been created in mcinfo_out_nexus (mcsiminfo_init) */ if (NXopengroup(f, "data", "NXdetector") == NX_OK) { /* create and open the data group */ /* this may fail when appending to list -> ignore/skip */ NXMDisableErrorReporting(); /* unactivate NeXus error messages, as creation may fail */ if (NXmakegroup(f, data_name, "NXdata") == NX_OK) if (NXopengroup(f, data_name, "NXdata") == NX_OK) { /* output metadata (as attributes) ======================================== */ nxprintattr(f, "Date", detector.date); nxprintattr(f, "type", detector.type); nxprintattr(f, "Source", detector.instrument); nxprintattr(f, "component", detector.component); nxprintattr(f, "position", detector.position); nxprintattr(f, "title", detector.title); nxprintattr(f, !mcget_run_num() || mcget_run_num() >= mcget_ncount() ? "Ncount" : "ratio", detector.ncount); if (strlen(detector.filename)) { nxprintattr(f, "filename", detector.filename); } nxprintattr(f, "statistics", detector.statistics); nxprintattr(f, "signal", detector.signal); nxprintattr(f, "values", detector.values); if (detector.rank >= 1) { nxprintattr(f, "xvar", detector.xvar); nxprintattr(f, "yvar", detector.yvar); nxprintattr(f, "xlabel", detector.xlabel); nxprintattr(f, "ylabel", detector.ylabel); if (detector.rank > 1) { nxprintattr(f, "zvar", detector.zvar); nxprintattr(f, "zlabel", detector.zlabel); } } nxprintattr(f, abs(detector.rank)==1 ? "xlimits" : "xylimits", detector.limits); nxprintattr(f, "variables", strcasestr(detector.format, "list") ? detector.ylabel : detector.variables); nxprintf(f, "distance", detector.position); nxprintf(f, "acquisition_mode", strcasestr(detector.format, "list") ? "event" : "summed"); NXclosegroup(f); } /* NXdata (filename) */ NXMEnableErrorReporting(); /* re-enable NeXus error messages */ NXclosegroup(f); } /* NXdetector (data) */ } /* mcdatainfo_out_nexus */ /******************************************************************************* * mcdetector_out_axis_nexus: write detector axis into current NXdata * requires: NXdata to be opened *******************************************************************************/ int mcdetector_out_axis_nexus(NXhandle f, char *label, char *var, int rank, long length, double min, double max) { if (!f || length <= 1 || mcdisable_output_files || max == min) return(NX_OK); else { double axis[length]; char valid[CHAR_BUF_LENGTH]; int dim=(int)length; int i; int nprimary=1; /* create an axis from [min:max] */ for(i = 0; i < length; i++) axis[i] = min+(max-min)*(i+0.5)/length; /* create the data set */ strcpy_valid(valid, label); NXcompmakedata(f, valid, NX_FLOAT64, 1, &dim, NX_COMPRESSION, &dim); /* open it */ if (NXopendata(f, valid) != NX_OK) { fprintf(stderr, "Warning: could not open axis rank %i '%s' (NeXus)\n", rank, valid); return(NX_ERROR); } /* put the axis and its attributes */ NXputdata (f, axis); nxprintattr(f, "long_name", label); nxprintattr(f, "short_name", var); NXputattr (f, "axis", &rank, 1, NX_INT32); nxprintattr(f, "units", var); NXputattr (f, "primary", &nprimary, 1, NX_INT32); NXclosedata(f); return(NX_OK); } } /* mcdetector_out_axis_nexus */ /******************************************************************************* * mcdetector_out_array_nexus: write detector array into current NXdata (1D,2D) * requires: NXdata to be opened *******************************************************************************/ int mcdetector_out_array_nexus(NXhandle f, char *part, double *data, MCDETECTOR detector) { int dims[3]={detector.m,detector.n,detector.p}; /* number of elements to write */ int fulldims[3]={detector.m,detector.n,detector.p}; int signal=1; int exists=0; int current_dims[3]={0,0,0}; int ret=NX_OK; if (!f || !data || !detector.m || mcdisable_output_files) return(NX_OK); /* when this is a list, we set 1st dimension to NX_UNLIMITED for creation */ if (strcasestr(detector.format, "list")) fulldims[0] = NX_UNLIMITED; /* create the data set in NXdata group */ NXMDisableErrorReporting(); /* unactivate NeXus error messages, as creation may fail */ ret = NXcompmakedata(f, part, NX_FLOAT64, detector.rank, fulldims, NX_COMPRESSION, dims); if (ret != NX_OK) { /* failed: data set already exists */ int datatype=0; int rank=0; exists=1; /* inquire current size of data set (nb of events stored) */ NXopendata(f, part); NXgetinfo(f, &rank, current_dims, &datatype); NXclosedata(f); } NXMEnableErrorReporting(); /* re-enable NeXus error messages */ /* open the data set */ if (NXopendata(f, part) == NX_ERROR) { fprintf(stderr, "Warning: could not open DataSet %s '%s' (NeXus)\n", part, detector.title); return(NX_ERROR); } if (strcasestr(detector.format, "list")) { current_dims[1] = current_dims[2] = 0; /* set starting location for writing slab */ NXputslab(f, data, current_dims, dims); if (!exists) printf("Events: \"%s\"\n", strlen(detector.filename) ? detector.filename : detector.component); else printf("Append: \"%s\"\n", strlen(detector.filename) ? detector.filename : detector.component); } else { NXputdata (f, data); } if (strstr(part,"data") || strstr(part, "events")) { NXputattr(f, "signal", &signal, 1, NX_INT32); nxprintattr(f, "short_name", detector.filename && strlen(detector.filename) ? detector.filename : detector.component); } nxprintattr(f, "long_name", "%s '%s'", part, detector.title); NXclosedata(f); return(NX_OK); } /* mcdetector_out_array_nexus */ /******************************************************************************* * mcdetector_out_data_nexus: write detector axes+data into current NXdata * The data:NXdetector is opened, then filename:NXdata * requires: NXentry to be opened *******************************************************************************/ int mcdetector_out_data_nexus(NXhandle f, MCDETECTOR detector) { char data_name[CHAR_BUF_LENGTH]; if (!f || !detector.m || mcdisable_output_files) return(NX_OK); strcpy_valid(data_name, detector.filename && strlen(detector.filename) ? detector.filename : detector.component); /* the NXdetector group has been created in mcinfo_out_nexus (mcsiminfo_init) */ if (NXopengroup(f, "data", "NXdetector") == NX_OK) { /* the NXdata group has been created in mcdatainfo_out_nexus */ if (NXopengroup(f, data_name, "NXdata") == NX_OK) { /* write axes, for histogram data sets, not for lists */ if (!strcasestr(detector.format, "list")) { mcdetector_out_axis_nexus(f, detector.xlabel, detector.xvar, 1, detector.m, detector.xmin, detector.xmax); mcdetector_out_axis_nexus(f, detector.ylabel, detector.yvar, 2, detector.n, detector.ymin, detector.ymax); mcdetector_out_axis_nexus(f, detector.zlabel, detector.zvar, 3, detector.p, detector.zmin, detector.zmax); } /* !list */ /* write the actual data (appended if already exists) */ if (!strcasestr(detector.format, "list")) { mcdetector_out_array_nexus(f, "data", detector.p1, detector); mcdetector_out_array_nexus(f, "errors", detector.p2, detector); mcdetector_out_array_nexus(f, "ncount", detector.p0, detector); } else mcdetector_out_array_nexus( f, "events", detector.p1, detector); NXclosegroup(f); } /* NXdata */ NXclosegroup(f); } /* NXdetector */ return(NX_OK); } /* mcdetector_out_array_nexus */ #ifdef USE_MPI /******************************************************************************* * mcdetector_out_list_slaves: slaves send their list data to master which writes * requires: NXentry to be opened * WARNING: this method has a flaw: it requires all nodes to flush the lists * the same number of times. In case one node is just below the buffer size * when finishing (e.g. monitor_nd), it may not trigger save but others may. * Then the number of recv/send is not constant along nodes, and simulation stalls. *******************************************************************************/ MCDETECTOR mcdetector_out_list_slaves(MCDETECTOR detector) { int node_i=0; MPI_MASTER( printf("\n** MPI master gathering slave node list data ** \n"); ); if (mpi_node_rank != mpi_node_root) { /* MPI slave: slaves send their data to master: 2 MPI_Send calls */ /* m, n, p must be sent first, since all slaves do not have the same number of events */ int mnp[3]={detector.m,detector.n,detector.p}; if (mc_MPI_Send(mnp, 3, MPI_INT, mpi_node_root)!= MPI_SUCCESS) fprintf(stderr, "Warning: proc %i to master: MPI_Send mnp list error (mcdetector_out_list_slaves)\n", mpi_node_rank); if (!detector.p1 || mc_MPI_Send(detector.p1, mnp[0]*mnp[1]*mnp[2], MPI_DOUBLE, mpi_node_root) != MPI_SUCCESS) fprintf(stderr, "Warning: proc %i to master: MPI_Send p1 list error: mnp=%i (mcdetector_out_list_slaves)\n", mpi_node_rank, abs(mnp[0]*mnp[1]*mnp[2])); /* slaves are done: sent mnp and p1 */ return (detector); } /* end slaves */ /* MPI master: receive data from slaves sequentially: 2 MPI_Recv calls */ if (mpi_node_rank == mpi_node_root) { for(node_i=0; node_i 1) { mcdetector_out_list_slaves(detector); } #endif /* USE_MPI */ return(detector); } /* mcdetector_out_2D_nexus */ #endif /* USE_NEXUS*/ /* ========================================================================== */ /* Main input functions */ /* DETECTOR_OUT_xD function calls -> ascii or NeXus */ /* ========================================================================== */ /******************************************************************************* * mcsiminfo_init: open SIM and write header *******************************************************************************/ FILE *mcsiminfo_init(FILE *f) { int exists=0; int index; /* check format */ if (!mcformat || !strlen(mcformat) || !strcasecmp(mcformat, "MCSTAS") || !strcasecmp(mcformat, "MCXTRACE") || !strcasecmp(mcformat, "PGPLOT") || !strcasecmp(mcformat, "GNUPLOT") || !strcasecmp(mcformat, "MCCODE") || !strcasecmp(mcformat, "MATLAB")) { mcformat="McCode"; #ifdef USE_NEXUS } else if (strcasestr(mcformat, "NeXus")) { /* Do nothing */ #endif } else { fprintf(stderr, "Warning: You have requested the output format %s which is unsupported by this binary. Resetting to standard %s format.\n",mcformat ,"McCode"); mcformat="McCode"; } /* open the SIM file if not defined yet */ if (mcsiminfo_file || mcdisable_output_files) return (mcsiminfo_file); #ifdef USE_NEXUS /* only master writes NeXus header: calls NXopen(nxhandle) */ if (mcformat && strcasestr(mcformat, "NeXus")) { MPI_MASTER( mcsiminfo_file = mcnew_file(mcsiminfo_name, "h5", &exists); if(!mcsiminfo_file) fprintf(stderr, "Warning: could not open simulation description file '%s'\n", mcsiminfo_name); else mcinfo_out_nexus(nxhandle); ); return(mcsiminfo_file); /* points to nxhandle */ } #endif /* write main description file (only MASTER) */ MPI_MASTER( mcsiminfo_file = mcnew_file(mcsiminfo_name, "sim", &exists); if(!mcsiminfo_file) fprintf(stderr, "Warning: could not open simulation description file '%s'\n", mcsiminfo_name); else { /* write SIM header */ time_t t=time(NULL); mcsiminfo_out("%s simulation description file for %s.\n", MCCODE_NAME, mcinstrument_name); mcsiminfo_out("Date: %s", ctime(&t)); /* includes \n */ mcsiminfo_out("Program: %s\n\n", MCCODE_STRING); mcsiminfo_out("begin instrument: %s\n", mcinstrument_name); mcinfo_out( " ", mcsiminfo_file); mcsiminfo_out("end instrument\n"); mcsiminfo_out("\nbegin simulation: %s\n", mcdirname); mcruninfo_out(" ", mcsiminfo_file); mcsiminfo_out("end simulation\n"); } return (mcsiminfo_file); ); /* MPI_MASTER */ } /* mcsiminfo_init */ /******************************************************************************* * mcsiminfo_close: close SIM *******************************************************************************/ void mcsiminfo_close() { MPI_MASTER( if(mcsiminfo_file && !mcdisable_output_files) { #ifdef USE_NEXUS if (mcformat && strcasestr(mcformat, "NeXus")) { time_t t=time(NULL); nxprintf(nxhandle, "end_time", ctime(&t)); nxprintf(nxhandle, "duration", "%li", (long)t-mcstartdate); NXclosegroup(nxhandle); /* NXentry */ NXclose(&nxhandle); } else #endif fclose(mcsiminfo_file); ); mcsiminfo_file = NULL; } } /* mcsiminfo_close */ /******************************************************************************* * mcdetector_out_0D: wrapper for 0D (single value). * Output single detector/monitor data (p0, p1, p2). * Title is t, component name is c. *******************************************************************************/ MCDETECTOR mcdetector_out_0D(char *t, double p0, double p1, double p2, char *c, Coords posa) { /* import and perform basic detector analysis (and handle MPI reduce) */ MCDETECTOR detector = mcdetector_import(mcformat, c, (t ? t : MCCODE_STRING " data"), 1, 1, 1, "I", "", "", "I", "", "", 0, 0, 0, 0, 0, 0, "", &p0, &p1, &p2, posa); /* write Detector: line */ #ifdef USE_NEXUS if (strcasestr(detector.format, "NeXus")) return(mcdetector_out_0D_nexus(detector)); else #endif return(mcdetector_out_0D_ascii(detector)); } /* mcdetector_out_0D */ /******************************************************************************* * mcdetector_out_1D: wrapper for 1D. * Output 1d detector data (p0, p1, p2) for n bins linearly * distributed across the range x1..x2 (x1 is lower limit of first * bin, x2 is upper limit of last bin). Title is t, axis labels are xl * and yl. File name is f, component name is c. *******************************************************************************/ MCDETECTOR mcdetector_out_1D(char *t, char *xl, char *yl, char *xvar, double x1, double x2, long n, double *p0, double *p1, double *p2, char *f, char *c, Coords posa) { /* import and perform basic detector analysis (and handle MPI_Reduce) */ MCDETECTOR detector = mcdetector_import(mcformat, c, (t ? t : MCCODE_STRING " 1D data"), n, 1, 1, xl, yl, (n > 1 ? "Signal per bin" : " Signal"), xvar, "(I,I_err)", "I", x1, x2, 0, 0, 0, 0, f, p0, p1, p2, posa); /* write Detector: line */ if (!detector.p1 || !detector.m) return(detector); #ifdef USE_NEXUS if (strcasestr(detector.format, "NeXus")) return(mcdetector_out_1D_nexus(detector)); else #endif return(mcdetector_out_1D_ascii(detector)); } /* mcdetector_out_1D */ /******************************************************************************* * mcdetector_out_2D: wrapper for 2D. * special case for list: master creates file first, then slaves append their blocks without header *******************************************************************************/ MCDETECTOR mcdetector_out_2D(char *t, char *xl, char *yl, double x1, double x2, double y1, double y2, long m, long n, double *p0, double *p1, double *p2, char *f, char *c, Coords posa) { char xvar[CHAR_BUF_LENGTH]; char yvar[CHAR_BUF_LENGTH]; /* create short axes labels */ if (xl && strlen(xl)) { strncpy(xvar, xl, CHAR_BUF_LENGTH); xvar[strcspn(xvar,"\n\r ")]='\0'; } else strcpy(xvar, "x"); if (yl && strlen(yl)) { strncpy(yvar, yl, CHAR_BUF_LENGTH); yvar[strcspn(yvar,"\n\r ")]='\0'; } else strcpy(yvar, "y"); MCDETECTOR detector; /* import and perform basic detector analysis (and handle MPI_Reduce) */ if (labs(m) == 1) {/* n>1 on Y, m==1 on X: 1D, no X axis*/ detector = mcdetector_import(mcformat, c, (t ? t : MCCODE_STRING " 1D data"), n, 1, 1, yl, "", "Signal per bin", yvar, "(I,Ierr)", "I", y1, y2, x1, x2, 0, 0, f, p0, p1, p2, posa); /* write Detector: line */ } else if (labs(n)==1) {/* m>1 on X, n==1 on Y: 1D, no Y axis*/ detector = mcdetector_import(mcformat, c, (t ? t : MCCODE_STRING " 1D data"), m, 1, 1, xl, "", "Signal per bin", xvar, "(I,Ierr)", "I", x1, x2, y1, y2, 0, 0, f, p0, p1, p2, posa); /* write Detector: line */ }else { detector = mcdetector_import(mcformat, c, (t ? t : MCCODE_STRING " 2D data"), m, n, 1, xl, yl, "Signal per bin", xvar, yvar, "I", x1, x2, y1, y2, 0, 0, f, p0, p1, p2, posa); /* write Detector: line */ } if (!detector.p1 || !detector.m) return(detector); #ifdef USE_NEXUS if (strcasestr(detector.format, "NeXus")) return(mcdetector_out_2D_nexus(detector)); else #endif return(mcdetector_out_2D_ascii(detector)); } /* mcdetector_out_2D */ /******************************************************************************* * mcdetector_out_list: wrapper for list output (calls out_2D with mcformat+"list"). * m=number of events, n=size of each event *******************************************************************************/ MCDETECTOR mcdetector_out_list(char *t, char *xl, char *yl, long m, long n, double *p1, char *f, char *c, Coords posa) { char format_new[CHAR_BUF_LENGTH]; char *format_org; MCDETECTOR detector; format_org = mcformat; strcpy(format_new, mcformat); strcat(format_new, " list"); mcformat = format_new; detector = mcdetector_out_2D(t, xl, yl, 1,labs(m),1,labs(n), m,n, NULL, p1, NULL, f, c, posa); mcformat = format_org; return(detector); } /******************************************************************************* * mcuse_dir: set data/sim storage directory and create it, * or exit with error if exists ******************************************************************************/ static void mcuse_dir(char *dir) { if (!dir || !strlen(dir)) return; #ifdef MC_PORTABLE fprintf(stderr, "Error: " "Directory output cannot be used with portable simulation (mcuse_dir)\n"); exit(1); #else /* !MC_PORTABLE */ /* handle file://directory URL type */ if (strncmp(dir, "file://", strlen("file://"))) mcdirname = dir; else mcdirname = dir+strlen("file://"); MPI_MASTER( if(mkdir(mcdirname, 0777)) { #ifndef DANSE fprintf(stderr, "Error: unable to create directory '%s' (mcuse_dir)\n", dir); fprintf(stderr, "(Maybe the directory already exists?)\n"); #endif #ifdef USE_MPI MPI_Abort(MPI_COMM_WORLD, -1); #endif exit(-1); } ); /* MPI_MASTER */ /* remove trailing PATHSEP (if any) */ while (strlen(mcdirname) && mcdirname[strlen(mcdirname) - 1] == MC_PATHSEP_C) mcdirname[strlen(mcdirname) - 1]='\0'; #endif /* !MC_PORTABLE */ } /* mcuse_dir */ /******************************************************************************* * mcinfo: display instrument simulation info to stdout and exit *******************************************************************************/ static void mcinfo(void) { fprintf(stdout, "begin instrument: %s\n", mcinstrument_name); mcinfo_out(" ", stdout); fprintf(stdout, "end instrument\n"); fprintf(stdout, "begin simulation: %s\n", mcdirname ? mcdirname : "."); mcruninfo_out_backend(" ", stdout,1); fprintf(stdout, "end simulation\n"); exit(0); /* includes MPI_Finalize in MPI mode */ } /* mcinfo */ #endif /* ndef MCCODE_R_IO_C */ /* end of the I/O section =================================================== */ /******************************************************************************* * mcset_ncount: set total number of rays to generate *******************************************************************************/ void mcset_ncount(unsigned long long int count) { mcncount = count; } /* mcget_ncount: get total number of rays to generate */ unsigned long long int mcget_ncount(void) { return mcncount; } /* mcget_run_num: get curent number of rays in TRACE */ unsigned long long int mcget_run_num(void) { return mcrun_num; } /* mcsetn_arg: get ncount from a string argument */ static void mcsetn_arg(char *arg) { mcset_ncount((long long int) strtod(arg, NULL)); } /* mcsetseed: set the random generator seed from a string argument */ static void mcsetseed(char *arg) { mcseed = atol(arg); if(mcseed) { srandom(mcseed); } else { fprintf(stderr, "Error: seed must not be zero (mcsetseed)\n"); exit(1); } } /* Following part is only embedded when not redundent with mccode-r.h ========= */ #ifndef MCCODE_H /* SECTION: MCDISPLAY support. =============================================== */ /******************************************************************************* * Just output MCDISPLAY keywords to be caught by an external plotter client. *******************************************************************************/ void mcdis_magnify(char *what){ // Do nothing here, better use interactive zoom from the tools } void mcdis_line(double x1, double y1, double z1, double x2, double y2, double z2){ printf("MCDISPLAY: multiline(2,%g,%g,%g,%g,%g,%g)\n", x1,y1,z1,x2,y2,z2); } void mcdis_dashed_line(double x1, double y1, double z1, double x2, double y2, double z2, int n){ int i; const double dx = (x2-x1)/(2*n+1); const double dy = (y2-y1)/(2*n+1); const double dz = (z2-z1)/(2*n+1); for(i = 0; i < n+1; i++) mcdis_line(x1 + 2*i*dx, y1 + 2*i*dy, z1 + 2*i*dz, x1 + (2*i+1)*dx, y1 + (2*i+1)*dy, z1 + (2*i+1)*dz); } void mcdis_multiline(int count, ...){ va_list ap; double x,y,z; printf("MCDISPLAY: multiline(%d", count); va_start(ap, count); while(count--) { x = va_arg(ap, double); y = va_arg(ap, double); z = va_arg(ap, double); printf(",%g,%g,%g", x, y, z); } va_end(ap); printf(")\n"); } void mcdis_rectangle(char* plane, double x, double y, double z, double width, double height){ /* draws a rectangle in the plane */ /* x is ALWAYS width and y is ALWAYS height */ if (strcmp("xy", plane)==0) { mcdis_multiline(5, x - width/2, y - height/2, z, x + width/2, y - height/2, z, x + width/2, y + height/2, z, x - width/2, y + height/2, z, x - width/2, y - height/2, z); } else if (strcmp("xz", plane)==0) { mcdis_multiline(5, x - width/2, y, z - height/2, x + width/2, y, z - height/2, x + width/2, y, z + height/2, x - width/2, y, z + height/2, x - width/2, y, z - height/2); } else if (strcmp("yz", plane)==0) { mcdis_multiline(5, x, y - height/2, z - width/2, x, y - height/2, z + width/2, x, y + height/2, z + width/2, x, y + height/2, z - width/2, x, y - height/2, z - width/2); } else { fprintf(stderr, "Error: Definition of plane %s unknown\n", plane); exit(1); } } /* draws a box with center at (x, y, z) and width (deltax), height (deltay), length (deltaz) */ void mcdis_box(double x, double y, double z, double width, double height, double length){ mcdis_rectangle("xy", x, y, z-length/2, width, height); mcdis_rectangle("xy", x, y, z+length/2, width, height); mcdis_line(x-width/2, y-height/2, z-length/2, x-width/2, y-height/2, z+length/2); mcdis_line(x-width/2, y+height/2, z-length/2, x-width/2, y+height/2, z+length/2); mcdis_line(x+width/2, y-height/2, z-length/2, x+width/2, y-height/2, z+length/2); mcdis_line(x+width/2, y+height/2, z-length/2, x+width/2, y+height/2, z+length/2); } void mcdis_circle(char *plane, double x, double y, double z, double r){ printf("MCDISPLAY: circle('%s',%g,%g,%g,%g)\n", plane, x, y, z, r); } /* Draws a circle with center (x,y,z), radius (r), and in the plane * with normal (nx,ny,nz)*/ void mcdis_Circle(double x, double y, double z, double r, double nx, double ny, double nz){ int i; if(nx==0 && ny && nz==0){ for (i=0;i<24; i++){ mcdis_line(x+r*sin(i*2*M_PI/24),y,z+r*cos(i*2*M_PI/24), x+r*sin((i+1)*2*M_PI/24),y,z+r*cos((i+1)*2*M_PI/24)); } }else{ double mx,my,mz; /*generate perpendicular vector using (nx,ny,nz) and (0,1,0)*/ vec_prod(mx,my,mz, 0,1,0, nx,ny,nz); NORM(mx,my,mz); /*draw circle*/ for (i=0;i<24; i++){ double ux,uy,uz; double wx,wy,wz; rotate(ux,uy,uz, mx,my,mz, i*2*M_PI/24, nx,ny,nz); rotate(wx,wy,wz, mx,my,mz, (i+1)*2*M_PI/24, nx,ny,nz); mcdis_line(x+ux*r,y+uy*r,z+uz*r, x+wx*r,y+wy*r,z+wz*r); } } } /* Draws a cylinder with center at (x,y,z) with extent (r,height). * The cylinder axis is along the vector nx,ny,nz. * N determines how many vertical lines are drawn.*/ void mcdis_cylinder( double x, double y, double z, double r, double height, int N, double nx, double ny, double nz){ int i; /*no lines make little sense - so trigger the default*/ if(N<=0) N=5; NORM(nx,ny,nz); double h_2=height/2.0; mcdis_Circle(x+nx*h_2,y+ny*h_2,z+nz*h_2,r,nx,ny,nz); mcdis_Circle(x-nx*h_2,y-ny*h_2,z-nz*h_2,r,nx,ny,nz); double mx,my,mz; /*generate perpendicular vector using (nx,ny,nz) and (0,1,0)*/ if(nx==0 && ny && nz==0){ mx=my=0;mz=1; }else{ vec_prod(mx,my,mz, 0,1,0, nx,ny,nz); NORM(mx,my,mz); } /*draw circle*/ for (i=0; i<24; i++){ double ux,uy,uz; rotate(ux,uy,uz, mx,my,mz, i*2*M_PI/24, nx,ny,nz); mcdis_line(x+nx*h_2+ux*r, y+ny*h_2+uy*r, z+nz*h_2+uz*r, x-nx*h_2+ux*r, y-ny*h_2+uy*r, z-nz*h_2+uz*r); } } /* draws a sphere with center at (x,y,z) with extent (r) * The sphere is drawn using N longitudes and N latitudes.*/ void mcdis_sphere(double x, double y, double z, double r, int N){ double nx,ny,nz; int i; /*no lines make little sense - so trigger the default*/ if(N<=0) N=5; nx=0;ny=0;nz=1; mcdis_Circle(x,y,z,r,nx,ny,nz); for (i=1;ix /= temp; c->y /= temp; c->z /= temp; } /******************************************************************************* * The Rotation type implements a rotation transformation of a coordinate * system in the form of a double[3][3] matrix. * * Contrary to the Coords type in coords.c, rotations are passed by * reference. Functions that yield new rotations do so by writing to an * explicit result parameter; rotations are not returned from functions. The * reason for this is that arrays cannot by returned from functions (though * structures can; thus an alternative would have been to wrap the * double[3][3] array up in a struct). Such are the ways of C programming. * * A rotation represents the tranformation of the coordinates of a vector when * changing between coordinate systems that are rotated with respect to each * other. For example, suppose that coordinate system Q is rotated 45 degrees * around the Z axis with respect to coordinate system P. Let T be the * rotation transformation representing a 45 degree rotation around Z. Then to * get the coordinates of a vector r in system Q, apply T to the coordinates * of r in P. If r=(1,0,0) in P, it will be (sqrt(1/2),-sqrt(1/2),0) in * Q. Thus we should be careful when interpreting the sign of rotation angles: * they represent the rotation of the coordinate systems, not of the * coordinates (which has opposite sign). *******************************************************************************/ /******************************************************************************* * rot_set_rotation: Get transformation for rotation first phx around x axis, * then phy around y, then phz around z. *******************************************************************************/ void rot_set_rotation(Rotation t, double phx, double phy, double phz) { if ((phx == 0) && (phy == 0) && (phz == 0)) { t[0][0] = 1.0; t[0][1] = 0.0; t[0][2] = 0.0; t[1][0] = 0.0; t[1][1] = 1.0; t[1][2] = 0.0; t[2][0] = 0.0; t[2][1] = 0.0; t[2][2] = 1.0; } else { double cx = cos(phx); double sx = sin(phx); double cy = cos(phy); double sy = sin(phy); double cz = cos(phz); double sz = sin(phz); t[0][0] = cy*cz; t[0][1] = sx*sy*cz + cx*sz; t[0][2] = sx*sz - cx*sy*cz; t[1][0] = -cy*sz; t[1][1] = cx*cz - sx*sy*sz; t[1][2] = sx*cz + cx*sy*sz; t[2][0] = sy; t[2][1] = -sx*cy; t[2][2] = cx*cy; } } /******************************************************************************* * rot_test_identity: Test if rotation is identity *******************************************************************************/ int rot_test_identity(Rotation t) { return (t[0][0] + t[1][1] + t[2][2] == 3); } /******************************************************************************* * rot_mul: Matrix multiplication of transformations (this corresponds to * combining transformations). After rot_mul(T1, T2, T3), doing T3 is * equal to doing first T2, then T1. * Note that T3 must not alias (use the same array as) T1 or T2. *******************************************************************************/ void rot_mul(Rotation t1, Rotation t2, Rotation t3) { if (rot_test_identity(t1)) { rot_copy(t3, t2); } else if (rot_test_identity(t2)) { rot_copy(t3, t1); } else { int i,j; for(i = 0; i < 3; i++) for(j = 0; j < 3; j++) t3[i][j] = t1[i][0]*t2[0][j] + t1[i][1]*t2[1][j] + t1[i][2]*t2[2][j]; } } /******************************************************************************* * rot_copy: Copy a rotation transformation (arrays cannot be assigned in C). *******************************************************************************/ void rot_copy(Rotation dest, Rotation src) { int i,j; for(i = 0; i < 3; i++) for(j = 0; j < 3; j++) dest[i][j] = src[i][j]; } /******************************************************************************* * rot_transpose: Matrix transposition, which is inversion for Rotation matrices *******************************************************************************/ void rot_transpose(Rotation src, Rotation dst) { dst[0][0] = src[0][0]; dst[0][1] = src[1][0]; dst[0][2] = src[2][0]; dst[1][0] = src[0][1]; dst[1][1] = src[1][1]; dst[1][2] = src[2][1]; dst[2][0] = src[0][2]; dst[2][1] = src[1][2]; dst[2][2] = src[2][2]; } /******************************************************************************* * rot_invert: Matrix inversion, in case a Rotatoin is used to represent a * a general non-orthonormal matrix. *******************************************************************************/ void rot_invert(Rotation t1, Rotation t2) { Rotation cofactors; int r,c; double det=0; for (r=0;r<3;r++){ for (c=0;c<3;c++){ /*this algorithm automatically takes care of the sign changes in computing cofactors*/ cofactors[r][c]=t1[(r+1) % 3][(c+1) % 3]*t1[(r+2) % 3][(c+2) % 3] - t1[(r+2) % 3][(c+1) % 3]*t1[(r+1) % 3][(c +2) % 3] ; } } det=t1[0][0]*cofactors[0][0] + t1[0][1]*cofactors[0][1] + t1[0][2]*cofactors[0][2]; if(det==0){ fprintf(stderr,"Warning: matrix not invertable\n"); } rot_transpose(cofactors,t2); /*the adjoint matrix should now be scaled by 1/det to get the inverse*/ for (r=0;r<3;r++){ for (c=0;c<3;c++){ t2[r][c]=t2[r][c]/det; } } } /******************************************************************************* * rot_apply: returns t*a *******************************************************************************/ Coords rot_apply(Rotation t, Coords a) { Coords b; if (rot_test_identity(t)) { return a; } else { b.x = t[0][0]*a.x + t[0][1]*a.y + t[0][2]*a.z; b.y = t[1][0]*a.x + t[1][1]*a.y + t[1][2]*a.z; b.z = t[2][0]*a.x + t[2][1]*a.y + t[2][2]*a.z; return b; } } /** * Pretty-printing of rotation matrices. */ void rot_print(Rotation rot) { printf("[ %4.2f %4.2f %4.2f ]\n", rot[0][0], rot[0][1], rot[0][2]); printf("[ %4.2f %4.2f %4.2f ]\n", rot[1][0], rot[1][1], rot[1][2]); printf("[ %4.2f %4.2f %4.2f ]\n\n", rot[2][0], rot[2][1], rot[2][2]); } /** * Vector product: used by vec_prod (mccode-r.h). Use coords_xp for Coords. */ mcstatic void vec_prod_func(double *x, double *y, double *z, double x1, double y1, double z1, double x2, double y2, double z2) { *x = (y1)*(z2) - (y2)*(z1); *y = (z1)*(x2) - (z2)*(x1); *z = (x1)*(y2) - (x2)*(y1); } /** * Scalar product: use coords_sp for Coords. */ mcstatic double scalar_prod( double x1, double y1, double z1, double x2, double y2, double z2) { return ((x1 * x2) + (y1 * y2) + (z1 * z2)); } /******************************************************************************* * mccoordschange: applies rotation to (x y z) and (vx vy vz) and Spin (sx,sy,sz) *******************************************************************************/ void mccoordschange(Coords a, Rotation t, double *x, double *y, double *z, double *vx, double *vy, double *vz, double *sx, double *sy, double *sz) { Coords b, c; b.x = *x; b.y = *y; b.z = *z; c = rot_apply(t, b); b = coords_add(c, a); *x = b.x; *y = b.y; *z = b.z; if ( (vz && vy && vx) && (*vz != 0.0 || *vx != 0.0 || *vy != 0.0) ) mccoordschange_polarisation(t, vx, vy, vz); if ( (sz && sy && sx) && (*sz != 0.0 || *sx != 0.0 || *sy != 0.0) ) mccoordschange_polarisation(t, sx, sy, sz); } /******************************************************************************* * mccoordschange_polarisation: applies rotation to vector (sx sy sz) *******************************************************************************/ void mccoordschange_polarisation(Rotation t, double *sx, double *sy, double *sz) { Coords b, c; b.x = *sx; b.y = *sy; b.z = *sz; c = rot_apply(t, b); *sx = c.x; *sy = c.y; *sz = c.z; } /* SECTION: vector math ==================================================== */ /* normal_vec_func: Compute normal vector to (x,y,z). */ mcstatic void normal_vec_func(double *nx, double *ny, double *nz, double x, double y, double z) { double ax = fabs(x); double ay = fabs(y); double az = fabs(z); double l; if(x == 0 && y == 0 && z == 0) { *nx = 0; *ny = 0; *nz = 0; return; } if(ax < ay) { if(ax < az) { /* Use X axis */ l = sqrt(z*z + y*y); *nx = 0; *ny = z/l; *nz = -y/l; return; } } else { if(ay < az) { /* Use Y axis */ l = sqrt(z*z + x*x); *nx = z/l; *ny = 0; *nz = -x/l; return; } } /* Use Z axis */ l = sqrt(y*y + x*x); *nx = y/l; *ny = -x/l; *nz = 0; } /* normal_vec */ /******************************************************************************* * solve_2nd_order: second order equation solve: A*t^2 + B*t + C = 0 * solve_2nd_order(&t1, NULL, A,B,C) * returns 0 if no solution was found, or set 't1' to the smallest positive * solution. * solve_2nd_order(&t1, &t2, A,B,C) * same as with &t2=NULL, but also returns the second solution. * EXAMPLE usage for intersection of a trajectory with a plane in gravitation * field (gx,gy,gz): * The neutron starts at point r=(x,y,z) with velocityv=(vx vy vz). The plane * has a normal vector n=(nx,ny,nz) and contains the point W=(wx,wy,wz). * The problem consists in solving the 2nd order equation: * 1/2.n.g.t^2 + n.v.t + n.(r-W) = 0 * so that A = 0.5 n.g; B = n.v; C = n.(r-W); * Without acceleration, t=-n.(r-W)/n.v ******************************************************************************/ int solve_2nd_order_old(double *t1, double *t2, double A, double B, double C) { int ret=0; if (!t1) return 0; *t1 = 0; if (t2) *t2=0; if (fabs(A) < 1E-10) /* approximate to linear equation: A ~ 0 */ { if (B) { *t1 = -C/B; ret=1; if (t2) *t2=*t1; } /* else no intersection: A=B=0 ret=0 */ } else { double D; D = B*B - 4*A*C; if (D >= 0) /* Delta > 0: two solutions */ { double sD, dt1, dt2; sD = sqrt(D); dt1 = (-B + sD)/2/A; dt2 = (-B - sD)/2/A; /* we identify very small values with zero */ if (fabs(dt1) < 1e-10) dt1=0.0; if (fabs(dt2) < 1e-10) dt2=0.0; /* now we choose the smallest positive solution */ if (dt1<=0.0 && dt2>0.0) ret=2; /* dt2 positive */ else if (dt2<=0.0 && dt1>0.0) ret=1; /* dt1 positive */ else if (dt1> 0.0 && dt2>0.0) { if (dt1 < dt2) ret=1; else ret=2; } /* all positive: min(dt1,dt2) */ /* else two solutions are negative. ret=-1 */ if (ret==1) { *t1 = dt1; if (t2) *t2=dt2; } else { *t1 = dt2; if (t2) *t2=dt1; } ret=2; /* found 2 solutions and t1 is the positive one */ } /* else Delta <0: no intersection. ret=0 */ } return(ret); } /* solve_2nd_order */ int solve_2nd_order(double *t0, double *t1, double A, double B, double C){ int retval=0; double sign=copysign(1.0,B); double dt0,dt1; dt0=0; dt1=0; *t0; if(t1){ *t1=0;} /*protect against rounding errors by locally equating DBL_EPSILON with 0*/ if (fabs(A)=0){ dt0=(-B - sign*sqrt(B*B-4*A*C))/(2*A); dt1=C/(A*dt0); retval=2; }else{ /*no real roots*/ retval=0; } } /*sort the solutions*/ if (retval==1){ /*put both solutions in t0 and t1*/ *t0=dt0; if(t1) *t1=dt1; }else{ /*we have two solutions*/ /*swap if both are positive and t1 smaller than t0 or t1 the only positive*/ int swap=0; if(dt1>0 && ( dt1) * * If height or width is zero, choose random direction in full 4PI, no target. * * Traditionally, this routine had the name randvec_target_rect - this is now a * a define (see mcstas-r.h) pointing here. If you use the old rouine, you are NOT * taking the local emmission coordinate into account. *******************************************************************************/ void randvec_target_rect_real(double *xo, double *yo, double *zo, double *solid_angle, double xi, double yi, double zi, double width, double height, Rotation A, double lx, double ly, double lz, int order) { double dx, dy, dist, dist_p, nx, ny, nz, mx, my, mz, n_norm, m_norm; double cos_theta; Coords tmp; Rotation Ainverse; rot_transpose(A, Ainverse); if(height == 0.0 || width == 0.0) { randvec_target_circle(xo, yo, zo, solid_angle, xi, yi, zi, 0); return; } else { /* Now choose point uniformly on rectangle within width x height */ dx = width*randpm1()/2.0; dy = height*randpm1()/2.0; /* Determine distance to target plane*/ dist = sqrt(xi*xi + yi*yi + zi*zi); /* Go to global coordinate system */ tmp = coords_set(xi, yi, zi); tmp = rot_apply(Ainverse, tmp); coords_get(tmp, &xi, &yi, &zi); /* Determine vector normal to trajectory axis (z) and gravity [0 1 0] */ vec_prod(nx, ny, nz, xi, yi, zi, 0, 1, 0); /* This now defines the x-axis, normalize: */ n_norm=sqrt(nx*nx + ny*ny + nz*nz); nx = nx/n_norm; ny = ny/n_norm; nz = nz/n_norm; /* Now, determine our y-axis (vertical in many cases...) */ vec_prod(mx, my, mz, xi, yi, zi, nx, ny, nz); m_norm=sqrt(mx*mx + my*my + mz*mz); mx = mx/m_norm; my = my/m_norm; mz = mz/m_norm; /* Our output, random vector can now be defined by linear combination: */ *xo = xi + dx * nx + dy * mx; *yo = yi + dx * ny + dy * my; *zo = zi + dx * nz + dy * mz; /* Go back to local coordinate system */ tmp = coords_set(*xo, *yo, *zo); tmp = rot_apply(A, tmp); coords_get(tmp, &*xo, &*yo, &*zo); /* Go back to local coordinate system */ tmp = coords_set(xi, yi, zi); tmp = rot_apply(A, tmp); coords_get(tmp, &xi, &yi, &zi); if (solid_angle) { /* Calculate vector from local point to remote random point */ lx = *xo - lx; ly = *yo - ly; lz = *zo - lz; dist_p = sqrt(lx*lx + ly*ly + lz*lz); /* Adjust the 'solid angle' */ /* 1/r^2 to the chosen point times cos(\theta) between the normal */ /* vector of the target rectangle and direction vector of the chosen point. */ cos_theta = (xi * lx + yi * ly + zi * lz) / (dist * dist_p); *solid_angle = width * height / (dist_p * dist_p); int counter; for (counter = 0; counter < order; counter++) { *solid_angle = *solid_angle * cos_theta; } } } } /* randvec_target_rect_real */ /* SECTION: random numbers ================================================== */ /* * Copyright (c) 1983 Regents of the University of California. * All rights reserved. * * Redistribution and use in source and binary forms are permitted * provided that the above copyright notice and this paragraph are * duplicated in all such forms and that any documentation, * advertising materials, and other materials related to such * distribution and use acknowledge that the software was developed * by the University of California, Berkeley. The name of the * University may not be used to endorse or promote products derived * from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ /* * This is derived from the Berkeley source: * @(#)random.c 5.5 (Berkeley) 7/6/88 * It was reworked for the GNU C Library by Roland McGrath. * Rewritten to use reentrant functions by Ulrich Drepper, 1995. */ /******************************************************************************* * Modified for McStas from glibc 2.0.7pre1 stdlib/random.c and * stdlib/random_r.c. * * This way random() is more than four times faster compared to calling * standard glibc random() on ix86 Linux, probably due to multithread support, * ELF shared library overhead, etc. It also makes McStas generated * simulations more portable (more likely to behave identically across * platforms, important for parrallel computations). *******************************************************************************/ #define TYPE_3 3 #define BREAK_3 128 #define DEG_3 31 #define SEP_3 3 static mc_int32_t randtbl[DEG_3 + 1] = { TYPE_3, -1726662223, 379960547, 1735697613, 1040273694, 1313901226, 1627687941, -179304937, -2073333483, 1780058412, -1989503057, -615974602, 344556628, 939512070, -1249116260, 1507946756, -812545463, 154635395, 1388815473, -1926676823, 525320961, -1009028674, 968117788, -123449607, 1284210865, 435012392, -2017506339, -911064859, -370259173, 1132637927, 1398500161, -205601318, }; static mc_int32_t *fptr = &randtbl[SEP_3 + 1]; static mc_int32_t *rptr = &randtbl[1]; static mc_int32_t *state = &randtbl[1]; #define rand_deg DEG_3 #define rand_sep SEP_3 static mc_int32_t *end_ptr = &randtbl[sizeof (randtbl) / sizeof (randtbl[0])]; mc_int32_t mc_random (void) { mc_int32_t result; *fptr += *rptr; /* Chucking least random bit. */ result = (*fptr >> 1) & 0x7fffffff; ++fptr; if (fptr >= end_ptr) { fptr = state; ++rptr; } else { ++rptr; if (rptr >= end_ptr) rptr = state; } return result; } void mc_srandom (unsigned int x) { /* We must make sure the seed is not 0. Take arbitrarily 1 in this case. */ state[0] = x ? x : 1; { long int i; for (i = 1; i < rand_deg; ++i) { /* This does: state[i] = (16807 * state[i - 1]) % 2147483647; but avoids overflowing 31 bits. */ long int hi = state[i - 1] / 127773; long int lo = state[i - 1] % 127773; long int test = 16807 * lo - 2836 * hi; state[i] = test + (test < 0 ? 2147483647 : 0); } fptr = &state[rand_sep]; rptr = &state[0]; for (i = 0; i < 10 * rand_deg; ++i) random (); } } /* "Mersenne Twister", by Makoto Matsumoto and Takuji Nishimura. */ /* See http://www.math.keio.ac.jp/~matumoto/emt.html for original source. */ /* A C-program for MT19937, with initialization improved 2002/1/26. Coded by Takuji Nishimura and Makoto Matsumoto. Before using, initialize the state by using mt_srandom(seed) or init_by_array(init_key, key_length). Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. The names of its contributors may not be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Any feedback is very welcome. http://www.math.keio.ac.jp/matumoto/emt.html email: matumoto@math.keio.ac.jp */ #include /* Period parameters */ #define N 624 #define M 397 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ #define UPPER_MASK 0x80000000UL /* most significant w-r bits */ #define LOWER_MASK 0x7fffffffUL /* least significant r bits */ static unsigned long mt[N]; /* the array for the state vector */ static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ /* initializes mt[N] with a seed */ void mt_srandom(unsigned long s) { mt[0]= s & 0xffffffffUL; for (mti=1; mti> 30)) + mti); /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ /* In the previous versions, MSBs of the seed affect */ /* only MSBs of the array mt[]. */ /* 2002/01/09 modified by Makoto Matsumoto */ mt[mti] &= 0xffffffffUL; /* for >32 bit machines */ } } /* initialize by an array with array-length */ /* init_key is the array for initializing keys */ /* key_length is its length */ void init_by_array(unsigned long init_key[], unsigned long key_length) { int i, j, k; mt_srandom(19650218UL); i=1; j=0; k = (N>key_length ? N : key_length); for (; k; k--) { mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + init_key[j] + j; /* non linear */ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ i++; j++; if (i>=N) { mt[0] = mt[N-1]; i=1; } if (j>=key_length) j=0; } for (k=N-1; k; k--) { mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i; /* non linear */ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ i++; if (i>=N) { mt[0] = mt[N-1]; i=1; } } mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ } /* generates a random number on [0,0xffffffff]-interval */ unsigned long mt_random(void) { unsigned long y; static unsigned long mag01[2]={0x0UL, MATRIX_A}; /* mag01[x] = x * MATRIX_A for x=0,1 */ if (mti >= N) { /* generate N words at one time */ int kk; if (mti == N+1) /* if mt_srandom() has not been called, */ mt_srandom(5489UL); /* a default initial seed is used */ for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; } for (;kk> 1) ^ mag01[y & 0x1UL]; } y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; mti = 0; } y = mt[mti++]; /* Tempering */ y ^= (y >> 11); y ^= (y << 7) & 0x9d2c5680UL; y ^= (y << 15) & 0xefc60000UL; y ^= (y >> 18); return y; } #undef N #undef M #undef MATRIX_A #undef UPPER_MASK #undef LOWER_MASK /* End of "Mersenne Twister". */ /* End of McCode random number routine. */ /* randnorm: generate a random number from normal law */ double randnorm(void) { static double v1, v2, s; static int phase = 0; double X, u1, u2; if(phase == 0) { do { u1 = rand01(); u2 = rand01(); v1 = 2*u1 - 1; v2 = 2*u2 - 1; s = v1*v1 + v2*v2; } while(s >= 1 || s == 0); X = v1*sqrt(-2*log(s)/s); } else { X = v2*sqrt(-2*log(s)/s); } phase = 1 - phase; return X; } /** * Generate a random number from -1 to 1 with triangle distribution */ double randtriangle(void) { double randnum = rand01(); if (randnum>0.5) return(1-sqrt(2*(randnum-0.5))); else return(sqrt(2*randnum)-1); } /** * Random number between 0.0 and 1.0 (including?) */ double rand01() { double randnum; randnum = (double) random(); randnum /= (double) MC_RAND_MAX + 1; return randnum; } /** * Return a random number between 1 and -1 */ double randpm1() { double randnum; randnum = (double) random(); randnum /= ((double) MC_RAND_MAX + 1) / 2; randnum -= 1; return randnum; } /** * Return a random number between 0 and max. */ double rand0max(double max) { double randnum; randnum = (double) random(); randnum /= ((double) MC_RAND_MAX + 1) / max; return randnum; } /** * Return a random number between min and max. */ double randminmax(double min, double max) { return rand0max(max - min) + max; } /* SECTION: main and signal handlers ======================================== */ /******************************************************************************* * mchelp: displays instrument executable help with possible options *******************************************************************************/ static void mchelp(char *pgmname) { int i; fprintf(stderr, "%s (%s) instrument simulation, generated with " MCCODE_STRING " (" MCCODE_DATE ")\n", mcinstrument_name, mcinstrument_source); fprintf(stderr, "Usage: %s [options] [parm=value ...]\n", pgmname); fprintf(stderr, "Options are:\n" " -s SEED --seed=SEED Set random seed (must be != 0)\n" " -n COUNT --ncount=COUNT Set number of " MCCODE_PARTICLE "s to simulate.\n" " -d DIR --dir=DIR Put all data files in directory DIR.\n" " -t --trace Enable trace of " MCCODE_PARTICLE "s through instrument.\n" " -g --gravitation Enable gravitation for all trajectories.\n" " --no-output-files Do not write any data files.\n" " -h --help Show this help message.\n" " -i --info Detailed instrument information.\n" " --format=FORMAT Output data files using FORMAT=" FLAVOR_UPPER #ifdef USE_NEXUS " NEXUS" #endif "\n\n" ); #ifdef USE_MPI fprintf(stderr, "This instrument has been compiled with MPI support.\n Use 'mpirun %s [options] [parm=value ...]'.\n", pgmname); #endif if(mcnumipar > 0) { fprintf(stderr, "Instrument parameters are:\n"); for(i = 0; i < mcnumipar; i++) if (mcinputtable[i].val && strlen(mcinputtable[i].val)) fprintf(stderr, " %-16s(%s) [default='%s']\n", mcinputtable[i].name, (*mcinputtypes[mcinputtable[i].type].parminfo)(mcinputtable[i].name), mcinputtable[i].val); else fprintf(stderr, " %-16s(%s)\n", mcinputtable[i].name, (*mcinputtypes[mcinputtable[i].type].parminfo)(mcinputtable[i].name)); } #ifndef NOSIGNALS fprintf(stderr, "Known signals are: " #ifdef SIGUSR1 "USR1 (status) " #endif #ifdef SIGUSR2 "USR2 (save) " #endif #ifdef SIGBREAK "BREAK (save) " #endif #ifdef SIGTERM "TERM (save and exit)" #endif "\n"); #endif /* !NOSIGNALS */ } /* mchelp */ /* mcshowhelp: show help and exit with 0 */ static void mcshowhelp(char *pgmname) { mchelp(pgmname); exit(0); } /* mcusage: display usage when error in input arguments and exit with 1 */ static void mcusage(char *pgmname) { fprintf(stderr, "Error: incorrect command line arguments\n"); mchelp(pgmname); exit(1); } /* mcenabletrace: enable trace/mcdisplay or error if requires recompile */ static void mcenabletrace(void) { if(mctraceenabled) mcdotrace = 1; else { fprintf(stderr, "Error: trace not enabled (mcenabletrace)\n" "Please re-run the " MCCODE_NAME " compiler " "with the --trace option, or rerun the\n" "C compiler with the MC_TRACE_ENABLED macro defined.\n"); exit(1); } } /******************************************************************************* * mcreadparams: request parameters from the prompt (or use default) *******************************************************************************/ void mcreadparams(void) { int i,j,status; static char buf[CHAR_BUF_LENGTH]; char *p; int len; MPI_MASTER(printf("Instrument parameters for %s (%s)\n", mcinstrument_name, mcinstrument_source)); for(i = 0; mcinputtable[i].name != 0; i++) { do { MPI_MASTER( if (mcinputtable[i].val && strlen(mcinputtable[i].val)) printf("Set value of instrument parameter %s (%s) [default='%s']:\n", mcinputtable[i].name, (*mcinputtypes[mcinputtable[i].type].parminfo) (mcinputtable[i].name), mcinputtable[i].val); else printf("Set value of instrument parameter %s (%s):\n", mcinputtable[i].name, (*mcinputtypes[mcinputtable[i].type].parminfo) (mcinputtable[i].name)); fflush(stdout); ); #ifdef USE_MPI if(mpi_node_rank == mpi_node_root) { p = fgets(buf, CHAR_BUF_LENGTH, stdin); if(p == NULL) { fprintf(stderr, "Error: empty input for paramater %s (mcreadparams)\n", mcinputtable[i].name); exit(1); } } else p = buf; MPI_Bcast(buf, CHAR_BUF_LENGTH, MPI_CHAR, mpi_node_root, MPI_COMM_WORLD); #else /* !USE_MPI */ p = fgets(buf, CHAR_BUF_LENGTH, stdin); if(p == NULL) { fprintf(stderr, "Error: empty input for paramater %s (mcreadparams)\n", mcinputtable[i].name); exit(1); } #endif /* USE_MPI */ len = strlen(buf); if (!len || (len == 1 && (buf[0] == '\n' || buf[0] == '\r'))) { if (mcinputtable[i].val && strlen(mcinputtable[i].val)) { strncpy(buf, mcinputtable[i].val, CHAR_BUF_LENGTH); /* use default value */ len = strlen(buf); } } for(j = 0; j < 2; j++) { if(len > 0 && (buf[len - 1] == '\n' || buf[len - 1] == '\r')) { len--; buf[len] = '\0'; } } status = (*mcinputtypes[mcinputtable[i].type].getparm) (buf, mcinputtable[i].par); if(!status) { (*mcinputtypes[mcinputtable[i].type].error)(mcinputtable[i].name, buf); if (!mcinputtable[i].val || strlen(mcinputtable[i].val)) { fprintf(stderr, " Change %s default value in instrument definition.\n", mcinputtable[i].name); exit(1); } } } while(!status); } } /* mcreadparams */ /******************************************************************************* * mcparseoptions: parse command line arguments (options, parameters) *******************************************************************************/ void mcparseoptions(int argc, char *argv[]) { int i, j; char *p; int paramset = 0, *paramsetarray; char *usedir=NULL; /* Add one to mcnumipar to avoid allocating zero size memory block. */ paramsetarray = (int*)malloc((mcnumipar + 1)*sizeof(*paramsetarray)); if(paramsetarray == NULL) { fprintf(stderr, "Error: insufficient memory (mcparseoptions)\n"); exit(1); } for(j = 0; j < mcnumipar; j++) { paramsetarray[j] = 0; if (mcinputtable[j].val != NULL && strlen(mcinputtable[j].val)) { int status; char buf[CHAR_BUF_LENGTH]; strncpy(buf, mcinputtable[j].val, CHAR_BUF_LENGTH); status = (*mcinputtypes[mcinputtable[j].type].getparm) (buf, mcinputtable[j].par); if(!status) fprintf(stderr, "Invalid '%s' default value %s in instrument definition (mcparseoptions)\n", mcinputtable[j].name, buf); else paramsetarray[j] = 1; } else { (*mcinputtypes[mcinputtable[j].type].getparm) (NULL, mcinputtable[j].par); paramsetarray[j] = 0; } } for(i = 1; i < argc; i++) { if(!strcmp("-s", argv[i]) && (i + 1) < argc) mcsetseed(argv[++i]); else if(!strncmp("-s", argv[i], 2)) mcsetseed(&argv[i][2]); else if(!strcmp("--seed", argv[i]) && (i + 1) < argc) mcsetseed(argv[++i]); else if(!strncmp("--seed=", argv[i], 7)) mcsetseed(&argv[i][7]); else if(!strcmp("-n", argv[i]) && (i + 1) < argc) mcsetn_arg(argv[++i]); else if(!strncmp("-n", argv[i], 2)) mcsetn_arg(&argv[i][2]); else if(!strcmp("--ncount", argv[i]) && (i + 1) < argc) mcsetn_arg(argv[++i]); else if(!strncmp("--ncount=", argv[i], 9)) mcsetn_arg(&argv[i][9]); else if(!strcmp("-d", argv[i]) && (i + 1) < argc) usedir=argv[++i]; /* will create directory after parsing all arguments (end of this function) */ else if(!strncmp("-d", argv[i], 2)) usedir=&argv[i][2]; else if(!strcmp("--dir", argv[i]) && (i + 1) < argc) usedir=argv[++i]; else if(!strncmp("--dir=", argv[i], 6)) usedir=&argv[i][6]; else if(!strcmp("-h", argv[i])) mcshowhelp(argv[0]); else if(!strcmp("--help", argv[i])) mcshowhelp(argv[0]); else if(!strcmp("-i", argv[i])) { mcformat=FLAVOR_UPPER; mcinfo(); } else if(!strcmp("--info", argv[i])) mcinfo(); else if(!strcmp("-t", argv[i])) mcenabletrace(); else if(!strcmp("--trace", argv[i])) mcenabletrace(); else if(!strcmp("--gravitation", argv[i])) mcgravitation = 1; else if(!strcmp("-g", argv[i])) mcgravitation = 1; else if(!strncmp("--format=", argv[i], 9)) { mcformat=&argv[i][9]; } else if(!strcmp("--format", argv[i]) && (i + 1) < argc) { mcformat=argv[++i]; } else if(!strcmp("--no-output-files", argv[i])) mcdisable_output_files = 1; else if(argv[i][0] != '-' && (p = strchr(argv[i], '=')) != NULL) { *p++ = '\0'; for(j = 0; j < mcnumipar; j++) if(!strcmp(mcinputtable[j].name, argv[i])) { int status; status = (*mcinputtypes[mcinputtable[j].type].getparm)(p, mcinputtable[j].par); if(!status || !strlen(p)) { (*mcinputtypes[mcinputtable[j].type].error) (mcinputtable[j].name, p); exit(1); } paramsetarray[j] = 1; paramset = 1; break; } if(j == mcnumipar) { /* Unrecognized parameter name */ fprintf(stderr, "Error: unrecognized parameter %s (mcparseoptions)\n", argv[i]); exit(1); } } else if(argv[i][0] == '-') { fprintf(stderr, "Error: unrecognized option argument %s (mcparseoptions). Ignored.\n", argv[i++]); } else { fprintf(stderr, "Error: unrecognized argument %s (mcparseoptions). Aborting.\n", argv[i]); mcusage(argv[0]); } } if(!paramset) mcreadparams(); /* Prompt for parameters if not specified. */ else { for(j = 0; j < mcnumipar; j++) if(!paramsetarray[j]) { fprintf(stderr, "Error: Instrument parameter %s left unset (mcparseoptions)\n", mcinputtable[j].name); exit(1); } } free(paramsetarray); #ifdef USE_MPI if (mcdotrace) mpi_node_count=1; /* disable threading when in trace mode */ #endif if (usedir && strlen(usedir) && !mcdisable_output_files) mcuse_dir(usedir); } /* mcparseoptions */ #ifndef NOSIGNALS mcstatic char mcsig_message[256]; /******************************************************************************* * sighandler: signal handler that makes simulation stop, and save results *******************************************************************************/ void sighandler(int sig) { /* MOD: E. Farhi, Sep 20th 2001: give more info */ time_t t1, t0; #define SIG_SAVE 0 #define SIG_TERM 1 #define SIG_STAT 2 #define SIG_ABRT 3 printf("\n# " MCCODE_STRING ": [pid %i] Signal %i detected", getpid(), sig); #ifdef USE_MPI printf(" [proc %i]", mpi_node_rank); #endif #if defined(SIGUSR1) && defined(SIGUSR2) && defined(SIGKILL) if (!strcmp(mcsig_message, "sighandler") && (sig != SIGUSR1) && (sig != SIGUSR2)) { printf("\n# Fatal : unrecoverable loop ! Suicide (naughty boy).\n"); kill(0, SIGKILL); /* kill myself if error occurs within sighandler: loops */ } #endif switch (sig) { #ifdef SIGINT case SIGINT : printf(" SIGINT (interrupt from terminal, Ctrl-C)"); sig = SIG_TERM; break; #endif #ifdef SIGILL case SIGILL : printf(" SIGILL (Illegal instruction)"); sig = SIG_ABRT; break; #endif #ifdef SIGFPE case SIGFPE : printf(" SIGFPE (Math Error)"); sig = SIG_ABRT; break; #endif #ifdef SIGSEGV case SIGSEGV : printf(" SIGSEGV (Mem Error)"); sig = SIG_ABRT; break; #endif #ifdef SIGTERM case SIGTERM : printf(" SIGTERM (Termination)"); sig = SIG_TERM; break; #endif #ifdef SIGABRT case SIGABRT : printf(" SIGABRT (Abort)"); sig = SIG_ABRT; break; #endif #ifdef SIGQUIT case SIGQUIT : printf(" SIGQUIT (Quit from terminal)"); sig = SIG_TERM; break; #endif #ifdef SIGTRAP case SIGTRAP : printf(" SIGTRAP (Trace trap)"); sig = SIG_ABRT; break; #endif #ifdef SIGPIPE case SIGPIPE : printf(" SIGPIPE (Broken pipe)"); sig = SIG_ABRT; break; #endif #ifdef SIGUSR1 case SIGUSR1 : printf(" SIGUSR1 (Display info)"); sig = SIG_STAT; break; #endif #ifdef SIGUSR2 case SIGUSR2 : printf(" SIGUSR2 (Save simulation)"); sig = SIG_SAVE; break; #endif #ifdef SIGHUP case SIGHUP : printf(" SIGHUP (Hangup/update)"); sig = SIG_SAVE; break; #endif #ifdef SIGBUS case SIGBUS : printf(" SIGBUS (Bus error)"); sig = SIG_ABRT; break; #endif #ifdef SIGURG case SIGURG : printf(" SIGURG (Urgent socket condition)"); sig = SIG_ABRT; break; #endif #ifdef SIGBREAK case SIGBREAK: printf(" SIGBREAK (Break signal, Ctrl-Break)"); sig = SIG_SAVE; break; #endif default : printf(" (look at signal list for signification)"); sig = SIG_ABRT; break; } printf("\n"); printf("# Simulation: %s (%s) \n", mcinstrument_name, mcinstrument_source); printf("# Breakpoint: %s ", mcsig_message); if (strstr(mcsig_message, "Save") && (sig == SIG_SAVE)) sig = SIG_STAT; SIG_MESSAGE("sighandler"); if (mcget_ncount() == 0) printf("(0 %%)\n" ); else { printf("%.2f %% (%10.1f/%10.1f)\n", 100.0*mcget_run_num()/mcget_ncount(), 1.0*mcget_run_num(), 1.0*mcget_ncount()); } t0 = (time_t)mcstartdate; t1 = time(NULL); printf("# Date: %s", ctime(&t1)); printf("# Started: %s", ctime(&t0)); if (sig == SIG_STAT) { printf("# " MCCODE_STRING ": Resuming simulation (continue)\n"); fflush(stdout); return; } else if (sig == SIG_SAVE) { printf("# " MCCODE_STRING ": Saving data and resume simulation (continue)\n"); mcsave(NULL); fflush(stdout); return; } else if (sig == SIG_TERM) { printf("# " MCCODE_STRING ": Finishing simulation (save results and exit)\n"); mcfinally(); exit(0); } else { fflush(stdout); perror("# Last I/O Error"); printf("# " MCCODE_STRING ": Simulation stop (abort).\n"); // This portion of the signal handling only works on UNIX #if defined(__unix__) || defined(__APPLE__) signal(sig, SIG_DFL); /* force to use default sighandler now */ kill(getpid(), sig); /* and trigger it with the current signal */ #endif exit(-1); } #undef SIG_SAVE #undef SIG_TERM #undef SIG_STAT #undef SIG_ABRT } /* sighandler */ #endif /* !NOSIGNALS */ /******************************************************************************* * mccode_main: McCode main() function. *******************************************************************************/ int mccode_main(int argc, char *argv[]) { /* double run_num = 0; */ time_t t; #ifdef USE_MPI char mpi_node_name[MPI_MAX_PROCESSOR_NAME]; int mpi_node_name_len; #endif /* USE_MPI */ #ifdef MAC argc = ccommand(&argv); #endif #ifdef USE_MPI MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD, &mpi_node_count); /* get number of nodes */ MPI_Comm_rank(MPI_COMM_WORLD, &mpi_node_rank); MPI_Comm_set_name(MPI_COMM_WORLD, mcinstrument_name); MPI_Get_processor_name(mpi_node_name, &mpi_node_name_len); #endif /* USE_MPI */ t = time(NULL); mcseed = (long)t+(long)getpid(); #ifdef USE_MPI /* *** print number of nodes *********************************************** */ if (mpi_node_count > 1) { MPI_MASTER( printf("Simulation '%s' (%s): running on %i nodes (master is '%s', MPI version %i.%i).\n", mcinstrument_name, mcinstrument_source, mpi_node_count, mpi_node_name, MPI_VERSION, MPI_SUBVERSION); ); } #endif /* USE_MPI */ mcstartdate = (long)t; /* set start date before parsing options and creating sim file */ /* *** parse options ******************************************************* */ SIG_MESSAGE("main (Start)"); mcformat=getenv(FLAVOR_UPPER "_FORMAT") ? getenv(FLAVOR_UPPER "_FORMAT") : FLAVOR_UPPER; mcinstrument_exe = argv[0]; /* store the executable path */ /* read simulation parameters and options */ mcparseoptions(argc, argv); /* sets output dir and format */ #ifdef USE_MPI if (mpi_node_count > 1) { /* share the same seed, then adapt random seed for each node */ MPI_Bcast(&mcseed, 1, MPI_LONG, 0, MPI_COMM_WORLD); /* root sends its seed to slaves */ mcseed += mpi_node_rank; /* make sure we use different seeds per node */ } #endif srandom(mcseed); /* *** install sig handler, but only once !! after parameters parsing ******* */ #ifndef NOSIGNALS #ifdef SIGQUIT if (signal( SIGQUIT ,sighandler) == SIG_IGN) signal( SIGQUIT,SIG_IGN); /* quit (ASCII FS) */ #endif #ifdef SIGABRT if (signal( SIGABRT ,sighandler) == SIG_IGN) signal( SIGABRT,SIG_IGN); /* used by abort, replace SIGIOT in the future */ #endif #ifdef SIGTERM if (signal( SIGTERM ,sighandler) == SIG_IGN) signal( SIGTERM,SIG_IGN); /* software termination signal from kill */ #endif #ifdef SIGUSR1 if (signal( SIGUSR1 ,sighandler) == SIG_IGN) signal( SIGUSR1,SIG_IGN); /* display simulation status */ #endif #ifdef SIGUSR2 if (signal( SIGUSR2 ,sighandler) == SIG_IGN) signal( SIGUSR2,SIG_IGN); #endif #ifdef SIGHUP if (signal( SIGHUP ,sighandler) == SIG_IGN) signal( SIGHUP,SIG_IGN); #endif #ifdef SIGILL if (signal( SIGILL ,sighandler) == SIG_IGN) signal( SIGILL,SIG_IGN); /* illegal instruction (not reset when caught) */ #endif #ifdef SIGFPE if (signal( SIGFPE ,sighandler) == SIG_IGN) signal( SIGSEGV,SIG_IGN); /* floating point exception */ #endif #ifdef SIGBUS if (signal( SIGBUS ,sighandler) == SIG_IGN) signal( SIGSEGV,SIG_IGN); /* bus error */ #endif #ifdef SIGSEGV if (signal( SIGSEGV ,sighandler) == SIG_IGN) signal( SIGSEGV,SIG_IGN); /* segmentation violation */ #endif #endif /* !NOSIGNALS */ mcsiminfo_init(NULL); /* open SIM */ SIG_MESSAGE("main (Init)"); mcinit(); #ifndef NOSIGNALS #ifdef SIGINT if (signal( SIGINT ,sighandler) == SIG_IGN) signal( SIGINT,SIG_IGN); /* interrupt (rubout) only after INIT */ #endif #endif /* !NOSIGNALS */ /* ================ main particle generation/propagation loop ================ */ #if defined (USE_MPI) /* sliced Ncount on each MPI node */ mcncount = mpi_node_count > 1 ? floor(mcncount / mpi_node_count) : mcncount; /* number of rays per node */ #endif /* main particle event loop */ while(mcrun_num < mcncount || mcrun_num < mcget_ncount()) { #ifndef NEUTRONICS mcgenstate(); #endif /* old init: mcsetstate(0, 0, 0, 0, 0, 1, 0, sx=0, sy=1, sz=0, 1); */ mcraytrace(); mcrun_num++; } #ifdef USE_MPI /* merge run_num from MPI nodes */ if (mpi_node_count > 1) { double mcrun_num_double = (double)mcrun_num; mc_MPI_Sum(&mcrun_num_double, 1); mcrun_num = (unsigned long long)mcrun_num_double; } #endif /* save/finally executed by master node/thread */ mcfinally(); #ifdef USE_MPI MPI_Finalize(); #endif /* USE_MPI */ return 0; } /* mccode_main */ #ifdef NEUTRONICS /*Main neutronics function steers the McStas calls, initializes parameters etc */ /* Only called in case NEUTRONICS = TRUE */ void neutronics_main_(float *inx, float *iny, float *inz, float *invx, float *invy, float *invz, float *intime, float *insx, float *insy, float *insz, float *inw, float *outx, float *outy, float *outz, float *outvx, float *outvy, float *outvz, float *outtime, float *outsx, float *outsy, float *outsz, float *outwgt) { extern double mcnx, mcny, mcnz, mcnvx, mcnvy, mcnvz; extern double mcnt, mcnsx, mcnsy, mcnsz, mcnp; /* External code governs iteration - McStas is iterated once per call to neutronics_main. I.e. below counter must be initiancated for each call to neutronics_main*/ mcrun_num=0; time_t t; t = (time_t)mcstartdate; mcstartdate = t; /* set start date before parsing options and creating sim file */ mcinit(); /* *** parse options *** */ SIG_MESSAGE("main (Start)"); mcformat=getenv(FLAVOR_UPPER "_FORMAT") ? getenv(FLAVOR_UPPER "_FORMAT") : FLAVOR_UPPER; /* Set neutron state based on input from neutronics code */ mcsetstate(*inx,*iny,*inz,*invx,*invy,*invz,*intime,*insx,*insy,*insz,*inw); /* main neutron event loop - runs only one iteration */ //mcstas_raytrace(&mcncount); /* prior to McStas 1.12 */ mcallowbackprop = 1; //avoid absorbtion from negative dt int argc=1; char *argv[0]; int dummy = mccode_main(argc, argv); *outx = mcnx; *outy = mcny; *outz = mcnz; *outvx = mcnvx; *outvy = mcnvy; *outvz = mcnvz; *outtime = mcnt; *outsx = mcnsx; *outsy = mcnsy; *outsz = mcnsz; *outwgt = mcnp; return; } /* neutronics_main */ #endif /*NEUTRONICS*/ #endif /* !MCCODE_H */ /* End of file "mccode-r.c". */