First stab at cogen support for N_trains et al

Author: willend

common/lib/share/mccode-r.c

@@ -60,6 +60,7 @@ 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 */
+mcstatic int  NTOF                   = 0; /* Number of TOF "sub-particles" in a TOF_TRAIN */
 mcstatic int  mcusedefaults          = 0; /* assume default value for all parameters */
 mcstatic int  mcdotrace              = 0; /* flag for --trace and messages for DISPLAY */
 mcstatic int  mcnexus_embed_idf      = 0; /* flag to embed xml-formatted IDF file for Mantid */
@@ -1234,6 +1235,9 @@ static void mcruninfo_out(char *pre, FILE *f)
   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");
+  #ifdef TOF_TRAIN
+  fprintf(f, "%sTOF_TRAIN: %d\n", pre, NTOF);
+  #endif
   snprintf(Parameters, CHAR_BUF_LENGTH, "%ld", mcseed);
   fprintf(f, "%sSeed: %s\n",        pre, Parameters);
   fprintf(f, "%sDirectory: %s\n",        pre, dirname ? dirname : ".");
@@ -4800,6 +4804,10 @@ mcparseoptions(int argc, char *argv[])
     else if(!strcmp("--format", argv[i]) && (i + 1) < argc) {
       mcformat=argv[++i];
     }
+#ifdef TOF_TRAIN
+    else if(!strncmp("--tof-trains=", argv[i], 13)) {
+      NTOF=atoi(&argv[i][13]);
+#endif
 #ifdef USE_NEXUS
     else if(!strcmp("--IDF", argv[i])) {
       mcnexus_embed_idf = 1;

mccode/nlib/share/mcstas-r.c

@@ -51,7 +51,7 @@ int (*mcMagneticField) (double, double, double, double,
 * mcsetstate: transfer parameters into global McStas variables
 *******************************************************************************/
 _class_particle mcsetstate(double x, double y, double z, double vx, double vy, double vz,
-			   double t, double sx, double sy, double sz, double p, int mcgravitation, void *mcMagnet, int mcallowbackprop)
+			   double t, double sx, double sy, double sz, double p, int mcgravitation, void *mcMagnet, int mcallowbackprop, int NTOF)
 {
   _class_particle mcneutron;
 
@@ -81,6 +81,14 @@ _class_particle mcsetstate(double x, double y, double z, double vx, double vy, d
   // what about mcneutron._logic ?
   mcneutron._logic.dummy=1;
   // init uservars via cogen'd-function
+
+  #ifdef TOF_TRAIN
+  mcneutron.N_trains=NTOF;
+  mcneutron.t_offset=malloc(NTOF*sizeof(double));
+  mcneutron.p_trains=malloc(NTOF*sizeof(double));
+  mcneutron.alive_trains=malloc(NTOF*sizeof(int));
+  #endif
+
   particle_uservar_init(&mcneutron);
 
   return(mcneutron);

mccode/src/cogen.c.in

@@ -2464,10 +2464,14 @@ cogen_header(struct instr_def *instr, char *output_name)
   cout("  randstate_t randstate[7];"); // for the KISS generator
   cout("  double t, p;     /* time, event weight */");
   cout("  #ifdef TOF_TRAIN");
-  cout("  int Ntof_train; /* initialised like e.g. ncount, seed from cmdline */");
-  cout("  double *tof_train;");
+  cout("  int N_trains; /* initialised like e.g. ncount, seed from cmdline */");
+  cout("  double *t_offset;");
+  cout("  double *p_trains;");
+  cout("  int *alive_trains;");
+  cout("  #pragma acc shape(t_offset[0:N_trains]) init_needed(N_trains)");
+  cout("  #pragma acc shape(p_trains[0:N_trains]) init_needed(N_trains)");
+  cout("  #pragma acc shape(alive_trains[0:N_trains]) init_needed(N_trains)");
   cout("  int Ntof_init;  /* 0/1 - set to 1 once \"times are filled in\" */");
-  cout("  #pragma acc shape(tof_train[0:Ntof_train]) init_needed(Ntof_train)");
   cout("  #endif /* TOF_TRAIN */");
   cout("  long long _uid;  /* Unique event ID */");
   cout("  long _index;     /* component index where to send this event */");
@@ -2542,29 +2546,33 @@ cogen_header(struct instr_def *instr, char *output_name)
 #if MCCODE_PROJECT == 1     /* neutron */
   cout("#pragma acc routine");
   cout("_class_particle mcsetstate(double x, double y, double z, double vx, double vy, double vz,");
-  cout("			   double t, double sx, double sy, double sz, double p, int mcgravitation, void *mcMagnet, int mcallowbackprop);");
+  cout("			   double t, double sx, double sy, double sz, double p, int mcgravitation, void *mcMagnet, int mcallowbackprop, int NTOF);");
   cout("#pragma acc routine");
   cout("_class_particle mcgetstate(_class_particle mcneutron, double *x, double *y, double *z,");
   cout("                           double *vx, double *vy, double *vz, double *t,");
   cout("                           double *sx, double *sy, double *sz, double *p);");
 #elif MCCODE_PROJECT == 2   /* xray */
   cout("#pragma acc routine");
   cout("_class_particle mcsetstate(double x, double y, double z, double kx, double ky, double kz,");
-  cout("			   double phi, double t, double Ex, double Ey, double Ez, double p, int mcgravitation, void *mcMagnet, int mcallowbackprop);");
+  cout("			   double phi, double t, double Ex, double Ey, double Ez, double p, int mcgravitation, void *mcMagnet, int mcallowbackprop, int NTOF);");
   cout("#pragma acc routine");
   cout("_class_particle mcgetstate(_class_particle mcphoton, double *x, double *y, double *z, double *kx, double *ky, double *kz,");
   cout("			   double *phi, double *t, double *Ex, double *Ey, double *Ez, double *p);");
 #endif
   cout("");
   cout("extern int mcgravitation;      /* flag to enable gravitation */");
   cout("#pragma acc declare create ( mcgravitation )");
+  cout("");
+
+  cout("extern int NTOF;      /* TOF_train suppport */");
+  cout("#pragma acc declare create ( NTOF )");
 
   cout("");
   cout("_class_particle mcgenstate(void) {");
 #if MCCODE_PROJECT == 1     /* neutron */
-  cout("  _class_particle particle = mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, mcgravitation, NULL, 0);");
+  cout("  _class_particle particle = mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, mcgravitation, NULL, 0, NTOF);");
 #elif MCCODE_PROJECT == 2   /* xray */
-  cout("  _class_particle particle = mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, mcgravitation, NULL, 0);");
+  cout("  _class_particle particle = mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, mcgravitation, NULL);");
 #endif
   cout("  return(particle);");
   cout("}");