STARSIM / user defined cylindrical decay region

agusumlif/agusumlif.F

@@ -0,0 +1,89 @@
+      subroutine agUsumlif( sumlifx, tlifex, vx, vy, vz )
+
+#include "geant321/gcbank.inc"
+#include "geant321/gckine.inc"
+#include "geant321/gcnum.inc"
+#include "geant321/gconsp.inc"
+#include "geant321/gcphys.inc"
+#include "geant321/gcstak.inc"
+#include "geant321/gctmed.inc"
+#include "geant321/gctrak.inc"
+#include "geant321/gcvolu.inc"
+
+      real :: vx, vy, vz     ! initial track position
+      real :: rmin =   46.6; ! TPC Inner field cage Rmin
+      real :: rmax =  200.0; ! Outer radius of the sensitive volume
+      real :: zmin = -210.0; ! Length of TPC
+      real :: zmax = +210.0; ! Length of TPC
+      real            :: sumlifx
+      real            :: tlifex
+      real            :: rndm(1)
+      real            :: dist, dist_p, dist_l
+
+      integer, parameter        :: nmaxp = 10
+      integer, dimension(nmaxp) :: particles=(/11,12,0,0,0,0,0,0,0,0 /)
+      integer                   :: nparticles = 2
+      integer                   :: i
+
+      character*10 :: command
+      integer      :: npar
+
+!     common /aguslctrl / nparticles, particles
+
+      sumlifx = 0
+
+ 1000 CALL GRNDM(RNDM,1)      
+      sumlifx = -clight * tlifex * log( rndm(1) ) 
+
+      do i=1,nparticles
+      if ( ipart==particles(i) ) then
+
+         ! Predicted path length to the decay point
+         dist = sumlifx * vect(7) / amass
+
+         dist_p = sumlifx * vect(7) * sqrt( vect(4)*vect(4)+vect(5)*vect(5) ) / amass
+         dist_l = sumlifx * vect(7) * vect(6) / amass
+
+         ! Resample if predicted decay point is outside of the bounds of the TPC.
+         ! n.b. Not good for low energy secondaries (loopers)
+         ! n.b. We ignore the initial radial position of the track
+
+         if ( dist_p .lt. rmin .or. dist_p .gt. rmax ) goto 1000
+
+         ! Correct for vertex z position
+         if ( dist_l .lt. zmin - vz .or. dist_l .gt. zmax - vz ) goto 1000
+
+         !write (*,*) ipart, clight, tlifex, sumlifx, r
+         write(*,*) dist_p, dist_l, vert(1:3)
+
+      endif
+      enddo
+
+      return
+c//____________________________________________________________________________________________
+      entry agusl_set
+      call kugetr( rmin )
+      call kugetr( rmax )
+      call kugetr( zmax )
+      zmin = -zmax
+      write(*,*) 'Particles ', particles(1:nparticles)
+      write(*,*) '... will be decayed w/in cylinder:'
+      write(*,*) '... zmin=', zmin, ' zmax=', zmax
+      write(*,*) '... rmin=', rmin, ' zmax=', rmax
+      return
+c//____________________________________________________________________________________________
+      entry agusl_part
+      call kupatl( command, npar )
+      nparticles = 0
+      do nparticles=1,npar
+         call kugeti( particles(nparticles) )
+         if (particles(nparticles).eq.0) exit
+      enddo
+      write(*,*) 'Particles ', particles(1:nparticles)
+      write(*,*) '... will be decayed w/in cylinder:'
+      write(*,*) '... zmin=', zmin, ' zmax=', zmax
+      write(*,*) '... rmin=', rmin, ' zmax=', rmax   
+      return
+c//____________________________________________________________________________________________
+      end subroutine agUsumlif
+

asps/Simulation/starsim/geant/gltrac.F

@@ -62,6 +62,9 @@ SUBROUTINE GLTRAC
 #endif
       PARAMETER (ONE=1)
       COMMON /GCKINE_CONT/ ITRO,TOFO,MECATO
+
+      real agsumlif
+
 C.
 C.    ------------------------------------------------------------------
 *
@@ -185,18 +188,18 @@ SUBROUTINE GLTRAC
       ELSE IF (ITRTYP.EQ.3) THEN
 *       Neutral hadrons
          CALL GRNDM(RNDM,2)
