.gitignore

@@ -55,3 +55,5 @@ __pycache__/
 *.py[cod]
 *$py.class
 .ipynb_checkpoints
+
+.cache/*

.gitlab-ci.yml

-image: eicweb.phy.anl.gov:4567/whit/image_recipes/root_base:latest
-#image: whit2333/root_base:latest
+image: eicweb.phy.anl.gov:4567/containers/image_recipes/ubuntu_dind:latest
 #
 
 simc_singularity:
-  tags:
-     - singularity
   script:
      - cp containers/singularity/Singularity Singularity.simc
-     - /bin/bash .gitlabci/setup.sh
      - /bin/bash .gitlabci/build.sh Singularity.simc
      - mkdir -p build 
      - cp Singularity.simc.simg build/.

CMakeLists.txt

@@ -5,7 +5,7 @@
 ################################################################################
 cmake_minimum_required (VERSION 3.8)
 
-project (simc VERSION 1.0.0 LANGUAGES CXX )
+project (simc VERSION 1.0.0 LANGUAGES C CXX Fortran)
 
 ################################################################################
 ## CMAKE Settings 

CTP/CMakeLists.txt

@@ -9,7 +9,7 @@ set (PROGRAM "${PROJECT_NAME}_makereg")
 
 set (TARGETS ${TARGETS} ${LIBRARY} ${CLIENT_LIBRARY} ${PROGRAM} PARENT_SCOPE)
 
-## for RPC Command ${RPCCOM}
+## for RPC Command ${RPCCOM} ${RPCGEN_FLAGS}
 include (compiler.cmake)
 
 ################################################################################
@@ -42,23 +42,23 @@ ADD_CUSTOM_COMMAND(
   DEPENDS ${RPC_DEP_HEADER_SRC})
 ADD_CUSTOM_COMMAND(
 	OUTPUT ${RPC_XDR_SOURCE}
-	COMMAND ${RPCCOM} -o ${RPC_XDR_SOURCE} -c ${RPC_RAW}
+	COMMAND ${RPCCOM} ${RPCGEN_FLAGS} -o ${RPC_XDR_SOURCE} -c ${RPC_RAW}
 	DEPENDS ${RPC_RAW})
 ADD_CUSTOM_COMMAND(
 	OUTPUT ${RPC_SVC_SOURCE}
-	COMMAND ${RPCCOM} -o ${RPC_SVC_SOURCE} -m ${RPC_RAW}
+	COMMAND ${RPCCOM} ${RPCGEN_FLAGS} -o ${RPC_SVC_SOURCE} -m ${RPC_RAW}
 	DEPENDS ${RPC_RAW})
 ADD_CUSTOM_COMMAND(
 	OUTPUT ${RPC_SVC_PATCH_SOURCE}
 	COMMAND sed -e 's/extern/static/g' ${RPC_SVC_SOURCE} > ${RPC_SVC_PATCH_SOURCE}
-	DEPENDS ${RCP_SVC_SOURCE})
+	DEPENDS ${RPC_SVC_SOURCE})
 ADD_CUSTOM_COMMAND(
 	OUTPUT ${RPC_CLNT_SOURCE}
-	COMMAND ${RPCCOM} -o ${RPC_CLNT_SOURCE} -l ${RPC_RAW}
+	COMMAND ${RPCCOM} ${RPCGEN_FLAGS} -o ${RPC_CLNT_SOURCE} -l ${RPC_RAW}
 	DEPENDS ${RPC_RAW})
 ADD_CUSTOM_COMMAND(
 	OUTPUT ${RPC_HEADER}
-	COMMAND ${RPCCOM} -o ${RPC_HEADER} -h ${RPC_RAW}
+	COMMAND ${RPCCOM} ${RPCGEN_FLAGS} -o ${RPC_HEADER} -h ${RPC_RAW}
 	DEPENDS ${RPC_RAW})
 
 ################################################################################
@@ -66,6 +66,7 @@ ADD_CUSTOM_COMMAND(
 ################################################################################
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include)
 include_directories(${CMAKE_CURRENT_BINARY_DIR})
+include_directories(SYSTEM /usr/include/tirpc)
 
 ################################################################################
 ## Compile and Link
@@ -76,15 +77,18 @@ ${RPC_SVC_SOURCE} ${RPC_XDR_SOURCE} ${RPC_HEADER})
 set_target_properties(${LIBRARY} PROPERTIES 
   VERSION ${PROJECT_VERSION} 
   SOVERSION ${PROJECT_VERSION_MAJOR})
+set_property(TARGET ${LIBRARY} PROPERTY C_STANDARD 11)
 
 add_library(${CLIENT_LIBRARY} STATIC 
 	${CLIENT_SOURCES} 
 	${RPC_CLNT_SOURCE} ${RPC_XDR_SOURCE} ${RPC_HEADER})
-set_target_properties(${LIBRARY} PROPERTIES 
+set_target_properties(${CLIENT_LIBRARY} PROPERTIES 
   VERSION ${PROJECT_VERSION} 
   SOVERSION ${PROJECT_VERSION_MAJOR})
+set_property(TARGET ${CLIENT_LIBRARY} PROPERTY C_STANDARD 11)
 
 add_executable(${PROGRAM} ${MAKEREG_SOURCES})
+set_property(TARGET ${PROGRAM} PROPERTY C_STANDARD 11)
 target_link_libraries(${PROGRAM} 
 	${CLIENT_LIBRARY} ${LIBRARY} 
 	${CLIENT_LIBRARY} ${LIBRARY} 

CTP/CTP/daVarRpcProc.c

@@ -154,7 +154,7 @@ RVALLIST *davar_readmultiple_1(NAMELIST *argp, CLIENT *clnt)
     }
     free(result.RVALLIST_val);
 #else
-    xdr_free(xdr_RVALLIST, (void *) &result);
+    xdr_free((xdrproc_t)xdr_RVALLIST, (void *) &result);
 #endif
   } else need_to_free = 1;
 /*  minfo = mallinfo();
@@ -180,7 +180,7 @@ ERRLIST *davar_writemultiple_1(WVALLIST *argp, CLIENT *clnt)
   static ERRLIST result;
   int i;
 
-  xdr_free(xdr_ERRLIST, (void *) &result);
+  xdr_free((xdrproc_t)xdr_ERRLIST, (void *) &result);
 
   result.ERRLIST_len = argp->WVALLIST_len;
   result.ERRLIST_val = (int *) malloc(result.ERRLIST_len*sizeof(any));
@@ -285,7 +285,7 @@ davar_readmultiple_test_cb_1(argp, clnt)
 #endif
 
 	memset((char *)&clnt_res, 0, sizeof(clnt_res));
-	if ((clnt_stat = clnt_call(clnt, DAVAR_READMULTIPLE_TEST_CB, xdr_RVALLIST, argp, xdr_int, &clnt_res, TIMEOUT)) != RPC_SUCCESS) {
+	if ((clnt_stat = clnt_call(clnt, DAVAR_READMULTIPLE_TEST_CB, (xdrproc_t)xdr_RVALLIST, argp, (xdrproc_t)xdr_int, &clnt_res, TIMEOUT)) != RPC_SUCCESS) {
 	  printf("clnt call failed, clnt_stat = %d\n",clnt_stat);
 	  clnt_perrno(clnt_stat);
 		return (NULL);
@@ -319,7 +319,7 @@ WVALLIST *davar_readpatternmatch_1(char **argp, CLIENT *clnt)
     }
     free(result.WVALLIST_val);
 #else
-    xdr_free(xdr_WVALLIST, (void *) &result);
+    xdr_free((xdrproc_t)xdr_WVALLIST, (void *) &result);
 #endif
   } else need_to_free = 1;
 /*  minfo = mallinfo();
@@ -439,9 +439,9 @@ int daVarCallBack()
       free(argp->test_condition);
 #else
       if(testresult) {
-	xdr_free(xdr_RVALLIST, (void *) &rpc);
+	xdr_free((xdrproc_t)xdr_RVALLIST, (void *) &rpc);
       }
-      xdr_free(xdr_TESTNAMELIST, (void *) argp);
+      xdr_free((xdrproc_t)xdr_TESTNAMELIST, (void *) argp);
 #endif
       free(argp);
       free(this->sock_in);

CTP/client/thClient.c

@@ -84,7 +84,7 @@ FCALLSCFUN1(INT,thPrintList,THPRTLIST,thprtlist,INT);
 /* Don't really understand the following.  What about ultrix?
    This is probably because of the _ in clnt_create */
 #ifndef __osf__
