Added muon stop position to ntuple output; x0,y0 offsets,

beamspot size scalable and exponential relaxation as additional
option.

geant3/src/lemsr/common.inc

@@ -43,7 +43,7 @@ c
 c      PARFILE: read input parameter from this file
 c
       real*4    MU_MASS, TWO_PI, POS_MASS, SQR_POS_MASS, 
-     1          MAX_TOTAL_ENERGY, TAU_MUON, MU_GYRO, R_MCP
+     1          MAX_TOTAL_ENERGY, TAU_MUON, MU_GYRO, R_MCP, C_LIGHT
       parameter (MU_MASS   = .105658389 )       !muon mass, GeV/c^2
       parameter (POS_MASS  =  0.51099906 )      ! e+ mass, MeV/c^2
       parameter (SQR_POS_MASS =  0.26112 )      ! squared positron mass
@@ -52,15 +52,21 @@ c
       parameter (TAU_MUON     = 2.19703E-6 )    ! muon life time in s
       parameter (MU_GYRO      = 13554.)         ! Hz/G
 c 
-      parameter (R_MCP        = 1.00 )          ! radius of MCP in cm
+      parameter (R_MCP        = 2.0 )          ! radius of MCP in cm
       parameter (TWO_PI = 6.283185307)
+      parameter (C_LIGHT = 2.99792458E10)      ! speed of light in cm/s 
 c
 c-------------------------------------------------------------------
 c
-      real*4      bfield(5)   ! new data card with B-field components
+      real*4      bfield(5)         ! new data card with B-field components
-c                             ! bfield(4) = polarization (%)
+c                                   ! bfield(4) = polarization (%)
-c			      ! bfield(5) = polarization (%) in sample
+c			            ! bfield(5) = polarization (%) in sample
-      integer*4 ix1           !random generator seed
+      real*4      beam_parameter(4) ! new data card with beam parameters
+                                    ! (1): x0
+                                    ! (2): y0
+                                    ! (3): r-fraction  r*R_MCP2 defines size
+                                    ! (4): exponential relaxation rate 1/ns
+      integer*4 ix1                 !random generator seed
 c
 c         count successfull and rejected MICHEL-energy throws
 c
@@ -114,6 +120,7 @@ c-------------------------------------------------------------------------
 c
 c---  the common blocks
 c 
+      common /beam_offset/ beam_parameter
       common /magnetic_field/ bfield
       common /seeds/ ix1
       common /counter/ cou, rejcou

geant3/src/lemsr/cwn.inc

@@ -12,6 +12,7 @@ c
       real*4      nt_x0, nt_y0, nt_z0, nt_px0, nt_py0, nt_pz0, nt_p0
       integer*4   nt_ipart
       real*4      nt_xe, nt_ye, nt_ze, nt_eposend
+      real*4      nt_xm, nt_ym, nt_zm, nt_emend
       real*4      nt_tmcp, nt_tsci    ! times are integer, use 1ns bins
       real*4      nt_xsci, nt_ysci, nt_zsci
       real*4      nt_desci, nt_descipos, nt_descigam, nt_desciele
@@ -32,6 +33,7 @@ c
      1                     nt_p0
       common /partcode/ nt_ipart
       common /end_pos/ nt_xe, nt_ye, nt_ze, nt_eposend
+      common /end_muon/ nt_xm, nt_ym, nt_zm, nt_emend
       common /times/ nt_tmcp, nt_tsci
       common /sci_pos/ nt_xsci, nt_ysci, nt_zsci
       common /sci_elos/ nt_desci, nt_descipos, nt_descigam, nt_desciele

geant3/src/lemsr/geant_lemsr_main.F

@@ -149,6 +149,7 @@ c
       logical*1   l_init_par  /.false./ 
       logical*1   l_code      /.false./
       logical*1   l_endp      /.false./
+      logical*1   l_endm      /.false./
       logical*1   l_time      /.false./
       logical*1   l_scip      /.false./
       logical*1   l_scie      /.false./
@@ -267,6 +268,10 @@ c     define new data card for magnetic field components
 c
       call ffkey('BFIE', bfield, 5, 'R')
 c
+c     define new data card for beam spot scaling and offset
+c
+      call ffkey('BPAR', beam_parameter, 4, 'R')
+c
 c     define new data card for NTUPLE variables
 c
       call ffkey('NTVA', ntvar, 20, 'M')
@@ -437,6 +442,12 @@ c
           l_endp = .true.
           write(6,*) 'NT variable: End position and end energy of e+.'
       endif
+c
+      call glook('endm', ntvar,20,ipres)
+      if ( ipres .ne. 0 ) then
+          l_endm = .true.
+          write(6,*) 'NT variable: End position and end energy of muon.'
+      endif
 c       
       call glook('time', ntvar,20,ipres)
       if ( ipres .ne. 0 ) then
@@ -670,9 +681,14 @@ c
             call hbname(111,'partcode', nt_ipart, 'Ipart:u*4:9')
       endif
       if ( l_endp .or. l_allv ) then