-         SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1))
+         sumlif = agsumlif(tlife) !SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1))
          ZINTHA = -LOG(RNDM(2))
       ELSE IF (ITRTYP.EQ.4) THEN
 *       Charged hadrons
          CALL GRNDM(RNDM,3)
-         SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1))
+         sumlif = agsumlif(tlife) !SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1))
          ZINTHA = -LOG(RNDM(2))
          ZINTDR = -LOG(RNDM(3))
       ELSE IF (ITRTYP.EQ.5) THEN
 *       Muons
          CALL GRNDM(RNDM,5)
-         SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1))
+         sumlif = agsumlif(tlife) !SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1))
          ZINTBR = -LOG(RNDM(2))
          ZINTPA = -LOG(RNDM(3))
          ZINTDR = -LOG(RNDM(4))
@@ -208,7 +211,7 @@ SUBROUTINE GLTRAC
       ELSE IF (ITRTYP.EQ.8) THEN
 *       Ions
          CALL GRNDM(RNDM,3)
-         SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1)) !-- Added lifetime for unstable nuclei (e.g. hypernuclei)
+         sumlif = agsumlif(tlife) !SUMLIF = -CLIGHT*TLIFE*LOG(RNDM(1)) !-- Added lifetime for unstable nuclei (e.g. hypernuclei)
          ZINTHA = -LOG(RNDM(2))
          ZINTDR = -LOG(RNDM(3))
       ENDIF
@@ -223,4 +226,35 @@ SUBROUTINE GLTRAC
 *                                                             END GLTRAC
       END
 
+      real function agsumlif( tlife, vertex )
+        integer, save :: uaddress = 0
+        logical, save :: first = .true.
+        integer       :: csaddr
+        real          :: rndm(1)
+        real          :: vertex(3)
+!#include "geant321/gconsp.inc"
+        double precision, PARAMETER :: CLIGHT=29979245800.
+
+        if ( first ) then
+           first = .false.
+           uaddress = CsADDR( 'agusumlif' )
+        endif
+
+        if ( uaddress .gt. 0 ) then
+
+           call CsJCAL(uaddress,5, agsumlif,tlife,
+     >                 vertex(1),vertex(2),vertex(3), 0,0,0,0,0)
+
+        else
+
+           CALL GRNDM(RNDM,1)          
+           agsumlif = -clight * tlife * log( rndm(1) )
+
+        endif
+
+
+
+        return
+
+      end function agsumlif
 

pams/sim/agusumlif/example.kumac

@@ -0,0 +1,12 @@
+*** Begin of history file: Mon Sep 19 12:58:28 2022
+detp geom y2019a
+gexe $STAR_LIB/StarAgmlUtil.so
+gexe $STAR_LIB/xgeometry.so
+gexe .$STAR_HOST_SYS/lib/libagusumlif.so
+
+* AGUSER/GKINE NTRACK ID [ PTLOW PTHIGH YLOW YHIGH PHILOW PHIHIGH ZLOW ZHIGH option ]
+gkine 10 11 0.1 1.0 -1 1
+trig 1
+gprint vert
+exit
+*** End   of history file: Mon Sep 19 13:00:57 2022

pams/sim/userlifetime/agusumlif.F