-FCALLSCFUN4(INT,clnt_create,CLNT_CREATE,clnt_create,STRING,INT,INT,STRING);
+FCALLSCFUN4(INT,myClntCreate,CLNT_CREATE,clnt_create,STRING,INT,INT,STRING);
 #else
 #ifdef BIT64
 FCALLSCFUN4(INT,myClntCreate,CLNT_CREATE,clnt_create,STRING,INT,INT,STRING);

CTP/compiler.cmake

@@ -5,6 +5,8 @@ set(RPCCOM rpcgen -b)
 ## C Compiler Settings 
 ################################################################################
 enable_language (C)
+set(CMAKE_C_STANDARD 89)
+set(CMAKE_C_STANDARD_REQUIRED yes)
 
 ## set special compiler flags
 get_filename_component(C_COMPILER_NAME ${CMAKE_C_COMPILER} NAME)
@@ -12,7 +14,7 @@ get_filename_component(C_COMPILER_NAME ${CMAKE_C_COMPILER} NAME)
 if (C_COMPILER_NAME MATCHES "cc.*" OR 
     C_COMPILER_NAME MATCHES "gcc.*" OR 
     C_COMPILER_NAME MATCHES "clang.*")
-  set (CC_EXTRA_FLAGS "-Wall -W")
+  set (CC_EXTRA_FLAGS "-Wall -W -fomit-frame-pointer -Wno-error=return-type -Wno-unused-but-set-variable -Wno-implicit-function-declaration -Wno-unused-parameter -Wno-unused-variable -Wno-sign-compare -Wno-parentheses")
   set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${CC_EXTRA_FLAGS}")
   set (CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} ${CC_EXTRA_FLAGS}")
   set (CMAKE_C_FLAGS_DEBUG   "${CMAKE_C_FLAGS_DEBUG} ${CC_EXTRA_FLAGS}")
