diff --git a/geant3/src/lemsr/common.inc b/geant3/src/lemsr/common.inc index 8c61f43..2dbea41 100644 --- a/geant3/src/lemsr/common.inc +++ b/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 -c ! bfield(4) = polarization (%) -c ! bfield(5) = polarization (%) in sample - integer*4 ix1 !random generator seed + real*4 bfield(5) ! new data card with B-field components +c ! bfield(4) = polarization (%) +c ! bfield(5) = polarization (%) in sample + 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 diff --git a/geant3/src/lemsr/cwn.inc b/geant3/src/lemsr/cwn.inc index 1845c43..db3f5fd 100644 --- a/geant3/src/lemsr/cwn.inc +++ b/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 diff --git a/geant3/src/lemsr/geant_lemsr_main.F b/geant3/src/lemsr/geant_lemsr_main.F index faa5809..1a9c221 100644 --- a/geant3/src/lemsr/geant_lemsr_main.F +++ b/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 diff --git a/geant3/src/lemsr/get_direction.F b/geant3/src/lemsr/get_direction.F index 2e39ebd..bba26ba 100644 --- a/geant3/src/lemsr/get_direction.F +++ b/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) diff --git a/geant3/src/lemsr/gudcay.F b/geant3/src/lemsr/gudcay.F index 174c43c..c9650dc 100644 --- a/geant3/src/lemsr/gudcay.F +++ b/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 diff --git a/geant3/src/lemsr/gukine.F b/geant3/src/lemsr/gukine.F index 839d93e..8f5c015 100644 --- a/geant3/src/lemsr/gukine.F +++ b/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) - y0 = radius * sin(phi0) + x0 = radius * cos(phi0) + beam_parameter(1) + 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) - y0 = radius * sin(phi0)+0.25 + x0 = radius * cos(phi0) + beam_parameter(1) + y0 = radius * sin(phi0) + beam_parameter(2) c c if sets(8) > 0 throw gaussian momentum distribution c diff --git a/geant3/src/lemsr/guout.F b/geant3/src/lemsr/guout.F index 2639c29..e21c08b 100644 --- a/geant3/src/lemsr/guout.F +++ b/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. diff --git a/geant3/src/lemsr/gustep.F b/geant3/src/lemsr/gustep.F index fbc5123..9a8a8e4 100644 --- a/geant3/src/lemsr/gustep.F +++ b/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