@@ -0,0 +1,87 @@
+      subroutine agUsumlif( sumlifx, tlifex, vx, vy, vz )
+
+#include "geant321/gcbank.inc"
+#include "geant321/gckine.inc"
+#include "geant321/gcnum.inc"
+#include "geant321/gconsp.inc"
+#include "geant321/gcphys.inc"
+#include "geant321/gcstak.inc"
+#include "geant321/gctmed.inc"
+#include "geant321/gctrak.inc"
+#include "geant321/gcvolu.inc"
+
+      real :: vx, vy, vz     ! initial track position
+      real :: rmin =   46.6; ! TPC Inner field cage Rmin
+      real :: rmax =  200.0; ! Outer radius of the sensitive volume
+      real :: zmin = -210.0; ! Length of TPC
+      real :: zmax = +210.0; ! Length of TPC
+      real            :: sumlifx
+      real            :: tlifex
+      real            :: rndm(1)
+      real            :: dist, dist_p, dist_l
+
+      integer, parameter        :: nmaxp = 10
+      integer, dimension(nmaxp) :: particles=(/0,0,0,0,0,0,0,0,0,0 /)
+      integer                   :: nparticles = 0
+      integer                   :: i
+
+      character*10 :: command
+      integer      :: npar
+
+!     common /aguslctrl / nparticles, particles
+
+      sumlifx = 0
+
+ 1000 CALL GRNDM(RNDM,1)      
+      sumlifx = -clight * tlifex * log( rndm(1) ) 
+
+      do i=1,nparticles
+      if ( ipart==particles(i) ) then
+
+         ! Predicted path length to the decay point
+         dist = sumlifx * vect(7) / amass
+
+         dist_p = sumlifx * vect(7) * sqrt( vect(4)*vect(4)+vect(5)*vect(5) ) / amass
+         dist_l = sumlifx * vect(7) * vect(6) / amass
+
+         ! Resample if predicted decay point is outside of the bounds of the TPC.
+         ! n.b. Not good for low energy secondaries (loopers)
+         ! n.b. We ignore the initial radial position of the track
+
+         if ( dist_p .lt. rmin .or. dist_p .gt. rmax ) goto 1000
+
+         ! Correct for vertex z position
+         if ( dist_l .lt. zmin - vz .or. dist_l .gt. zmax - vz ) goto 1000
+
+         !write (*,*) ipart, clight, tlifex, sumlifx, r
+         write(*,*) dist_p, dist_l, vert(1:3)
+
+      endif
+      enddo
+
+      return
+c//____________________________________________________________________________________________
+      entry agusl_set
+      call kugetr( rmin )
+      call kugetr( rmax )
+      call kugetr( zmax )
+      zmin = -zmax
+      write(*,*) 'Particles ', particles(1:nparticles)
+      write(*,*) '... will be decayed w/in cylinder:'
+      write(*,*) '... zmin=', zmin, ' zmax=', zmax
+      write(*,*) '... rmin=', rmin, ' zmax=', rmax
+      return
+c//____________________________________________________________________________________________
+      entry agusl_part
+      if ( nparticles .le. 10 ) then
+         nparticles=nparticles+1
+         call kugeti( particles(nparticles) )
+      endif      
+      write(*,*) 'Particles ', particles(1:nparticles)
+      write(*,*) '... will be decayed w/in cylinder:'
+      write(*,*) '... zmin=', zmin, ' zmax=', zmax
+      write(*,*) '... rmin=', rmin, ' zmax=', rmax   
+      return
+c//____________________________________________________________________________________________
+      end subroutine agUsumlif
+

pams/sim/userlifetime/example.kumac

@@ -0,0 +1,16 @@
+*** Begin of history file: Mon Sep 19 12:58:28 2022
+detp geom y2019a
+gexe $STAR_LIB/StarAgmlUtil.so
+gexe $STAR_LIB/xgeometry.so
+gexe .$STAR_HOST_SYS/lib/userlifetime.so
+
+addpart 11
+addpart 12
+
+* AGUSER/GKINE NTRACK ID [ PTLOW PTHIGH YLOW YHIGH PHILOW PHIHIGH ZLOW ZHIGH option ]
+gkine 1 11 0.1 1.0 -1 1
+trig 1
+gprint vert 1
+gprint vert 2
+
+*** End   of history file: Mon Sep 19 13:00:57 2022

pams/sim/userlifetime/userlifetime.cdf

@@ -0,0 +1,26 @@
+*
+* Steering for agusumlif package
+*
+>Name userlifetime
+>MENU userlife
+
+>Guidance
+Configuration of the user lifetime package, allowing users
+to specify a cylindrical region for particle decays
+
+>Command cylinder
+>Guidance
+Specifies the dimensions of the cylindrical region
+>Parameters
+rmin 'Inner radius [cm]' R D=50.0
+rmax 'Outer radius [cm]' R D=200.0
+dz   'Half z [cm]'       R D=200.0
+>action agusl_set
+
+>Command addpart
+>Guidance 
+Add a particle (G3 id) to the list of particles which will be
+decayed within the specified cylinder.  Max 10.
+>Parameters
+ipart 'G3 id' I D=11
+>action agusl_part