@@ -33,3 +35,27 @@ endif ()
 ################################################################################
 ## f2c settings
 add_definitions("-Df2cFortran")
+
+################################################################################
+## RPCGEN Configuration (macOS Homebrew vs. default)
+################################################################################
+## On macOS, prefer Homebrew rpcsvc-proto rpcgen to avoid broken system rpcgen.
+## On other Unix platforms, use the system-default rpcgen without changes.
+set(RPCCOM "rpcgen")  # default fallback
+set(RPCGEN_FLAGS "")
+
+if (APPLE)
+  find_program(HOMEBREW_RPCGEN NAMES rpcgen PATHS /opt/homebrew/opt/rpcsvc-proto/bin NO_DEFAULT_PATH)
+  if (HOMEBREW_RPCGEN)
+    message(STATUS "Using Homebrew rpcgen: ${HOMEBREW_RPCGEN}")
+    set(RPCCOM ${HOMEBREW_RPCGEN})
+    set(RPCGEN_FLAGS -b)
+  else ()
+    message(WARNING "Homebrew rpcgen not found. Using system default.")
+    set(RPCCOM rpcgen)
+    set(RPCGEN_FLAGS -b)
+  endif ()
+else ()
+  set(RPCCOM rpcgen)
+  set(RPCGEN_FLAGS -b)
+endif ()

cmake/compiler.cmake

@@ -17,7 +17,7 @@ get_filename_component(C_COMPILER_NAME ${CMAKE_C_COMPILER} NAME)
 if (C_COMPILER_NAME MATCHES "cc.*" OR 
     C_COMPILER_NAME MATCHES "gcc.*" OR 
     C_COMPILER_NAME MATCHES "clang.*")
-  set (CC_EXTRA_FLAGS "-fomit-frame-pointer -Wno-error=return-type")
+  set (CC_EXTRA_FLAGS "-fomit-frame-pointer -Wno-error=return-type -std=gnu89 -Wno-unused-but-set-variable -Wno-implicit-function-declaration -Wno-unused-parameter -Wno-unused-variable -Wno-sign-compare -Wno-parentheses")
   set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${CC_EXTRA_FLAGS}")
   set (CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} ${CC_EXTRA_FLAGS}")
   set (CMAKE_C_FLAGS_DEBUG   "${CMAKE_C_FLAGS_DEBUG} ${CC_EXTRA_FLAGS}")
