! ********************************************************* SUBROUTINE AVAN(IP1,Z1,A1,XIA2P,T2,T4) ! ! Auxiliary routine for SAMCS ! C ------------------------------------------------------------------------------------ C T2*XIC2 - < tet**2> due to MCS C T4*XIC2/PLAB**2 - < tet**4> due to MCS C ------------------------------------------------------------------------------------ C IP1 - projectile type C ------------------------------------------------------------------------------------- C Z1,A1 - charge and atomic mass of target C ------------------------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON/CSF/Z,A,IP,MOM DIMENSION S2(100,21),S4(100,21) DIMENSION IPS(21),ISP(8) DIMENSION ATBL(100) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOGICAL FIRSTCALL/.TRUE./ SAVE FIRSTCALL,S2,S4 PARAMETER (SC=5.06765D0) EXTERNAL SEXT DATA IPS/1,0,2,2,3,3,4,4,0,4,4,1,0,5,6,7,8,4*0/ DATA ISP/1,3,5,7,14,15,16,17/ DATA ATBL/1.00797,4.0026,6.939,9.0122,10.811,12.01115,14.0067, 15 *.9994,18.9984,20.183,22.9898,24.312,26.9815,28.088,30.9738, 32.064 *,35.453,39.948,39.102,40.08,44.956,47.90,50.942,51.998, 54.9380,55 *.847,58.9332,58.71,63.54,65.37,69.72,72.59,74.9216, 78.96,79.808,8 *3.80,85.47,87.62,88.905,91.22,92.906,95.94,99.0, 101.07,102.905,10 *6.4,107.87,112.4,114.82,118.69,121.75,127.60, 126.9044,131.30,132. *905,137.34,138.91, 140.12,140.907,144.24,147.,150.35,151.98,157.25 *,158.924,162.50, 164.930,167.26,168.934,173.04,174.97,178.49,180.9 *48,183.85, 186.2,190.2,192.2,195.08,196.987,200.59,204.37,207.19,2 *08.980, 210.,210.,222.,223.,226.,227.,232.036,231.,238.03,237.,242 *., 243.,247.,247.,248.,254.,253./ ! IF ( FIRSTCALL ) THEN EPS=1.D-4 NZMAX=100 IPMAX=8 DO NZ=1,NZMAX Z=1.D0*NZ A=ATBL(NZ) DO NIP=1,IPMAX IP=ISP(NIP) R2=MAX(R2P(IP),RAV2(Z,A)) ! fm**2 Q1=SQRT(1.D-1/R2)/SC ! lower limit GeV QMAX=SQRT(1.D3/R2)/SC ! upper limit GeV MOM=2 CALL SGINT1(SEXT,Q1,QMAX,EPS,S2(NZ,NIP)) MOM=4 CALL SGINT1(SEXT,Q1,QMAX,EPS,S4(NZ,NIP)) ENDDO ENDDO FIRSTCALL = .FALSE. IP=IP1 Z=Z1 A=A1 ENDIF R2=RAV2(Z1,A1) R2M=MAX(R2P(IP1),RAV2(Z1,A1)) ! fm**2 Q12=1.D-1/R2M/SC**2 ! GeV**2 IZ=Z1 IPP=IPS(IP1) R2T=R2+R2P(IP1) PAR=Q12*R2T*SC**2/3.D0 T2=LOG(Q12/XIA2P)-1.D0-PAR+S2(IZ,IPP) T4=Q12*(1.D0-PAR/2.d0)+S4(IZ,IPP) RETURN END c ********************************************************* subroutine BETHE_BLOCH(iel,pavg,dens,dedxm) ! ! Get mean energy loss rate from Bethe-Bloch formula ! ! iel: index for element arrays ! pavg: momentum ! dens: density effect parameter ! dedxm: energy loss rate in GeV / (g/cm^2) ! implicit double precision (a-h,o-z) include 'icommon.inc' real*8 ip,mel ! data bbcon/0.307075d-3/ ! [GeV/g cm^2] c mel = mass(1) e = SQRT(pavg**2 + pmass**2) gam = e / pmass bet = pavg / pmass / gam bg = bet * gam bet2 = bet * bet zm = zelcomp(iel,krad) am = aelcomp(iel,krad) ip = ipelcomp(iel,krad) fac = bbcon * (ABS(pcharge)/ee/bet)**2 * zm/am ! tmax = 2d0*mel*bg*bg/(1d0+2d0*gam*mel/pmass+(mel/pmass)**2) arg = 2d0*mel*bg**2*tmax*1d18 / ip**2 if( arg .le. 0. ) then iflag = -25 go to 900 end if dedxm = fac * ( 0.5d0*LOG(arg) - bet2 - dens/2d0 ) ! 900 continue end c ********************************************************* function CS_INEL(p,iptyp) ! ! Get total inelastic cross section [mb] for hadron interactions ! in the range from 0 - 1 GeV/c (0 - 2 GeV/c for protons) ! ! Cross section data from Rev. Part. Prop. (2000) ! implicit double precision (a-h,o-z) ! if( iptyp .eq. 3 ) then ! pi+ if( p .le. 0.425 ) then cs_inel = 0. else if( p.gt.0.425 .and. p.lt.1.000 ) then cs_inel = -117.59 + 793.79*p - 1958.95*p**2 + 2083.40*p**3 & - 789.68*p**4 else cs_inel = 10 end if ! else if( iptyp .eq. -3 ) then ! pi- if( p .le. 0.100 ) then cs_inel = 7 else if( p.gt.0.100 .and. p.lt.1.000 ) then cs_inel = 241.88 - 5790.80*p + 51512.1*p**2 - 213391.*p**3 & + 466937.*p**4 - 557466.*p**5 + 343528.*p**6 - 85537.9*p**7 else cs_inel = 33 end if ! else if( iptyp .eq. 4 ) then ! K+ if( p .le. 0.265 ) then cs_inel = 0. else if( p.gt.0.265 .and. p.lt.0.550 ) then cs_inel = 1 else if( p.gt.0.550 .and. p.lt.1.000 ) then cs_inel = 3.51 - 13.1985*p + 13.2155*p**2 else cs_inel = 5 end if ! else if( iptyp .eq. -4 ) then ! K- if( p .le. 0.245 ) then cs_inel = 60 else if( p.gt.0.245 .and. p.lt.0.290 ) then cs_inel = 61 - 555.6*(p-0.245) else if( p.gt.0.290 .and. p.lt.0.365 ) then cs_inel = 36 else if( p.gt.0.365 .and. p.lt.0.405 ) then cs_inel = 32 + 28*EXP(-2e4*(p-0.385)**2) else if( p.gt.0.405 .and. p.lt.1.000 ) then cs_inel = 4437.18 - 58983.5*p + 325312*p**2 - 955527*p**3 & + 1.61593e6*p**4 - 1.5786e6*p**5 + 827978*p**6 & - 180518*p**7 else cs_inel = 31 end if ! else if( iptyp .eq. 5 ) then ! proton if( p .le. 0.800 ) then cs_inel = 0. else if( p.gt.0.800 .and. p.lt.0.975 ) then cs_inel = 2 else if( p.gt.0.975 .and. p.lt.2.000 ) then cs_inel = -75.5644 + 109.224*p - 29.3942*p**2 else cs_inel = 25 end if ! else if( iptyp .eq. -5 ) then ! pbar if( p .le. 0.240 ) then cs_inel = 172 else if( p.gt.0.240 .and. p.lt.1.000 ) then cs_inel = 254.486 - 390.950*p + 209.053*p**2 else cs_inel = 72 end if end if ! end c ********************************************************* subroutine DECAY(pm) c c Process particle decays c implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 polp(3),jackson real*8 ctau(6),ppd(4),ppn(4) real*8 pcm(4,3),bet(4),bet0(4) ! data ctau/0d0, 658.654d0, 7.8045d0, 3.713d0, 0d0,0d0/ ! TRUE VALUES ! data ctau/0.,1e12,1e-12,3.713,0.,0./ ! fast pi, slow mu decay <<<<<< c c Return immediately for e and p if( ABS(iptyp).eq.1 .or. ABS(iptyp).ge.5 ) go to 900 c if( fastdecay ) then ! force immediate particle decays do i=2,4 ctau(i) = 1.d-12 end do end if ! c Does this particle decay during this step? bg = pm / pmass beta = bg / SQRT( 1d0 + bg**2 ) gamma = bg / beta arg = sdid / ( ctau(ABS(iptyp)) * bg ) ! if( ABS(iptyp).eq.2 .and. declev.eq.5 ) then ! A patch for Bill to suppress muon decay ! Save accumulated normalized proper "time" as "SPINx" pol(1) = pol(1) + arg go to 900 end if ! probdecay = 1d0 - EXP(-arg) rn = RAN1(irnarg) if( rn .gt. probdecay ) go to 900 ! ! --------------------------------------------------------- c Particle decayed if( .not. dectrk ) then ! stop tracking iflag = -52 go to 900 end if ! if( ABS(iptyp).eq.2 .and. (declev.eq.3 .or. declev.eq.4) ) then ! A patch for Bob to suppress muon decay go to 900 end if ! ! Scale momentum vector if PM doesn't match current PP vector call VMAG( pp,pmp ) do i=1,3 pp(i) = pp(i) * pm / pmp end do ! ! Save the parent-particle-in-LAB-frame information pmp = pm do i=1,3 ppd(i) = pp(i) end do call PANG( ppd,cthp,sthp,cphp,sphp ) ! ! Follow the decay product( daughter particle) ! ipnum = ipnum + 1 ! ! The equations in subr GLOREN require beta to be negative bet0(1) = 0d0 bet0(2) = 0d0 bet0(3) = -beta bet0(4) = gamma ! ! --------------------------------------------------------- if( ABS(iptyp) .eq. 2 ) then ! Muon decay mu -> e nu nu ! ! Get the electron 4-vector in the muon rest frame call MUON_DECAY(pcm) if( iflag .lt. 0 ) go to 900 ! DO I=1,4 ! PPN(I) = PCM(I,1) ! END DO ! CALL OUTPUT_NU(1,PPN,1.d0) ! IPTYP = 1 ! GO TO 900 ! Boost electron to frame with z axis aligned with muon in the LAB call GLOREN( bet0,pcm(1,1),ppd ) ! Rotate electron momentum into frame with z axis aligned to ! the problem geometry call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(1,iptyp) ! n.b. the electron polarization is not implemented ! ! Neutrino production data if( neutrino.gt.19 .and. neutrino.lt.100 .and. & nnudk.gt.0 ) then wt = evtwt / nnudk ! do 50 in=1,nnudk call MUON_DECAY(pcm) if( iflag .lt. 0 ) go to 50 ! call GLOREN( bet0,pcm(1,2),ppn ) ! boosted nu ! call GDROT( ppn,cthp,sthp,cphp,sphp ) if( iptyp .gt. 0 ) then ! mu+ decay id = 27 else ! mu- decay id = 26 end if call OUTPUT_NU(id,ppn,wt) ! call GLOREN( bet0,pcm(1,3),ppn ) ! boosted anti-nu call GDROT( ppn,cthp,sthp,cphp,sphp ) if( iptyp .gt. 0 ) then ! mu+ decay id = -26 else ! mu- decay id = -27 end if call OUTPUT_NU(id,ppn,wt) 50 continue end if ! ! --------------------------------------------------------- else if( ABS(iptyp) .eq. 3 ) then ! Pion decay pi -> mu nu ! call GDECA2( mass(3),mass(2),0.d0,pcm ) ! if( declev.eq.2 .or. declev.eq.4) then ! A patch for Bob to force production at 90 degrees in pion frame pmr = 0.02978d0 phi = RANV(0.d0,twopi) pcm(1,1) = pmr * COS(phi) pcm(2,1) = pmr * SIN(phi) pcm(3,1) = 0d0 pcm(4,1) = SQRT( pcm(1,1)**2 + pcm(2,1)**2 + mass(2)**2 ) end if ! epilab = SQRT( pmp**2 + pmass**2 ) ! refers to parent particle do i=1,3 bet(i) = -pp(i) / epilab end do bet(4) = bet0(4) ! call GLOREN( bet,pcm(1,1),ppd ) ! ! Don't do the following rotation since we did the boost along ! the pion direction in the lab and not along the z axis in the lab ! call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(2,iptyp) ! ! Neutrino production data if( neutrino.gt.19 .and. neutrino.lt.100 .and. & nnudk.gt.0 ) then wt = evtwt / nnudk if( iptyp .gt. 0 ) then ! pi+ decay id = 36 else ! pi- decay id = -36 end if ! do in=1,nnudk call GDECA2( mass(3),mass(2),0.d0,pcm ) call GLOREN( bet0,pcm(1,2),ppn ) call GDROT( ppn,cthp,sthp,cphp,sphp ) call OUTPUT_NU(id,ppn,wt) end do end if ! ......................................................... ! Compute muon polarization ! Including many contributions from Y. Fukui, February 2000 ! if( spin ) then ! Ref. S. Hayakawa, Phys. Rev. 108:1533, 1957, eqn. 11 pmucm = SQRT( pcm(1,1)**2 + pcm(2,1)**2 + pcm(3,1)**2 ) pmulab = SQRT( ppd(1)**2 + ppd(2)**2 + ppd(3)**2 ) denom = pmulab * pmucm emucm = pcm(4,1) emulab = ppd(4) hayakawa = emulab*emucm/denom - gamma*mass(2)**2/denom hayakawa = -hayakawa * SIGN(1,iptyp) ! lab frame ! ! Use Wigner angle method (JCG 25 Oct 99) ppilab = pmp ! cthcm = pcm(3,1) / pmucm cthcm = pp(1)*pcm(1,1) + pp(2)*pcm(2,1) + pp(3)*pcm(3,1) cthcm = cthcm / pmucm / ppilab coswig = epilab*pmucm + cthcm*emucm*ppilab coswig = coswig / mass(3) / pmulab wigner = -coswig * SIGN(1,iptyp) ! if( spintrk .eq. 0 ) then pol(1) = 0d0 pol(2) = 0d0 pol(3) = hayakawa if( jpar .le. nprnt ) then write(2,100) jpar,hayakawa,wigner 100 format(' ip,helicity : ',i3,1x,2f10.4) end if go to 800 ! else ! If spin tracking is to be done, use instead a normalized ! 3-component spin vector. ! ! Consider first the PION rest frame where the helicity is -1 for mu+ ! Let POLP be the unit vector along the momentum direction do i=1,3 ! polp(i) = -pcm(i,1) / pmucm * SIGN(1,iptyp) polp(i) = pcm(i,1) / pmucm end do ! ! Boost the spin from the PION frame to the LAB do i=1,3 bet(i) = pp(i) / epilab ! pion beta in lab end do bds = polp(1)*bet(1)+polp(2)*bet(2)+polp(3)*bet(3) do i=1,3 pol(i) = emucm*polp(i) & + (gamma-1d0)*emucm/beta**2*bds*bet(i) & + gamma*pmucm*bet(i) pol(i) = pol(i) / mass(2) ! Hayakawa (8) end do call VMAG(pol,smag) if( smag .lt. 1.d0 ) then iflag = -62 go to 900 end if s0 = SQRT( smag**2 - 1.d0 ) ! Hayakawa (4) if( emulab .lt. gamma*mass(2)**2/emucm ) s0 = -s0 hlab = s0 * mass(2) / pmulab ! Hayakawa (9) hlab = -hlab * SIGN(1,iptyp) ! if( jpar .le. nprnt ) then write(2,110) jpar,hayakawa,wigner,hlab 110 format(' ip,helicity : ',i3,1x,3f10.4) write(2,115) polp,pol 115 format(' Spi,Slab:',2(3f9.4,3x)) end if end if ! end of test on spin tracking ! ! Now let POLP be the spin vector in the muon rest frame ! Use Jackson, 3rd, ed, Eq. 11.158 do i=1,3 bet(i) = ppd(i) / emulab ! muon beta in lab end do gmulab = emulab / mass(2) bds = pol(1)*bet(1)+pol(2)*bet(2)+pol(3)*bet(3) term = gmulab/(gmulab+1d0)*bds do i=1,3 polp(i) = pol(i) - term*bet(i) end do ! ! Ensure proper normalization and load final spin into POL call VMAG(polp,smag) do i=1,3 pol(i) = -polp(i) / smag * SIGN(1,iptyp) ! pol(i) = polp(i) / smag end do ! ! Get helicity from definition jackson = pol(1)*ppd(1) + pol(2)*ppd(2) + pol(3)*ppd(3) ! jackson = -jackson / pmulab * SIGN(1,iptyp) jackson = jackson / pmulab if( jpar .le. nprnt ) then write(2,120) pol,jackson 120 format(' Smu,helicity:',3f9.4,5x,f9.4) end if end if ! end of test on spin ! --------------------------------------------------------- else if( ABS(iptyp) .eq. 4 ) then ! kaon decay ! ! n.b. the muon (and electron) polarization in K decay is ignored c Which K decay mode? rn = RAN1(irnarg) ! if( rn .le. 0.6351d0 ) then ! K -> mu nu call GDECA2( mass(4),mass(2),0.d0,pcm ) call GLOREN( bet0,pcm(1,1),ppd ) call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(2,iptyp) ! ! Neutrino production data if( neutrino.gt.19 .and. neutrino.lt.100 .and. & nnudk.gt.0 ) then wt = evtwt / nnudk if( iptyp .gt. 0 ) then ! K+ decay id = 46 else ! K- decay id = -46 end if ! do in=1,nnudk call GDECA2( mass(4),mass(2),0.d0,pcm ) call GLOREN( bet0,pcm(1,2),ppn ) call GDROT( ppn,cthp,sthp,cphp,sphp ) call OUTPUT_NU(id,ppn,wt) end do end if ! ......................................................... else if( rn .gt. 0.6351d0 .and. rn .le. 0.8467d0 ) then ! K -> pi pi0 call GDECA2( mass(4),mass(2),0.134976d0,pcm ) call GLOREN( bet0,pcm(1,1),ppd ) call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(3,iptyp) ! ......................................................... else if( rn .gt. 0.8467d0 .and. rn .le. 0.8949d0 ) then ! K -> e pi0 nu call GDECA3( mass(4),mass(1),0.134976d0,0.d0,pcm ) call GLOREN( bet0,pcm(1,1),ppd ) call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(1,iptyp) ! ! Neutrino production data if( neutrino.gt.19 .and. neutrino.lt.100 .and. & nnudk.gt.0 ) then wt = evtwt / nnudk if( iptyp .gt. 0 ) then ! K+ decay id = 47 else ! K- decay id = -47 end if ! do in=1,nnudk call GDECA3( mass(4),mass(1),0.134976d0,0.d0,pcm ) call GLOREN( bet0,pcm(1,3),ppn ) call GDROT( ppn,cthp,sthp,cphp,sphp ) call OUTPUT_NU(id,ppn,wt) end do end if ! ......................................................... else if( rn .gt. 0.8949d0 .and. rn .le. 0.9267d0 ) then ! K -> mu pi0 nu call GDECA3( mass(4),mass(2),0.134976d0,0.d0,pcm ) call GLOREN( bet0,pcm(1,1),ppd ) call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(2,iptyp) ! ! Neutrino production data if( neutrino.gt.19 .and. neutrino.lt.100 .and. & nnudk.gt.0 ) then wt = evtwt / nnudk if( iptyp .gt. 0 ) then ! K+ decay id = 46 else ! K- decay id = -46 end if ! do in=1,nnudk call GDECA3( mass(4),mass(2),0.134976d0,0.d0,pcm ) call GLOREN( bet0,pcm(1,3),ppn ) call GDROT( ppn,cthp,sthp,cphp,sphp ) call OUTPUT_NU(id,ppn,wt) end do end if ! ......................................................... else if( rn .gt. 0.9267d0 .and. rn .le. 0.9440d0 ) then ! K -> pi pi0 pi0 call GDECA3( mass(4),mass(3),0.134976d0,0.134976d0,pcm ) call GLOREN( bet0,pcm(1,1),ppd ) call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(3,iptyp) ! ......................................................... else if( rn .gt. 0.9440d0 ) then ! K -> pi pi pi ! n.b. only the first pion is tracked; the other two are ignored call GDECA3( mass(4),mass(3),mass(3),mass(3),pcm ) call GLOREN( bet0,pcm(1,1),ppd ) call GDROT( ppd,cthp,sthp,cphp,sphp ) iptyp = SIGN(3,iptyp) end if end if ! end of test on particle type ! --------------------------------------------------------- ! 800 continue do i=1,3 pp(i) = ppd(i) end do call VMAG(pp,pm) ! return new momentum value ! 900 continue return end ! ********************************************************* subroutine DEDX(hstp,pavg) ! ! Ionization energy loss ! This routine treats continuous energy loss and its fluctuations ! implicit double precision (a-h,o-z) include 'icommon.inc' c real*8 mel,kappa,landaur c The following function gives the maximum value for the c Landau lambda function; cf. GEANT manual, PHYS332-4 PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753, + P5=.74442,P6=1.1934) FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5) ! data c2/0.1571d-4/ data c3/0.1534d-1/ ! [GeV] data gam1/-0.422785/ ! ! Find mean (or most probable) energy loss call DEDXMEAN(hstp,pavg,demean) dedxm = demean / hstp ! ! --------------------------------------------------------- ! de = 0. if( lstrag ) then ! Determine fluctuation of energy loss around mean value ! mel = mass(1) e = SQRT(pavg**2 + pmass**2) gam = e / pmass bet = pavg / pmass / gam bg = bet * gam bet2 = bet * bet zm = mparm(1,krad) am = aelcomp(1,krad) dm = mparm(2,krad) densel = zm * densat ! [# cm^-3] ! if( straglev .eq. 1 ) then ! Gaussian distribution ! Ref. S. Ahlen, RMP 52:142, 1980 arg = c2 * (ABS(pcharge)/ee)**2 * zm*dm/wtmol & * gam**2 * (1.-0.5*bet2) * hstp sige = SQRT( arg ) de = GAUSS( 0.d0,sige ) go to 800 end if ! ! Straggling chi parameter Geant3 manual, p. 247 xi = c3 * (ABS(pcharge)/ee/bet)**2 * zm*dm/wtmol * hstp ! [GeV] eta = bg ratio = mel / pmass emax = (2.*mel*eta**2) / (1.+2.*ratio*gam+ratio**2) kappa = xi / emax if( kappa .lt. 0.d0 ) then iflag = -25 go to 900 else if( kappa .eq. 0.d0 ) then ! allow zero pressure go to 800 end if xmean = -bet2 - LOG(kappa) + gam1 xlamx = FLAND( xmean ) ! if( straglev .eq. 2 ) then ! Landau distribution 100 continue call LANDAU(xlam) if( xlam .gt. xlamx ) go to 100 de = xi * (xlam - xmean ) ! else if( straglev .eq. 3 ) then ! Alternative routine for the Landau distribution 110 continue rn = RAN1(irnarg) ! This test on rn exists in a commented-out section of GEANT routine GFLUCT ! if( rn .gt. 0.980 ) go to 110 ! ??? ! LANDAUR is a minimal repackaging of the GEANT routine GLANDR ! The comments at the head of that routine say it is a copy of the ! CERN library routine RANLAN and that it is "for restricted ! Landau distribution". However, there is no mention in the CERNLIB ! write up that RANLAN has any type of TCUT restrictions. ! xlam = LANDAUR(rn) if( xlam .gt. xlamx ) go to 110 de = xi * (xlam - xmean ) ! else if( straglev .eq. 4 ) then ! Vavilov distribution if( kappa .lt. 0.01 ) then ! Landau limit 120 continue call LANDAU(xlam) if( xlam .gt. xlamx ) go to 120 de = xi * (xlam - xmean ) c else if( kappa .gt. 10. ) then ! Gausssian limit arg = xi * emax * (1. - 0.5*bet2) sige = SQRT( arg ) de = GAUSS( 0.d0,sige ) c else rn = RAN1(irnarg) xlam = VAVILOV(kappa,bet2,rn) de = xi * (xlam - xmean ) end if ! end of test on kappa ! else if( straglev .eq. 5 ) then ! restricted fluctuations ! This takes into account the total energy loss ! for continuous processes for energy transfers less than DCUTx ! RESTRICT_LOSS is different from the other fluctuation routines ! because it returns DE as the TOTAL loss and not the deviation ! from call RESTRICT_LOSS(zm,hstp,pavg,e,dedxm,de) if( iflag .lt. 0 ) go to 900 eloss = de go to 805 end if ! end of test on straglev end if ! end of test on lstrag ! ! --------------------------------------------------------- c Correct the particle 4-momentum c 800 continue eloss = demean + de ! total kinetic energy loss ! 805 continue ! if( spinmatter.eq.1 .and. ABS(iptyp).eq.2 ) then ! Compute spin flip probability, cf. Rossmanith pflip = bet*bet * mel/mass(2) * eloss/e rn = RAN1(irnarg) if( rn .lt. pflip ) then do i=1,3 pol(i) = -pol(i) end do end if end if ! ! Apply the computed energy loss to the actual particle 4-vector call VMAG(pp,pm) e = SQRT(pm**2 + pmass**2) ekin = e - pmass - eloss ! remaining kinetic energy engmin = SQRT(pzmintrk**2 + pmass**2) - pmass ! if( ekin .gt. engmin ) then ! apply loss and continue e = e - eloss ! else ! bring particle to rest e = pmass iflag = -43 go to 900 end if ! pmnew = SQRT(e**2 - pmass**2) c do i=1,3 pp(i) = pp(i) * pmnew / pm end do c if( dedxm .gt. 0. ) then hnextde = MAX( epsf*ekin,engmin ) / dedxm else hnextde = zstep * scalestep end if pavg = pmnew ! 900 continue end ! ********************************************************* subroutine DEDXMEAN(hstp,pavg,demean) ! ! Find mean (or most probable) ionization energy loss ! implicit double precision (a-h,o-z) include 'icommon.inc' real*8 mel,ip data bbcon/0.307075d-3/ ! [GeV/g cm^2] ! mel = mass(1) e = SQRT(pavg**2 + pmass**2) gam = e / pmass bet = pavg / pmass / gam bet2 = bet * bet bg = bet * gam ! ! Get density effect factor if( delev .gt. 1 ) then call DENSITY(bg,dens) else dens = 0d0 end if ! ! Determine mean value of energy loss ! Use Geant3 manual, p. 289, eq. 6 dedxm = 0d0 ! do k=1,nelcomp(krad) am = aelcomp(k,krad) wt = felcomp(k,krad)*am ! if( delev.eq.1 .or. delev.eq.2 )then call BETHE_BLOCH(k,pavg,dens,dedxk) ! else if( delev .eq. 3 ) then ! Restricted Bethe-Bloch loss ! cf. Review of Particle Properties, 2006, eq. 27.7 ! This takes into account the mean energy loss for continuous processes ! for energy transfers less than DCUTx zm = zelcomp(k,krad) ip = ipelcomp(k,krad) fac = bbcon * (ABS(pcharge)/ee/bet)**2 * zm/am if( ABS(iptyp) .eq. 1 ) then tcut = dcute else tcut = dcutm end if tmax = 2d0*mel*bg*bg/(1d0+2d0*gam*mel/pmass+(mel/pmass)**2) tupper = MIN(tcut,tmax) ! arg = 2.*mel*bg**2*tupper*1e18 / ip**2 if( arg .le. 0. ) then iflag = -25 go to 900 end if dedxk = fac * (0.5*LOG(arg)-0.5*(1.+tupper/tmax)*bet2-dens/2.) ! else if( delev .eq. 4 ) then ! most probable energy loss ! cf. Review of Particle Properties, 2010, eq. 27.10 zm = zelcomp(k,krad) ip = ipelcomp(k,krad) * 1d-6 ! [MeV] fac = bbcon/2d0 * (ABS(pcharge)/ee/bet)**2 * zm/am ! [GeV/g cm^2] arg = 2.*mel*1d3*bg**2 / ip psi = fac*1d3 * mparm(2,krad) * hstp*1d2 ! [MeV] dedxk = fac*(LOG(arg)+LOG(psi/ip)+0.2-bet2-dens) ! [GeV/g cm^2] ! else if( delev .eq. 5 ) then ! test mode with dE = const * dz ! Use stepsize along z, not s call BETHE_BLOCH(k,pdelev5,dens,dedxk) end if ! dedxm = dedxm + wt*dedxk end do ! end of loop on nelcomp ! dedxm = dedxm / wtmol * mparm(2,krad) * 1d2 ! [GeV/m] ! if( delev .lt. 5 ) then demean = dedxm * hstp ! [GeV] else demean = dedxm * hdid end if ! 900 continue end ! ********************************************************* subroutine DELTA_RAY(hstp,pm) ! ! Simulate elastic scattering from an atomic electron ! (production of delta ray) ! ! Adapted from GEANT 3.16/21 routines GDRSGA and GDRAY ! implicit double precision (a-h,o-z) include 'icommon.inc' ! emass = mass(1) tkin = pmass * ( SQRT(1.+(pm/pmass)**2) - 1. ) XE=TKIN+ pMASS GAM=XE/ pMASS GAM2=GAM*GAM BET2=1.-1./GAM2 T=GAM-1. ndr = 0 ! number of delta rays ! ! --------------------------------------------------------- if( iptyp .eq. -1 ) then ! electron IF(TKIN.GT.DCUTE)THEN C X=DCUTE/(T*EMASS) X1=1.-X C C Moller scattering C IF(T.GT.2.*DCUTE/EMASS)THEN B1=T*T/GAM2 B2=(2.*GAM-1.)/GAM2 SIG=(B1*(0.5-X)+1./X-1./X1-B2*LOG(X1/X))/BET2 SIG= drnorm *SIG/T IF(SIG.GT.0.) then xint = 1./SIG / densat / 100. ! [m] xm = hstp / xint ndr = POIDEV( xm,irnarg ) if( ndr .le. 0 ) go to 900 else go to 900 end if ! rn = RAN1(idum) ! prob = 1. - EXP(-hstp/xint) ! if( rn .gt. prob ) go to 900 ENDIF end if ! end of test on tkin > dcute ! IF(X.GE.0.5) GO TO 900 CCc=1.-2.*X C do j=1,ndr 10 rn = RAN1(irnarg) E=X/(1.-CCc * rn) C B1=4./(9.*GAM2-10.*GAM+5.) B2=T*T*B1 B3=(2.*GAM2+2.*GAM-1.)*B1 E1=1.-E C SCREJ=B2*E*E-B3*E/E1+B1*GAM2/(E1*E1) C IF(RAN1(irnarg) .GT. SCREJ) GOTO 10 ! call DRVEC( pm,xe,t,e ) end do ! ! --------------------------------------------------------- else if( iptyp .eq. 1 ) then ! positron X=DCUTE/(T*EMASS) X1=1.-X IF(TKIN.GT.DCUTE)THEN C Bhabha scattering C IF(T.GT.DCUTE/EMASS)THEN Y=1./(GAM+1.) Y2=Y*Y C=1.-2.*Y B1=2.-Y2 B2=C*(3.+Y2) C2=C*C B4=C2*C B3=C2+B4 SIG=(1./X-1.)/BET2+B1*LOG(X) + +B2*(1.-X)-B3*(1.-X*X)/2.+B4*(1.-X**3)/3. SIG= drnorm *SIG/T IF(SIG.GT.0.) then xint = 1./SIG / densat / 100. ! [m] xm = hstp / xint ndr = POIDEV( xm,irnarg ) if( ndr .le. 0 ) go to 900 else go to 900 end if ENDIF end if ! end of test on tkin > dcute ! IF(X.GE.1.) GO TO 900 X1=1.-X ! do j=1,ndr 20 rn = RAN1(irnarg) E=X/(1.-X1*rn) C Y=1./(GAM+1.) Y2=Y*Y C=1.-2.*Y B1=2.-Y2 B2=C*(3.+Y2) C2=C*C B4=C2*C B3=C2+B4 B0=GAM2/(GAM2-1.) C SCREJ=(((B4*E-B3)*E+B2)*E-B1)*E+B0 SCREJ=SCREJ/((((B4*X-B3)*X+B2)*X-B1)*X+B0) IF( RAN1(irnarg) .GT.SCREJ) GOTO 20 ! call DRVEC( pm,xe,t,e ) end do ! ! --------------------------------------------------------- else if( ABS(iptyp).eq.2 .or. ABS(iptyp).eq.5 ) then ! spin 1/2 particles, here muon or proton ! IF(TKIN.GT.DCUTM)THEN emmu = pmass TMAX=2.*EMASS*(GAM2-1.)/ + (1.+2.*GAM*EMASS/EMMU+(EMASS/EMMU)**2) ! IF(TMAX.GT.DCUTM)THEN Y=DCUTM/TMAX SIG=(1.-Y+BET2*Y*LOG(Y))/DCUTM SIG=SIG+(TMAX-DCUTM)/(2.*XE**2) SIG=SIG* drnorm *EMASS/BET2 IF(SIG.GT.0.) then xint = 1./SIG / densat / 100. ! [m] xm = hstp / xint ndr = POIDEV( xm,irnarg ) if( ndr .le. 0 ) go to 900 else go to 900 end if end if ! end of test on tmax > dcutm end if ! end of test on tkin > dcutm ! do j=1,ndr 40 rn = RAN1(irnarg) E=1./DCUTM+rn*(1./TMAX-1./DCUTM) E=1./E SCREJ=1.-BET2*(E/TMAX) SCREJ=SCREJ+0.5*(E/ xe)**2 ! extra term for spin 1/2 rn = RAN1(irnarg) IF( rn .GT.SCREJ) GO TO 40 ! E=E/(T*EMASS) call DRVEC( pm,xe,t,e ) end do ! ! --------------------------------------------------------- else if( ABS(iptyp).eq.3 .or. ABS(iptyp).eq.4 ) then ! spin 0 particles, here pion or kaon ! cf. GEANT 3.16 manual, p.290 ! IF(TKIN.GT.DCUTM)THEN emmu = pmass TMAX=2.*EMASS*(GAM2-1.)/ + (1.+2.*GAM*EMASS/EMMU+(EMASS/EMMU)**2) ! IF(TMAX.GT.DCUTM)THEN Y=DCUTM/TMAX SIG=(1.-Y+BET2*Y*LOG(Y))/DCUTM SIG=SIG* drnorm *EMASS/BET2 IF(SIG.GT.0.) then xint = 1./SIG / densat / 100. ! [m] xm = hstp / xint ndr = POIDEV( xm,irnarg ) if( ndr .le. 0 ) go to 900 else go to 900 end if end if ! end of test on tmax > dcutm end if ! end of test on tkin > dcutm ! do j=1,ndr 50 rn = RAN1(irnarg) E=1./DCUTM+rn*(1./TMAX-1./DCUTM) E=1./E SCREJ=1.-BET2*(E/TMAX) IF(RAN1(irnarg) .GT.SCREJ) GO TO 50 ! E=E/(T*EMASS) call DRVEC( pm,xe,t,e ) end do ! end if ! end of test on particle type ! 900 continue return end ! ********************************************************* subroutine DENSITY(bg,dens) ! ! Get density effect parameter ! implicit double precision (a-h,o-z) include 'icommon.inc' ! x0 = mparm(5,krad) x1 = mparm(6,krad) sa = mparm(7,krad) sc = mparm(8,krad) sm = mparm(9,krad) dens = 0.d0 ! if( x0 .eq. 0d0 ) go to 900 ! undefined density parameters ! if( bg .le. 0. ) then iflag = -25 go to 900 end if ! xs = LOG10( bg ) if( xs.ge.x0 .and. xs.le.x1 )then dens = 4.6052 * xs + sa * (x1 - xs)**sm + sc else if( xs .gt. x1 )then dens = 4.6052 * xs + sc end if ! 900 continue end ! ********************************************************* REAL*8 FUNCTION DFL(AM2,AM3) ! ! Auxiliary routine for SAMCS ! IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) common/rans/iseed ! REAL RNDM PARAMETER (PI=3.141592653589793227D+00) C=AM3/AM2**(1.5) U1=SQRT(1.D0+C**2/4.D0) ZU=C/2.D0-U1 IF(ZU.GE.0.D0) THEN U=(C/2.D0+U1)**(1.D0/3.D0)+ZU**(1.D0/3.D0) ELSE U=(C/2.D0+U1)**(1.D0/3.D0)-(-ZU)**(1.D0/3.D0) ENDIF ALP=SQRT(LOG(1.D0+U**2)) ! R1=DBLE(RNDM(-1.)) ! R2=DBLE(RNDM(-1.)) R1 = RAN1(iseed) R2 = RAN1(iseed) R12=SIN(2.D0*PI*R1)*SQRT(-2.D0*LOG(R2)) DFL=SQRT(AM2)*(EXP(ALP*R12)/SQRT(1.D0+U**2)-1.D0)/U RETURN END ! ********************************************************* real*8 function DOT(a,b) ! ! 4-vector product of A and B ! Adapted from J. Gallardo 10 April 2000 ! implicit real*8(a-h,o-z) dimension a(4),b(4) ! prod=0 ! do i=1,4 if(i.eq.1)then prod=prod+a(i)*b(i) else prod=prod-a(i)*b(i) ! METRIC (+,-,-,-) endif end do ! dot=prod ! end ! ********************************************************* subroutine DRVEC(pm,xe,t,e) ! ! Update particle momentum vector for delta ray energy loss ! ! Adapted from GEANT 3.16/21 routine GDRAY ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 p(3) real*8 pels(3),gkin(3),p4(3),me,mu ! ! ndeltas = ndeltas + 1 ! emass = mass(1) EEL=(T*E+1.)*EMASS ! recoil electron total energy TEL=EEL-EMASS ! kinetic energy PEL=SQRT(ABS((EEL+EMASS)*TEL)) COSTH=(XE*EEL+EMASS*(TEL-XE))/(Pm*PEL) IF(COSTH.GE.1.) THEN COSTH=1. SINTH=0. ELSEIF(COSTH.LE.-1.) THEN COSTH=-1. SINTH=0. ELSE SINTH=SQRT((1.+COSTH)*(1.-COSTH)) ENDIF PHI = TWOPI* RAN1(irnarg) COSPHI = COS(PHI) SINPHI = SIN(PHI) ! C Polar co-ordinates to momentum components. GKIN(1)=PEL*SINTH*COSPHI GKIN(2)=PEL*SINTH*SINPHI GKIN(3)=PEL*COSTH C PELS(1)=-GKIN(1) PELS(2)=-GKIN(2) PELS(3)=Pm -GKIN(3) C do i=1,3 p4(i) = pp(i) ! at this point pp is still incident momentum end do ! Get angle of incident particle in LAB call PANG(p4,cthet,sthet,cphi,sphi) ! ! Rotate scattered particle to correct angle in LAB call GDROT(pels,cthet,sthet,cphi,sphi) ! if( spintrk .gt. 0 ) then do i=1,3 p(i) = pels(i) end do ! Rotate the spin such that helicity is conserved call ROTATE_SPIN(p) ! if( spinmatter.eq.2 .and. ABS(iptyp).eq.2 ) then ! Simulate spin flip in mu-e elastic scattering me = mass(1) mu = mass(2) e1 = SQRT( pm**2 + mu**2 ) s = mu**2 + me**2 + 2.*e1*me ecm = SQRT( s ) ! Get CM scattering angle from Jackson, 1st ed, Eq. 12.48 pms = SQRT( pels(1)**2+pels(2)**2+pels(3)**2 ) e3 = SQRT( pms**2 + mu**2 ) cthcm = e3-0.5*(e1+me)*(1.+(mu**2-me**2)/ecm**2) arg = (1.-((mu+me)/ecm)**2)*(1.-((mu-me)/ecm)**2) ! if( arg .ge. 0d0 ) then ! write(*,*) 'Hello from DRVEC' denom = SQRT(arg) ! else ! stop 'Negative square root in DRVEC' ! end if denom = 0.5*pm * denom cthcm = cthcm / denom ! if( cthcm .ge. 1.d0 ) then cthcm = 1.d0 else if( cthcm .le. -1.d0 ) then cthcm = -1.d0 end if ! call FLIP_ELASTIC( s,cthcm,pr ) ! rn = RAN1(irnarg) if( rn .lt. pr ) then do i=1,3 ! flip spin direction pol(i) = -pol(i) end do ! nflips = nflips + 1 end if end if end if ! do i=1,3 ! load pp with new, rotated scattered momentum pp(i) = pels(i) end do call VMAG(pp,pm) ! return end c ********************************************************* subroutine ELMS(hstp,pavg) c c Ionization energy loss c c This routine treats continuous energy loss and its fluctuations c ! Based on a routine by Simon Holmes, 22 NOV 2004 ! implicit double precision (a-h,o-z) include 'icommon.inc' c ! real demean ! real*8 mel,ip,kappa,landaur c elms variables real*8 dspran ! integer jint c REAL*8 pigreek, c, density PARAMETER (pigreek = 3.1415926535d0) PARAMETER (c = 2.997924d8) ! PARAMETER (elmsmass=105.658) PARAMETER (density=0.0708d0) !g cm2 c CHARACTER (100) trajectoryfolder, logname, stagenum, runnumber !!! NOW IN ICOMMON ! REAL sran ! REAL*8 dran ! INTEGER*2 ihr,imin,isec,ihun ! REAL*8 sigamin(8) REAL*8 meanin(8) c DB variables recorded in the directory (in icommon) c c incoming vectors etc ! INTEGER stage, rstage !cooling stage and number of stages, stage to be analysed INTEGER nstage !cooling stage and number of stages, stage to be analysed INTEGER maxtraj PARAMETER (maxtraj = 100001) !max number of trajectories INTEGER ntraj !actual number of trajectories ! REAL*8 px1(maxtraj), py1(maxtraj), pz1(maxtraj) REAL*8 pzoffset ! REAL*8 x1(maxtraj), y1(maxtraj), z1(maxtraj), t1(maxtraj) c input 6D phase space trajectories, and range (px, py, pz, x, y, z, t) REAL*8 slabthk, zend !slab thickness in z (metres), end of this absorber ! REAL rfboost, rfdelay1, rffreq !pz peak boost per absorber at time rfdelay*stage. ! REAL rfdelay2 ! REAL betam ! REAL drift !drift distance between absorbers ! REAL forward !output 6D phase space trajectories and range(z2) ! REAL*8 px2(maxtraj), py2(maxtraj), pz2(maxtraj) ! REAL*8 x2(maxtraj), y2(maxtraj), z2(maxtraj), t2(maxtraj) c ! working vars ! INTEGER I, J, K, L, N, Ijunk, Kjunk, Itraj, LL, IL, Jmom, ! & KP, KPU, KPD ! REAL*8 a1, a2, a3, a4, a5, a6, px, py, pz, ! & x, y, z, t, p, ke, pnew, dv ! REAL*8 deltapl, deltapt, deltae, steplen, vel REAL*8 UCpl1,UCpl2,UCpt1,UCpt2,UCang1,UCang2 REAL*8 UCen1,UCen2 ! store two pt,pl for uncorrelated INTEGER UCINT REAL*8 llenwt, cphi, sphi, ctheta, stheta REAL*8 oldi, oldj, oldk, cphir, sphir, cthetas, sthetas REAL*8 v(7) EQUIVALENCE (v(1), px), (v(2), py), (v(3), pz), (v(4), x), & (v(5), y), (v(6), z), (v(7), t) !for trajectory concerned ! REAL tthk, pmom ! external ranlux ! --------------------------------------------------------- slabthk=hstp !*** get from icool meanin(3)=pp(3)*1000 !*** get from icool pzoffset=meanin(3) nstage = 1 ntraj = 1 drift = 0 betam = meanin(3) / dsqrt (meanin(3)**2 + elmsmass**2) px=pp(1)*1000 py=pp(2)*1000 pz=meanin(3) !*** get from icool c zend = stage * slabthk + (stage - 1) * drift zend=slabthk meanout = 0.0 matrixout = 0.0 meanboost = 0.0 meanloss = 0.0 Jmom = 2 x=0.0 y=0.0 z=0.0 t=0.0 ! K = Itraj / 10000 ! not used c IF (Itraj - 10000*K .eq. 0) write (6,*) Itraj x=0.0 y=0.0 z=0.0 t=0.0 en1 = dsqrt(px**2 + py**2 + pz**2 + elmsmass**2) DO 30000 LL=1,10000 p = dSQRT(px**2 + py**2 + pz**2) ke = dSQRT(p**2 + elmsmass**2) - elmsmass IF( zend - z .LE. 1.0E-6 .OR. p .le. 5.0.OR.pz.le.5.0) & GOTO 30010 !unit vector in old direction oldi = px / p oldj = py / p oldk = pz / p !find the closest sample set in momentum IF (LL .EQ. 1 .OR. 2*P .GT. mtmdb(KP) + mtmdb(KPU) & .OR. 2*P .LT. mtmdb(KP) + mtmdb(KPD)) THEN N = npdone Jmom = 1 DO 30100 L = 1, 8 N = min0((N+1)/2, npdone-Jmom-1) IF (P .GT. mtmdb(pdone(Jmom+N))) THEN Jmom = Jmom + N IF (Jmom .GE. npdone-1) GOTO 30110 ELSEIF (N .EQ. 1) THEN GOTO 30110 ENDIF 30100 CONTINUE 30110 a1 = 0.5 * (mtmdb(pdone(Jmom))+mtmdb(pdone(Jmom+1))) IF (P .GT. a1) THEN Jmom = Jmom + 1 !take nearest a2 = a1 a1 = 0.5 * (mtmdb(pdone(Jmom))+mtmdb(pdone(Jmom+1))) ELSE a2 = 0.5 * (mtmdb(pdone(Jmom-1))+mtmdb(pdone(Jmom))) ENDIF KP = pdone(Jmom) KPU = pdone(Jmom+1) KPD = pdone(Jmom-1) IF (P .GT. a1 .OR. P .LT. a2) THEN c write(6,5721)P, mtmdb(KPD), mtmdb(KP), mtmdb(KPU) write(2,5721)P, mtmdb(KPD), mtmdb(KP), mtmdb(KPU) 5721 Format(' Momentum interpol fail ',4F10.3) endif ENDIF IF (ELMSCOR.eq.1)THEN !now sample momentum Jmom !generate integer between 0 and norm(Jmom) ! call ranlux(randvw,2) randvw(1) = RAN1(irnarg) randvw(2) = RAN1(irnarg) c dspran=(randvw(1)*1000000000+randvw(2))/1000000000 old technique (flawed) dspran=(IDint(randvw(1)*32768.0D0)+randvw(2))/32768.0D0 ! combine two sp to get one dp a6 = normelms(KP) * dspran !find it I = firstpoint(KP) N = lastpoint(KP) - I + 1 DO 30200 L = 1, 20 N = min0((N+1)/2, lastpoint(KP) - I) IF (a6 .GT. probarray(I+N)) THEN I = I + N IF (I.eq.lastpoint(KP)) GOTO 30210 ELSEIF (N .EQ. 1) THEN GOTO 30210 ENDIF IF (I .GE. lastpoint(KP) .OR. N .EQ. 0) GOTO 30210 !N=0 shouldnt happen, misordered? 30200 CONTINUE !check algorithm 30210 IF (I .GT. lastpoint(KP)) THEN c write (6,3461) a6, probarray(lastpoint(KP)), I, N write (8,3461) a6, probarray(lastpoint(KP)), I, N 3461 Format(' Select fail1', 2D15.8, 2I6) endif IF (L.GT.19.OR.N.EQ.0)THEN c write (6,3462) a6, L,N,I !L=20 would imply a choice of more than 10**6 entries write (8,3462) a6, L,N,I 3462 Format(' Select fail0', D15.8, 3I6) endif IF ((a6.LT.probarray(I-1).AND.I.GE.firstpoint(KP)) .OR. & (a6 .GT. probarray(I+1) .AND. I .NE. lastpoint(KP))) THEN c write(6,3463)a6,probarray(I),probarray(I+1), L, N, I write(2,3463)a6,probarray(I),probarray(I+1), L, N, I 3463 Format( 'select fail3', 3E15.8, 3I6) ENDIF deltapt = dptelms(I) * (1.0 + pslope(KP) * (P - mtmdb(KP))) !interpolation deltae = elmsde(I) * (1.0 + eslope(KP) * (P - mtmdb(KP))) deltapl = p - dsqrt(dabs(p**2 - deltapt**2 + deltae**2 & - 2.0 * deltae * (ke + elmsmass) )) !interpolation check a2 = abs(pslope(KP) * (P - mtmdb(KP))) a6 = abs(eslope(KP) * (P - mtmdb(KP))) IF (P.lt.100.)THEN a3 = 0.10 !check at 10% for low momentum ELSE a3 = 0.03 !otherwise 3% ENDIF IF(a2.gt.a3.or.a6.gt.a3)then c write(6,3464)P,a2,a6 write(2,3464)P,a2,a6 3464 Format(' interpolation too big: mtm, pt2, de', F8.1, 2F8.3) endif ELSE !Uncorrolated version 30300 CONTINUE DO 30400 UCINT=1,2 !now sample momentum Jmom !generate integer between 0 and norm(Jmom) ! call ranlux(randvw,2) randvw(1) = RAN1(irnarg) randvw(2) = RAN1(irnarg) dspran=(IDint(randvw(1)*32768.0D0)+randvw(2))/32768.0D0 ! combine two sp to get one dp a6 = normelms(KP) * dspran !find it I = firstpoint(KP) N = lastpoint(KP) - I + 1 DO 30500 L = 1, 20 N = min0((N+1)/2, lastpoint(KP) - I) IF (a6 .GT. probarray(I+N)) THEN I = I + N IF (I.eq.lastpoint(KP)) GOTO 30510 ELSEIF (N .EQ. 1) THEN GOTO 30510 ENDIF IF (I .GE. lastpoint(KP) .OR. N .EQ. 0) GOTO 30510 !N=0 shouldnt happen, misordered? 30500 CONTINUE !check algorithm 30510 IF (I .GT. lastpoint(KP)) THEN write (6,3461) a6, probarray(lastpoint(KP)), I, N write (8,3461) a6, probarray(lastpoint(KP)), I, N endif IF (L.GT.19.OR.N.EQ.0)THEN write (6,3462) a6, L,N,I !L=20 would imply a choice of more than 10**6 entries write (8,3462) a6, L,N,I endif IF ((a6.LT.probarray(I-1).AND.I.GE.firstpoint(KP)) .OR. & (a6 .GT. probarray(I+1) .AND. I .NE. lastpoint(KP))) THEN write(2,3463)a6,probarray(I),probarray(I+1), L, N, I ! write(8,3463)a6,probarray(I),probarray(I+1), L, N, I ENDIF deltapt = dptelms(I) * (1.0 + pslope(KP) * (P - mtmdb(KP))) !interpolation deltae = elmsde(I) * (1.0 + eslope(KP) * (P - mtmdb(KP))) deltapl = p - dsqrt(dabs(p**2 - deltapt**2 + deltae**2 & - 2.0 * deltae * (ke + elmsmass) )) !interpolation check a2 = abs(pslope(KP) * (P - mtmdb(KP))) a6 = abs(eslope(KP) * (P - mtmdb(KP))) IF (P.lt.100.)THEN a3 = 0.10 !check at 10% for low momentum ELSE a3 = 0.03 !otherwise 3% ENDIF IF(a2.gt.a3.or.a6.gt.a3)then ! write(6,3464)P,a2,a6 write(2,3464)P,a2,a6 endif IF (UCINT.EQ.1) then UCPL1=deltapl UCPT1=deltapt UCen1=deltae UCang1=atan(deltapt/p) end if IF (UCINT.EQ.2) then UCPL2=deltapl UCPT2=deltapt UCen2=deltae UCang2=atan(deltapt/p) end if 30400 CONTINUE c make uncorrelated - check!!! deltapt= P*tan(UCang2) deltapl= p - dsqrt(dabs(p**2 - deltapt**2 + UCen1**2 & - 2.0 * UCen1 * (ke + elmsmass) )) pnew = dsqrt((p - deltapl)**2 + deltapt**2) IF(Pnew .lt.p*(1.0001))then GOTO 30310 end if write(2,*)'increase p',pnew,p GOTO 30300 30310 CONTINUE ENDIF !step shorten if it exits steplen = thkdb(KP) llenwt = (zend - z) * p / (pz * steplen) IF(llenwt.le.0.0)then write(2,*)'negative lenwt',a2,llenwt,mtmdb(Jmom) goto 30010 endif IF (llenwt .LT.1.0 .AND. pz .GT. 0.0) THEN deltapl = deltapl * llenwt deltapt = deltapt * llenwt steplen = steplen * llenwt ENDIF !combine this step, assuming that the energy loss and scattering occurs at the END of the step x = x + steplen * oldi y = y + steplen * oldj z = z + steplen * oldk ! t = t + steplen * dsqrt(ivel2) / c ! not used forward = sign(1d0, p - deltapl) pnew = dsqrt((p - deltapl)**2 + deltapt**2) IF(Pnew .gt.p*(1.0001))then write(2,*)'increase p',pnew,p endif p = pnew IF (p .LE. 5.0) THEN px = 0.0 py = 0.0 pz = 0.0 GOTO 30010 ELSE !unit vector in new direction ! call ranlux(randvw,1) randvw(1) = RAN1(irnarg) sran = randvw(1) * 2.0 * pigreek cphir = cos(sran) sphir = sin(sran) sthetas = min (1.0, deltapt / p) cthetas = dsqrt(1.0 - sthetas**2) * forward ctheta = oldk stheta = dsqrt(1.0 - ctheta**2) IF (stheta .GT. 1.0E-10) THEN sphi = oldi / stheta cphi = oldj / stheta ELSE sphi = 0.0 cphi = 1.0 ENDIF px = p * (cphi * sthetas * sphir + sphi * ctheta & * sthetas * cphir + sphi * stheta * cthetas) py = p * (- sphi * sthetas * sphir + cphi * ctheta & * sthetas * cphir + cphi * stheta * cthetas) pz = p * (- stheta * sthetas * cphir + & ctheta * cthetas) ke = dsqrt(p**2+elmsmass**2)-elmsmass !jan 2004*********************************** IF(dabs((px**2+py**2+pz**2)/p**2 - 1.0) .gt. 5.0E-9) & THEN write(2,*) "fail 1", dabs((px**2+py**2+pz**2)/p**2 - 1.0) ENDIF !REcheck at E-10 level IF(abs((oldi*px+oldj*py+oldk*pz)/p - cthetas) .gt. & 5.0E-9) stop "fail 2" !Rechecks at E-10 level ENDIF 30000 CONTINUE ! --------------------------------------------------------- 30010 pp(1)= px/1000 pp(2)= py/1000 pp(3)= pz/1000 call VMAG(pp,pm) pavg=pm e = SQRT(pm**2 + pmass**2) ekin=ke 800 continue 805 continue 900 continue return end ! ********************************************************* subroutine ELMS_INIT ! ! Set up arrays required by ELMS ! ! Based on a routine by Simon Holmes, 22 NOV 2004 ! implicit double precision (a-h,o-z) include 'icommon.inc' integer elb,plb,ne3,np3 logical done character*100 runnumber character*50 elmsdir,elmsroot ! c read in DB on first run only c elmsmass = 105.658 ! ! myownifrun = 0 OPEN(UNIT=12, file='elmscom.txt',status='unknown',iostat=ioc) if( ioc .ne. 0 ) then stop 'Couldnt open ELMSCOM.TXT' end if read(12,'(a50)') elmsdir read(12,'(a50)') elmsroot ! read (12,1368) ELMSCOR !correlations on (1) or off (0) ! read (12,1369) jint !random seed CLOSE (12) ! jint=234561 ! call rluxgo(4,jint,0,0) ! OPEN(UNIT=12, file='ELMSDB/RunDirectory.txt') open(unit=12,file=elmsdir(1:LASTNB(elmsdir)), & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(*,'(a,a)') 'ELMSDIR = ',elmsdir stop 'Couldnt open ELMS run directory' end if ! elmsfirst = .true. ! DO 30600 nmtm = 1,80 calculated(nmtm) = 0 iF (nmtm .EQ. 1)THEN firstpoint(1) = 1 ELSE firstpoint(nmtm) = lastpoint(nmtm-1) + 1 lastpoint(nmtm) = lastpoint(nmtm-1) ENDIF READ (12, 1335) mrun, mtmdb(nmtm), thkdb(nmtm), & colprobperm(nmtm), meanex(nmtm), meanp2x(nmtm) L = nmtm / 10 K = nmtm - 10 * L J = L/10 L = L - 10 * J ! runnumber = 'ELMSDB/ELMSFv3run' // char(48 + J) // ! & char(48 + L) // char(48+K) // '.elm' runnumber = elmsroot(1:LASTNB(elmsroot)) // char(48 + J) // & char(48 + L) // char(48+K) // '.elm' INQUIRE(FILE = runnumber, EXIST = done ) IF (.NOT. done) GOTO 30600 calculated(nmtm) = 1 OPEN (UNIT=13, ACCESS = 'SEQUENTIAL',file=runnumber,iostat=ioc) if( ioc .ne. 0 ) then write(*,'(a,a)') 'ELMS file = ',runnumber stop 'Couldnt open ELMS data file' end if LL = firstpoint(nmtm) READ(13,7832)IL, a1, a2, factorp, factore, pscale(2), & escale(2),ELB, ne3, PLB, np3 IF (IL.NE.nmtm) STOP 'ELMS db check fail' DO I = 2, maxe-1 escale(I+1) = escale(I) * factore pscale(I+1) = pscale(I) * factorp ENDDO pscale(PLB) = 0.0 escale(ELB) = 0.0 energy = dsqrt(mtmdb(nmtm)**2 + elmsmass**2) a1 = 0.0 a2 = 0.0 DO J=ELB,ne3 !read actual map. Store nonzero values as *10**8 integers, integrated !store mean pl (a1), mean de (a2), prob lostas being unphysical READ(13,2335)(map3(I),I=PLB,np3) DO 30700 I=PLB,np3 IF (map3(I) .LE. 1.0D-10) GOTO 30700 IF (LL .EQ. firstpoint(nmtm)) THEN probarray(LL) = map3(I) ELSE probarray(LL) = map3(I) + probarray(LL-1) ENDIF elmsde(LL) = 0.5*(escale(J)+escale(J+1)) dptelms(LL) = 0.5 * (pscale(I) + pscale(I+1)) lastpoint(nmtm) = LL !check physical a4 = (energy - elmsde(LL))**2 - dptelms(LL)**2 - elmsmass**2 !pl squared IF(a4.le.0.0)THEN !unphysical combination of pt and de. Note if sizeable IF(map3(I).GT.1.0E-9) write(2,*) 'unphysical', elmsde(LL), & dptelms(LL), a4, map3(I) GOTO 30700 ENDIF !accumulate checks a2 = a2 + elmsde(LL) * map3(I) a1 = a1 + dptelms(LL)**2 * map3(I) LL = LL + 1 30700 CONTINUE ENDDO CLOSE (13) normelms(nmtm) = probarray(lastpoint(nmtm)) a1 = a1 / (meanp2x(nmtm) * thkdb(nmtm)) a2 = a2 / (meanex(nmtm) * thkdb(nmtm)) N = lastpoint(nmtm)-firstpoint(nmtm)+1 write(2,7833) nmtm, mtmdb(nmtm), lastpoint(nmtm), & normelms(nmtm), a2, a1 30600 CONTINUE nmtm = 138 CLOSE(12) 1335 Format(I4,F13.3,E13.6,3E17.9) 7832 Format(I5,2E12.5,2F6.3,2E10.2,4I4) 7833 Format(I4, ' Mom', F8.0, I8, ' Norm', F10.7, & ' Ratio E, Pt2', 2F8.5) 2335 Format(10F13.10) 1334 Format(12E15.8) 1368 Format(I1) 1369 Format(I6) ! !Find the value of "eslope" and "pslope" for each momentum for interpolation J = 0 DO I = 1, nmtm IF (calculated(I).NE.1) CYCLE IF (J .NE. 0) THEN eslope(I) = (meanex(I)/meanex(J) -1.0)/(mtmdb(I) - mtmdb(J)) pslope(I) = (dsqrt(meanp2x(I)/meanp2x(J)) -1.0) & /(mtmdb(I) - mtmdb(J)) ENDIF J = I ENDDO !vector of momenta actually run DO I = 1, nmtm IF (calculated(I).NE.1)CYCLE npdone = npdone + 1 pdone(npdone) = I ENDDO !find max interpolations DO I = 1, npdone IF(I.EQ.1)CYCLE a1=0.5*eslope(pdone(I))*(mtmdb(pdone(I))-mtmdb(pdone(I-1))) a2=0.5*pslope(pdone(I))*(mtmdb(pdone(I))-mtmdb(pdone(I-1))) write(2,7145)mtmdb(pdone(I)),a1,a2 ! write(8,7145)mtmdb(pdone(I)),a1,a2 7145 Format(' Max mtm interpol at ',F6.0,' MeV in de,dpt',2F6.3) ENDDO ! end c ********************************************************* subroutine FANO(omega,chic2,pavg,th) c c Determine polar scattering angle from Fano distribution ! ! Ref. U. Fano, Phys. Rev. 93:117-120, 1954 ! Derived from GEANT program GMOLIE c C. * * C. * THRED(NA)=reduced angles of Moliere theory * C. * * C. * F0I(NA),F1I(NA),F2I(NA)= integrals of Moliere functions * C. * * implicit double precision (a-h,o-z) include 'icommon.inc' ! PARAMETER (ENEPER = 2.7182818d0) DIMENSION TINT(40),ARG(4),VAL(4),THRED(40),F0I(40),F1I(40),F2I(40) DATA THRED/ + 0.00, 0.10, 0.20, 0.30 +, 0.40, 0.50, 0.60, 0.70 +, 0.80, 0.90, 1.00, 1.10 +, 1.20, 1.30, 1.40, 1.50 +, 1.60, 1.70, 1.80, 1.90 +, 2.00, 2.20, 2.40, 2.60 +, 2.80, 3.00, 3.20, 3.40 +, 3.60, 3.80, 4.00, 5.00 +, 6.00, 7.00, 8.00, 9.00 +, 10.00,11.00,12.00,13.00/ DATA F0I/ + 0.000000E+00 ,0.995016E-02 ,0.392106E-01 ,0.860688E-01 + ,0.147856E+00 ,0.221199E+00 ,0.302324E+00 ,0.387374E+00 + ,0.472708E+00 ,0.555142E+00 ,0.632121E+00 ,0.701803E+00 + ,0.763072E+00 ,0.815480E+00 ,0.859142E+00 ,0.894601E+00 + ,0.922695E+00 ,0.944424E+00 ,0.960836E+00 ,0.972948E+00 + ,0.981684E+00 ,0.992093E+00 ,0.996849E+00 ,0.998841E+00 + ,0.999606E+00 ,0.999877E+00 ,0.999964E+00 ,0.999990E+00 + ,0.999998E+00 ,0.999999E+00 ,0.100000E+01 ,0.100000E+01 + ,0.100000E+01 ,0.100000E+01 ,0.100000E+01 ,0.100000E+01 + ,1.,1.,1.,1./ DATA F1I/ + 0.000000E+00,0.414985E-02,0.154894E-01,0.310312E-01 + ,0.464438E-01,0.569008E-01,0.580763E-01,0.468264E-01 + ,0.217924E-01,-0.163419E-01,-0.651205E-01,-0.120503E+00 + ,-0.178272E+00,-0.233580E+00,-0.282442E+00,-0.321901E+00 + ,-0.350115E+00,-0.366534E+00,-0.371831E+00,-0.367378E+00 + ,-0.354994E+00,-0.314803E+00,-0.266539E+00,-0.220551E+00 + ,-0.181546E+00,-0.150427E+00,-0.126404E+00,-0.107830E+00 + ,-0.933106E-01,-0.817375E-01,-0.723389E-01,-0.436650E-01 + ,-0.294700E-01,-0.212940E-01,-0.161406E-01,-0.126604E-01 + ,-0.102042E-01,-0.840465E-02,-0.704261E-02,-0.598886E-02/ DATA F2I/ + 0.0,0.121500E-01,0.454999E-01,0.913000E-01 + ,0.137300E+00,0.171400E+00,0.183900E+00,0.170300E+00 + ,0.132200E+00,0.763000E-01,0.126500E-01,-0.473500E-01 + ,-0.936000E-01,-0.119750E+00,-0.123450E+00,-0.106300E+00 + ,-0.732800E-01,-0.312400E-01,0.128450E-01,0.528800E-01 + ,0.844100E-01,0.114710E+00,0.106200E+00,0.765830E-01 + ,0.435800E-01,0.173950E-01,0.695001E-03,-0.809500E-02 + ,-0.117355E-01,-0.125449E-01,-0.120280E-01,-0.686530E-02 + ,-0.385275E-02,-0.231115E-02,-0.147056E-02,-0.982480E-03 + ,-0.682440E-03,-0.489715E-03,-0.361190E-03,-0.272582E-03/ * * ------------------------------------------------------------------ * IF(OMEGA.LE.ENEPER)GO TO 90 CHIC = SQRT(chic2) CNST=LOG(OMEGA) ! Moliere b parameter bet = pavg/pmass / SQRT(1.+(pavg/pmass)**2) z = mparm(1,krad) if( NINT(z) .eq. 1 ) then uin = 3.6 ! Fano p.118 else if( NINT(z) .eq. 2 ) then uin = 4.2 ! numerical integration of Hubbell S(x) data else if( NINT(z) .eq. 3 ) then uin = 4.6 ! Fano p.118 else if( NINT(z) .eq. 4 ) then uin = 4.3 ! numerical integration of Hubbell S(x) data else ! Z > 4 uin = 5. end if if( mt .eq. 'LIH' ) then uin = 4.1 ! use average of Li and H end if ! if( ABS(iptyp) .eq. 1 ) then ! electron argm = 0.160 * z**(-0.67) * (1.+3.87e-3*z/bet) db = 1./(z+1.)*(LOG(argm)+uin) else ! heavy particle argm = 1130.*z**(-1.33)/(1./bet**2-1.) db = 1./z*(LOG(argm)+uin-0.5*bet*bet) end if cnst = cnst + db ! --------------------------------------------------------- B=5. c DO 10 L=1,10 IF(ABS(B).LT.1.E-10)THEN B=1.E-10 ENDIF DB=-(B-LOG(ABS(B))-CNST)/(1.-1./B) B=B+DB IF(ABS(DB).LE.0.0001)GO TO 20 10 CONTINUE c GO TO 90 c 20 CONTINUE IF(B.LE.0.)GO TO 90 BINV = 1./B TINT(1) = 0. c DO 30 JA=2,4 TINT(JA)=F0I(JA)+(F1I(JA)+F2I(JA)*BINV)*BINV 30 CONTINUE c NMAX = 4 40 CONTINUE XINT = RAN1(irnarg) c DO 50 NA=3,40 IF(NA.GT.NMAX) THEN TINT(NA)=F0I(NA)+(F1I(NA)+F2I(NA)*BINV)*BINV NMAX=NA ENDIF IF(XINT.LE.TINT(NA-1)) GO TO 60 50 CONTINUE c IF(XINT.LE.TINT(40)) THEN NA=40 GOTO 60 ELSE TMP=1.-(1.-B*(1.-XINT))**5 IF(TMP.LE.0.)GO TO 40 THRI=5./TMP GO TO 80 ENDIF c 60 CONTINUE NA = MAX(NA-1,3) NA3 = NA-3 c DO 70 M=1,4 NA3M=NA3+M ARG(M)=TINT(NA3M) VAL(M)=THRED(NA3M)**2 70 CONTINUE c call POLINT(arg,val,4,xint,thri,err) ! Num Rec c 80 CONTINUE TH=CHIC*SQRT(ABS(B*THRI)) IF(TH.GT.PI)GO TO 40 SINTH=SIN(TH) TEST=TH*(RAN1(irnarg))**2 IF(TEST.GT.SINTH)GO TO 40 GOTO 900 ! 90 CONTINUE iflag = -24 c 900 continue END ! ********************************************************* subroutine FLIP_COULOMB(p,h,pr) ! ! Compute probability of muon spin flip in mu-Z Coulomb scattering ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 me,mu ! data con/4.375e-8/,alp/7.299e-3/ ! me = mass(1) ! [GeV] mu = mass(2) e = SQRT( p*p + mu*mu ) b = p / e z = mparm(1,krad) a = aelcomp(1,krad) r = mparm(2,krad) ! [g cm^-3] thmin = alp*z**(0.333333)*me/p thmax = 2./alp/a**(0.333333)*me/p arg = (1.-COS(thmax)) / (1.-COS(thmin)) ! pr = con*h*r/a*(z/b/p/e)**2 * LOG( arg ) hel = HELICITY( pol ) pr = pr * ABS( hel ) ! end ! ********************************************************* subroutine FLIP_ELASTIC(s,z,pr) ! ! Compute probability of muon spin flip in mu-e elastic scattering ! Adapted from J. Gallardo 10 April 2000 ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension s1(4),s2(4),q1(4),q2(4),p1(4),p2(4) ! ! xme=0.51099907d-3 ! mass in GeV ! xmu=105.658389d-3 xme = mass(1) xmu = mass(2) ! if( ABS(z) .le. 1.d0 ) then sinz=sqrt((1.d0-z**2)) else iflag = -64 go to 900 end if ! c cm of momentum frame of reference emu1cm=0.5*(s+xmu**2-xme**2)/sqrt(s) qmu1cm=sqrt(emu1cm**2-xmu**2) t=-2.*qmu1cm**2*(1.d0-z) ! q1(1)=emu1cm q1(2)= 0.d0 q1(3)= 0.d0 q1(4)=qmu1cm p1(1)=0.5*(s+xme**2-xmu**2)/sqrt(s) p1(2)= 0.d0 p1(3)= 0.d0 p1(4)=-qmu1cm c q2(1)=emu1cm q2(2)=sinz*qmu1cm q2(3)= 0.d0 q2(4)=z*qmu1cm p2(1)=0.5*(s+xme**2-xmu**2)/sqrt(s) p2(2)=-sinz*qmu1cm p2(3)= 0.d0 p2(4)=-z*qmu1cm ! s1(1)=qmu1cm/xmu s1(2)= 0.d0 s1(3)= 0.d0 s1(4)=q1(1)/xmu s2(1)=qmu1cm/xmu s2(2)=emu1cm*sinz/xmu s2(3)= 0.d0 s2(4)=emu1cm*z/xmu ! denominator = 32*((s-(xme**2+xmu**2))**2+0.5*t**2+s*t) c xnumerator = -32*dot(s1,s2)*((s-(xmu**2+xme**2))**2+s*t)+ + 64*s*(dot(p1,s2)*dot(q2,s1)+dot(p2,s1)*dot(q1,s2)+ + dot(q1,s2)*dot(q2,s1))+ + 64*t*(dot(p1,s2)*dot(p2,s1)+dot(p1,s2)*dot(q2,s1)+ + dot(p2,s1)*dot(q1,s2)+dot(q1,s2)*dot(q2,s1))- + 64*(dot(p1,s2)*dot(q2,s1)*(xmu**2+xme**2)+dot(p2,s1)* + dot(q1,s2)*(xmu**2+xme**2)+ + (xmu**2+2*xme**2)*dot(q1,s2)*dot(q2,s1)) ! polar=xnumerator/denominator ! polarization pr =0.5*(1.0d0-polar) ! Prob(up --->down) f = 0.5d0*(1.0d0 - z) ! spin flip vs helicity flip pr = f - pr hel = HELICITY( pol ) pr = pr * ABS( hel ) ! 900 continue end ! ********************************************************* REAL*8 FUNCTION FNEAR(C) ! ! Auxiliary routine for SAMCS ! IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) common/rans/iseed ! REAL RNDM ! 1 RNF=DBLE(RNDM(-1.)) 1 RNF = RAN1(iseed) F1=1.D0-C/(1.D0+2.D0*C) IF(RNF.LE.F1) THEN ! RNF=DBLE(RNDM(-1.)) RNF = RAN1(iseed) X=-LOG(RNF) ELSE ! RNF1=DBLE(RNDM(-1.)) ! RNF2=DBLE(RNDM(-1.)) ! RNF3=DBLE(RNDM(-1.)) RNF1 = RAN1(iseed) RNF2 = RAN1(iseed) RNF3 = RAN1(iseed) X=-LOG(RNF1*RNF2*RNF3) ENDIF F2=1.-2.D0*C*X/(1.D0+C+C*X**2/2.D0) ! RNF=DBLE(RNDM(-1.)) RNF = RAN1(iseed) IF(RNF-F2) 2,1,1 2 FNEAR=X RETURN END ! ********************************************************* DOUBLE PRECISION FUNCTION FORMN(QG,Z,A) ! ! Auxiliary routine for SAMCS ! C ------------------------------------------------------------------------------------ C Electromagnetic form factor of target atom squared C ------------------------------------------------------------------------------------ C QG - momentum transfer in GeV/c C ------------------------------------------------------------------------------------- C Z,A - charge and atomic mass of target C ------------------------------------------------------------------------------------- C Data from: C H. de Vries, C.W. de Jager and C. de Vries C Atomic Data and Nuclear Data Tables 36 (1987) 495-536 C C V.V. Burov et al. Phys. Atom. Nucl. 61 (1998) 526-532 C ----------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) DIMENSION AHO(8) PARAMETER (SC=5.06765D0) PARAMETER (PI=3.141592653589793227D+00) DATA AHO/0.85d0,0.,1.77d0,1.78d0,1.69d0,0.,1.745d0,1.833d0/ ! Q=QG*SC C Gaussian formfactor - for comparizon with analitical calculation c if(z.le.300.d0) then c r2=rav2(z,a) c formn=exp(-q**2*r2/3.d0) c return c endif IF(Z.EQ.1.D0) THEN FORMN = 1.D0/(1.D0+(Q*AHO(1))**2/12.D0)**4 return endif IF(Z.eq.2.D0) then R=1.238 B=0.38 FORMN = ZFORM(Q,R,B)**2 return endif IF(Z.eq.6.D0) then R=2.214 B=0.488 FORMN = ZFORM(Q,R,B)**2 return endif IF(Z.LE.8.D0) then IZ=Z B2=AHO(IZ)**2 A0=(B2-.4824)*A/(A-1.) Z1=1.-(Z-2.)*A0*Q**2/6./Z FORMN=(Z1*EXP(-B2*Q**2/4.))**2 ELSE A13=A**(1./3.) R=1.141*A13-1.047/A13 R2=(0.898*A13+1.198/A13)**2 B=R*SQRT(3.*(5.*R2/3./R**2-1.)/7.)/PI FORMN = ZFORM(Q,R,B)**2 ENDIF RETURN END ! ********************************************************* DOUBLE PRECISION FUNCTION FORMP(QG,IP) ! ! Auxiliary routine for SAMCS ! C ------------------------------------------------------------------------------------ C Electromagnetic form factor of projectile squared C ------------------------------------------------------------------------------------ IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) C --------------- GeV to 1/fermi Q=QG*5.06765D0 IF(IP.LE.13) THEN C --------------- Dipole form factor for hadrons FORMP = 1.D0/(1.D0+Q**2*R2P(IP)/12.D0)**4 ELSE C --------------- Gaussian form factor for light nuclei FORMP=EXP(-Q**2*R2P(IP)/3.D0) ENDIF RETURN END ! ********************************************************* SUBROUTINE GDECA2(XM0,XM1,XM2,PCM) C. C. ****************************************************************** C. * * C. * Simulates two body decay process with isotropic angular * C. * distribution in CMS. * C. * * C. * Author G.Patrick ********* * C. * * C. ****************************************************************** C. ! Adapted from Geant 3.2.1 ! implicit double precision (a-h,o-z) DIMENSION PCM(4,2) DIMENSION RNDM(2) C. data twopi/6.2831853d0/ C. ------------------------------------------------------------------ C. C Generate first decay product in CMS. C E1=(XM0*XM0+XM1*XM1-XM2*XM2)/(2d0*XM0) P1=SQRT(ABS((E1-XM1)*(E1+XM1))) C C Isotropic decay angular distribution. C CALL GRNDM(RNDM,2) COSTH=2d0*RNDM(1)-1d0 IF(ABS(COSTH).GE.1d0) THEN COSTH=SIGN(1.d0,COSTH) SINTH=0d0 ELSE SINTH=SQRT((1d0-COSTH)*(1d0+COSTH)) ENDIF PHI=TWOPI*RNDM(2) ! sinth=0.173648 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! costh=0.984808 ! phi=0. C C Polar co-ordinates to momentum components. C PCM(1,1)=P1*SINTH*COS(PHI) PCM(2,1)=P1*SINTH*SIN(PHI) PCM(3,1)=P1*COSTH PCM(4,1)=E1 C C Generate second decay product. C PCM(1,2)=-PCM(1,1) PCM(2,2)=-PCM(2,1) PCM(3,2)=-PCM(3,1) PCM(4,2)=SQRT(PCM(1,2)**2+PCM(2,2)**2+PCM(3,2)**2+XM2*XM2) C END ! ********************************************************* SUBROUTINE GDECA3(XM0,XM1,XM2,XM3,PCM) C. C. ****************************************************************** C. * * C. * This routine simulates three-body decays in center-of-mass. * C. * Written by Jurgen Schukraft and Vincent Hedberg for R807. * C. * Adapted for Geant3 by Sverker Johansson (Aug-1987). * C. * * C. * As it is, the decay is done according to pure phase-space. * C. * There is code in the routine to properly treat K decays (V-A)* C. * but it is not used. (Requires an additional parameter * C. * in the call, to tell which type of decay.) * C. * * C. * Author S. Johansson ********* * C. * * C. ****************************************************************** C. ! Adapted from Geant 3.2.1 ! implicit double precision (a-h,o-z) DIMENSION P(20,20) DIMENSION HM(4,7) DIMENSION PCM(4,3) DIMENSION RNDM(3) * **************************************************** * Translation from new to old input parameters : * ( S.J. ) **************************************************** ID = 1 HM(1,ID) = XM1 HM(2,ID) = XM2 HM(3,ID) = XM3 HM(4,ID) = XM0 K1 = 1 K2 = 2 K3 = 3 PI2 = 2d0 * 3.141592654d0 *********************** * Original code : C-------------------------------------------------------- C SIMULATION OF SECONDARY PARTICLE MOMENTA C IN THE REST FRAME OF PARENT PARTICLE C-------------------------------------------------------- T0=HM(4,ID)-HM(1,ID)-HM(2,ID)-HM(3,ID) 10 CALL GRNDM(RNDM,2) G1=RNDM(1) G2=RNDM(2) GMI=MIN(G1,G2) GMA=MAX(G1,G2) TE=GMI*T0 TM=(1d0-GMA)*T0 TN=(GMA-GMI)*T0 PA1=SQRT(TE**2+2d0*HM(1,ID)*TE) PA2=SQRT(TM**2+2d0*HM(2,ID)*TM) PA3=SQRT(TN**2+2d0*HM(3,ID)*TN) C-------------------------------------------------------- C CHECK OF THE MOMENTUM CONSERVATION C-------------------------------------------------------- PMMM=MAX(PA1,PA2) PMAX=MAX(PA3,PMMM) SP=PA1+PA2+PA3 IF(PMAX.GE.SP-PMAX) GOTO 10 P(4,K1)=TE+HM(1,ID) P(4,K2)=TM+HM(2,ID) P(4,K3)=TN+HM(3,ID) C-------------------------------------------------------- C ACCOUNTING FOR MATRIX ELEMENT (FOR ID=4 AND 5) * (** Calculates correct (V-A) decay for kaons. **) * (** Now inactive. Change ID to activate it. (S.J) **) C-------------------------------------------------------- IF(ID.NE.4.AND.ID.NE.5) GOTO 20 EEM=((HM(4,ID)-HM(2,ID))**2+HM(1,ID)**2- * HM(3,ID)**2)/(2d0*(HM(4,ID)-HM(2,ID))) HMMAX=EEM*(HM(4,ID)-HM(2,ID)-EEM)+ * 0.5d0*(EEM**2-HM(1,ID)**2) HMA=P(4,K1)*P(4,K3)+0.25d0*(PA1**2+PA3**2-PA2**2) CALL GRNDM(RNDM,1) GG=RNDM(1) IF(HMA.LT.GG*HMMAX) GOTO 10 C-------------------------------------------------------- C CALCULATIONS OF MOMENTA C-------------------------------------------------------- 20 CONTINUE CALL GRNDM(RNDM,3) CTE=2d0*RNDM(1)-1d0 STE=SQRT(ABS(1d0-CTE**2)) FE =PI2*RNDM(2) CFE=COS(FE) SFE=SIN(FE) P(1,K1)=PA1*STE*CFE P(2,K1)=PA1*STE*SFE P(3,K1)=PA1*CTE CTEN=(PA2**2-PA1**2-PA3**2)/(2d0*PA1*PA3) STEN=SQRT(ABS(1d0-CTEN**2)) FEN =PI2*RNDM(3) CFEN=COS(FEN) SFEN=SIN(FEN) P(1,K3)=PA3*(STEN*CFEN*CTE*CFE-STEN*SFEN*SFE+CTEN*STE*CFE) P(2,K3)=PA3*(STEN*CFEN*CTE*SFE+STEN*SFEN*CFE+CTEN*STE*SFE) P(3,K3)=PA3*(-STEN*CFEN*STE+CTEN*CTE) P(1,K2)=-P(1,K1)-P(1,K3) P(2,K2)=-P(2,K1)-P(2,K3) P(3,K2)=-P(3,K1)-P(3,K3) C------------------------------------------------------- *** Re-translation from old to new output arrays : *** *** (S.J.) *** ***--------------------------------------------------*** DO 37 I=1,4 DO 38 J=1,3 PCM(I,J) = P(I,J) 38 CONTINUE 37 CONTINUE END c ********************************************************* SUBROUTINE GDROT(P,COSTH,SINTH,COSPH,SINPH) C. C. ****************************************************************** C. * * C. * Rotates vector from one reference system into another. * C. * THETA and PHI are anti-clockwise Eulerian angles between the * C. * two systems. * C. * * C. * Author M.Hansroul, G.Patrick ********* * C. * * C. * * C. ****************************************************************** C. implicit double precision (a-h,o-z) DIMENSION P(3) C. C. ------------------------------------------------------------------ C. P1=P(1) P2=P(2) P3=P(3) P(1)=P1*COSTH*COSPH - P2*SINPH + P3*SINTH*COSPH P(2)=P1*COSTH*SINPH + P2*COSPH + P3*SINTH*SINPH P(3)=-P1*SINTH + P3*COSTH C END ! ********************************************************* SUBROUTINE GLOREN(BETA,PA,PB) C. C. ****************************************************************** C. * * C * Routine to transform momentum and energy from the * C * Lorentz frame A to the Lorentz frame B * C * * C * PA(1) * C * PA(2) Momentum components in frame A * C * PA(3) * C * PA(4) Energy * C * PB(..) same quantities in frame B * C * * C * BETA(1) Components of velocity of frame B * C * BETA(2) as seen from frame A * C * BETA(3) * C * BETA(4) 1./SQRT(1.-BETA**2) * C. * * C. * Author M.Hansroul ********* * C. * * C. ****************************************************************** C. implicit double precision (a-h,o-z) DIMENSION BETA(4),PA(4),PB(4) C. C. ------------------------------------------------------------------ C. ! gam = beta(4) ! bg = beta(3)*beta(4) ! pb(1) = pa(1) ! pb(2) = pa(2) ! pb(3) = gam*pa(3) - bg*pa(4) ! pb(4) = gam*pa(4) - bg*pa(3) ! return ! BETPA = BETA(1)*PA(1) + BETA(2)*PA(2) + BETA(3)*PA(3) BPGAM = (BETPA * BETA(4)/(BETA(4) + 1d0) - PA(4)) * BETA(4) PB(1) = PA(1) + BPGAM * BETA(1) PB(2) = PA(2) + BPGAM * BETA(2) PB(3) = PA(3) + BPGAM * BETA(3) PB(4) =(PA(4) - BETPA) * BETA(4) END c ********************************************************* SUBROUTINE GPOISS(AMVEC,NPVEC,LEN) ! ! Poisson random number interface for GEANT routines ! C. C. ****************************************************************** C. * * C. * Generates a vector NPVEC of LEN random numbers * C. * POISSON distribued with mean values AMVEC * C. * * C. * If the mean value A greater than PLIM, N is calculated * C. * according to the Gaussian approximation of the Poisson * C. * distribution. * C. * * C. * * C. * Author : L.Urban * C. * Date : 28.04.1988 Last update : 1.02.1990 * C. * * C. ****************************************************************** C. ! implicit double precision (a-h,o-z) include 'icommon.inc' ! REAL*8 AMVEC(*),RNDM(2), N INTEGER NPVEC(*) PARAMETER (PLIM=16.,HMXINT=2E+9) ! ! do i=1,len ! xm = amvec(i) ! npvec(i) = POIDEV(xm,idum) ! end do ! DO 30 I=1,LEN * Protection against negative mean values N=0. IF(AMVEC(I).GT.0.) THEN IF(AMVEC(I).LE.PLIM) THEN CALL GRNDM(RNDM,1) R=RNDM(1) P=EXP(-AMVEC(I)) S=P IF(R.LE.S) GOTO 20 10 N=N+1. P=P*AMVEC(I)/N S=S+P IF(S.LT.R.AND.P.GT.1.E-30) GOTO 10 ELSE CALL GRNDM(RNDM,2) RR=SQRT(-2.*LOG(RNDM(1))) PHI=TWOPI*RNDM(2) X=RR*COS(PHI) N=MIN(MAX(AMVEC(I)+X*SQRT(AMVEC(I)),0.d0),HMXINT) ENDIF ENDIF * 20 NPVEC(I) = N 30 CONTINUE * END c ********************************************************* subroutine GRNDM(rvec,len) ! ! Uniform random number interface for GEANT routines ! implicit double precision (a-h,o-z) real*8 rvec(*) ! do i=1,len rvec(i) = RANV(0.d0,1.d0) end do ! return end c ********************************************************* real*8 function HELICITY(s) ! ! Compute current muon helicity in the LAB frame ! s is the current value of particle spin vector in muon REST frame ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 s(3),p1(3) ! call VMAG(pp,pm) ! do i=1,3 p1(i) = pp(i) / pm ! unit vector along p end do ! helicity = DOT3(s,p1) ! return end ! ********************************************************* function ILVSDMI(a,n,inc) ! ! Find index of smallest element in array A(N) ! Replaces CERNLIB routine LVSDMI ! Used by SAMCS ! implicit double precision (a-h,o-z) dimension a(n) ! ilvsdmi = 1 amax = 1d12 do i=1,n,inc if( a(i) .lt. amax ) then amax = a(i) ilvsdmi = i end if end do ! end ! ********************************************************* function ILVSDMX(a,n,inc) ! ! Find index of largest element in array A(N) ! Replaces CERNLIB routine LVSDMX ! Used by SAMCS ! implicit double precision (a-h,o-z) dimension a(n) ! ilvsdmx = 1 amin = -1d12 do i=1,n,inc if( a(i) .gt. amin ) then amin = a(i) ilvsdmx = i end if end do ! end c ********************************************************* subroutine INTERACTION(hstp,pavg) ! ! Check for inelastic nuclear interactions of hadrons ! implicit double precision (a-h,o-z) include 'icommon.inc' ! if( ABS(iptyp) .lt. 3 ) go to 900 ! skip e and mu p4 = pavg sig = CS_INEL(p4,iptyp) * 1e-27 ! get inelastic cross section if( sig .eq. 0. ) go to 900 a = aelcomp(1,krad) xint = 0.01 / ( densat * sig * a**(0.71) ) ! [m] xm = hstp / xint nin = POIDEV( xm,irnarg ) if( nin .ge. 1 ) then if( intlev.eq.2 .and. iptyp.eq.5 ) then call WANG else iflag = -73 end if end if ! 900 continue end c ********************************************************* logical function IN_ASPWEDGE(hstp) ! ! Determine if particle step goes thru an ! Azimuthally Symmetric Polynomial shaped wedge ! ! Describe polynomial shape in 1st quadrant of coordinate system ! centered on the wedge: r(z) = a0+a1*z+a2*z**2+a3*z**3 ! ! The wedge body has reflection symmetry thru x-y plane ! implicit double precision (a-h,o-z) include 'icommon.inc' ! logical in0,in1 ! in_aspwedge = .false. if( mgt .ne. 'ASPW' ) go to 900 ! in0 = .false. in1 = .false. dr0 = 0. dr1 = 0. ! zv = sr0 + gparm(1,krad) ! z of wedge center zcut = gparm(2,krad) ! z offset from center to edge a0 = gparm(3,krad) a1 = gparm(4,krad) a2 = gparm(5,krad) a3 = gparm(6,krad) ! rp = SQRT( xp(1)**2 + xp(2)**2 ) zmin = zv - zcut zmax = zv + zcut ! ! Check end point ! n.b. xp contains coordinates at end of the time step if( xp(3).gt.zmin .and. xp(3).lt.zmax ) then dz = ABS( xp(3) - zv ) rsurf = a0 + a1*dz + a2*dz**2 + a3*dz**3 dr1 = rp - rsurf if( dr1 .le. 0. ) in1 = .true. end if ! ! Check starting point x1 = xp(1) - pp(1)/pp(3)*hdid y1 = xp(2) - pp(2)/pp(3)*hdid z1 = xp(3) - hdid rp = SQRT( x1*x1 + y1*y1 ) ! if( z1.gt.zmin .and. z1.lt.zmax ) then dz = ABS( z1 - zv ) rsurf = a0 + a1*dz + a2*dz**2 + a3*dz**3 dr0 = rp - rsurf if( dr0 .le. 0. ) in0 = .true. end if ! if( in0 .and. in1 ) then in_aspwedge = .true. ! else if( in0 .and. (.not. in1) ) then in_aspwedge = .true. hstp = hstp * ABS(dr0) / ( ABS(dr0)+ABS(dr1) ) ! else if( in1 .and. (.not. in0) ) then in_aspwedge = .true. hstp = hstp * ABS(dr1) / ( ABS(dr0)+ABS(dr1) ) end if ! 900 continue return end c ********************************************************* logical function IN_ASRWEDGE(hstp) ! ! Determine if particle step goes thru an ! Azimuthally Symmetric, Radial Wedge ! (defines thickness as a function of R) ! ! Wedge is symmetric about xy plane at z=zv ! Describe polynomial shape in terms of half-thickness ! (distance of edge from zv) ! dz(r) = b0+b1*r+b2*r**2+b3*r**3 ! dz < 0 means no wedge at that radius ! implicit double precision (a-h,o-z) include 'icommon.inc' ! logical in0,in1 ! in_asrwedge = .false. if( mgt .ne. 'ASRW' ) go to 900 ! in0 = .false. in1 = .false. dr0 = 0. dr1 = 0. ! zv = sr0 + gparm(1,krad) ! z of center of wedge zcut = gparm(2,krad) ! maximum allowed half-thickness b0 = gparm(3,krad) b1 = gparm(4,krad) b2 = gparm(5,krad) b3 = gparm(6,krad) ! ! Check end point ! n.b. xp contains coordinates at end of the time step rp = SQRT( xp(1)**2 + xp(2)**2 ) dzsurf = b0 + b1*rp + b2*rp**2 + b3*rp**3 ! dzsurf <0 means no wedge at this radius if (dzsurf .lt. 0.) dzsurf = 0. ! max thickness zcut if (dzsurf .gt. zcut) dzsurf = zcut zmin = zv - dzsurf zmax = zv + dzsurf ! if( xp(3).gt.zmin .and. xp(3).lt.zmax ) then ! if inside, find out how far from left edge (in case crossing occurred) in1 = .true. dz1 = xp(3) - zmin else ! if outside, find out how far from right edge (in case crossing occurred) dz1 = zmax - xp(3) end if ! ! Check starting point x1 = xp(1) - pp(1)/pp(3)*hdid y1 = xp(2) - pp(2)/pp(3)*hdid z1 = xp(3) - hdid rp = SQRT( x1*x1 + y1*y1 ) dzsurf = b0 + b1*rp + b2*rp**2 + b3*rp**3 ! dzsurf <0 means no wedge at this radius if (dzsurf .lt. 0.) dzsurf = 0. ! max thickness zcut if (dzsurf .gt. zcut) dzsurf = zcut zmin = zv - dzsurf zmax = zv + dzsurf ! if( z1.gt.zmin .and. z1.lt.zmax ) then ! if inside, find out how far from right edge (in case crossing occurred) in0 = .true. dz0 = zmax - z1 else ! if outside, find out how far from left edge (in case crossing occurred) dz0 = z1 - zmin end if ! if( in0 .and. in1 ) then in_asrwedge = .true. ! else if( in0 .and. (.not. in1) ) then in_asrwedge = .true. hstp = hstp * ABS(dz0) / ( ABS(dz0)+ABS(dz1) ) ! else if( in1 .and. (.not. in0) ) then in_asrwedge = .true. hstp = hstp * ABS(dz1) / ( ABS(dz0)+ABS(dz1) ) end if ! 900 continue return end c ********************************************************* logical function IN_HEMIWIN() ! ! Determine how a particle step proceeds thru a ! hemispherical end region of an absorber ! ! IHWIN returned as current ICOOL "radial" region ! Returns IFLAG=1 if the step crosses a region boundary ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 dp(3),dz(3),x0(3),x(3) ! data eps/2.e-6/ ! in_hemiwin = .false. if( mgt .ne. 'HWIN' ) go to 900 ! i0 = 0 ! subregion of initial step position i1 = 0 ! subregion of final step position ! ......................................................... ! Set geometry parameters rs1 = gparm(2,1) ! inner spherical radius rs2 = rs1 + gparm(3,1) ! outer spherical radius dzc = gparm(4,1) ! offset of center of shell r3 = rhigh(1) ! maximum cylindrical radius dz(1) = rs1 - dzc ! absorber thickness dz(2) = gparm(3,1) ! window thickness dz(3) = slen - rs2 + dzc ! exterior thickness ! if( gparm(1,1) .lt. 0. ) then ! ENTRANCE end zoff = slen else ! EXIT end zoff = 0. end if ! ! ......................................................... ! Check location of starting point ! Use coordinate system centered on the spherical shell x0(1) = xp0(1) x0(2) = xp0(2) x0(3) = zoff + dzc + gparm(1,1)*(xp0(3)-sr0) rsp0 = SQRT( x0(1)**2 + x0(2)**2 + x0(3)**2 ) ! if( rsp0 .lt. rs1 ) then i0 = 1 else if( rsp0.ge.rs1 .and. rsp0.le.rs2 ) then i0 = 2 else rp = SQRT( xp0(1)**2 + xp0(2)**2 ) if( rp .lt. r3 ) then i0 = 3 else go to 900 end if end if ! ......................................................... ! Check location of end point x(1) = xp(1) x(2) = xp(2) x(3) = zoff + dzc + gparm(1,1)*(xp(3)-sr0) rsp = SQRT( x(1)**2 + x(2)**2 + x(3)**2 ) ! if( rsp .lt. rs1 ) then i1 = 1 else if( rsp.ge.rs1 .and. rsp.le.rs2 ) then i1 = 2 else rp = SQRT( xp(1)**2 + xp(2)**2 ) if( rp .lt. r3 ) then i1 = 3 else go to 900 end if end if ! ......................................................... if( iflag .eq. 1 ) then ! this is the correction step ! Prepare stepping for the NEXT region ihwin = i0 + NINT( gparm(1,1) ) hregn = dz(ihwin) / 10 iflag = 0 go to 800 end if ! ihwin = i0 ! default values for current region hregn = dz(i0) / 10. iflag = 0 ! ! --------------------------------------------------------- ! Step handling logic ! There are 18 cases to consider, including jumping over windows altogether ! if( gparm(1,1) .lt. 0. ) then ! ENTRANCE end if( i0.eq.3 .and. i1.eq.3 ) then diff = ABS( rsp - rs2 ) if( diff .le. eps ) then ! step ends on 2nd sphere boundary ihwin = 2 hregn = dz(2) / 10. ! step size for NEXT step end if go to 800 ! good step ! else if( i0.eq.3 .and. i1.eq.2 ) then diff0 = ABS( rsp0 - rs2 ) diff = ABS( rsp - rs2 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 2nd sphere boundary ihwin = 2 hregn = dz(2) / 10. if( diff0 .lt. eps ) i0 = 2 ! get right material go to 800 ! good step end if r2 = rs2**2 ! set radius for step adjustment logic ! else if( i0.eq.3 .and. i1.eq.1 ) then diff0 = ABS( rsp0 - rs1 ) if( diff0 .le. eps ) then ! step starts on 2nd sphere boundary ihwin = 2 hregn = dz(2) / 10. iflag = 1 go to 900 ! redo step end if r2 = rs2**2 ! set radius for step adjustment logic ! else if( i0.eq.2 .and. i1.eq.2 ) then diff = ABS( rsp - rs1 ) if( diff .le. eps ) then ! step ends on 1st sphere boundary ihwin = 1 hregn = dz(1) / 10. ! step size for NEXT step end if go to 800 ! good step ! else if( i0.eq.2 .and. i1.eq.1 ) then diff0 = ABS( rsp0 - rs1 ) diff = ABS( rsp - rs1 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 1st sphere boundary ihwin = 1 hregn = dz(1) / 10. if( diff0 .lt. eps ) i0 = 1 ! get right material go to 800 ! good step end if r2 = rs1**2 ! set radius for step adjustment logic ! else if( i0.eq.1 .and. i1.eq.1 ) then go to 800 ! good step ! ! It's also possible for trajectories to cross subregions backwards else if( i0.eq.1 .and. i1.eq.2 ) then diff0 = ABS( rsp0 - rs1 ) diff = ABS( rsp - rs1 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 1st sphere boundary ! Use the small stepping appropriate to region 2 ihwin = 2 hregn = dz(2) / 10. if( diff0 .lt. eps ) i0 = 2 ! use right material go to 800 ! good step end if r2 = rs1**2 ! set radius for step adjustment logic ! else if( i0.eq.2 .and. i1.eq.3 ) then diff0 = ABS( rsp0 - rs2 ) diff = ABS( rsp - rs2 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 2nd sphere boundary ihwin = 2 hregn = dz(2) / 10. if( diff0 .lt. eps ) i0 = 3 ! use right material go to 800 ! good step end if r2 = rs2**2 ! set radius for step adjustment logic ! else if( i0.eq.1 .and. i1.eq.3 ) then diff0 = ABS( rsp0 - rs1 ) if( diff0 .le. eps ) then ! step starts on 1st sphere boundary ihwin = 2 hregn = dz(2) / 10. iflag = 1 go to 900 ! redo step end if r2 = rs1**2 ! set radius for step adjustment logic end if ! end of test on i0 and i1 ! ! else ! EXIT end if( i0.eq.1 .and. i1.eq.1 ) then diff = ABS( rsp - rs1 ) if( diff .le. eps ) then ! step ends on 1st sphere boundary ihwin = 2 hregn = dz(2) / 10. ! step size for NEXT step end if go to 800 ! good step ! else if( i0.eq.1 .and. i1.eq.2 ) then diff0 = ABS( rsp0 - rs1 ) diff = ABS( rsp - rs1 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 1st sphere boundary ihwin = 2 hregn = dz(2) / 10. if( diff0 .lt. eps ) i0 = 2 ! get right material go to 800 ! good step end if r2 = rs1**2 ! set radius for step adjustment logic ! else if( i0.eq.1 .and. i1.eq.3 ) then diff0 = ABS( rsp0 - rs1 ) if( diff0 .le. eps ) then ! step starts on 1st sphere boundary ihwin = 2 hregn = dz(2) / 10. iflag = 1 go to 900 ! redo step end if r2 = rs1**2 ! set radius for step adjustment logic ! else if( i0.eq.2 .and. i1.eq.2 ) then diff = ABS( rsp - rs2 ) if( diff .le. eps ) then ! step ends on 2nd sphere boundary ihwin = 3 hregn = dz(3) / 10. ! step size for NEXT step end if go to 800 ! good step ! else if( i0.eq.2 .and. i1.eq.3 ) then diff0 = ABS( rsp0 - rs2 ) diff = ABS( rsp - rs2 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 2nd sphere boundary ihwin = 3 hregn = dz(3) / 10. if( diff0 .lt. eps ) i0 = 3 ! get right material go to 800 ! good step end if r2 = rs2**2 ! set radius for step adjustment logic ! else if( i0.eq.3 .and. i1.eq.3 ) then go to 800 ! good step ! ! It's also possible for trajectories to cross subregions backwards else if( i0.eq.2 .and. i1.eq.1 ) then diff0 = ABS( rsp0 - rs1 ) diff = ABS( rsp - rs1 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 1st sphere boundary ! Leave the small stepping appropriate to region 2 if( diff0 .lt. eps ) i0 = 1 ! but use right material go to 800 ! good step end if r2 = rs1**2 ! set radius for step adjustment logic ! else if( i0.eq.3 .and. i1.eq.2 ) then diff0 = ABS( rsp0 - rs2 ) diff = ABS( rsp - rs2 ) if( diff0.le.eps .or. diff.le.eps ) then ! One end of the step ends on the 2nd sphere boundary ihwin = 2 hregn = dz(2) / 10. if( diff0 .lt. eps ) i0 = 2 ! use right material go to 800 ! good step end if r2 = rs2**2 ! set radius for step adjustment logic ! else if( i0.eq.3 .and. i1.eq.1 ) then diff0 = ABS( rsp0 - rs2 ) if( diff0 .le. eps ) then ! step starts on 2nd sphere boundary ihwin = 2 hregn = dz(2) / 10. iflag = 1 go to 900 ! redo step end if r2 = rs2**2 ! set radius for step adjustment logic ! end if ! end of test on i0 and i1 end if ! end of test on gparm(1,1) ! --------------------------------------------------------- ! Step crosses significantly over the region boundary ! Readjust the step size to the remaining path length ! in the incident region do i=1,3 dp(i) = x(i) - x0(i) end do ! dp2 = dp(1)**2 + dp(2)**2 + dp(3)**2 ! if( dp2 .eq. 0. ) then iflag = -58 go to 900 end if ! p12 = x0(1)**2 + x0(2)**2 + x0(3)**2 aqf = dp2 bqf = 2. * DOT3(x0,dp) cqf = p12 - r2 u = -bqf + gparm(1,1)*SQRT( bqf**2 - 4.*aqf*cqf ) u = u / (2.*aqf) if( u.lt.-eps .or. u.gt.1.+eps ) then iflag = -58 go to 900 end if hregn = u * (xp(3)-xp0(3)) ! fractional step size iflag = 1 ! signal to redo this step go to 900 ! ......................................................... 800 continue ! good step mt = mtag(i0) if( mt .ne. 'VAC' ) then call LOAD_MATERIAL in_hemiwin = .true. end if ! ......................................................... 900 continue return end c ********************************************************* logical function IN_NIA(hstp) ! ! Determine if particle step goes thru a Non-Isoceles Absorber ! ! Ref. J.S. Berg, MUC-NOTE-COOL_THEORY-261, October 2002 ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 p0(3),p1(3),u0(3),u1(3),x1(3),t(3),v(3) logical in0,in1 data rad/0.01745329/ ! pi / 180 ! dum = hstp ! in_nia = .false. if( mgt .ne. 'NIA' ) go to 900 ! in0 = .false. in1 = .false. ! zv = sr0 + gparm(1,krad) dz0 = gparm(2,krad) dz1 = gparm(3,krad) th0 = gparm(4,krad) * rad ph0 = gparm(5,krad) * rad th1 = gparm(6,krad) * rad ph1 = gparm(7,krad) * rad ! do i=1,2 p0(i) = 0.d0 p1(i) = 0.d0 end do p0(3) = zv - dz0 p1(3) = zv + dz1 ! sth0 = SIN(th0) u0(1) = -sth0 * COS(ph0) u0(2) = -sth0 * SIN(ph0) u0(3) = -COS(th0) ! sth1 = SIN(th1) u1(1) = -sth1 * COS(ph1) u1(2) = -sth1 * SIN(ph1) u1(3) = COS(th1) ! call VMAG(pp,pm) t(1) = pp(1) / pm t(2) = pp(2) / pm t(3) = pp(3) / pm ! ! Check end point do i=1,3 v(i) = xp(i) - p0(i) end do val0 = DOT3( v,u0 ) do i=1,3 v(i) = xp(i) - p1(i) end do val1 = DOT3( v,u1 ) if( val0.lt.0. .and. val1.lt.0. ) in1 = .true. ! ! Check starting point x1(1) = xp(1) - pp(1)/pp(3)*hdid x1(2) = xp(2) - pp(2)/pp(3)*hdid x1(3) = xp(3) - hdid do i=1,3 v(i) = x1(i) - p0(i) end do val0 = DOT3( v,u0 ) do i=1,3 v(i) = x1(i) - p1(i) end do val1 = DOT3( v,u1 ) if( val0.lt.0. .and. val1.lt.0. ) in0 = .true. ! if( in0 .and. in1 ) then in_nia = .true. ! ! else if( in0 .and. (.not. in1) ) then ! in_nia = .true. ! hstp = DOT3( (p1-x1),t ) ! ! else if( in1 .and. (.not. in0) ) then ! in_nia = .true. ! hstp = DOT3( (xp-p0),t ) end if ! 900 continue return end c ********************************************************* logical function IN_PWEDGE(hstp) ! ! Determine if particle step goes thru an asymmetric, polynomial ! wedge oriented at an angle phi from the x-axis ! ! Describe polynomial shape in upper half-plane of x-z coord system ! dz(x) = a0+a1*x+a2*x**2+a3*x**3 ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 x(2),p(2) logical in0,in1 data rad/0.01745329/ ! pi / 180 ! in_pwedge = .false. if( mgt .ne. 'PWEDGE' ) go to 900 ! in0 = .false. in1 = .false. dz0 = 0. dz1 = 0. ! xv = gparm(2,krad) zv = sr0 + gparm(3,krad) phi = gparm(4,krad) * rad w = gparm(5,krad) h = gparm(6,krad) a0 = gparm(7,krad) a1 = gparm(8,krad) a2 = gparm(9,krad) a3 = gparm(10,krad) ! ! Get limits of problem box xmin = xv - w xmax = xv ymin = -h / 2. ymax = h / 2. dz = a0 + a1*xmin + a2*xmin**2 + a3*xmin**3 zmin = zv - dz zmax = zv + dz ! ! Rotate particle coordinates into frame where wedge vertex lies ! along the +ve x-axis ca = COS(phi) sa = SIN(phi) x(1) = xp(1)*ca + xp(2)*sa x(2) = -xp(1)*sa + xp(2)*ca p(1) = pp(1)*ca + pp(2)*sa p(2) = -pp(1)*sa + pp(2)*ca ! ! Check end point if( x(1).gt.xmin .and. x(1).lt.xmax .and. & x(2).gt.ymin .and. x(2).lt.ymax .and. & xp(3).gt.zmin .and. xp(3).lt.zmax ) then ! By construction zsurf is a positive quantity in a c.s. ! centered at zv zsurf = a0 + a1*x(1) + a2*x(1)**2 + a3*x(1)**3 dz = xp(3) - zv if( dz .ge. 0. ) then dz1 = xp(3) - (zv + zsurf) if( dz1 .le. 0. ) in1 = .true. else ! xp(3) < zv dz1 = xp(3) - (zv - zsurf) if( dz1 .ge. 0. ) in1 = .true. end if end if ! ! Check starting point x1 = x(1) - p(1)/pp(3)*hdid y1 = x(2) - p(2)/pp(3)*hdid z1 = xp(3) - hdid ! if( x1.gt.xmin .and. x1.lt.xmax .and. y1.gt.ymin .and. & y1.lt.ymax .and. z1.gt.zmin .and. z1.lt.zmax ) then zsurf = a0 + a1*x1 + a2*x1**2 + a3*x1**3 dz = z1 - zv if( dz .ge. 0. ) then dz0 = z1 - (zv + zsurf) if( dz0 .le. 0. ) in0 = .true. else ! z1 < zv dz0 = z1 - (zv - zsurf) if( dz1 .ge. 0. ) in0 = .true. end if end if ! if( in0 .and. in1 ) then in_pwedge = .true. ! else if( in0 .and. (.not. in1) ) then in_pwedge = .true. hstp = hstp * ABS(dz0) / ( ABS(dz0)+ABS(dz1) ) ! else if( in1 .and. (.not. in0) ) then in_pwedge = .true. hstp = hstp * ABS(dz1) / ( ABS(dz0)+ABS(dz1) ) end if ! 900 continue return end c ********************************************************* logical function IN_RING(hstp) ! ! Determine how much of particle step goes thru an annular ring of material ! implicit double precision (a-h,o-z) include 'icommon.inc' ! logical in1,in2 ! in_ring = .false. stp = hstp if( mgt .ne. 'RING' ) go to 900 ! in1 = .false. in2 = .false. ! ! Set geometrical parameters a1 = gparm(1,krad) ! inner radius of ring a2 = gparm(2,krad) ! outer radius of ring ! ! Check location of end of step r2 = SQRT( xp(1)**2 + xp(2)**2 ) if( r2.ge.a1 .and. r2.le.a2 ) in2 = .true. ! ! Check location of start of step xpr = pp(1) / pp(3) ypr = pp(2) / pp(3) x1 = xp(1) - xpr*hdid y1 = xp(2) - ypr*hdid r1 = SQRT( x1*x1 + y1*y1 ) if( r1.ge.a1 .and. r1.le.a2 ) in1 = .true. ! ! --------------------------------------------------------- ! Step handling logic (6 cases of step logic) if( in1 .and. in2 ) then ! all of hstp is in material in_ring = .true. stp = hstp go to 900 end if ! if( .not.in1 .and. .not.in2 ) then ! none of hstp in the material stp = 0. go to 900 end if ! ! Step crossed a radial boundary z1 = xp(3) - hdid x2 = xp(1) y2 = xp(2) z2 = xp(3) denom = xpr*xpr + ypr*ypr f1 = x1*xpr + y1*ypr ! if( in1 .and. r2.lt.a1 ) then ! case 3 arg = f1*f1 - denom*(x1*x1+y1*y1-a1*a1) if( f1 .ge. 0. ) then zb = (-f1 + SQRT(arg)) / denom else zb = (-f1 - SQRT(arg)) / denom end if xb = xpr*zb yb = ypr*zb stp = SQRT( xb*xb + yb*yb + zb*zb ) ! else if( in1 .and. r2.gt.a2 ) then ! case 4 arg = f1*f1 - denom*(x1*x1+y1*y1-a2*a2) if( f1 .ge. 0. ) then zb = (-f1 + SQRT(arg)) / denom else zb = (-f1 - SQRT(arg)) / denom end if xb = xpr*zb yb = ypr*zb stp = SQRT( xb*xb + yb*yb + zb*zb ) ! else if( in2 .and. r1.lt.a1 ) then ! case 5 arg = f1*f1 - denom*(x1*x1+y1*y1-a1*a1) if( f1 .ge. 0. ) then zb = (-f1 + SQRT(arg)) / denom else zb = (-f1 - SQRT(arg)) / denom end if xb = xpr*zb yb = ypr*zb stp = hstp - SQRT( xb*xb + yb*yb + zb*zb ) ! else if( in2 .and. r1.gt.a2 ) then ! case 6 arg = f1*f1 - denom*(x1*x1+y1*y1-a2*a2) if( f1 .ge. 0. ) then zb = (-f1 + SQRT(arg)) / denom else zb = (-f1 - SQRT(arg)) / denom end if xb = xpr*zb yb = ypr*zb stp = hstp - SQRT( xb*xb + yb*yb + zb*zb ) end if in_ring = .true. ! --------------------------------------------------------- 900 continue ! if( stp.lt.0. .or. stp.gt.hstp ) then ! write(*,'(a,2i5,2e12.4)') 'RING ',j7rec,jpar,xp(3),stp ! stop ! end if hstp = stp end c ********************************************************* logical function IN_WEDGE(hstp) ! ! Determine if particle step goes thru an asymmetric wedge ! oriented at an angle phi from the x-axis ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 x(2),p(2) logical in1,in2 data rad/0.01745329d0/ ! pi / 180 ! in_wedge = .false. if( mgt .ne. 'WEDGE' ) go to 900 ! in1 = .false. in2 = .false. ! alp = gparm(1,krad) * rad xv = gparm(2,krad) zv = sr0 + gparm(3,krad) phi = gparm(4,krad) * rad w = gparm(5,krad) h = gparm(6,krad) ta = TAN(alp/2.d0) dz = w * ta ! xmin = xv - w xmax = xv ymin = -h / 2.d0 ymax = h / 2.d0 zmin = zv - dz zmax = zv + dz ! ! Rotate particle coordinates into frame where wedge vertex lies ! along the +ve x-axis ca = COS(phi) sa = SIN(phi) x(1) = xp(1)*ca + xp(2)*sa x(2) = -xp(1)*sa + xp(2)*ca p(1) = pp(1)*ca + pp(2)*sa p(2) = -pp(1)*sa + pp(2)*ca ! ! Check end point if( x(1).gt.xmin .and. x(1).lt.xmax .and. & x(2).gt.ymin .and. x(2).lt.ymax .and. & xp(3).gt.zmin .and. xp(3).lt.zmax ) then dx = ABS(x(1)-xv) dzs = dx * ta ! dz for surface dzp = ABS(xp(3)-zv) ! dz for particle if( dzp .le. dzs ) in2 = .true. end if ! ! Check starting point x1 = x(1) - p(1)/pp(3)*hdid y1 = x(2) - p(2)/pp(3)*hdid z1 = xp(3) - hdid ! if( x1.gt.xmin .and. x1.lt.xmax .and. y1.gt.ymin .and. & y1.lt.ymax .and. z1.gt.zmin .and. z1.lt.zmax ) then dx = ABS(x1-xv) dzs = dx * ta ! dz for surface dzp = ABS(z1-zv) ! dz for particle if( dzp .le. dzs ) in1 = .true. end if ! if( in1 .and. in2 ) then in_wedge = .true. ! else if( in1 .and. (.not. in2) ) then in_wedge = .true. dx1 = ABS(x1-xv) dx2 = ABS(xp(1)-xv) a = zv - z1 + dx1*ta b = xp(3) - dx2*ta hstp = hstp * ABS(a) / ( ABS(a)+ABS(b) ) ! else if( in2 .and. (.not. in1) ) then in_wedge = .true. dx1 = ABS(x1-xv) dx2 = ABS(xp(1)-xv) a = -dx1*ta - z1 b = xp(3) - zv + dx2*ta hstp = hstp * ABS(b) / ( ABS(a)+ABS(b) ) end if ! 900 continue end c ********************************************************* subroutine LANDAU(yran) c c Determine fluctuation in ionization energy loss c from Landau distribution c c Derived from GEANT routine GLANDG C. C. ****************************************************************** C. * * C. * Copy of the CERN library routine GENLAN * C. * Generation of LANDAU-distributed random numbers by 4-point * C. * interpolation in the previously-tabulated inverse cumulative * C. * distribution * C. * * C. ****************************************************************** ! implicit double precision (a-h,o-z) include 'icommon.inc' ! PARAMETER (RANGE1=0.807069,RANGE2=0.994869) PARAMETER (X1BOT =0.0 ,X2BOT =0.791240) PARAMETER (GP1INV=118.836 ,GP2INV=486.178 ) C DIMENSION XCUM1(100) , XCUM2(100) DIMENSION RNDM(1) save xcum1,xcum2 ! C====> 1ST TABLE OF INVERSE CUMULATIVE LANDAU POINTS. ( 0.0 2ND TABLE OF INVERSE CUMULATIVE LANDAU POINTS. ( 0.791 !) * C. * * C. * Author : L.Urban * C. * Date : 28.04.1988 Last update : 1.02.90 * C. * * C. ****************************************************************** C. ! implicit double precision (a-h,o-z) include 'icommon.inc' ! PARAMETER (MAXRND=100) DIMENSION RNDM(MAXRND),APOIS(3),NPOIS(3),FPOIS(3) PARAMETER (HMXINT=2**30) ! PARAMETER (HMXINT=2**45) PARAMETER (ONE=1,HALF=ONE/2,ZERO=0) PARAMETER (RCD=0.40, RCD1=1.-RCD, PROBLM=0.01) PARAMETER( C1=4, C2=16) * * ------------------------------------------------------------------ * amass = pmass emass = mass(1) ! POTI=1.6E-8*Z**0.9 POTIL=LOG(POTI) * IF(Z.GT.2.) THEN F2=2./Z ELSE F2=0. ENDIF F1=1.-F2 * E2 = 1.E-8*Z*Z E2L= LOG(E2) E1L= (POTIL-F2*E2L)/F1 E1 = EXP(E1L) * * P2=P*P B2=P2/(E*E) BG2=P2/(AMASS*AMASS) ! IF(ITRTYP.EQ.2) THEN if( ABS(iptyp) .eq. 1 ) then TM=P2/(E+EMASS) ! IF(CHARGE.LT.0.) TM=TM/2. IF(pCHARGE.LT.0.) TM=TM/2. TM=TM-POTI IF(TM.GT.DCUTE) TM=DCUTE ELSE TM=EMASS*P2/(0.5*AMASS*AMASS+EMASS*E) TM=TM-POTI IF(TM.GT.DCUTM) TM=DCUTM ENDIF * * * *** Protection against negative TM --------------------- * TM can be negative only for heavy particles with a very * low kinetic energy (e.g. for proton with T 100-300 kev) TM=MAX(TM,ZERO) * W = TM+POTI WW = W/POTI WWW= 2.*EMASS*BG2 WL = LOG(WWW) CSB=STEP*RCD1*DEDX/(WL-POTIL-B2) APOIS(1)=CSB*F1*(WL-E1L-B2)/E1 APOIS(2)=CSB*F2*(WL-E2L-B2)/E2 * IF(TM.GT.0.) THEN APOIS(3)=RCD*DEDX*STEP*TM/(POTI*W*LOG(WW)) ELSE APOIS(1)=APOIS(1)/RCD1 APOIS(2)=APOIS(2)/RCD1 APOIS(3)=0. ENDIF * * calculate the probability of the zero energy loss * APSUM=APOIS(1)+APOIS(2)+APOIS(3) if( apsum .lt. 0d0 ) then iflag = -29 go to 900 end if IF(APSUM.LT.50. ) THEN PROB=EXP(-APSUM) ELSE PROB=0. ENDIF * * * do it differently if prob > problm <==================== IF(PROB.GT.PROBLM) THEN E0=1.E-8 EMEAN=DEDX*STEP IF(TM.LE.0.) THEN * excitation only .... APOIS(1)=EMEAN/E0 * CALL GPOISS(APOIS,NPOIS,1) FPOIS(1)=NPOIS(1) DE=FPOIS(1)*E0 * ELSE * ionization only .... EM=TM+E0 APOIS(1)=EMEAN*(EM-E0)/(EM*E0*LOG(EM/E0)) CALL GPOISS(APOIS,NPOIS,1) NN=NPOIS(1) DE=0. * IF(NN.GT.0) THEN RCORR=1. IF(NN.GT.MAXRND) THEN RCORR=FLOAT(NN)/MAXRND NN=MAXRND * ENDIF W=(EM-E0)/EM CALL GRNDM(RNDM,NN) DO 10 I=1,NN DE=DE+E0/(1.-W*RNDM(I)) 10 CONTINUE DE=RCORR*DE * ENDIF ENDIF GOTO 999 ENDIF * IF(MAX(APOIS(1),APOIS(2),APOIS(3)).LT.HMXINT) THEN CALL GPOISS(APOIS,NPOIS,3) FPOIS(1)=NPOIS(1) FPOIS(2)=NPOIS(2) FPOIS(3)=NPOIS(3) ELSE DO 20 JPOIS=1, 3 IF(APOIS(JPOIS).LT.HMXINT) THEN CALL GPOISS(APOIS(JPOIS),NPOIS(JPOIS),1) FPOIS(JPOIS)=NPOIS(JPOIS) ELSE CALL GRNDM(RNDM,2) FPOIS(JPOIS)=ABS(SQRT(-2.*LOG(RNDM(1)*ONE) + *APOIS(JPOIS))*SIN(TWOPI*RNDM(2)*ONE)+APOIS(JPOIS)) ENDIF 20 CONTINUE ENDIF * * Now we have all three numbers in REAL/DOUBLE * variables. REALK is actually an INTEGER that now may * exceed the machine representation limit for integers. * DE=FPOIS(1)*E1+FPOIS(2)*E2 * * smear to avoid peaks in the energy loss (note: E1< th0, use Tollestrup eqs. 19-23 zsmcs = 0. do i=1,nelcomp(krad) zsmcs = zsmcs + felcomp(i,krad) * zelcomp(i,krad)**2 end do pc2 = 0.1571d12*dens/a*(hstp*100.)*zsmcs/(bet*bet) chicsq = pc2 / (pavg*1d9)**2 if( n .le. 4 ) then bel = tab3(n) + LOG(th0*(pavg*1d9)/1d5) cnst = LOG(pc2) - 2.*tab1el(n) - 0.154 + 2.*bel/z else ! use standard Moliere method with Z*Z ! Recompute compound sums with Z*Z zemcs = 0. zxmcs = 0. do i=1,nelcomp(krad) zemcs = zemcs + felcomp(i,krad) * & zelcomp(i,krad)**2 * & LOG( zelcomp(i,krad)**(-0.6667) ) zxmcs = zxmcs + felcomp(i,krad) * & zelcomp(i,krad)**2 * & LOG( 1. + 3.34 * (zelcomp(i,krad)/137.036/bet)**2 ) end do chia2 = 2.016e-5 * facmms / & (pavg*1e3)**2 * EXP( (zxmcs-zemcs)/zsmcs ) cnst = LOG(chicsq/(1.167*chia2)) end if ! end of test on n (=Z) if( mt .eq. 'LIH' ) then ! use eq. 21 pc2H = 0.1571d12*0.0708/1.008*(hstp*100.)*1./(bet*bet) ! [eV^2] pc2Li = 0.1571d12*0.534/6.941*(hstp*100.)*9./(bet*bet) ! [eV^2] const = LOG(pc2H)-2.*tab1el(1) + LOG(pc2Li)-2.*tab1el(3) & - 0.154 + 2.*tab3(1) + 2.*tab3(3)/3. & + LOG(th0*(pavg*1d9)/1d5) end if end if ! end of test on thmsw ! --------------------------------------------------------- chic = SQRT(chicsq) ! ! Only CHIC and CONST affect the rest of this routine B=5. c DO 10 L=1,10 IF(ABS(B).LT.1.E-10)THEN B=1.E-10 ENDIF DB=-(B-LOG(ABS(B))-CNST)/(1.-1./B) B=B+DB IF(ABS(DB).LE.0.0001)GO TO 20 10 CONTINUE c GO TO 90 c 20 CONTINUE IF(B.LE.0.)GO TO 90 BINV = 1./B TINT(1) = 0. c DO 30 JA=2,4 TINT(JA)=F0I(JA)+(F1I(JA)+F2I(JA)*BINV)*BINV 30 CONTINUE c NMAX = 4 40 CONTINUE XINT = RAN1(irnarg) c DO 50 NA=3,40 IF(NA.GT.NMAX) THEN TINT(NA)=F0I(NA)+(F1I(NA)+F2I(NA)*BINV)*BINV NMAX=NA ENDIF IF(XINT.LE.TINT(NA-1)) GO TO 60 50 CONTINUE c IF(XINT.LE.TINT(40)) THEN NA=40 GOTO 60 ELSE TMP=1.-(1.-B*(1.-XINT))**5 IF(TMP.LE.0.)GO TO 40 THRI=5./TMP GO TO 80 ENDIF c 60 CONTINUE NA = MAX(NA-1,3) NA3 = NA-3 c DO 70 M=1,4 NA3M=NA3+M ARG(M)=TINT(NA3M) VAL(M)=THRED(NA3M)**2 70 CONTINUE c call POLINT(arg,val,4,xint,thri,err) ! Num Rec c 80 CONTINUE TH=CHIC*SQRT(ABS(B*THRI)) IF(TH.GT.PI)GO TO 40 SINTH=SIN(TH) TEST=TH*(RAN1(irnarg))**2 IF(TEST.GT.SINTH)GO TO 40 ! ! Check if this angle is in the tail for thmsw < th0 case if( thmsw .le. th0 ) then if( th .gt. th0 ) then test = z / (z+1.) if( RAN1(irnarg) .gt. test ) then go to 40 ! deplete tail, find some other angle for this step end if end if end if ! GOTO 900 ! 90 CONTINUE iflag = -26 c 900 continue END ! ********************************************************* SUBROUTINE T2ION(XMIN,XMAX,P,PM,SPIN,RES) ! ! Auxiliary routine for SAMCS ! IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) PARAMETER (PME=0.51099906D-03) E=SQRT(P**2+PM**2) ARA=2.*PME*E/(PME**2+PM**2) C2=ARA/(1.D0+ARA) BET2=(P/E)**2 a1aa=-1.+bet2*(1.+ara)**2/ara**2 ! -bet2**2*(2.+ara)/ara*spin a2new=(-1.+bet2/c2-bet2**2*spin)*(1.-c2) a2t=(a2new/(1.-c2*xmin)/(1.-c2*xmax)-bet2**2*spin) ! *(xmax-xmin) res=a1aa*log((1.-c2*xmax)/(1.-c2*xmin))+a2t return end ! ********************************************************* SUBROUTINE T4ION(YMIN,YMAX,P,PM,SPIN,RES) ! ! Auxiliary routine for SAMCS ! IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) PARAMETER (PME=0.51099906D-03) E=SQRT(P**2+PM**2) R=2.*PME*E/(PME**2+PM**2) C2=R/(1.D0+R) BET2=(P/E)**2 C3=BET2**2*C2**2*SPIN d1=-2.-bet2 d2=1.+2.*bet2+c3 d3=-bet2-2.*c3 zmax=1.-c2*Ymax zmin=1.-c2*Ymin f1=1.+d1/c2+ d2/c2**2+ d3/c2**3+ c3/c2**4 f2= -d1/c2-2.*d2/c2**2-3.*d3/c2**3-4.*c3/c2**4 f3= d2/c2**2+3.*d3/c2**3+6.*c3/c2**4 f4= -d3/c2**3-4.*c3/c2**4 f5= c3/c2**4 s1=-f1*(-1./zmax**3+1./zmin**3)/3./c2 s2=-f2*(-1./zmax**2+1./zmin**2)/2./c2 s3=-f3*(-1./zmax+1./zmin)/c2 s4=-f4*log(zmax/zmin)/c2 s5=-f5*(zmax-zmin)/c2 RES=s1+s2+s3+s4+s5 RETURN END c ********************************************************* function VAVILOV(rkappa,beta2,ran) c c Determine ionization energy loss fluctuation c from Vavilov distribution c c Derived from GEANT routine GVAVIV C. C. ****************************************************************** C. * * C. * Initializes the parameters for a Vavilov distribution * C. * * C. * BETA2 is the particle velocity squared * C. * RKAPPA the usual straggling parameter * C. * * C. * This routine has been extracted from a proposed CERNLIB * C. * set of routines (VAVCOE,VAVFSM). * C. * * C. * Alberto Rotondi, Paolo Montagna * C. * Fast calculation of Vavilov distribution NIM B47:215, 1990 * C. * * C. * Author K.S.Koelbig ******** * C. * * C. ****************************************************************** C. implicit double precision (a-h,o-z) DIMENSION AC(0:12),HC(0:8),H(9) PARAMETER 1(BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05, 2 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1 , 3 BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1 ) PARAMETER 1(FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2), 2 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1), 3 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3)) DIMENSION EDGEC(2:7),FNINV(5),DRK(5),DSIGM(5),ALFA(2:5) DIMENSION U1(13),U2(13),U3(13),U4(12),U5(13),U6(13),U7( 8),U8(13) DIMENSION V1(12),V2(12),V3(13),V4(12),V5(13),V6(13),V7(11),V8(11) DIMENSION W1(13),W2(11),W3(13),W4(13),W5(13),W6(13), W8( 8) DATA FNINV /1, 0.5, 0.33333333, 0.25, 0.2/ DATA (EDGEC(J),J=2,7) 1/ 0.16666667E+0, 0.41666667E-1, 0.83333333E-2, 2 0.13888889E-1, 0.69444444E-2, 0.77160493E-3/ DATA (U1(K),K=1,13) 1/ 0.25850868E+0, 0.32477982E-1, -0.59020496E-2, 2 0. , 0.24880692E-1, 0.47404356E-2, 3 -0.74445130E-3, 0.73225731E-2, 0. , 4 0.11668284E-2, 0. , -0.15727318E-2,-0.11210142E-2/ DATA (U2(K),K=1,13) 1/ 0.43142611E+0, 0.40797543E-1, -0.91490215E-2, 2 0. , 0.42127077E-1, 0.73167928E-2, 3 -0.14026047E-2, 0.16195241E-1, 0.24714789E-2, 4 0.20751278E-2, 0. , -0.25141668E-2,-0.14064022E-2/ DATA (U3(K),K=1,13) 1/ 0.25225955E+0, 0.64820468E-1, -0.23615759E-1, 2 0. , 0.23834176E-1, 0.21624675E-2, 3 -0.26865597E-2, -0.54891384E-2, 0.39800522E-2, 4 0.48447456E-2, -0.89439554E-2, -0.62756944E-2,-0.24655436E-2/ DATA (U4(K),K=1,12) 1/ 0.12593231E+1, -0.20374501E+0, 0.95055662E-1, 2 -0.20771531E-1, -0.46865180E-1, -0.77222986E-2, 3 0.32241039E-2, 0.89882920E-2, -0.67167236E-2, 4 -0.13049241E-1, 0.18786468E-1, 0.14484097E-1/ DATA (U5(K),K=1,13) 1/-0.24864376E-1, -0.10368495E-2, 0.14330117E-2, 2 0.20052730E-3, 0.18751903E-2, 0.12668869E-2, 3 0.48736023E-3, 0.34850854E-2, 0. , 4 -0.36597173E-3, 0.19372124E-2, 0.70761825E-3, 0.46898375E-3/ DATA (U6(K),K=1,13) 1/ 0.35855696E-1, -0.27542114E-1, 0.12631023E-1, 2 -0.30188807E-2, -0.84479939E-3, 0. , 3 0.45675843E-3, -0.69836141E-2, 0.39876546E-2, 4 -0.36055679E-2, 0. , 0.15298434E-2, 0.19247256E-2/ DATA (U7(K),K=1,8) 1/ 0.10234691E+2, -0.35619655E+1, 0.69387764E+0, 2 -0.14047599E+0, -0.19952390E+1, -0.45679694E+0, 3 0. , 0.50505298E+0/ DATA (U8(K),K=1,13) 1/ 0.21487518E+2, -0.11825253E+2, 0.43133087E+1, 2 -0.14500543E+1, -0.34343169E+1, -0.11063164E+1, 3 -0.21000819E+0, 0.17891643E+1, -0.89601916E+0, 4 0.39120793E+0, 0.73410606E+0, 0. ,-0.32454506E+0/ DATA (V1(K),K=1,12) 1/ 0.27827257E+0, -0.14227603E-2, 0.24848327E-2, 2 0. , 0.45091424E-1, 0.80559636E-2, 3 -0.38974523E-2, 0. , -0.30634124E-2, 4 0.75633702E-3, 0.54730726E-2, 0.19792507E-2/ DATA (V2(K),K=1,12) 1/ 0.41421789E+0, -0.30061649E-1, 0.52249697E-2, 2 0. , 0.12693873E+0, 0.22999801E-1, 3 -0.86792801E-2, 0.31875584E-1, -0.61757928E-2, 4 0. , 0.19716857E-1, 0.32596742E-2/ DATA (V3(K),K=1,13) 1/ 0.20191056E+0, -0.46831422E-1, 0.96777473E-2, 2 -0.17995317E-2, 0.53921588E-1, 0.35068740E-2, 3 -0.12621494E-1, -0.54996531E-2, -0.90029985E-2, 4 0.34958743E-2, 0.18513506E-1, 0.68332334E-2,-0.12940502E-2/ DATA (V4(K),K=1,12) 1/ 0.13206081E+1, 0.10036618E+0, -0.22015201E-1, 2 0.61667091E-2, -0.14986093E+0, -0.12720568E-1, 3 0.24972042E-1, -0.97751962E-2, 0.26087455E-1, 4 -0.11399062E-1, -0.48282515E-1, -0.98552378E-2/ DATA (V5(K),K=1,13) 1/ 0.16435243E-1, 0.36051400E-1, 0.23036520E-2, 2 -0.61666343E-3, -0.10775802E-1, 0.51476061E-2, 3 0.56856517E-2, -0.13438433E-1, 0. , 4 0. , -0.25421507E-2, 0.20169108E-2,-0.15144931E-2/ DATA (V6(K),K=1,13) 1/ 0.33432405E-1, 0.60583916E-2, -0.23381379E-2, 2 0.83846081E-3, -0.13346861E-1, -0.17402116E-2, 3 0.21052496E-2, 0.15528195E-2, 0.21900670E-2, 4 -0.13202847E-2, -0.45124157E-2, -0.15629454E-2, 0.22499176E-3/ DATA (V7(K),K=1,11) 1/ 0.54529572E+1, -0.90906096E+0, 0.86122438E-1, 2 0. , -0.12218009E+1, -0.32324120E+0, 3 -0.27373591E-1, 0.12173464E+0, 0. , 4 0. , 0.40917471E-1/ DATA (V8(K),K=1,11) 1/ 0.93841352E+1, -0.16276904E+1, 0.16571423E+0, 2 0. , -0.18160479E+1, -0.50919193E+0, 3 -0.51384654E-1, 0.21413992E+0, 0. , 4 0. , 0.66596366E-1/ DATA (W1(K),K=1,13) 1/ 0.29712951E+0, 0.97572934E-2, 0. , 2 -0.15291686E-2, 0.35707399E-1, 0.96221631E-2, 3 -0.18402821E-2, -0.49821585E-2, 0.18831112E-2, 4 0.43541673E-2, 0.20301312E-2, -0.18723311E-2,-0.73403108E-3/ DATA (W2(K),K=1,11) 1/ 0.40882635E+0, 0.14474912E-1, 0.25023704E-2, 2 -0.37707379E-2, 0.18719727E+0, 0.56954987E-1, 3 0. , 0.23020158E-1, 0.50574313E-2, 4 0.94550140E-2, 0.19300232E-1/ DATA (W3(K),K=1,13) 1/ 0.16861629E+0, 0. , 0.36317285E-2, 2 -0.43657818E-2, 0.30144338E-1, 0.13891826E-1, 3 -0.58030495E-2, -0.38717547E-2, 0.85359607E-2, 4 0.14507659E-1, 0.82387775E-2, -0.10116105E-1,-0.55135670E-2/ DATA (W4(K),K=1,13) 1/ 0.13493891E+1, -0.26863185E-2, -0.35216040E-2, 2 0.24434909E-1, -0.83447911E-1, -0.48061360E-1, 3 0.76473951E-2, 0.24494430E-1, -0.16209200E-1, 4 -0.37768479E-1, -0.47890063E-1, 0.17778596E-1, 0.13179324E-1/ DATA (W5(K),K=1,13) 1/ 0.10264945E+0, 0.32738857E-1, 0. , 2 0.43608779E-2, -0.43097757E-1, -0.22647176E-2, 3 0.94531290E-2, -0.12442571E-1, -0.32283517E-2, 4 -0.75640352E-2, -0.88293329E-2, 0.52537299E-2, 0.13340546E-2/ DATA (W6(K),K=1,13) 1/ 0.29568177E-1, -0.16300060E-2, -0.21119745E-3, 2 0.23599053E-2, -0.48515387E-2, -0.40797531E-2, 3 0.40403265E-3, 0.18200105E-2, -0.14346306E-2, 4 -0.39165276E-2, -0.37432073E-2, 0.19950380E-2, 0.12222675E-2/ DATA (W8(K),K=1,8) 1/ 0.66184645E+1, -0.73866379E+0, 0.44693973E-1, 2 0. , -0.14540925E+1, -0.39529833E+0, 3 -0.44293243E-1, 0.88741049E-1/ vavilov=0 IF(RKAPPA .LT. 0.01 .OR. RKAPPA .GT. 12) RETURN IF(RKAPPA .GE. 0.29) THEN ITYPE=1 NPT=100 WK=1/SQRT(RKAPPA) AC(0)=(-0.032227*BETA2-0.074275)*RKAPPA+ 1 (0.24533*BETA2+0.070152)*WK+(-0.55610*BETA2-3.1579) AC(8)=(-0.013483*BETA2-0.048801)*RKAPPA+ 1 (-1.6921*BETA2+8.3656)*WK+(-0.73275*BETA2-3.5226) DRK(1)=WK**2 DSIGM(1)=SQRT(RKAPPA/(1-0.5*BETA2)) DO 1 J = 1,4 DRK(J+1)=DRK(1)*DRK(J) DSIGM(J+1)=DSIGM(1)*DSIGM(J) 1 ALFA(J+1)=(FNINV(J)-BETA2*FNINV(J+1))*DRK(J) HC(0)=LOG(RKAPPA)+BETA2+0.42278434 HC(1)=DSIGM(1) HC(2)=ALFA(3)*DSIGM(3) HC(3)=(3*ALFA(2)**2+ALFA(4))*DSIGM(4)-3 HC(4)=(10*ALFA(2)*ALFA(3)+ALFA(5))*DSIGM(5)-10*HC(2) HC(5)=HC(2)**2 HC(6)=HC(2)*HC(3) HC(7)=HC(2)*HC(5) DO 2 J = 2,7 2 HC(J)=EDGEC(J)*HC(J) HC(8)=0.39894228*HC(1) ELSEIF(RKAPPA .GE. 0.22) THEN ITYPE=2 NPT=150 X=1+(RKAPPA-BKMXX3)*FBKX3 Y=1+(SQRT(BETA2)-BKMXY3)*FBKY3 XX=2*X YY=2*Y X2=XX*X-1 X3=XX*X2-X Y2=YY*Y-1 Y3=YY*Y2-Y XY=X*Y P2=X2*Y P3=X3*Y Q2=Y2*X Q3=Y3*X PQ=X2*Y2 AC(1)=W1(1)+W1(2)*X+W1(4)*X3+W1(5)*Y+W1(6)*Y2+W1(7)*Y3+ 1 W1(8)*XY+W1(9)*P2+W1(10)*P3+W1(11)*Q2+W1(12)*Q3+W1(13)*PQ AC(2)=W2(1)+W2(2)*X+W2(3)*X2+W2(4)*X3+W2(5)*Y+W2(6)*Y2+ 1 W2(8)*XY+W2(9)*P2+W2(10)*P3+W2(11)*Q2 AC(3)=W3(1)+W3(3)*X2+W3(4)*X3+W3(5)*Y+W3(6)*Y2+W3(7)*Y3+ 1 W3(8)*XY+W3(9)*P2+W3(10)*P3+W3(11)*Q2+W3(12)*Q3+W3(13)*PQ AC(4)=W4(1)+W4(2)*X+W4(3)*X2+W4(4)*X3+W4(5)*Y+W4(6)*Y2+W4(7)*Y3+ 1 W4(8)*XY+W4(9)*P2+W4(10)*P3+W4(11)*Q2+W4(12)*Q3+W4(13)*PQ AC(5)=W5(1)+W5(2)*X+W5(4)*X3+W5(5)*Y+W5(6)*Y2+W5(7)*Y3+ 1 W5(8)*XY+W5(9)*P2+W5(10)*P3+W5(11)*Q2+W5(12)*Q3+W5(13)*PQ AC(6)=W6(1)+W6(2)*X+W6(3)*X2+W6(4)*X3+W6(5)*Y+W6(6)*Y2+W6(7)*Y3+ 1 W6(8)*XY+W6(9)*P2+W6(10)*P3+W6(11)*Q2+W6(12)*Q3+W6(13)*PQ AC(8)=W8(1)+W8(2)*X+W8(3)*X2+W8(5)*Y+W8(6)*Y2+W8(7)*Y3+W8(8)*XY AC(0)=-3.05 ELSEIF(RKAPPA .GE. 0.12) THEN ITYPE=3 NPT=200 X=1+(RKAPPA-BKMXX2)*FBKX2 Y=1+(SQRT(BETA2)-BKMXY2)*FBKY2 XX=2*X YY=2*Y X2=XX*X-1 X3=XX*X2-X Y2=YY*Y-1 Y3=YY*Y2-Y XY=X*Y P2=X2*Y P3=X3*Y Q2=Y2*X Q3=Y3*X PQ=X2*Y2 AC(1)=V1(1)+V1(2)*X+V1(3)*X2+V1(5)*Y+V1(6)*Y2+V1(7)*Y3+ 1 V1(9)*P2+V1(10)*P3+V1(11)*Q2+V1(12)*Q3 AC(2)=V2(1)+V2(2)*X+V2(3)*X2+V2(5)*Y+V2(6)*Y2+V2(7)*Y3+ 1 V2(8)*XY+V2(9)*P2+V2(11)*Q2+V2(12)*Q3 AC(3)=V3(1)+V3(2)*X+V3(3)*X2+V3(4)*X3+V3(5)*Y+V3(6)*Y2+V3(7)*Y3+ 1 V3(8)*XY+V3(9)*P2+V3(10)*P3+V3(11)*Q2+V3(12)*Q3+V3(13)*PQ AC(4)=V4(1)+V4(2)*X+V4(3)*X2+V4(4)*X3+V4(5)*Y+V4(6)*Y2+V4(7)*Y3+ 1 V4(8)*XY+V4(9)*P2+V4(10)*P3+V4(11)*Q2+V4(12)*Q3 AC(5)=V5(1)+V5(2)*X+V5(3)*X2+V5(4)*X3+V5(5)*Y+V5(6)*Y2+V5(7)*Y3+ 1 V5(8)*XY+V5(11)*Q2+V5(12)*Q3+V5(13)*PQ AC(6)=V6(1)+V6(2)*X+V6(3)*X2+V6(4)*X3+V6(5)*Y+V6(6)*Y2+V6(7)*Y3+ 1 V6(8)*XY+V6(9)*P2+V6(10)*P3+V6(11)*Q2+V6(12)*Q3+V6(13)*PQ AC(7)=V7(1)+V7(2)*X+V7(3)*X2+V7(5)*Y+V7(6)*Y2+V7(7)*Y3+ 1 V7(8)*XY+V7(11)*Q2 AC(8)=V8(1)+V8(2)*X+V8(3)*X2+V8(5)*Y+V8(6)*Y2+V8(7)*Y3+ 1 V8(8)*XY+V8(11)*Q2 AC(0)=-3.04 ELSE ITYPE=4 IF(RKAPPA .GE. 0.02) ITYPE=3 NPT=200 X=1+(RKAPPA-BKMXX1)*FBKX1 Y=1+(SQRT(BETA2)-BKMXY1)*FBKY1 XX=2*X YY=2*Y X2=XX*X-1 X3=XX*X2-X Y2=YY*Y-1 Y3=YY*Y2-Y XY=X*Y P2=X2*Y P3=X3*Y Q2=Y2*X Q3=Y3*X PQ=X2*Y2 IF(ITYPE .EQ. 4) GO TO 4 AC(1)=U1(1)+U1(2)*X+U1(3)*X2+U1(5)*Y+U1(6)*Y2+U1(7)*Y3+ 1 U1(8)*XY+U1(10)*P3+U1(12)*Q3+U1(13)*PQ AC(2)=U2(1)+U2(2)*X+U2(3)*X2+U2(5)*Y+U2(6)*Y2+U2(7)*Y3+ 1 U2(8)*XY+U2(9)*P2+U2(10)*P3+U2(12)*Q3+U2(13)*PQ AC(3)=U3(1)+U3(2)*X+U3(3)*X2+U3(5)*Y+U3(6)*Y2+U3(7)*Y3+ 1 U3(8)*XY+U3(9)*P2+U3(10)*P3+U3(11)*Q2+U3(12)*Q3+U3(13)*PQ AC(4)=U4(1)+U4(2)*X+U4(3)*X2+U4(4)*X3+U4(5)*Y+U4(6)*Y2+U4(7)*Y3+ 1 U4(8)*XY+U4(9)*P2+U4(10)*P3+U4(11)*Q2+U4(12)*Q3 AC(5)=U5(1)+U5(2)*X+U5(3)*X2+U5(4)*X3+U5(5)*Y+U5(6)*Y2+U5(7)*Y3+ 1 U5(8)*XY+U5(10)*P3+U5(11)*Q2+U5(12)*Q3+U5(13)*PQ AC(6)=U6(1)+U6(2)*X+U6(3)*X2+U6(4)*X3+U6(5)*Y+U6(7)*Y3+ 1 U6(8)*XY+U6(9)*P2+U6(10)*P3+U6(12)*Q3+U6(13)*PQ 4 AC(7)=U7(1)+U7(2)*X+U7(3)*X2+U7(4)*X3+U7(5)*Y+U7(6)*Y2+U7(8)*XY AC(8)=U8(1)+U8(2)*X+U8(3)*X2+U8(4)*X3+U8(5)*Y+U8(6)*Y2+U8(7)*Y3+ 1 U8(8)*XY+U8(9)*P2+U8(10)*P3+U8(11)*Q2+U8(13)*PQ AC(0)=-3.03 ENDIF AC(9)=(AC(8)-AC(0))/NPT IF(ITYPE .EQ. 3) THEN X=(AC(7)-AC(8))/(AC(7)*AC(8)) Y=1/LOG(AC(8)/AC(7)) P2=AC(7)**2 AC(11)=P2*(AC(1)*EXP(-AC(2)*(AC(7)+AC(5)*P2)- 1 AC(3)*EXP(-AC(4)*(AC(7)+AC(6)*P2)))-0.045*Y/AC(7))/ 2 (1+X*Y*AC(7)) AC(12)=(0.045+X*AC(11))*Y ENDIF IF(ITYPE .EQ. 4) AC(10)=0.995/vavilov_s(AC(8)) T = 2 * RAN / AC(9) RLAM=AC(0) FL=0 S=0 DO 21 N = 1,NPT RLAM=RLAM+AC(9) IF(ITYPE .EQ. 1) THEN FN=1 X=(RLAM+HC(0))*HC(1) H(1)=X H(2)=X**2-1 DO 31 K = 2,8 FN=FN+1 31 H(K+1)=X*H(K)-FN*H(K-1) Y=1+HC(7)*H(9) DO 32 K = 2,6 32 Y=Y+HC(K)*H(K+1) FU=HC(8)*EXP(-0.5*X**2)*MAX(Y,0.d0) ELSEIF(ITYPE .EQ. 2) THEN X=RLAM**2 FU=AC(1)*EXP(-AC(2)*(RLAM+AC(5)*X)- 1 AC(3)*EXP(-AC(4)*(RLAM+AC(6)*X))) ELSEIF(ITYPE .EQ. 3) THEN IF(RLAM .LT. AC(7)) THEN X=RLAM**2 FU=AC(1)*EXP(-AC(2)*(RLAM+AC(5)*X)- 1 AC(3)*EXP(-AC(4)*(RLAM+AC(6)*X))) ELSE X=1/RLAM FU=(AC(11)*X+AC(12))*X ENDIF ELSE FU=AC(10)*vavilov_e(RLAM) ENDIF S=S+FL+FU IF(S .GT. T) GO TO 22 21 FL=FU 22 S0=S-FL-FU vavilov=RLAM-AC(9) IF(S .GT. S0) vavilov=vavilov+AC(9)*(T-S0)/(S-S0) RETURN END c ********************************************************* function VAVILOV_E(x) c c Utility routine for calculating Vavilov distribution c c Derived from GEANT routine GLANDE c C. ****************************************************************** C. * * C. * Copy of the CERN library routine DENLAN (G110) * C. * * C. ****************************************************************** C. implicit double precision (a-h,o-z) DIMENSION P1(0:4),P2(0:4),P3(0:4),P4(0:4),P5(0:4),P6(0:4) DIMENSION Q1(0:4),Q2(0:4),Q3(0:4),Q4(0:4),Q5(0:4),Q6(0:4) DIMENSION A1(1:3),A2(1:2) C DATA (P1(I),I=0,4),(Q1(J),J=0,4) 1/ 0.42598 94875E+0,-0.12497 62550E+0, 0.39842 43700E-1, 2 -0.62982 87635E-2, 0.15111 62253E-2, 3 1.0 ,-0.33882 60629E+0, 0.95943 93323E-1, 4 -0.16080 42283E-1, 0.37789 42063E-2/ C DATA (P2(I),I=0,4),(Q2(J),J=0,4) 1/ 0.17885 41609E+0, 0.11739 57403E+0, 0.14888 50518E-1, 2 -0.13949 89411E-2, 0.12836 17211E-3, 3 1.0 , 0.74287 95082E+0, 0.31539 32961E+0, 4 0.66942 19548E-1, 0.87906 09714E-2/ C DATA (P3(I),I=0,4),(Q3(J),J=0,4) 1/ 0.17885 44503E+0, 0.93591 61662E-1, 0.63253 87654E-2, 2 0.66116 67319E-4,-0.20310 49101E-5, 3 1.0 , 0.60978 09921E+0, 0.25606 16665E+0, 4 0.47467 22384E-1, 0.69573 01675E-2/ C DATA (P4(I),I=0,4),(Q4(J),J=0,4) 1/ 0.98740 54407E+0, 0.11867 23273E+3, 0.84927 94360E+3, 2 -0.74377 92444E+3, 0.42702 62186E+3, 3 1.0 , 0.10686 15961E+3, 0.33764 96214E+3, 4 0.20167 12389E+4, 0.15970 63511E+4/ C DATA (P5(I),I=0,4),(Q5(J),J=0,4) 1/ 0.10036 75074E+1, 0.16757 02434E+3, 0.47897 11289E+4, 2 0.21217 86767E+5,-0.22324 94910E+5, 3 1.0 , 0.15694 24537E+3, 0.37453 10488E+4, 4 0.98346 98876E+4, 0.66924 28357E+5/ C DATA (P6(I),I=0,4),(Q6(J),J=0,4) 1/ 0.10008 27619E+1, 0.66491 43136E+3, 0.62972 92665E+5, 2 0.47555 46998E+6,-0.57436 09109E+7, 3 1.0 , 0.65141 01098E+3, 0.56974 73333E+5, 4 0.16591 74725E+6,-0.28157 59939E+7/ C DATA (A1(I),I=1,3) 1/ 0.41666 66667E-1,-0.19965 27778E-1, 0.27095 38966E-1/ C DATA (A2(I),I=1,2) 1/-0.18455 68670E+1,-0.42846 40743E+1/ C V=X IF(V .LT. -5.5) THEN U=EXP(V+1.0) vavilov_e=0.3989422803*(EXP(-1.0/U)/SQRT(U))* 1 (1.0+(A1(1)+(A1(2)+A1(3)*U)*U)*U) ELSE IF(V .LT. -1.0) THEN U=EXP(-V-1.0) vavilov_e=EXP(-U)*SQRT(U)* 1 (P1(0)+(P1(1)+(P1(2)+(P1(3)+P1(4)*V)*V)*V)*V)/ 2 (Q1(0)+(Q1(1)+(Q1(2)+(Q1(3)+Q1(4)*V)*V)*V)*V) ELSE IF(V .LT. 1.0) THEN vavilov_e=(P2(0)+(P2(1)+(P2(2)+(P2(3)+P2(4)*V)*V)*V)*V)/ 1 (Q2(0)+(Q2(1)+(Q2(2)+(Q2(3)+Q2(4)*V)*V)*V)*V) ELSE IF(V .LT. 5.0) THEN vavilov_e=(P3(0)+(P3(1)+(P3(2)+(P3(3)+P3(4)*V)*V)*V)*V)/ 1 (Q3(0)+(Q3(1)+(Q3(2)+(Q3(3)+Q3(4)*V)*V)*V)*V) ELSE IF(V .LT. 12.0) THEN U=1.0/V vavilov_e=U**2*(P4(0)+(P4(1)+(P4(2)+(P4(3)+P4(4)*U)*U)*U)*U)/ 1 (Q4(0)+(Q4(1)+(Q4(2)+(Q4(3)+Q4(4)*U)*U)*U)*U) ELSE IF(V .LT. 50.0) THEN U=1.0/V vavilov_e=U**2*(P5(0)+(P5(1)+(P5(2)+(P5(3)+P5(4)*U)*U)*U)*U)/ 1 (Q5(0)+(Q5(1)+(Q5(2)+(Q5(3)+Q5(4)*U)*U)*U)*U) ELSE IF(V .LT. 300.0) THEN U=1.0/V vavilov_e=U**2*(P6(0)+(P6(1)+(P6(2)+(P6(3)+P6(4)*U)*U)*U)*U)/ 1 (Q6(0)+(Q6(1)+(Q6(2)+(Q6(3)+Q6(4)*U)*U)*U)*U) ELSE U=1.0/(V-V*LOG(V)/(V+1.0)) vavilov_e=U**2*(1.0+(A2(1)+A2(2)*U)*U) END IF RETURN END c ********************************************************* function VAVILOV_S(x) c c Utility routine for calculating Vavilov distribution c c Derived from GEANT routine GLANDS c C. ****************************************************************** C. * * C. * Copy of the CERN library routine DSTLAN (G110) * C. * * C. ****************************************************************** C. implicit double precision (a-h,o-z) DIMENSION P1(0:4),P2(0:3),P3(0:3),P4(0:3),P5(0:3),P6(0:3) DIMENSION Q1(0:4),Q2(0:3),Q3(0:3),Q4(0:3),Q5(0:3),Q6(0:3) DIMENSION A1(1:3),A2(1:3) C DATA (P1(I),I=0,4),(Q1(J),J=0,4) 1/ 0.25140 91491E+0,-0.62505 80444E-1, 0.14583 81230E-1, 2 -0.21088 17737E-2, 0.74112 47290E-3, 3 1.0 ,-0.55711 75625E-2, 0.62253 10236E-1, 4 -0.31373 78427E-2, 0.19314 96439E-2/ C DATA (P2(I),I=0,3),(Q2(J),J=0,3) 1/ 0.28683 28584E+0, 0.35643 63231E+0, 0.15235 18695E+0, 2 0.22513 04883E-1, 3 1.0 , 0.61911 36137E+0, 0.17207 21448E+0, 4 0.22785 94771E-1/ C DATA (P3(I),I=0,3),(Q3(J),J=0,3) 1/ 0.28683 29066E+0, 0.30038 28436E+0, 0.99509 51941E-1, 2 0.87338 27185E-2, 3 1.0 , 0.42371 90502E+0, 0.10956 31512E+0, 4 0.86938 51567E-2/ C DATA (P4(I),I=0,3),(Q4(J),J=0,3) 1/ 0.10003 51630E+1, 0.45035 92498E+1, 0.10858 83880E+2, 2 0.75360 52269E+1, 3 1.0 , 0.55399 69678E+1, 0.19335 81111E+2, 4 0.27213 21508E+2/ C DATA (P5(I),I=0,3),(Q5(J),J=0,3) 1/ 0.10000 06517E+1, 0.49094 14111E+2, 0.85055 44753E+2, 2 0.15321 53455E+3, 3 1.0 , 0.50099 28881E+2, 0.13998 19104E+3, 4 0.42000 02909E+3/ C DATA (P6(I),I=0,3),(Q6(J),J=0,3) 1/ 0.10000 00983E+1, 0.13298 68456E+3, 0.91621 49244E+3, 2 -0.96050 54274E+3, 3 1.0 , 0.13398 87843E+3, 0.10559 90413E+4, 4 0.55322 24619E+3/ C DATA (A1(I),I=1,3) 1/-0.45833 33333E+0, 0.66753 47222E+0,-0.16417 41416E+1/ C DATA (A2(I),I=1,3) 1/ 1.0 ,-0.42278 43351E+0,-0.20434 03138E+1/ C V=X IF(V .LT. -5.5) THEN U=EXP(V+1.0) vavilov_s=0.3989422803*EXP(-1.0/U)*SQRT(U)* 1 (1.0+(A1(1)+(A1(2)+A1(3)*U)*U)*U) ELSE IF(V .LT. -1.0) THEN U=EXP(-V-1.0) vavilov_s=(EXP(-U)/SQRT(U))* 1 (P1(0)+(P1(1)+(P1(2)+(P1(3)+P1(4)*V)*V)*V)*V)/ 2 (Q1(0)+(Q1(1)+(Q1(2)+(Q1(3)+Q1(4)*V)*V)*V)*V) ELSE IF(V .LT. 1.0) THEN vavilov_s=(P2(0)+(P2(1)+(P2(2)+P2(3)*V)*V)*V)/ 1 (Q2(0)+(Q2(1)+(Q2(2)+Q2(3)*V)*V)*V) ELSE IF(V .LT. 4.0) THEN vavilov_s=(P3(0)+(P3(1)+(P3(2)+P3(3)*V)*V)*V)/ 1 (Q3(0)+(Q3(1)+(Q3(2)+Q3(3)*V)*V)*V) ELSE IF(V .LT. 12.0) THEN U=1.0/V vavilov_s=(P4(0)+(P4(1)+(P4(2)+P4(3)*U)*U)*U)/ 1 (Q4(0)+(Q4(1)+(Q4(2)+Q4(3)*U)*U)*U) ELSE IF(V .LT. 50.0) THEN U=1.0/V vavilov_s=(P5(0)+(P5(1)+(P5(2)+P5(3)*U)*U)*U)/ 1 (Q5(0)+(Q5(1)+(Q5(2)+Q5(3)*U)*U)*U) ELSE IF(V .LT. 300.0) THEN U=1.0/V vavilov_s=(P6(0)+(P6(1)+(P6(2)+P6(3)*U)*U)*U)/ 1 (Q6(0)+(Q6(1)+(Q6(2)+Q6(3)*U)*U)*U) ELSE U=1.0/(V-V*LOG(V)/(V+1.0)) vavilov_s=1.0-(A2(1)+(A2(2)+A2(3)*U)*U)*U END IF RETURN END ! ********************************************************* subroutine WANG ! ! Create a pion from WANG distribution ! implicit double precision (a-h,o-z) include 'icommon.inc' dimension pms(3) ! knt = 0 pmax = ABS( wangpmx ) ! 100 continue knt = knt + 1 if( knt .gt. 1000 ) then write(*,*) 'WANG sampling does not converge' write(2,*) 'WANG sampling does not converge' stop end if x = RANV(0d0,1d0) pt = RANV(0d0,pmax) arg = -wangb*x**wangc - wangd*pt cs = wanga * pmax * x * (1d0-x) * EXP(arg) v = RANV(0d0,wangfmx) if( v .gt. cs ) go to 100 ! ! Convert this proton track into a pion ! ipnum = ipnum + 1 if( wangpmx .gt. 0d0 ) then iptyp = 3 else iptyp = -3 end if phi = RANV(0d0,twopi) ! Pion in fixed coordinate system pms(1) = pt * COS(phi) pms(2) = pt * SIN(phi) pms(3) = x * pmax ! Angles of proton in LAB call VMAG(pp,pm) cth = pp(3) / pm pt = SQRT(pp(1)**2 + pp(2)**2) sth = pt / pm if( pt .gt. 0. ) then cph = pp(1) / pt sph = pp(2) / pt else cph = 1. sph = 0. end if ! Pion in the LAB frame pp(1) = pms(1)*cth*cph - pms(2)*sph + pms(3)*sth*cph pp(2) = pms(1)*cth*sph + pms(2)*cph + pms(3)*sth*sph pp(3) = -pms(1)*sth + pms(3)*cth ! end ! ********************************************************* DOUBLE PRECISION FUNCTION ZFORM(Q,R,B) ! ! Auxiliary routine for SAMCS ! IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) PARAMETER (PI=3.141592653589793227D+00) Z=PI*B/R Z2=Z*Z R2=3.*R**2*(1.+7.*Z2/3.)/5. R4=3.*R**4*(1.+6.*Z2+31.*Z2**2/3.)/7. Q2=Q**2 DD=0.02*R2/R4 IF(Q2.GT.DD) THEN EX=EXP(-PI*B*Q) EX2=EX*EX CT=(EX2+1.)/(1.-EX2) ST=1.-EX2 SK=-COS(Q*R) + Z*SIN(Q*R)*CT ZFORM=6.*PI*B*EX*SK/Q/ST/R/R/(1.+Z2) ELSE ZFORM=1.-Q2*R2/6.+Q2**2*R4/120. ENDIF RETURN END