-            call hbname(111,'endpos'  , nt_xe,
+            call hbname(111,'end_pos'  , nt_xe,
      1             'xe:r*4, ye:r*4, ze:r*4, epos:r*4')
       endif
+      if ( l_endm .or. l_allv ) then
+            call hbname(111,'end_muon'  , nt_xm,
+     1             'xm:r*4, ym:r*4, zm:r*4, em:r*4')
+      endif
+
       if ( l_time .or. l_allv ) then
             call hbname(111,'times'   , nt_tmcp, 'tmcp:r*4, tsci:r*4')
       endif

geant3/src/lemsr/get_direction.F

@@ -16,8 +16,10 @@ c
       implicit none
 c
 #define CERNLIB_TYPE
+#include "geant321/gclist.inc"
 #include "geant321/gctrak.inc"
 #include        "common.inc"
+#include        "cwn.inc"
 c
       real*4 pos_tot                  ! total e+ energy including mass
       real*4 p                        ! polarization == bfield(4)
@@ -28,6 +30,8 @@ c
       real*4 Dw                       ! asymmetry factor
       real*4 random(3)
       real*4 x,y,z,r
+      real*4 t_in_sample,tof,beta,gamma_sqr ! to add relaxation in sample
+      real*4 mass, mass_sqr, pz_sqr, z0
 c
 c-----------------------------------------------------------------------------
 c
@@ -35,10 +39,30 @@ c
         y = vect(2)
         z = vect(3)
         r = sqrt(x**2 + y**2)
+  
         if (1.4 .le. z .and. z .le. 1.41 .and. r .le. 3.5) then
-           p = bfield(5) / 100.
+          if ( lsets(2) .eq. 0 ) then ! muon created in MCP2 or sample
+            t_in_sample = tofg*1.e9   ! time in nsec
+          else
+            mass     = MU_MASS*1000. ! muon mass in MeV/c
+            mass_sqr = mass*mass
+            pz_sqr   = nt_pz0*nt_pz0
+            gamma_sqr = (mass_sqr + pz_sqr)/mass_sqr  
+            beta = sqrt(gamma_sqr - 1.)
+            if ( lsets(2) .lt. 0 ) z0 = -lsets(2)
+            tof         = z0/beta/C_LIGHT
+            t_in_sample = (tofg - tof)*1e9 ! time in nsec 
+          endif
+          p = bfield(5) / 100. * exp(-beam_parameter(4)*t_in_sample) ! expo relaxation
+c          print*,' lsets(2)               = ',lsets(2)
+c          print*,' nt_pz0                 = ',nt_pz0
+c          print*,' beta                   = ',beta
+c          print*,' tofg                   = ',tofg
+c          print*,' tof                    = ',tof
+c          print*,' t_in_sample            = ',t_in_sample
+c          print*,' Polarization           = ',p
         else
-           p = bfield(4) / 100.
+          p = bfield(4) / 100.
         endif
 c
 #if defined (OS_UNIX)

geant3/src/lemsr/gudcay.F

@@ -102,7 +102,7 @@ c
       gkin(3,1) = decay_lab(3)
       gkin(4,1) = decay_lab(4)
       gkin(5,1) = 2
-        gpos(1,1) = vect(1)
+      gpos(1,1) = vect(1)
       gpos(2,1) = vect(2)
       gpos(3,1) = vect(3)
 c

geant3/src/lemsr/gukine.F

@@ -74,11 +74,11 @@ c
 c                       ! z=1.40002 corresponds to an implatation depth of the
 c                       ! mu+ of 200 nm.
 c
-            radius = R_MCP * sqrt(random_1(1)) ! uniform distribution on MCP
+            radius = beam_parameter(3) * R_MCP * sqrt(random_1(1)) ! uniform distribution on MCP
 c           radius = 4.7
             phi0   = TWO_PI * random_1(2)
-            x0     = radius * cos(phi0) 
+            x0     = radius * cos(phi0) + beam_parameter(1)
-            y0     = radius * sin(phi0) 
+            y0     = radius * sin(phi0) + beam_parameter(2)
 c
 c---------------------------
 c
@@ -131,6 +131,8 @@ c
 c
             read(70) x0, y0, z0, px0, py0, pz0
             z0 = z0-0.3
+            x0 = x0 + beam_parameter(1)
+            y0 = y0 + beam_parameter(2)
 c
 c-------------------------
 c
@@ -149,10 +151,10 @@ c
 
 c           z0  = -74.              ! start in vacuum upstream from MCP
             z0  = float(lsets(2))   ! start in vacuum at z=z0;
-            radius = R_MCP * sqrt(random_2(1)) ! uniform distribution on MCP
+            radius = beam_parameter(3) * R_MCP * sqrt(random_2(1)) ! uniform distribution on MCP
             phi0   = TWO_PI * random_2(2)
-            x0     = radius * cos(phi0)
+            x0     = radius * cos(phi0) + beam_parameter(1)
-            y0     = radius * sin(phi0)+0.25
+            y0     = radius * sin(phi0) + beam_parameter(2)
 c
 c           if sets(8) > 0 throw gaussian momentum distribution
 c 

geant3/src/lemsr/guout.F

@@ -341,6 +341,10 @@ c
       nt_ye       = 0.
       nt_ze       = 0.
       nt_eposend  = 0.
+      nt_xm       = 0.
+      nt_ym       = 0.
+      nt_zm       = 0.
+      nt_emend    = 0.
       nt_tmcp     = 0
       nt_tsci     = 0
       nt_xsci     = 0.

geant3/src/lemsr/gustep.F

@@ -93,16 +93,23 @@ c
 c       call gsxyz
 c       call gpcxyz
 c
-c       look where the e+ stops
+c     look where the muon decays (user defined partcode 500)
 c
-        if (ipart.eq.2.and.(inwvol.eq.3.or.istop.ge.1)) then
+      if (ipart.eq.500.and.istop.ge.1) then
+            nt_xm      = vect(1)
+            nt_ym      = vect(2)
+            nt_zm      = vect(3)
+            nt_emend   = 1000. * gekin
+      endif
 c
+c     look where the e+ stops
+c
+      if (ipart.eq.2.and.(inwvol.eq.3.or.istop.ge.1)) then
             nt_xe   = vect(1)
             nt_ye   = vect(2)
             nt_ze   = vect(3)
             nt_eposend   = 1000. * gekin
-c
+      endif
-        endif
 c
 c       look now for detector hits
 c