@@ -73,6 +73,9 @@ get_filename_component (Fortran_COMPILER_NAME ${CMAKE_Fortran_COMPILER} NAME)
 if (Fortran_COMPILER_NAME MATCHES "gfortran.*")
 
   set (gfortran_EXTRA_FLAGS "-fomit-frame-pointer -fno-automatic -fno-second-underscore -ffixed-line-length-none -fno-range-check")
+  if (${CMAKE_Fortran_COMPILER_VERSION} VERSION_GREATER 10.0.0)
+     set (gfortran_EXTRA_FLAGS "${gfortran_EXTRA_FLAGS} -fallow-argument-mismatch -fallow-invalid-boz")
+  endif()
   set (CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${gfortran_EXTRA_FLAGS}")
   set (CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} ${gfortran_EXTRA_FLAGS}")
   set (CMAKE_Fortran_FLAGS_DEBUG   "${CMAKE_Fortran_FLAGS_DEBUG} ${gfortran_EXTRA_FLAGS}")

containers/singularity/Singularity

 Bootstrap: docker
-From: eicweb.phy.anl.gov:4567/whit/image_recipes/root_base:latest
+From: eicweb.phy.anl.gov:4567/containers/image_recipes/root_base:latest
 
 %help
   Container for file-input simc

launcher/debug-simc-launcher.in

@@ -51,7 +51,7 @@ tmp_data_path = os.path.basename(data_path)
 tmp_conf_path = os.path.basename(conf_path)
 tmp_conf_name, tmp_conf_ext = os.path.splitext(tmp_conf_path)
 ## make sure our tmp_conf_path ends on '.inp'
-if tmp_conf_ext is '.inp':
+if tmp_conf_ext == '.inp':
     pass
 else:
     tmp_conf_path = '{}.inp'.format(tmp_conf_name)

launcher/simc-launcher.in

@@ -49,7 +49,7 @@ tmp_data_path = os.path.basename(data_path)
 tmp_conf_path = os.path.basename(conf_path)
 tmp_conf_name, tmp_conf_ext = os.path.splitext(tmp_conf_path)
 ## make sure our tmp_conf_path ends on '.inp'
-if tmp_conf_ext is '.inp':
+if tmp_conf_ext == '.inp':
     pass
 else:
     tmp_conf_path = '{}.inp'.format(tmp_conf_name)

nanocernlib @ 55a2f9a4

-Subproject commit 13987d3179c837fd1e8a1c3e40848b8b31085fe5
+Subproject commit 55a2f9a4ebea3d2339f54082c20a50d5c9aa842d

simc/compiler.cmake

@@ -11,7 +11,8 @@ if (Fortran_COMPILER_NAME MATCHES "gfortran.*")
   set (CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} ${gfortran_EXTRA_FLAGS}")
   set (CMAKE_Fortran_FLAGS_DEBUG   "${CMAKE_Fortran_FLAGS_DEBUG} ${gfortran_EXTRA_FLAGS}")
   set (CMAKE_Fortran_FLAGS_RELWITHDEBINFO "${CMAKE_Fortran_FLAGS_RELWITHDEBINFO} ${gfortran_EXTRA_FLAGS}")
-
+## gcc10 optimization
+  set (gfortran_EXTRA_FLAGS "${gfortran_EXTRA_FLAGS} -fallow-argument-mismatch -fallow-invalid-boz")
 ## g77
 elseif (Fortran_COMPILER_NAME MATCHES "g77.*")
 

simc/constants.inc

@@ -9,12 +9,14 @@
 !	All roads lead to infinity
 
 ! ... some physical constants
-	real*8 Me, Me2, Mp, Mp2, Mn, Mn2, Mpi, Mpi2, Mpi0, Mpi02
+	real*8 Me, Me2, Mmu, Mmu2, Mp, Mp2, Mn, Mn2, Mpi, Mpi2, Mpi0, Mpi02
 	real*8 Mk, Mk2, Mrho, Mrho2,Md, Md2
 	real*8 Mlambda, Msigma0, Msigma_minus
 	real*8 amu, hbarc, pi, twopi, alpha, alpi, degrad, euler
 	parameter (Me=0.51099906)
 	parameter (Me2=Me**2)
+	parameter (Mmu=113.4289257)
+	parameter (Mmu2=Mmu**2)
 	parameter (Mp=938.27231)
 	parameter (Mp2=Mp**2)
 	parameter (Mn=939.56563)

simc/dbase.f

@@ -193,6 +193,9 @@ C DJG:
         else if (doing_positron) then
 	   Mh = Me
 	   doing_eep=.false.
+        else if (doing_muons) then
+	   Mh = Mmu
+	   doing_eep=.false.
 	else		!doing_eep if nothing else set.
 	  Mh=Mp
 	  doing_eep = .true.
@@ -304,7 +307,7 @@ C DJG:
 	      sign_hadron=-1.0
 	   endif
 
-	else if(doing_positron) then
+	else if(doing_positron.or.doing_muons) then
 	   targ%Mtar_struck = Mp
 	   targ%Mrec_struck = Mp
 
@@ -776,6 +779,8 @@ C DJG:
 	   endif
 	else if (doing_positron) then
 	  write(6,*) ' ****--------  Doing positron  --------****'
+	else if (doing_muons) then
+	  write(6,*) ' ****--------  Doing muon  --------****'
 	else if (doing_phsp) then
 	  write(6,*) ' ****--------  PHASE SPACE - NO physics, NO radiation  --------****'
 	else
@@ -894,6 +899,7 @@ C DJG:
 	ierr = regparmint('doing_2pi',doing_2pi,0)
 	ierr = regparmint('doing_rho',doing_rho,0)
 	ierr = regparmint('doing_positron',doing_positron,0)
+	ierr = regparmint('doing_muons',doing_muons,0)
 	ierr = regparmint('doing_decay',doing_decay,0)
 	ierr = regparmdouble('ctau',ctau,0)
 

simc/shared/musc.f

@@ -2,7 +2,7 @@
 C+_____________________________________________________________________
 !
 ! MUSC - Simulate multiple scattering of any particle.
-!
+! 
 ! ASSUMPTIONS: DTH and DPH given in milli-radians, RAD_LEN in radiation
 !   lengths. The formula used is due to Rossi and Greisen (See the book
 !   by Segre, NUCLEI AND PARTICLES, 1982, p. 48.) The formula assumes a
@@ -24,11 +24,13 @@ C+_____________________________________________________________________
 C-_____________________________________________________________________
 
 	implicit none
+	include '../simulate.inc'
 
 	real*8 musc_pdg, musc_with_tail
 
 	real*8 rad_len, dth, dph
 	real*8 beta, m2, p
+	logical heavy_tail
 
 	if (rad_len.eq.0) return
 	if (p.lt.25.) write(6,*)
@@ -38,8 +40,13 @@ C-_____________________________________________________________________
 
 	beta = p / sqrt(m2+p*p)
 
-	dth = dth + musc_with_tail(beta, p, rad_len)
-	dph = dph + musc_with_tail(beta, p, rad_len)
+	if (doing_positron) then
+	  dth = dth + musc_with_tail(beta, p, rad_len)
+	  dph = dph + musc_with_tail(beta, p, rad_len)
+	else
+	  dth = dth + musc_pdg(beta, p, rad_len)
+	  dph = dph + musc_pdg(beta, p, rad_len)
+	endif
 
 	return
 	end

simc/simc.f

@@ -795,6 +795,8 @@ c	  write(7,*) 'BP thingie in/out     ',shmsSTOP_BP_in,shmsSTOP_BP_out
 	  write(iun,*) '              ****--- PHASE SPACE - NO physics, NO radiation (may not work)---****'
 	else if (doing_positron) then
 	  write(iun,*) '              ****--- doing positron****'
+	else if (doing_muons) then
+	  write(iun,*) '              ****--- doing muons****'
 	else
 	  stop 'I don''t have ANY idea what we''re doing!!!'
 	endif
@@ -1319,6 +1321,7 @@ c	enddo
 	real*8 eloss_E_arm, eloss_P_arm, r, beta, dangles(2), dang_in(2)
 	logical success
 	logical ok_E_arm, ok_P_arm
+	logical ms_flag, wcs_flag !Mult. scat. and wire chamber smearing
 	type(event):: orig, recon
 	type (event_main)::	main
 
@@ -1350,6 +1353,11 @@ c	enddo
 	  frx = 0.0
 	endif
 
+! Disable multiple scattering for muons, but not the wire chamber
+! smearing
+      ms_flag = mc_smear .and. (.not. doing_muons)
+      wcs_flag = mc_smear
+
 !BEAM MULTIPLE SCATTERING:  NOT YET IMPLEMENTED, JUST GETTING IT READY!
 
 ! ... multiple scattering of the beam.  Generate angles for beam deflection,
@@ -1357,7 +1365,7 @@ c	enddo
 ! ... The yptar offset goes directly to yptar of both arms (small angle approx).
 ! ... xptar is multiplied by cos(theta) to get the xptar offsets.
 
-	if (mc_smear) then
+	if (ms_flag) then
 	  call target_musc(orig%Ein,beta_electron,main%target%teff(1),dang_in)
 	else
 	  dang_in(1)=0.0	!yp offset, goes directly to yp of both arms
@@ -1384,7 +1392,7 @@ c	enddo
 
 ! ... multiple scattering
 
-	if (mc_smear) then
+	if (ms_flag) then
 	  beta = orig%p%p/orig%p%E
 	  call target_musc(orig%p%p, beta, main%target%teff(3), dangles)
 	else
@@ -1401,7 +1409,7 @@ c	enddo
 ! ... change to P arm spectrometer coordinates (TRANSPORT system),
 
 	  if (abs(cos(spec%p%phi)).gt.0.0001) then  !phi not at +/- pi/2
-	    write(6,*) 'y_P_arm, z_P_arm will be incorrect if spec.p.phi <> pi/2 or 3*pi/2'
+	    write(6,*) 'y_P_arm, z_P_arm is incorrect if spec.p.phi <> pi/2 or 3pi/2'
 	    write(6,*) 'spec%p%phi=',spec%p%phi,'=',spec%p%phi*180/pi,'degrees'
 	  endif
 	  delta_P_arm = main%SP%p%delta
@@ -1466,27 +1474,27 @@ C DJG moved this to the last part of generate!!!
 	  if (hadron_arm.eq.1) then
 	    call mc_hms(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
      >		y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
-     >		m2, mc_smear, mc_smear, doing_decay,
+     >		m2, ms_flag, wcs_flag, doing_decay,
      >		ntup%resfac, fry, ok_P_arm, pathlen)
 	  else if (hadron_arm.eq.2) then
 	    call mc_sos(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
      >		y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
-     >		m2, mc_smear, mc_smear, doing_decay,
+     >		m2, ms_flag, wcs_flag, doing_decay,
      >		ntup%resfac, fry, ok_P_arm, pathlen)
 	  else if (hadron_arm.eq.3) then
 	    call mc_hrsr(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
      >		y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
-     >		m2, mc_smear, mc_smear, doing_decay,
+     >		m2, ms_flag, wcs_flag, doing_decay,
      >		ntup%resfac, fry, ok_P_arm, pathlen)
 	  else if (hadron_arm.eq.4) then
 	    call mc_hrsl(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
      >		y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
-     >		m2, mc_smear, mc_smear, doing_decay,
+     >		m2, ms_flag, wcs_flag, doing_decay,
      >		ntup%resfac, fry, ok_P_arm, pathlen)
 	  else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then
 	    call mc_shms(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
      >		y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
-     >		m2, mc_smear, mc_smear, doing_decay,
+     >		m2, ms_flag, wcs_flag, doing_decay,
      >		ntup%resfac, fry, ok_P_arm, pathlen, hadron_arm, use_first_cer)
 	  endif
 
@@ -1578,12 +1586,7 @@ C	  recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100.
 
 ! ... multiple scattering
 
-	if (mc_smear) then
-	  call target_musc(orig%e%p, beta_electron, main%target%teff(2), dangles)
-	else
-	  dangles(1)=0.0
-	  dangles(2)=0.0
-	endif
+	call target_musc(orig%e%p, beta_electron, main%target%teff(2), dangles)
 
 	main%SP%e%yptar = orig%e%yptar + dangles(1) + dang_in(1)
 	main%SP%e%xptar = orig%e%xptar + dangles(2) + dang_in(2)*spec%e%cos_th
@@ -1595,7 +1598,7 @@ C	  recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100.
 ! ... change to E arm spectrometer coordinates (TRANSPORT system),
 
 	  if (abs(cos(spec%e%phi)).gt.0.0001) then  !phi not at +/- pi/2
-	    write(6,*) 'y_E_arm, z_E_arm will be incorrect if spec.e.phi <> pi/2 or 3*pi/2'
+	    write(6,*) 'y_E_arm, z_E_arm is incorrect if spec.e.phi <> pi/2 or 3*pi/2'
 	    write(6,*) 'spec.e.phi=',spec%e%phi,'=',spec%e%phi*180/pi,'degrees'
 	  endif
 	  delta_E_arm = main%SP%e%delta
@@ -1645,31 +1648,34 @@ C	  recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100.
 !	ok flag
 
 	  m2 = me2
+          if (doing_muons) then
+             m2 = Mmu2
+          endif
 	  pathlen = 0.0
 	  if (electron_arm.eq.1) then
 	    call mc_hms(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
      >		y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
-     >		me2, mc_smear, mc_smear, .false.,
+     >		m2, ms_flag, wcs_flag, .false.,
      >		tmpfact, fry, ok_E_arm, pathlen)
 	  else if (electron_arm.eq.2) then
 	    call mc_sos(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
      >		y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
-     >		me2, mc_smear, mc_smear, .false.,
+     >		m2, ms_flag, wcs_flag, .false.,
      >		tmpfact, fry, ok_E_arm, pathlen)
 	  else if (electron_arm.eq.3) then
 	    call mc_hrsr(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
      >		y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
-     >		me2, mc_smear, mc_smear, .false.,
+     >		m2, ms_flag, wcs_flag, .false.,
      >		tmpfact, fry, ok_E_arm, pathlen)
 	  else if (electron_arm.eq.4) then
 	    call mc_hrsl(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
      >		y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
-     >		me2, mc_smear, mc_smear, .false.,
+     >		m2, ms_flag, wcs_flag, .false.,
      >		tmpfact, fry, ok_E_arm, pathlen)
 	  else if (electron_arm.eq.5 .or. electron_arm.eq.6) then
 	    call mc_shms(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
      >		y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
-     >		me2, mc_smear, mc_smear, .false.,
+     >		m2, ms_flag, wcs_flag, .false.,
      >		tmpfact, fry, ok_E_arm, pathlen, electron_arm, use_first_cer)
 	  else if (electron_arm.eq.7 .or. electron_arm .eq. 8) then
              if (abs(spec%p%phi-pi/2) .eq. 10.) then
@@ -1679,7 +1685,7 @@ C	  recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100.
 	     endif
 	    call mc_calo(spec%e%p, spec%e%theta, delta_e_arm, x_e_arm,
      >		y_e_arm, z_e_arm, dx_e_arm, dy_e_arm, xfp, dxfp, yfp, dyfp,
-     >		m2, mc_smear, mc_smear, doing_decay,
+     >		m2, ms_flag, wcs_flag, doing_decay,
      >		ntup%resfac, frx, fry, ok_e_arm, pathlen, using_tgt_field,
      >          zhadron,electron_arm,drift_to_cal)
 	  endif

simc/simulate.inc

@@ -68,8 +68,8 @@
 	logical doing_hyd_elast, doing_deuterium, doing_heavy
 	logical doing_eep, doing_pion, doing_delta, doing_2pi, doing_kaon, doing_rho
 	logical doing_semi, doing_semipi, doing_semika, doing_hplus
-	logical doing_positron
-	integer*4 which_kaon, which_pion
+	logical doing_positron, doing_muons
+	integer*4 padding, which_kaon, which_pion
 	logical using_cit_generation, using_Coulomb, using_Eloss
 	logical correct_Eloss, correct_raster,do_fermi
 	logical using_tgt_field
@@ -92,10 +92,11 @@
      >    extra_dbase_file,tgt_field_file,using_E_arm_montecarlo,using_P_arm_montecarlo,
      >    doing_phsp, using_rad, hard_cuts, doing_hyd_elast, doing_deuterium, doing_heavy,
      >    doing_eep, doing_pion, doing_delta, doing_2pi, doing_kaon, doing_rho, doing_semi,doing_positron,
+     >    doing_muons, 
      >    doing_semipi, doing_semika, doing_hplus, which_kaon,
      >    which_pion, using_cit_generation, using_Coulomb, using_Eloss,
      >    correct_Eloss, correct_raster, do_fermi,using_tgt_field,
-     >    electron_arm, hadron_arm, use_first_cer,
+     >    electron_arm, hadron_arm, use_first_cer, padding,
      >    transparency, use_benhar_sf, random_state_file, random_seed
 
 ! ........ note: make these first two parameters at least ONE BIGGER than the actual dimensions you want to read in