From c314fb01115e61da379d5661e01ceaf353d537cf Mon Sep 17 00:00:00 2001 From: Kurt Date: Thu, 8 Oct 2015 14:57:55 -0400 Subject: [PATCH] upgrade AMPT to v2.25t5 --- .../AMPTInterface/src/amptsub.f | 847 ++++++++------- GeneratorInterface/AMPTInterface/src/art1f.f | 976 ++++++++++++++---- .../AMPTInterface/src/hijing1.383_ampt.f | 5 +- .../AMPTInterface/src/hipyset1.35.f | 276 ++--- GeneratorInterface/AMPTInterface/src/linana.f | 727 +++++++++---- GeneratorInterface/AMPTInterface/src/main.f | 125 +-- GeneratorInterface/AMPTInterface/src/zpc.f | 194 ++-- 7 files changed, 2023 insertions(+), 1127 deletions(-) diff --git a/GeneratorInterface/AMPTInterface/src/amptsub.f b/GeneratorInterface/AMPTInterface/src/amptsub.f index 125544097db1c..ddb0b8c5afc3c 100755 --- a/GeneratorInterface/AMPTInterface/src/amptsub.f +++ b/GeneratorInterface/AMPTInterface/src/amptsub.f @@ -12,12 +12,9 @@ c======================================================================= c.....subroutine to set up ART parameters and analysis files c.....before looping different events -cms -cms dlw & gsfs 8/2009 commented out lots of output files -cms SUBROUTINE ARTSET c - PARAMETER (AMU= 0.9383) + PARAMETER (AMU= 0.9383,nxymax=10001) double precision dpcoal,drcoal,ecritl INTEGER ZTA, ZPR common /gg/ dx,dy,dz,dpx,dpy,dpz @@ -41,6 +38,9 @@ SUBROUTINE ARTSET common /coal/dpcoal,drcoal,ecritl common/anim/nevent,isoft,isflag,izpc common /para7/ ioscar,nsmbbbar,nsmmeson + common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd, + 1 psembd,tmaxembd,phidecomp + common/xyembed/nxyjet,xyjet(nxymax,2) SAVE clin-10/03/03 ecritl: local energy density below which a parton c will freeze out (in GeV/fm^3), for improvements on string melting, @@ -161,24 +161,45 @@ SUBROUTINE ARTSET clin-5b/2008 c if(ioscar.eq.2) then if(ioscar.eq.2.or.ioscar.eq.3) then -cms OPEN (92,FILE='ana/parton-initial-afterPropagation.dat', -cms 1 STATUS = 'UNKNOWN') +c OPEN (92,FILE='ana/parton-initial-afterPropagation.dat', +c 1 STATUS = 'UNKNOWN') endif if(ioscar.eq.3) then clin-6/2009 write out full parton collision history: -cms OPEN (95,FILE='ana/parton-collisionsHistory.dat', -cms 1 STATUS='UNKNOWN') +c OPEN (95,FILE='ana/parton-collisionsHistory.dat', +c 1 STATUS='UNKNOWN') clin-6/2009 write out initial minijet information: -cms OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat', -cms 1 STATUS='UNKNOWN') +c OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat', +c 1 STATUS='UNKNOWN') clin-6/2009 write out parton info after coalescence: if(isoft.eq.4.or.isoft.eq.5) then -cms OPEN (85,FILE='ana/parton-after-coalescence.dat', -cms 1 STATUS='UNKNOWN') +c OPEN (85,FILE='ana/parton-after-coalescence.dat', +c 1 STATUS='UNKNOWN') endif endif clin-6/2009 write out initial transverse positions of initial nucleons: -cms OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN') +c OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN') +c +clin-8/2009 In case that random positions are used to embed high-Pt jets: + if(iembed.eq.3.or.iembed.eq.4) then +c OPEN (97,FILE='embed-jet-xy.txt',STATUS = 'UNKNOWN') +c read(97,*) nxyjet +c Save positions in array to reuse when embedding more jet pairs +c than the number of entries in the position file: + if(nevent.gt.nxyjet) then + if(nxyjet.gt.nxymax) then +c print *, 'Too many lines in embed-jet-xy.txt: +c 1 increase value of the parameter nxymax' + stop + elseif(nxyjet.le.0) then +c print *, 'Check number of entries in embed-jet-xy.txt' + stop + endif + do ixy=1,nxyjet +c read(97,*) xyjet(ixy,1),xyjet(ixy,2) + enddo + endif + endif RETURN END @@ -252,15 +273,13 @@ SUBROUTINE ARINI1 COMMON/RNDF77/NSEED cc SAVE /RNDF77/ common /para8/ idpert,npertd,idxsec - - logical iwrite - data iwrite / .false. / SAVE + clin-5/2008 for perturbatively-produced hadrons (currently only deuterons): -cms OPEN (91, FILE = 'ana/deuteron_processes.dat', -cms 1 STATUS = 'UNKNOWN') +c OPEN (91, FILE = 'ana/deuteron_processes.dat', +c 1 STATUS = 'UNKNOWN') if(idpert.eq.1.or.idpert.eq.2) then -cms OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN') +c OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN') endif c.....generate formation time and position at formation time. TAU0 = ARPAR1(1) @@ -270,20 +289,28 @@ SUBROUTINE ARINI1 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then if(NP.le.nattzp) return do 1001 I = nattzp+1, NP - IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN - if ( iwrite ) then - PRINT *, ' IN ARINI1' - PRINT *, 'ABS(PZ) .GE. EE for particle ', I - PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), - & ' PY = ', PYAR(I) - PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I) - PRINT *, ' XM = ', XMAR(I) - endif - RAP = 1000000.0 - GOTO 50 - END IF - RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I))) - 50 CONTINUE +clin-9/2012 determine rapidity more generally +c to prevent overflow when Pt~=0 and E=|Pz|: +c IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN +c PRINT *, ' IN ARINI1' +c PRINT *, 'ABS(PZ) .GE. EE for particle ', I +c PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), +c & ' PY = ', PYAR(I) +c PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I) +c PRINT *, ' XM = ', XMAR(I) +c RAP = 1000000.0 +c GOTO 50 +c END IF +cc RAP=0.5*LOG((PEAR(I)+PZAR(I))/(PEAR(I)-PZAR(I))) +c RAP=0.5*LOG((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5)) +c 50 CONTINUE + if((XMAR(I)**2+PXAR(I)**2+PYAR(I)**2).gt.0.) then + RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2)) + else + PRINT *, ' IN ARINI1 mt=0' + RAP = 1000000.0*sign(1.,PZAR(I)) + endif + VX = PXAR(I) / PEAR(I) VY = PYAR(I) / PEAR(I) FTAR(I) = TAU0 * COSH(RAP) @@ -292,32 +319,43 @@ SUBROUTINE ARINI1 GZAR(I) = TAU0 * SINH(RAP) clin-5/2009 No formation time for spectator projectile or target nucleons: if(PXAR(I).eq.0.and.PYAR(I).eq.0 - 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99 2 .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then +clin-2/2013 for spectator target nucleons in LAB frame: +c 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99 + if((PEAR(I)/HINT1(6).gt.0.99.and.PEAR(I)/HINT1(6).lt.1.01) + 1 .or.(PEAR(I)/HINT1(7).gt.0.99.and.PEAR(I)/HINT1(7).lt.1.01)) then +c TAUI=1.E-20 FTAR(I)=TAUI*COSH(RAP) GZAR(I)=TAUI*SINH(RAP) + endif endif 1001 continue clin-7/10/01-end clin-3/2009 cleanup of program flow: else DO 1002 I = 1, NP - IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN - if ( iwrite ) then - PRINT *, ' IN ARINI1' - PRINT *, 'ABS(PZ) .GE. EE for particle ', I - PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), - & ' PY = ', PYAR(I) - PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I) - PRINT *, ' XM = ', XMAR(I) - endif - RAP = 1000000.0 - GOTO 100 -c STOP - END IF - RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I))) - 100 CONTINUE +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN +c PRINT *, ' IN ARINI1' +c PRINT *, 'ABS(PZ) .GE. EE for particle ', I +c PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), +c & ' PY = ', PYAR(I) +c PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I) +c PRINT *, ' XM = ', XMAR(I) +c RAP = 1000000.0 +c GOTO 100 +cc STOP +c END IF +c 100 CONTINUE +c RAP=0.5*LOG((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5)) + if((XMAR(I)**2+PXAR(I)**2+PYAR(I)**2).gt.0.) then + RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2)) + else + PRINT *, ' IN ARINI1 mt=0' + RAP = 1000000.0*sign(1.,PZAR(I)) + endif + VX = PXAR(I) / PEAR(I) VY = PYAR(I) / PEAR(I) c.....give initial formation time shift @@ -334,13 +372,18 @@ SUBROUTINE ARINI1 cbz1/28/99end c 10/05/01 no formation time for spectator projectile or target nucleons: if(PXAR(I).eq.0.and.PYAR(I).eq.0 - 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99 2 .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then +clin-2/2013 for spectator target nucleons in LAB frame: +c 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99 + if((PEAR(I)/HINT1(6).gt.0.99.and.PEAR(I)/HINT1(6).lt.1.01) + 1 .or.(PEAR(I)/HINT1(7).gt.0.99.and.PEAR(I)/HINT1(7).lt.1.01)) then +c clin-5/2008: c TAUI=0.00001 TAUI=1.E-20 FTAR(I)=TAUI*COSH(RAP) GZAR(I)=TAUI*SINH(RAP)+zsmear + endif endif 1002 CONTINUE clin-3/2009 cleanup of program flow: @@ -2118,7 +2161,6 @@ subroutine fstate(iseed,srt,dm3,dm4,px,py,pz,iflag) cc SAVE /RNDF77/ SAVE - iseed=iseed iflag=-1 c iflag=-1: fail to find momenta c = 1: success @@ -2280,14 +2322,6 @@ subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt, 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) SAVE - iseed=iseed - dt=dt - nchrg=nchrg - nt=nt - ratiok=ratiok - px(1)=px(1) - py(1)=py(1) - pz(1)=pz(1) px1cm=pcx py1cm=pcy pz1cm=pcz @@ -2506,10 +2540,6 @@ subroutine pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2, COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - irun=irun - iseed=iseed - nt=nt - dt=dt px1cm=pcx py1cm=pcy pz1cm=pcz @@ -2676,11 +2706,6 @@ subroutine kaonN(brel,brsgm,irun,iseed,dt,nt, COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - dt=dt - irun=irun - iseed=iseed - nchrg=nchrg - nt=nt px1cm=pcx py1cm=pcy pz1cm=pcz @@ -2830,7 +2855,7 @@ SUBROUTINE ARTAN1 PARAMETER (YMT1 = -1.0, YMT2 = 1.0) cbz3/17/99 end c.....bin width for mt spectrum and y spectrum -clin-9/26/03 no symmetrization in y (or eta) for ana.dat: +c.... clin-9/26/03 no symmetrization in y (or eta) for ana dat files c PARAMETER (BMT = 0.05, BY = 0.2) PARAMETER (BMT = 0.05, BY = 0.4) COMMON /RUN/ NUM @@ -2856,9 +2881,6 @@ SUBROUTINE ARTAN1 & DY1PI0(50), DY1PR(50), DY1PB(50) & ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50) cc SAVE /ARANA1/ - - logical iwrite - data iwrite / .false. / SAVE cbz3/17/99 end @@ -2873,31 +2895,49 @@ SUBROUTINE ARTAN1 c 2/24/03 leptons and photons: if(xm.lt.0.01) goto 200 ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2) - eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5)) +clin-9/2012 determine pseudo-rapidity more generally: +c eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5)) + if((PX**2+PY**2).gt.0.) then + eta=asinh(PZ/sqrt(PX**2+PY**2)) + else + eta = 1000000.0*sign(1.,PZ) +clin-2/2013 for spectator target nucleons in LAB frame, +c note that finite precision of HBOOST +c would give spectator target nucleons a small but non-zero pz: + if(abs(pz).le.1e-3) eta=0. + endif XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) - IF (ABS(PZ) .GE. EE) THEN - if ( iwrite ) then - PRINT *, 'IN ARTAN1' - PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR' -cbzdbg2/16/99 - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY -cbzdbg2/16/99 -cbzdbg2/15/99 - PRINT *, ' PZ = ', PZ, ' EE = ', EE -cbzdbg2/16/99 - PRINT *, ' XM = ', XM -cbzdbg2/16/99end - endif - GOTO 200 -c STOP -cbzdbg2/15/99end - END IF DXMT = XMT - XM - Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. EE) THEN +c PRINT *, 'IN ARTAN1' +c PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR' +ccbzdbg2/16/99 +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +ccbzdbg2/16/99 +ccbzdbg2/15/99 +c PRINT *, ' PZ = ', PZ, ' EE = ', EE +ccbzdbg2/16/99 +c PRINT *, ' XM = ', XM +ccbzdbg2/16/99end +c GOTO 200 +c else +cc Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) +c Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ +1e-5)) +cc STOP +ccbzdbg2/15/99end +c END IF + if(XMT.gt.0.) then + Y=asinh(PZ/XMT) + else + PRINT *, ' IN ARTAN1 mt=0' + Y = 1000000.0*sign(1.,PZ) + endif + c.....rapidity cut for the rapidity distribution IF (ABS(Y) .GE. 10.0) GOTO 100 -clin-9/26/03 no symmetrization in y (or eta) for ana.dat: +c.....clin-9/26/03 no symmetrization in y (or eta) for ana dat files c IY = 1 + int(ABS(Y) / BY) c Ieta = 1 + int(ABS(eta) / BY) IF (ABS(eta) .GE. 10.0) GOTO 100 @@ -3054,9 +3094,6 @@ SUBROUTINE ARTAN2 & ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50) cbz3/17/99 end cc SAVE /ARANA2/ - - logical iwrite - data iwrite / .false. / SAVE DO 1002 J = 1, NUM @@ -3070,28 +3107,42 @@ SUBROUTINE ARTAN2 XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) c 2/24/03 leptons and photons: if(xm.lt.0.01) goto 200 + DXMT = XMT - XM ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2) - eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5)) +clin-9/2012 determine pseudo-rapidity more generally: +c eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5)) + if((PX**2+PY**2).gt.0.) then + eta=asinh(PZ/sqrt(PX**2+PY**2)) + else + eta = 1000000.0*sign(1.,PZ) +clin-2/2013 for spectator target nucleons in LAB frame: + if(abs(pz).le.1e-3) eta=0. + endif + +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. EE) THEN +c PRINT *, 'IN ARTAN2' +c PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR' +ccbzdbg2/16/99 +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +ccbzdbg2/16/99 +ccbzdbg2/15/99 +c PRINT *, ' PZ = ', PZ, ' EE = ', EE +ccbzdbg2/16/99 +c PRINT *, ' XM = ', XM +ccbzdbg2/16/99end +c GOTO 200 +cc STOP +ccbzdbg2/15/99end +c END IF +c Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5)) + if(XMT.gt.0.) then + Y=asinh(PZ/XMT) + else + PRINT *, ' IN ARTAN2 mt=0' + Y = 1000000.0*sign(1.,PZ) + endif - IF (ABS(PZ) .GE. EE) THEN - if ( iwrite ) then - PRINT *, 'IN ARTAN2' - PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR' -cbzdbg2/16/99 - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY -cbzdbg2/16/99 -cbzdbg2/15/99 - PRINT *, ' PZ = ', PZ, ' EE = ', EE -cbzdbg2/16/99 - PRINT *, ' XM = ', XM -cbzdbg2/16/99end - endif - GOTO 200 -c STOP -cbzdbg2/15/99end - END IF - DXMT = XMT - XM - Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) c.....rapidity cut for the rapidity distribution IF (ABS(Y) .GE. 10.0) GOTO 100 c IY = 1 + int(ABS(Y) / BY) @@ -3261,157 +3312,157 @@ SUBROUTINE ARTOUT(NEVNT) cc SAVE /ARANA2/ SAVE cbz3/17/99 end -cms OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN') -cms OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN') -cms OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN') -cms OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN') -cms OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN') -cms OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN') -cms OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN') -cms OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN') -cms OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN') +cms OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN') +cms OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN') +cms OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN') +cms OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN') +cms OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN') +cms OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN') +cms OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN') +cms OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN') +cms OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN') clin-4/24/03 c OPEN (39, FILE = 'ana/dndy_k0s.dat', STATUS = 'UNKNOWN') -cms OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN') -cms OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN') -cms OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN') -cms OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN') +cms OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN') +cms OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN') +cms OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN') +cms OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN') cbz3/17/99 -cms OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN') -cms OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN') -cms OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN') -cms OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN') -cms OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN') -cms OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN') +cms OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN') +cms OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN') +cms OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN') +cms OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN') +cms OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN') +cms OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN') cbz3/17/99 end -cms OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN') -cms OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN') -cms OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN') -cms OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN') -cms OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN') -cms OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN') -cms OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN') -cms OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN') -cms OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN') -clin-9/26/03 no symmetrization in y (or eta) for ana.dat: +cms OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN') +cms OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN') +cms OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN') +cms OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN') +cms OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN') +cms OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN') +cms OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN') +cms OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN') +cms OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN') +clin-9/26/03 no symmetrization in y (or eta) for ana dat files c SCALE1 = 1. / REAL(NEVNT * NUM) / BY / 2.0 SCALE1 = 1. / REAL(NEVNT * NUM) / BY SCALE2 = 1. / REAL(NEVNT * NUM) / BMT / (YMT2 - YMT1) c DO 1001 I = 1, 50 ymid=-10.+BY * (I - 0.5) -cms WRITE (30, 333) ymid, SCALE1 * dy1ntb(I) -cms WRITE (31, 333) ymid, SCALE1 * dy1ntp(I) -cms WRITE (32, 333) ymid, SCALE1 * DY1HM(I) -cms WRITE (37, 333) ymid, SCALE1 * DY1KP(I) -cms WRITE (38, 333) ymid, SCALE1 * DY1KM(I) -cms WRITE (39, 333) ymid, SCALE1 * DY1K0S(I) -cms WRITE (40, 333) ymid, SCALE1 * DY1LA(I) -cms WRITE (41, 333) ymid, SCALE1 * DY1LB(I) -cms WRITE (42, 333) ymid, SCALE1 * DY1PHI(I) -cms WRITE (33, 333) ymid, SCALE1 * DY1NEG(I) -cms WRITE (34, 333) ymid, SCALE1 * DY1CH(I) -cms WRITE (35, 333) ymid, SCALE1 * DE1NEG(I) -cms WRITE (36, 333) ymid, SCALE1 * DE1CH(I) -cms WRITE (43, 333) ymid, SCALE1 * dy1msn(I) -cms WRITE (44, 333) ymid, SCALE1 * DY1PIP(I) -cms WRITE (45, 333) ymid, SCALE1 * DY1PIM(I) -cms WRITE (46, 333) ymid, SCALE1 * DY1PI0(I) -cms WRITE (47, 333) ymid, SCALE1 * DY1PR(I) -cms WRITE (48, 333) ymid, SCALE1 * DY1PB(I) +cms WRITE (30, 333) ymid, SCALE1 * dy1ntb(I) +cms WRITE (31, 333) ymid, SCALE1 * dy1ntp(I) +cms WRITE (32, 333) ymid, SCALE1 * DY1HM(I) +cms WRITE (37, 333) ymid, SCALE1 * DY1KP(I) +cms WRITE (38, 333) ymid, SCALE1 * DY1KM(I) +cms WRITE (39, 333) ymid, SCALE1 * DY1K0S(I) +cms WRITE (40, 333) ymid, SCALE1 * DY1LA(I) +cms WRITE (41, 333) ymid, SCALE1 * DY1LB(I) +cms WRITE (42, 333) ymid, SCALE1 * DY1PHI(I) +cms WRITE (33, 333) ymid, SCALE1 * DY1NEG(I) +cms WRITE (34, 333) ymid, SCALE1 * DY1CH(I) +cms WRITE (35, 333) ymid, SCALE1 * DE1NEG(I) +cms WRITE (36, 333) ymid, SCALE1 * DE1CH(I) +cms WRITE (43, 333) ymid, SCALE1 * dy1msn(I) +cms WRITE (44, 333) ymid, SCALE1 * DY1PIP(I) +cms WRITE (45, 333) ymid, SCALE1 * DY1PIM(I) +cms WRITE (46, 333) ymid, SCALE1 * DY1PI0(I) +cms WRITE (47, 333) ymid, SCALE1 * DY1PR(I) +cms WRITE (48, 333) ymid, SCALE1 * DY1PB(I) IF (dm1pip(I) .NE. 0.0) THEN -cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I) +c WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I) END IF IF (dm1pim(I) .NE. 0.0) THEN -cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * -cms & dm1pim(I) +c WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * +c & dm1pim(I) END IF IF (DMT1PR(I) .NE. 0.0) THEN -cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I) +c WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I) END IF IF (DMT1PB(I) .NE. 0.0) THEN -cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I) +c WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I) END IF IF (DMT1KP(I) .NE. 0.0) THEN -cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I) +c WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I) END IF IF (dm1km(I) .NE. 0.0) THEN -cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * -cms & dm1km(I) +c WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * +c & dm1km(I) END IF IF (dm1k0s(I) .NE. 0.0) THEN -cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I) +c WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I) END IF IF (DMT1LA(I) .NE. 0.0) THEN -cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I) +c WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I) END IF IF (DMT1LB(I) .NE. 0.0) THEN -cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I) +c WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I) END IF 1001 CONTINUE c DO 1002 I = 30, 48 -cms WRITE (I, *) 'after hadron evolution' +c WRITE (I, *) 'after hadron evolution' 1002 CONTINUE DO 1003 I = 50, 58 -cms WRITE (I, *) 'after hadron evolution' +c WRITE (I, *) 'after hadron evolution' 1003 CONTINUE DO 1004 I = 1, 50 ymid=-10.+BY * (I - 0.5) -cms WRITE (30, 333) ymid, SCALE1 * dy2ntb(I) -cms WRITE (31, 333) ymid, SCALE1 * dy2ntp(I) -cms WRITE (32, 333) ymid, SCALE1 * DY2HM(I) -cms WRITE (37, 333) ymid, SCALE1 * DY2KP(I) -cms WRITE (38, 333) ymid, SCALE1 * DY2KM(I) -cms WRITE (39, 333) ymid, SCALE1 * DY2K0S(I) -cms WRITE (40, 333) ymid, SCALE1 * DY2LA(I) -cms WRITE (41, 333) ymid, SCALE1 * DY2LB(I) -cms WRITE (42, 333) ymid, SCALE1 * DY2PHI(I) -cms WRITE (33, 333) ymid, SCALE1 * DY2NEG(I) -cms WRITE (34, 333) ymid, SCALE1 * DY2CH(I) -cms WRITE (35, 333) ymid, SCALE1 * DE2NEG(I) -cms WRITE (36, 333) ymid, SCALE1 * DE2CH(I) -cms WRITE (43, 333) ymid, SCALE1 * dy2msn(I) -cms WRITE (44, 333) ymid, SCALE1 * DY2PIP(I) -cms WRITE (45, 333) ymid, SCALE1 * DY2PIM(I) -cms WRITE (46, 333) ymid, SCALE1 * DY2PI0(I) -cms WRITE (47, 333) ymid, SCALE1 * DY2PR(I) -cms WRITE (48, 333) ymid, SCALE1 * DY2PB(I) +cms WRITE (30, 333) ymid, SCALE1 * dy2ntb(I) +cms WRITE (31, 333) ymid, SCALE1 * dy2ntp(I) +cms WRITE (32, 333) ymid, SCALE1 * DY2HM(I) +cms WRITE (37, 333) ymid, SCALE1 * DY2KP(I) +cms WRITE (38, 333) ymid, SCALE1 * DY2KM(I) +cms WRITE (39, 333) ymid, SCALE1 * DY2K0S(I) +cms WRITE (40, 333) ymid, SCALE1 * DY2LA(I) +cms WRITE (41, 333) ymid, SCALE1 * DY2LB(I) +cms WRITE (42, 333) ymid, SCALE1 * DY2PHI(I) +cms WRITE (33, 333) ymid, SCALE1 * DY2NEG(I) +cms WRITE (34, 333) ymid, SCALE1 * DY2CH(I) +cms WRITE (35, 333) ymid, SCALE1 * DE2NEG(I) +cms WRITE (36, 333) ymid, SCALE1 * DE2CH(I) +cms WRITE (43, 333) ymid, SCALE1 * dy2msn(I) +cms WRITE (44, 333) ymid, SCALE1 * DY2PIP(I) +cms WRITE (45, 333) ymid, SCALE1 * DY2PIM(I) +cms WRITE (46, 333) ymid, SCALE1 * DY2PI0(I) +cms WRITE (47, 333) ymid, SCALE1 * DY2PR(I) +cms WRITE (48, 333) ymid, SCALE1 * DY2PB(I) c IF (dm2pip(I) .NE. 0.0) THEN -cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I) +c WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I) END IF IF (dm2pim(I) .NE. 0.0) THEN -cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * -cms & dm2pim(I) +c WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * +c & dm2pim(I) END IF IF (DMT2PR(I) .NE. 0.0) THEN -cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I) +c WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I) END IF IF (DMT2PB(I) .NE. 0.0) THEN -cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I) +c WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I) END IF IF (DMT2KP(I) .NE. 0.0) THEN -cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I) +c WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I) END IF IF (dm2km(I) .NE. 0.0) THEN -cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * -cms & dm2km(I) +c WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * +c & dm2km(I) END IF IF (dm2k0s(I) .NE. 0.0) THEN -cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I) +c WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I) END IF IF (DMT2LA(I) .NE. 0.0) THEN -cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I) +c WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I) END IF IF (DMT2LB(I) .NE. 0.0) THEN -cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I) +c WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I) END IF 1004 CONTINUE -cms 333 format(2(f12.5,1x)) + 333 format(2(f12.5,1x)) RETURN END @@ -3461,9 +3512,6 @@ SUBROUTINE HJANA1 cc SAVE /AREVT/ COMMON /AROUT/ IOUT cc SAVE /AROUT/ - - logical iwrite - data iwrite / .false. / SAVE DATA IW/0/ @@ -3524,20 +3572,27 @@ SUBROUTINE HJANA1 PZ = PJPZ(I, J) PE = PJPE(I, J) PM = PJPM(I, J) - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA1, PROJ STR ', I, ' PART ', J - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 200 - END IF - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA1, PROJ STR ', I, ' PART ', J +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 200 +c END IF +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 100 +clin-8/2014 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 100 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 100 dyp1(IY) = dyp1(IY) + 1.0 DEYP1(IY) = DEYP1(IY) + XMT IF (ITYP .EQ. 21) THEN @@ -3564,20 +3619,28 @@ SUBROUTINE HJANA1 PZ = PJTZ(I, J) PE = PJTE(I, J) PM = PJTM(I, J) - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA1, TARG STR ', I, ' PART ', J - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 400 - END IF - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA1, TARG STR ', I, ' PART ', J +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 400 +c END IF +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA1 mt=0' + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 300 +clin-8/2014 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 300 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 300 dyp1(IY) = dyp1(IY) + 1.0 DEYP1(IY) = DEYP1(IY) + XMT IF (ITYP .EQ. 21) THEN @@ -3604,20 +3667,28 @@ SUBROUTINE HJANA1 PZ = PZSG(I, J) PE = PESG(I, J) PM = PMSG(I, J) - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA1, INDP STR ', I, ' PART ', J - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 600 - END IF - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA1, INDP STR ', I, ' PART ', J +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 600 +c END IF +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA1 mt=0' + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 500 +clin-8/2014 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 500 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 500 dyp1(IY) = dyp1(IY) + 1.0 DEYP1(IY) = DEYP1(IY) + XMT IF (ITYP .EQ. 21) THEN @@ -3674,20 +3745,28 @@ SUBROUTINE HJANA1 PZ = sngl(PZ0(I)) PE = sngl(E0(I)) PM = sngl(XMASS0(I)) - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA1, GLUON ', I - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 800 - END IF - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA1, GLUON ', I +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 800 +c END IF +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA1 mt=0' + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 700 +clin-8/2014 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 700 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 700 dyg1c(IY) = dyg1c(IY) + 1.0 deyg1c(IY) = deyg1c(IY) + XMT 700 CONTINUE @@ -3939,9 +4018,6 @@ SUBROUTINE HJANA2 & GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3), & K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR) cc SAVE /SOFT/ - - logical iwrite - data iwrite / .false. / SAVE DATA IW/0/ @@ -4014,23 +4090,31 @@ SUBROUTINE HJANA2 PZ = PJPZ(I, J) PE = PJPE(I, J) PM = PJPM(I, J) -cbzdbg2/16/99 -c IF (ABS(PZ) .GE. PE) GOTO 200 - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA2, PROJ STR ', I, ' PART ', J - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 200 - END IF -cbzdbg2/16/99end - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +clin-9/2012 determine rapidity more generally: +ccbzdbg2/16/99 +cc IF (ABS(PZ) .GE. PE) GOTO 200 +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA2, PROJ STR ', I, ' PART ', J +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 200 +c END IF +ccbzdbg2/16/99end +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA2 mt=0' + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 100 +clin-8/2014 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 100 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 100 dyp2(IY) = dyp2(IY) + 1.0 DEYP2(IY) = DEYP2(IY) + XMT IF (ITYP .EQ. 21) THEN @@ -4057,23 +4141,31 @@ SUBROUTINE HJANA2 PZ = PJTZ(I, J) PE = PJTE(I, J) PM = PJTM(I, J) -cbzdbg2/16/99 -c IF (ABS(PZ) .GE. PE) GOTO 400 - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA2, TARG STR ', I, ' PART ', J - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 400 - END IF -cbzdbg2/16/99end - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +clin-9/2012 determine rapidity more generally: +ccbzdbg2/16/99 +cc IF (ABS(PZ) .GE. PE) GOTO 400 +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA2, TARG STR ', I, ' PART ', J +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 400 +c END IF +ccbzdbg2/16/99end +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA2 mt=0' + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 300 +clin-8/2014 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 300 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 300 dyp2(IY) = dyp2(IY) + 1.0 DEYP2(IY) = DEYP2(IY) + XMT IF (ITYP .EQ. 21) THEN @@ -4120,23 +4212,31 @@ SUBROUTINE HJANA2 endif clin-4/25/01-end -cbzdbg2/16/99 -c IF (ABS(PZ) .GE. PE) GOTO 600 - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA2, INDP STR ', I, ' PART ', J - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 600 - END IF -cbzdbg2/16/99end - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) +clin-9/2012 determine rapidity more generally: XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +ccbzdbg2/16/99 +cc IF (ABS(PZ) .GE. PE) GOTO 600 +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA2, INDP STR ', I, ' PART ', J +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 600 +c END IF +ccbzdbg2/16/99end +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA2 mt=0' + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 500 +clin-8/2014 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 500 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 500 dyp2(IY) = dyp2(IY) + 1.0 DEYP2(IY) = DEYP2(IY) + XMT IF (ITYP .EQ. 21) THEN @@ -4216,25 +4316,31 @@ SUBROUTINE HJANA2 PZ = sngl(PZ5(I)) PE = sngl(E5(I)) PM = sngl(XMASS5(I)) -cbzdbg2/16/99 -c IF (ABS(PZ) .GE. PE) GOTO 800 - - IF (ABS(PZ) .GE. PE) THEN - if ( iwrite ) then - PRINT *, ' IN HJANA2, GLUON ', I - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', PE - PRINT *, ' XM = ', PM - endif - GOTO 800 - END IF - -cbzdbg2/16/99end - RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) +clin-9/2012 determine rapidity more generally: XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) DXMT = XMT - PM +ccbzdbg2/16/99 +cc IF (ABS(PZ) .GE. PE) GOTO 800 +c IF (ABS(PZ) .GE. PE) THEN +c PRINT *, ' IN HJANA2, GLUON ', I +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', PE +c PRINT *, ' XM = ', PM +c GOTO 800 +c END IF +ccbzdbg2/16/99end +c RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5)) + if(XMT.gt.0.) then + RAP=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA2 mt=0' + RAP = 1000000.0*sign(1.,PZ) + endif + IY = 1 + int(ABS(RAP) / DY) - IF (IY .GT. 50) GOTO 700 +clin-9/2012 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 700 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 700 dyg2c(IY) = dyg2c(IY) + 1.0 deyg2c(IY) = deyg2c(IY) + XMT 700 CONTINUE @@ -4552,9 +4658,6 @@ SUBROUTINE HJANA3 cc SAVE /AROUT/ COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult cc SAVE /iflow/ - - logical iwrite - data iwrite / .false. / SAVE DATA IW/0/ @@ -4569,22 +4672,32 @@ SUBROUTINE HJANA3 EE = EE1(I, J) XM = XM1(I, J) XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) - IF (ABS(PZ) .GE. EE) THEN - if ( iwrite) then - PRINT *, 'IN HJANA3' - PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR' - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', EE - PRINT *, ' XM = ', XM - endif - GOTO 200 - END IF DXMT = XMT - XM - Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. EE) THEN +c PRINT *, 'IN HJANA3' +c PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR' +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', EE +c PRINT *, ' XM = ', XM +c GOTO 200 +c END IF +c Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5)) + if(XMT.gt.0.) then + Y=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA3 mt=0' + Y = 1000000.0*sign(1.,PZ) + endif + c.....rapidity cut for the rapidity distribution c IY = 1 + int(ABS(Y) / DY) - IY = 1 + int((Y+10.) / DY) - IF (IY .GT. 50) GOTO 100 +clin-8/2014 no rapidity shift here: +c IY = 1 + int((Y+10.) / DY) + IY = 1 + int(Y/DY) +clin-9/2012 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 100 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 100 dndyh3(IY) = dndyh3(IY) + 1.0 DEYH3(IY) = DEYH3(IY) + XMT 100 CONTINUE @@ -4628,9 +4741,6 @@ SUBROUTINE HJANA4 cc SAVE /AROUT/ COMMON /fflow/ v2f,etf,xmultf,v2fpi,xmulpi cc SAVE /fflow/ - - logical iwrite - data iwrite / .false. / SAVE DATA IW/0/ @@ -4645,22 +4755,32 @@ SUBROUTINE HJANA4 EE = EE1(I, J) XM = XM1(I, J) XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) - IF (ABS(PZ) .GE. EE) THEN - if ( iwrite ) then - PRINT *, 'IN HJANA4' - PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR' - PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY - PRINT *, ' PZ = ', PZ, ' EE = ', EE - PRINT *, ' XM = ', XM - endif - GOTO 200 - END IF DXMT = XMT - XM - Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) +clin-9/2012 determine rapidity more generally: +c IF (ABS(PZ) .GE. EE) THEN +c PRINT *, 'IN HJANA4' +c PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR' +c PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY +c PRINT *, ' PZ = ', PZ, ' EE = ', EE +c PRINT *, ' XM = ', XM +c GOTO 200 +c END IF +c Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5)) + if(XMT.gt.0.) then + Y=asinh(PZ/XMT) + else + PRINT *, ' IN HJANA4 mt=0' + Y = 1000000.0*sign(1.,PZ) + endif + c.....rapidity cut for the rapidity distribution c IY = 1 + int(ABS(Y) / DY) - IY = 1 + int((Y+10.) / DY) - IF (IY .GT. 50) GOTO 100 +clin-8/2014 no rapidity shift here: +c IY = 1 + int((Y+10.) / DY) + IY = 1 + int(Y/DY) +clin-9/2012 prevent possible segmentation fault (due to IY<=0): +c IF (IY .GT. 50) GOTO 100 + IF (IY.lt.1 .or.IY .GT. 50) GOTO 100 dndyh4(IY) = dndyh4(IY) + 1.0 DEYH4(IY) = DEYH4(IY) + XMT 100 CONTINUE @@ -4822,6 +4942,17 @@ SUBROUTINE zpstrg RETURN END +clin-10/01/03 random number generator for f77: +c function RANART(NSEED) +c SAVE +clin-4/2008 ran(nseed) is renamed to avoid conflict with system functions: +c ran=rand() +c ranart=rand(0) +c one may also use the following random number generator in PYTHIA/JETSET: +c ranart=rlu(0) +c return +c end + clin-3/2009 c Initialize hadron weights; c Can add initial hadrons before the hadron cascade starts (but after ZPC). @@ -4840,6 +4971,7 @@ subroutine addhad COMMON/RNDF77/NSEED common /para8/ idpert,npertd,idxsec SAVE + c All hadrons at the start of hadron cascade have the weight of 1 c except those inserted by the user in this subroutine: np0=IAINT2(1) @@ -4851,7 +4983,9 @@ subroutine addhad tau0=ARPAR1(1) DO 100 i=np0+1,np0+nadd ITYPAR(I)=42 - dpertp(I)=1d0/dble(nadd) +clin-5/2012 fix type mismatch: +c dpertp(I)=1d0/dble(nadd) + dpertp(I)=1./float(nadd) GXAR(I)=5.*(1.-2.*RANART(NSEED)) GYAR(I)=5.*(1.-2.*RANART(NSEED)) GZAR(I)=2.*(1.-2.*RANART(NSEED)) @@ -4862,7 +4996,10 @@ subroutine addhad XMAR(I)=xmd c PEAR(I)=sqrt(PXAR(I)**2+PYAR(I)**2+PZAR(I)**2+XMAR(I)**2) - RAP=0.5*alog((PEAR(I)+PZAR(I))/(PEAR(I)-PZAR(I))) +clin-9/2012 determine rapidity more generally: +c RAP=0.5*alog((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5)) + RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2)) +c VX=PXAR(I)/PEAR(I) VY=PYAR(I)/PEAR(I) c.....give initial formation time shift and boost according to rapidity: @@ -4890,3 +5027,15 @@ subroutine addhad c return end + +clin-8/2014 define function asinh(): + FUNCTION asinh(x) + SAVE + if(x.gt.0) then + ASINH=alog(x+sqrt(x**2+1.)) + else +c a la suggestion de YP Liu: + ASINH=-alog(-x+sqrt(x**2+1.)) + endif + return + end diff --git a/GeneratorInterface/AMPTInterface/src/art1f.f b/GeneratorInterface/AMPTInterface/src/art1f.f index 449beb1030fb8..669c9e0244824 100755 --- a/GeneratorInterface/AMPTInterface/src/art1f.f +++ b/GeneratorInterface/AMPTInterface/src/art1f.f @@ -274,7 +274,7 @@ SUBROUTINE ARTMN cbz11/16/98end common /lastt/itimeh,bimp cc SAVE /lastt/ - common/snn/efrm,npart1,npart2 + common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG cc SAVE /snn/ COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast cc SAVE /hbt/ @@ -581,7 +581,6 @@ SUBROUTINE ARTMN cbz11/16/98end DO 10000 NT = 1,NTMAX - *TEMPORARY PARTICLE COUNTERS *4.2 PION COUNTERS : LP1,LP2 AND LP3 ARE THE NO. OF P+,P0 AND P- LP1=0 @@ -1624,7 +1623,7 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR) dimension ftpisv(MAXSTR,MAXR),fttemp(MAXSTR) common /dpi/em2,lb2 - common/phidcy/iphidcy,pttrig,ntrig,maxmiss + common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy clin-5/2008: DIMENSION dptemp(MAXSTR) common /para8/ idpert,npertd,idxsec @@ -1748,12 +1747,15 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, DO 800 J1 = J10,MASSR(IRUN) I1 = J1 + MSUM * E(I)=0 are for pions having been absorbed or photons which do not enter here: - IF(E(I1).EQ.0.)GO TO 800 - +clin-4/2012 option of pi0 decays: +c IF(E(I1).EQ.0.)GO TO 800 + IF(E(I1).EQ.0.)GO TO 798 c To include anti-(Delta,N*1440 and N*1535): c IF ((LB(I1) .LT. -13 .OR. LB(I1) .GT. 28) c 1 .and.iabs(LB(I1)) .ne. 30 ) GOTO 800 - IF (LB(I1) .LT. -45 .OR. LB(I1) .GT. 45) GOTO 800 +clin-4/2012 option of pi0 decays: +c IF (LB(I1) .LT. -45 .OR. LB(I1) .GT. 45) GOTO 800 + IF (LB(I1) .LT. -45 .OR. LB(I1) .GT. 45) GOTO 798 X1 = R(1,I1) Y1 = R(2,I1) Z1 = R(3,I1) @@ -1787,7 +1789,10 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, & .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30 & .or.(iabs(lb1).ge.6.and.iabs(lb1).le.13) & .or.(iksdcy.eq.1.and.lb1.eq.24) - & .or.iabs(lb1).eq.16) then + & .or.iabs(lb1).eq.16 + & .or.(ipi0dcy.eq.1.and.nt.eq.ntmax.and.lb1.eq.4)) then +clin-4/2012-above for option of pi0 decay: +c & .or.iabs(lb1).eq.16) then continue else goto 1 @@ -1820,6 +1825,9 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, WID=W1440(EM1) ELSEIF((iabs(LB1).EQ.12).OR.(iabs(LB1).EQ.13)) then WID=W1535(EM1) +clin-4/2012 for option of pi0 decay: + ELSEIF(ipi0dcy.eq.1.and.nt.eq.ntmax.and.lb1.eq.4) then + wid=7.85e-9 ENDIF * if it is the last time step, FORCE all resonance to strong-decay @@ -1828,6 +1836,10 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, pdecay=1.1 clin-5b/2008 forbid phi decay at the end of hadronic cascade: if(iphidcy.eq.0.and.iabs(LB1).eq.29) pdecay=0. +ctest off clin-9/2012 forbid long-time decays (eta,omega,K*,Sigma0) +c at the end of hadronic cascade to analyze freezeout time: +c if(LB1.eq.0.or.LB1.eq.28.or.iabs(LB1).eq.30 +c 1 .or.iabs(LB1).eq.16) pdecay=0. else T0=0.19733/WID GFACTR=E1/EM1 @@ -1859,7 +1871,10 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, & .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30 & .or.(iabs(lb1).ge.6.and.iabs(lb1).le.9) & .or.(iksdcy.eq.1.and.lb1.eq.24) - & .or.iabs(lb1).eq.16) then + & .or.iabs(lb1).eq.16 + & .or.(ipi0dcy.eq.1.and.nt.eq.ntmax.and.lb1.eq.4)) then +clin-4/2012 Above for option of pi0 decay: +c & .or.iabs(lb1).eq.16) then c previous rho decay performed in rhodecay(): c nnn=nnn+1 c call rhodecay(idecay,i1,nnn,iseed) @@ -1867,7 +1882,10 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, ctest off record decays of phi,K*,Lambda(1520) resonances: c if(lb1.eq.29.or.iabs(lb1).eq.30) c 1 write(18,112) 'decay',lb1,px1,py1,pz1,am1,nt - call resdec(i1,nt,nnn,wid,idecay) +c +clin-4/2012 option of pi0 decays: +c call resdec(i1,nt,nnn,wid,idecay) + call resdec(i1,nt,nnn,wid,idecay,0) p(1,i1)=px1n p(2,i1)=py1n p(3,i1)=pz1n @@ -1932,7 +1950,12 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, elseif(lb(i1).eq.0) then wid=1.18e-6 elseif(lb(i1).eq.24.and.iksdcy.eq.1) then - wid=7.36e-17 +clin-4/2012 corrected K0s decay width: +c wid=7.36e-17 + wid=7.36e-15 +clin-4/2012 option of pi0 decays: + elseif(ipi0dcy.eq.1.and.lb(i1).eq.4) then + wid=7.85e-9 else goto 9000 endif @@ -1942,7 +1965,9 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, PZ1=P(3,I1) EM1=E(I1) E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2) - call resdec(i1,nt,nnn,wid,idecay) +clin-4/2012 option of pi0 decays: +c call resdec(i1,nt,nnn,wid,idecay) + call resdec(i1,nt,nnn,wid,idecay,0) p(1,i1)=px1n p(2,i1)=py1n p(3,i1)=pz1n @@ -1954,12 +1979,38 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, dpertp(i1)=dp1n endif +c Decay daughter of the above decay in lb(i1) may be a pi0: + if(nt.eq.ntmax.and.ipi0dcy.eq.1.and.lb(i1).eq.4) then + wid=7.85e-9 + LB1=LB(I1) + PX1=P(1,I1) + PY1=P(2,I1) + PZ1=P(3,I1) + EM1=E(I1) + E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2) + call resdec(i1,nt,nnn,wid,idecay,0) + p(1,i1)=px1n + p(2,i1)=py1n + p(3,i1)=pz1n + R(1,i1)=xfnl + R(2,i1)=yfnl + R(3,i1)=zfnl + tfdcy(i1)=tfnl + dpertp(i1)=dp1n + endif + * negelecting the Pauli blocking at high energies - 9000 go to 800 +clin-4/2012 option of pi0 decays: +c 9000 go to 800 + 9000 go to 798 + ENDIF * LOOP OVER ALL PSEUDOPARTICLES 2 IN THE SAME RUN * SAVE ALL THE COORDINATES FOR POSSIBLE CHANGE IN THE FOLLOWING COLLISION - 1 if(nt.eq.ntmax)go to 800 +clin-4/2012 option of pi0 decays: +c 1 if(nt.eq.ntmax)go to 800 + 1 if(nt.eq.ntmax)go to 798 + X1 = R(1,I1) Y1 = R(2,I1) Z1 = R(3,I1) @@ -1970,6 +2021,7 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, IF(E(I2).EQ.0.) GO TO 600 clin-5/2008 in case the first particle is already destroyed: IF(E(I1).EQ.0.) GO TO 800 +clin-4/2012 option of pi0 decays: IF (LB(I2) .LT. -45 .OR. LB(I2) .GT. 45) GOTO 600 clin-7/26/03 improve speed X2=R(1,I2) @@ -2222,7 +2274,8 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, ipz2 = nint(pz2/dpz) * FIND MOMENTA OF PARTICLES IN THE CMS OF THE TWO COLLIDING PARTICLES * AND THE CMS ENERGY SRT - CALL CMS(I1,I2,PCX,PCY,PCZ,SRT) + CALL CMS(I1,I2,PCX,PCY,PCZ,SRT) + clin-7/26/03 improve speed drmax=dr0max call distc0(drmax,deltr0,DT, @@ -4325,6 +4378,28 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, c endif 600 CONTINUE + +clin-4/2012 option of pi0 decays: +c particles in lpion() may be a pi0, and when ipi0dcy=1 +c we need to decay them at nt=ntmax after all lb(i1) decays are done: + 798 if(nt.eq.ntmax.and.ipi0dcy.eq.1 + 1 .and.i1.eq.(MASSR(IRUN)+MSUM)) then + do ipion=1,NNN + if(LPION(ipion,IRUN).eq.4) then + wid=7.85e-9 + call resdec(i1,nt,nnn,wid,idecay,ipion) + endif + enddo + endif +ctest off +c if(nt.eq.ntmax.and.i1.eq.(MASSR(IRUN)+MSUM)) then +c do ip=1,i1 +c write(98,*) lb(ip),e(ip),ip +c enddo +c endif + +clin-4/2012 option of pi0 decays-end + 800 CONTINUE * RELABLE MESONS LEFT IN THIS RUN EXCLUDING THOSE BEING CREATED DURING * THIS TIME STEP AND COUNT THE TOTAL NO. OF PARTICLES IN THIS RUN @@ -4459,6 +4534,53 @@ SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk, c RETURN END + +clin-9/2012: use double precision for S in CMS(): to avoid crash +c (segmentation fault due to s<0, which happened at high energies +c such as LHC with large NTMAX for two almost-comoving hadrons +c that have small Pt but large |Pz|): +**************************************** +c SUBROUTINE CMS(I1,I2,PX1CM,PY1CM,PZ1CM,SRT) +* PURPOSE : FIND THE MOMENTA OF PARTICLES IN THE CMS OF THE +* TWO COLLIDING PARTICLES +* VARIABLES : +***************************************** +c PARAMETER (MAXSTR=150001) +c COMMON /AA/ R(3,MAXSTR) +ccc SAVE /AA/ +c COMMON /BB/ P(3,MAXSTR) +ccc SAVE /BB/ +c COMMON /CC/ E(MAXSTR) +ccc SAVE /CC/ +c COMMON /BG/ BETAX,BETAY,BETAZ,GAMMA +ccc SAVE /BG/ +c SAVE +c PX1=P(1,I1) +c PY1=P(2,I1) +c PZ1=P(3,I1) +c PX2=P(1,I2) +c PY2=P(2,I2) +c PZ2=P(3,I2) +c EM1=E(I1) +c EM2=E(I2) +c E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2) +c E2=SQRT(EM2**2 + PX2**2 + PY2**2 + PZ2**2 ) +c S=(E1+E2)**2-(PX1+PX2)**2-(PY1+PY2)**2-(PZ1+PZ2)**2 +c SRT=SQRT(S) +c*LORENTZ-TRANSFORMATION IN I1-I2-C.M. SYSTEM +c ETOTAL = E1 + E2 +c BETAX = (PX1+PX2) / ETOTAL +c BETAY = (PY1+PY2) / ETOTAL +c BETAZ = (PZ1+PZ2) / ETOTAL +c GAMMA = 1.0 / SQRT(1.0-BETAX**2-BETAY**2-BETAZ**2) +c*TRANSFORMATION OF MOMENTA (PX1CM = - PX2CM) +c P1BETA = PX1*BETAX + PY1*BETAY + PZ1 * BETAZ +c TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) - E1 ) +c PX1CM = BETAX * TRANSF + PX1 +c PY1CM = BETAY * TRANSF + PY1 +c PZ1CM = BETAZ * TRANSF + PZ1 +c RETURN +c END **************************************** SUBROUTINE CMS(I1,I2,PX1CM,PY1CM,PZ1CM,SRT) * PURPOSE : FIND THE MOMENTA OF PARTICLES IN THE CMS OF THE @@ -4466,41 +4588,51 @@ SUBROUTINE CMS(I1,I2,PX1CM,PY1CM,PZ1CM,SRT) * VARIABLES : ***************************************** PARAMETER (MAXSTR=150001) - COMMON /AA/ R(3,MAXSTR) -cc SAVE /AA/ + double precision px1,py1,pz1,px2,py2,pz2,em1,em2,e1,e2, + 1 s,ETOTAL,P1BETA,TRANSF,dBETAX,dBETAY,dBETAZ,dGAMMA,scheck COMMON /BB/ P(3,MAXSTR) -cc SAVE /BB/ COMMON /CC/ E(MAXSTR) -cc SAVE /CC/ COMMON /BG/ BETAX,BETAY,BETAZ,GAMMA -cc SAVE /BG/ - SAVE - PX1=P(1,I1) - PY1=P(2,I1) - PZ1=P(3,I1) - PX2=P(1,I2) - PY2=P(2,I2) - PZ2=P(3,I2) - EM1=E(I1) - EM2=E(I2) - E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2) - E2=SQRT(EM2**2 + PX2**2 + PY2**2 + PZ2**2 ) + SAVE + PX1=dble(P(1,I1)) + PY1=dble(P(2,I1)) + PZ1=dble(P(3,I1)) + PX2=dble(P(1,I2)) + PY2=dble(P(2,I2)) + PZ2=dble(P(3,I2)) + EM1=dble(E(I1)) + EM2=dble(E(I2)) + E1=dSQRT(EM1**2+PX1**2+PY1**2+PZ1**2) + E2=dSQRT(EM2**2+PX2**2+PY2**2+PZ2**2) S=(E1+E2)**2-(PX1+PX2)**2-(PY1+PY2)**2-(PZ1+PZ2)**2 - SRT=SQRT(S) + IF(S.LE.0) S=0d0 + SRT=sngl(dSQRT(S)) *LORENTZ-TRANSFORMATION IN I1-I2-C.M. SYSTEM - ETOTAL = E1 + E2 - BETAX = (PX1+PX2) / ETOTAL - BETAY = (PY1+PY2) / ETOTAL - BETAZ = (PZ1+PZ2) / ETOTAL - GAMMA = 1.0 / SQRT(1.0-BETAX**2-BETAY**2-BETAZ**2) + ETOTAL = E1 + E2 + dBETAX = (PX1+PX2) / ETOTAL + dBETAY = (PY1+PY2) / ETOTAL + dBETAZ = (PZ1+PZ2) / ETOTAL +clin-9/2012: check argument in sqrt(): + scheck=1.d0-dBETAX**2-dBETAY**2-dBETAZ**2 + if(scheck.le.0d0) then + write(99,*) 'scheck1: ', scheck + stop + endif + dGAMMA=1.d0/dSQRT(scheck) *TRANSFORMATION OF MOMENTA (PX1CM = - PX2CM) - P1BETA = PX1*BETAX + PY1*BETAY + PZ1 * BETAZ - TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) - E1 ) - PX1CM = BETAX * TRANSF + PX1 - PY1CM = BETAY * TRANSF + PY1 - PZ1CM = BETAZ * TRANSF + PZ1 - RETURN - END + P1BETA = PX1*dBETAX + PY1*dBETAY + PZ1 * dBETAZ + TRANSF = dGAMMA * ( dGAMMA * P1BETA / (dGAMMA + 1d0) - E1 ) + PX1CM = sngl(dBETAX * TRANSF + PX1) + PY1CM = sngl(dBETAY * TRANSF + PY1) + PZ1CM = sngl(dBETAZ * TRANSF + PZ1) + BETAX = sngl(dBETAX) + BETAY = sngl(dBETAY) + BETAZ = sngl(dBETAZ) + GAMMA = sngl(dGAMMA) + RETURN + END +clin-9/2012-end + *************************************** SUBROUTINE DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT 1 ,IC,PX1CM,PY1CM,PZ1CM) @@ -5338,7 +5470,15 @@ SUBROUTINE CRNN(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK, cc1=ptr(xptr,iseed) clin-10/25/02-end - c1=sqrt(pr**2-cc1**2)/pr +clin-9/2012: check argument in sqrt(): + scheck=pr**2-cc1**2 + if(scheck.lt.0) then + write(99,*) 'scheck2: ', scheck + scheck=0. + endif + c1=sqrt(scheck)/pr +c c1=sqrt(pr**2-cc1**2)/pr + endif T1 = 2.0 * PI * RANART(NSEED) if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then @@ -6406,7 +6546,16 @@ SUBROUTINE CRNN(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK, S1 = 1.0 - C1**2 IF(S1.LE.0)S1=0 S1=SQRT(S1) - S2 = SQRT( 1.0 - C2**2 ) + +clin-9/2012: check argument in sqrt(): + scheck=1.0 - C2**2 + if(scheck.lt.0) then + write(99,*) 'scheck3: ', scheck + scheck=0. + endif + S2=SQRT(scheck) +c S2 = SQRT( 1.0 - C2**2 ) + CT1 = COS(T1) ST1 = SIN(T1) CT2 = COS(T2) @@ -6578,7 +6727,7 @@ SUBROUTINE CRPP(PX,PY,PZ,SRT,I1,I2,IBLOCK, PR=SQRT(PR2)/(2.*SRT) C1 = 1.0 - 2.0 * RANART(NSEED) T1 = 2.0 * PI * RANART(NSEED) - S1 = SQRT( 1.0 - C1**2 ) + S1 = SQRT( 1.0 - C1**2 ) CT1 = COS(T1) ST1 = SIN(T1) PZ = PR * C1 @@ -7526,7 +7675,15 @@ SUBROUTINE CRND(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK, cc1=ptr(xptr,iseed) clin-10/25/02-end - c1=sqrt(pr**2-cc1**2)/pr +clin-9/2012: check argument in sqrt(): + scheck=pr**2-cc1**2 + if(scheck.lt.0) then + write(99,*) 'scheck4: ', scheck + scheck=0. + endif + c1=sqrt(scheck)/pr +c c1=sqrt(pr**2-cc1**2)/pr + endif T1 = 2.0 * PI * RANART(NSEED) IBLOCK=3 @@ -7543,8 +7700,25 @@ SUBROUTINE CRND(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK, ELSE T2=ATAN2(PY,PX) END IF - S1 = SQRT( 1.0 - C1**2 ) - S2 = SQRT( 1.0 - C2**2 ) + +clin-9/2012: check argument in sqrt(): + scheck=1.0 - C1**2 + if(scheck.lt.0) then + write(99,*) 'scheck5: ', scheck + scheck=0. + endif + S1=SQRT(scheck) +c S1 = SQRT( 1.0 - C1**2 ) + +clin-9/2012: check argument in sqrt(): + scheck=1.0 - C2**2 + if(scheck.lt.0) then + write(99,*) 'scheck6: ', scheck + scheck=0. + endif + S2=SQRT(scheck) +c S2 = SQRT( 1.0 - C2**2 ) + CT1 = COS(T1) ST1 = SIN(T1) CT2 = COS(T2) @@ -8815,7 +8989,15 @@ SUBROUTINE CRDD(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK, cc1=ptr(xptr,iseed) clin-10/25/02-end - c1=sqrt(pr**2-cc1**2)/pr +clin-9/2012: check argument in sqrt(): + scheck=pr**2-cc1**2 + if(scheck.lt.0) then + write(99,*) 'scheck7: ', scheck + scheck=0. + endif + c1=sqrt(scheck)/pr +c c1=sqrt(pr**2-cc1**2)/pr + endif T1 = 2.0 * PI * RANART(NSEED) if(ianti.eq.1 .and. lb(i1).ge.1 .and. lb(i2).ge.1)then @@ -8824,8 +9006,25 @@ SUBROUTINE CRDD(IRUN,PX,PY,PZ,SRT,I1,I2,IBLOCK, endif ENDIF *COM: SET THE NEW MOMENTUM COORDINATES -107 S1 = SQRT( 1.0 - C1**2 ) - S2 = SQRT( 1.0 - C2**2 ) + +clin-9/2012: check argument in sqrt(): + 107 scheck=1.0 - C1**2 + if(scheck.lt.0) then + write(99,*) 'scheck8: ', scheck + scheck=0. + endif + S1=SQRT(scheck) +c107 S1 = SQRT( 1.0 - C1**2 ) + +clin-9/2012: check argument in sqrt(): + scheck=1.0 - C2**2 + if(scheck.lt.0) then + write(99,*) 'scheck9: ', scheck + scheck=0. + endif + S2=SQRT(scheck) +c S2 = SQRT( 1.0 - C2**2 ) + CT1 = COS(T1) ST1 = SIN(T1) CT2 = COS(T2) @@ -9184,13 +9383,21 @@ SUBROUTINE INIT(MINNUM,MAXNUM,NUM,RADIUS,X0,Z0,P0, *---------------------------------------------------------------------- * PREPARATION FOR LORENTZ-TRANSFORMATIONS * - ISEED=ISEED IF (P0 .NE. 0.) THEN SIGN = P0 / ABS(P0) ELSE SIGN = 0. END IF - BETA = SIGN * SQRT(GAMMA**2-1.)/GAMMA + +clin-9/2012: check argument in sqrt(): + scheck=GAMMA**2-1. + if(scheck.lt.0) then + write(99,*) 'scheck10: ', scheck + scheck=0. + endif + BETA=SIGN*SQRT(scheck)/GAMMA +c BETA = SIGN * SQRT(GAMMA**2-1.)/GAMMA + *----------------------------------------------------------------------- * TARGET-ID = 1 AND PROJECTILE-ID = -1 * @@ -9885,12 +10092,30 @@ REAL FUNCTION FD5(DMASS,SRT,CON) FD=4.*(AM0**2)*W1535(DMASS)/((DMASS**2-1.535**2)**2 1 +AM0**2*W1535(DMASS)**2) IF(CON.EQ.1.)THEN - P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 - 1 /(4.*SRT**2)-DMASS**2) + +clin-9/2012: check argument in sqrt(): + scheck=(SRT**2+DMASS**2-AMN**2)**2/(4.*SRT**2)-DMASS**2 + if(scheck.lt.0) then + write(99,*) 'scheck11: ', scheck + scheck=0. + endif + P1=SQRT(scheck) +c P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 +c 1 /(4.*SRT**2)-DMASS**2) + ELSE DMASS=AMN+AVPI - P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 - 1 /(4.*SRT**2)-DMASS**2) + +clin-9/2012: check argument in sqrt(): + scheck=(SRT**2+DMASS**2-AMN**2)**2/(4.*SRT**2)-DMASS**2 + if(scheck.lt.0) then + write(99,*) 'scheck12: ', scheck + scheck=0. + endif + P1=SQRT(scheck) +c P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 +c 1 /(4.*SRT**2)-DMASS**2) + ENDIF FD5=FD*P1*DMASS RETURN @@ -9906,12 +10131,29 @@ REAL FUNCTION FNS(DMASS,SRT,CON) AN0=1.43 FN=4.*(AN0**2)*WIDTH/((DMASS**2-1.44**2)**2+AN0**2*WIDTH**2) IF(CON.EQ.1.)THEN - P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 - 1 /(4.*SRT**2)-DMASS**2) + +clin-9/2012: check argument in sqrt(): + scheck=(SRT**2+DMASS**2-AMN**2)**2/(4.*SRT**2)-DMASS**2 + if(scheck.lt.0) then + write(99,*) 'scheck13: ', scheck + scheck=0. + endif + P1=SQRT(scheck) +c P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 +c 1 /(4.*SRT**2)-DMASS**2) + ELSE DMASS=AMN+AVPI - P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 - 1 /(4.*SRT**2)-DMASS**2) +clin-9/2012: check argument in sqrt(): + scheck=(SRT**2+DMASS**2-AMN**2)**2/(4.*SRT**2)-DMASS**2 + if(scheck.lt.0) then + write(99,*) 'scheck14: ', scheck + scheck=0. + endif + P1=SQRT(scheck) +c P1=SQRT((SRT**2+DMASS**2-AMN**2)**2 +c 1 /(4.*SRT**2)-DMASS**2) + ENDIF FNS=FN*P1*DMASS RETURN @@ -10142,7 +10384,6 @@ SUBROUTINE DKINE(IRUN,I,NNN,NLAB,ISEED,wid,nt) 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) EXTERNAL IARFLV, INVFLV SAVE - ISEED=ISEED * READ IN THE COORDINATES OF DELTA OR N* UNDERGOING DECAY PX=P(1,I) PY=P(2,I) @@ -10235,9 +10476,9 @@ SUBROUTINE DKINE(IRUN,I,NNN,NLAB,ISEED,wid,nt) tfdpi(NNN,IRUN)=tfnl endif -cc 200 format(a30,2(1x,e10.4)) -cc 210 format(i6,5(1x,f8.3)) -cc 220 format(a2,i5,5(1x,f8.3)) + 200 format(a30,2(1x,e10.4)) + 210 format(i6,5(1x,f8.3)) + 220 format(a2,i5,5(1x,f8.3)) RETURN END @@ -10412,8 +10653,7 @@ SUBROUTINE DKINE2(IRUN,I,NNN,NLAB,ISEED,wid,nt) 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) SAVE - - ISEED=ISEED + * READ IN THE COORDINATES OF THE N*(1440) UNDERGOING DECAY PX=P(1,I) PY=P(2,I) @@ -10429,7 +10669,16 @@ SUBROUTINE DKINE2(IRUN,I,NNN,NLAB,ISEED,wid,nt) IF(NLAB.EQ.1)AM=AMP * THE MAXIMUM MOMENTUM OF THE NUCLEON FROM THE DECAY OF A N* PMAX2=(DM**2-(AM+PM1+PM2)**2)*(DM**2-(AM-PM1-PM2)**2)/4/DM**2 - PMAX=SQRT(PMAX2) + +clin-9/2012: check argument in sqrt(): + scheck=PMAX2 + if(scheck.lt.0) then + write(99,*) 'scheck15: ', scheck + scheck=0. + endif + PMAX=SQRT(scheck) +c PMAX=SQRT(PMAX2) + * GENERATE THE MOMENTUM OF THE NUCLEON IN THE N* REST FRAME CSS=1.-2.*RANART(NSEED) SSS=SQRT(1-CSS**2) @@ -10444,7 +10693,16 @@ SUBROUTINE DKINE2(IRUN,I,NNN,NLAB,ISEED,wid,nt) BETAX=-PX0/(DM-EP0) BETAY=-PY0/(DM-EP0) BETAZ=-PZ0/(DM-EP0) - GD1=1./SQRT(1-BETAX**2-BETAY**2-BETAZ**2) + +clin-9/2012: check argument in sqrt(): + scheck=1-BETAX**2-BETAY**2-BETAZ**2 + if(scheck.le.0) then + write(99,*) 'scheck16: ', scheck + stop + endif + GD1=1./SQRT(scheck) +c GD1=1./SQRT(1-BETAX**2-BETAY**2-BETAZ**2) + FGD1=GD1/(1+GD1) * GENERATE THE MOMENTA OF PIONS IN THE CMS OF PION+PION- Q2=((DM-EP0)/(2.*GD1))**2-PM1**2 @@ -10574,9 +10832,9 @@ SUBROUTINE DKINE2(IRUN,I,NNN,NLAB,ISEED,wid,nt) tfdpi(NNN+1,IRUN)=tfnl endif -cc 200 format(a30,2(1x,e10.4)) -cc 210 format(i6,5(1x,f8.3)) -cc 220 format(a2,i5,5(1x,f8.3)) + 200 format(a30,2(1x,e10.4)) + 210 format(i6,5(1x,f8.3)) + 220 format(a2,i5,5(1x,f8.3)) RETURN END @@ -10591,6 +10849,8 @@ SUBROUTINE DRESON(I1,I2) PARAMETER (MAXSTR=150001,MAXR=1, 1 AMN=0.939457,AMP=0.93828, 2 AP1=0.13496,AP2=0.13957,AM0=1.232,PI=3.1415926) +clin-9/2012: improve precision for argument in sqrt(): + double precision e10,e20,scheck,p1,p2,p3 COMMON /AA/ R(3,MAXSTR) cc SAVE /AA/ COMMON /BB/ P(3,MAXSTR) @@ -10611,8 +10871,17 @@ SUBROUTINE DRESON(I1,I2) cc SAVE /PD/ SAVE * 1. DETERMINE THE MOMENTUM COMPONENT OF DELTA/N* IN THE LAB. FRAME - E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) - E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) +clin-9/2012: improve precision for argument in sqrt(): +c E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) +c E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) + E10=dSQRT(dble(E(I1))**2+dble(P(1,I1))**2 + 1 +dble(P(2,I1))**2+dble(P(3,I1))**2) + E20=dSQRT(dble(E(I2))**2+dble(P(1,I2))**2 + 1 +dble(P(2,I2))**2+dble(P(3,I2))**2) + p1=dble(P(1,I1))+dble(P(1,I2)) + p2=dble(P(2,I1))+dble(P(2,I2)) + p3=dble(P(3,I1))+dble(P(3,I2)) + IF(iabs(LB(I2)) .EQ. 1 .OR. iabs(LB(I2)) .EQ. 2 .OR. & (iabs(LB(I2)) .GE. 6 .AND. iabs(LB(I2)) .LE. 17)) THEN E(I1)=0. @@ -10625,7 +10894,19 @@ SUBROUTINE DRESON(I1,I2) P(2,I)=P(2,I1)+P(2,I2) P(3,I)=P(3,I1)+P(3,I2) * 2. DETERMINE THE MASS OF DELTA/N* BY USING THE REACTION KINEMATICS - DM=SQRT((E10+E20)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2) + +clin-9/2012: check argument in sqrt(): + scheck=(E10+E20)**2-p1**2-p2**2-p3**2 + if(scheck.lt.0) then + write(99,*) 'scheck17: ', scheck + write(99,*) 'scheck17', scheck,E10,E20,P(1,I),P(2,I),P(3,I) + write(99,*) 'scheck17-1',E(I1),P(1,I1),P(2,I1),P(3,I1) + write(99,*) 'scheck17-2',E(I2),P(1,I2),P(2,I2),P(3,I2) + scheck=0.d0 + endif + DM=SQRT(sngl(scheck)) +c DM=SQRT((E10+E20)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2) + E(I)=DM RETURN END @@ -10637,6 +10918,8 @@ SUBROUTINE RHORES(I1,I2) PARAMETER (MAXSTR=150001,MAXR=1, 1 AMN=0.939457,AMP=0.93828, 2 AP1=0.13496,AP2=0.13957,AM0=1.232,PI=3.1415926) +clin-9/2012: improve precision for argument in sqrt(): + double precision e10,e20,scheck,p1,p2,p3 COMMON /AA/ R(3,MAXSTR) cc SAVE /AA/ COMMON /BB/ P(3,MAXSTR) @@ -10658,13 +10941,31 @@ SUBROUTINE RHORES(I1,I2) SAVE * 1. DETERMINE THE MOMENTUM COMPONENT OF THE RHO IN THE CMS OF NN FRAME * WE LET I1 TO BE THE RHO AND ABSORB I2 - E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) - E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) +clin-9/2012: improve precision for argument in sqrt(): +c E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) +c E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) + E10=dSQRT(dble(E(I1))**2+dble(P(1,I1))**2 + 1 +dble(P(2,I1))**2+dble(P(3,I1))**2) + E20=dSQRT(dble(E(I2))**2+dble(P(1,I2))**2 + 1 +dble(P(2,I2))**2+dble(P(3,I2))**2) + p1=dble(P(1,I1))+dble(P(1,I2)) + p2=dble(P(2,I1))+dble(P(2,I2)) + p3=dble(P(3,I1))+dble(P(3,I2)) + P(1,I1)=P(1,I1)+P(1,I2) P(2,I1)=P(2,I1)+P(2,I2) P(3,I1)=P(3,I1)+P(3,I2) * 2. DETERMINE THE MASS OF THE RHO BY USING THE REACTION KINEMATICS - DM=SQRT((E10+E20)**2-P(1,I1)**2-P(2,I1)**2-P(3,I1)**2) + +clin-9/2012: check argument in sqrt(): + scheck=(E10+E20)**2-p1**2-p2**2-p3**2 + if(scheck.lt.0) then + write(99,*) 'scheck18: ', scheck + scheck=0.d0 + endif + DM=SQRT(sngl(scheck)) +c DM=SQRT((E10+E20)**2-P(1,I1)**2-P(2,I1)**2-P(3,I1)**2) + E(I1)=DM E(I2)=0 RETURN @@ -10680,6 +10981,8 @@ REAL FUNCTION XNPI(I1,I2,LA,XMAX) PARAMETER (MAXSTR=150001,MAXR=1, 1 AMN=0.939457,AMP=0.93828, 2 AP1=0.13496,AP2=0.13957,AM0=1.232,PI=3.1415926) +clin-9/2012: improve precision for argument in sqrt(): + double precision e10,e20,scheck,p1,p2,p3 COMMON /AA/ R(3,MAXSTR) cc SAVE /AA/ COMMON /BB/ P(3,MAXSTR) @@ -10702,13 +11005,31 @@ REAL FUNCTION XNPI(I1,I2,LA,XMAX) AVMASS=0.5*(AMN+AMP) AVPI=(2.*AP2+AP1)/3. * 1. DETERMINE THE MOMENTUM COMPONENT OF DELTA IN THE LAB. FRAME - E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) - E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) - P1=P(1,I1)+P(1,I2) - P2=P(2,I1)+P(2,I2) - P3=P(3,I1)+P(3,I2) +clin-9/2012: improve precision for argument in sqrt(): +c E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) +c E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) + E10=dSQRT(dble(E(I1))**2+dble(P(1,I1))**2 + 1 +dble(P(2,I1))**2+dble(P(3,I1))**2) + E20=dSQRT(dble(E(I2))**2+dble(P(1,I2))**2 + 1 +dble(P(2,I2))**2+dble(P(3,I2))**2) +c P1=P(1,I1)+P(1,I2) +c P2=P(2,I1)+P(2,I2) +c P3=P(3,I1)+P(3,I2) + p1=dble(P(1,I1))+dble(P(1,I2)) + p2=dble(P(2,I1))+dble(P(2,I2)) + p3=dble(P(3,I1))+dble(P(3,I2)) + * 2. DETERMINE THE MASS OF DELTA BY USING OF THE REACTION KINEMATICS - DM=SQRT((E10+E20)**2-P1**2-P2**2-P3**2) + +clin-9/2012: check argument in sqrt(): + scheck=(E10+E20)**2-p1**2-p2**2-p3**2 + if(scheck.lt.0) then + write(99,*) 'scheck19: ', scheck + scheck=0.d0 + endif + DM=SQRT(sngl(scheck)) +c DM=SQRT((E10+E20)**2-P1**2-P2**2-P3**2) + IF(DM.LE.1.1) THEN XNPI=1.e-09 RETURN @@ -10810,7 +11131,16 @@ REAL FUNCTION SIGMA(SRT,ID,IOI,IOF) ENDIF SS0=AM0**2 Q02=(SS0-(AMU-AMP)**2)*(SS0-(AMU+AMP)**2)/(4.*SS0) - Q0=SQRT(Q02) + +clin-9/2012: check argument in sqrt(): + scheck=Q02 + if(scheck.lt.0) then + write(99,*) 'scheck20: ', scheck + scheck=0. + endif + Q0=SQRT(scheck) +c Q0=SQRT(Q02) + SIGMA=PI*(HC)**2/(2.*P2)*ALFA*(PR/P0)**BETA*AM0**2*T**2 1 *(Q/Q0)**3/((SS-AM0**2)**2+AM0**2*T**2) SIGMA=SIGMA*10. @@ -10886,7 +11216,6 @@ real function ang(srt,iseed) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED c if(srt.le.2.14)then c b1s=0.5 c b2s=0. @@ -11189,6 +11518,8 @@ REAL FUNCTION XN1535(I1,I2,LA) PARAMETER (MAXSTR=150001,MAXR=1, 1 AMN=0.939457,AMP=0.93828,ETAM=0.5475, 2 AP1=0.13496,AP2=0.13957,AM0=1.232,PI=3.1415926) +clin-9/2012: improve precision for argument in sqrt(): + double precision e10,e20,scheck,p1,p2,p3 COMMON /AA/ R(3,MAXSTR) cc SAVE /AA/ COMMON /BB/ P(3,MAXSTR) @@ -11211,13 +11542,31 @@ REAL FUNCTION XN1535(I1,I2,LA) AVMASS=0.5*(AMN+AMP) AVPI=(2.*AP2+AP1)/3. * 1. DETERMINE THE MOMENTUM COMPONENT OF N*(1535) IN THE LAB. FRAME - E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) - E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) - P1=P(1,I1)+P(1,I2) - P2=P(2,I1)+P(2,I2) - P3=P(3,I1)+P(3,I2) +clin-9/2012: improve precision for argument in sqrt(): +c E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) +c E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) + E10=dSQRT(dble(E(I1))**2+dble(P(1,I1))**2 + 1 +dble(P(2,I1))**2+dble(P(3,I1))**2) + E20=dSQRT(dble(E(I2))**2+dble(P(1,I2))**2 + 1 +dble(P(2,I2))**2+dble(P(3,I2))**2) +c P1=P(1,I1)+P(1,I2) +c P2=P(2,I1)+P(2,I2) +c P3=P(3,I1)+P(3,I2) + p1=dble(P(1,I1))+dble(P(1,I2)) + p2=dble(P(2,I1))+dble(P(2,I2)) + p3=dble(P(3,I1))+dble(P(3,I2)) + * 2. DETERMINE THE MASS OF DELTA BY USING OF THE REACTION KINEMATICS - DM=SQRT((E10+E20)**2-P1**2-P2**2-P3**2) + +clin-9/2012: check argument in sqrt(): + scheck=(E10+E20)**2-p1**2-p2**2-p3**2 + if(scheck.lt.0) then + write(99,*) 'scheck21: ', scheck + scheck=0.d0 + endif + DM=SQRT(sngl(scheck)) +c DM=SQRT((E10+E20)**2-P1**2-P2**2-P3**2) + IF(DM.LE.1.1) THEN XN1535=1.E-06 RETURN @@ -11348,7 +11697,16 @@ SUBROUTINE ddp2(SRT,ISEED,PX,PY,PZ,DM1,PNX, bx=-Px/eln by=-Py/eln bz=-Pz/eln - ga=1./sqrt(1.-bx**2-by**2-bz**2) + +clin-9/2012: check argument in sqrt(): + scheck=1.-bx**2-by**2-bz**2 + if(scheck.le.0) then + write(99,*) 'scheck22: ', scheck + stop + endif + ga=1./sqrt(scheck) +c ga=1./sqrt(1.-bx**2-by**2-bz**2) + * the momentum of delta2 and pion in their cms frame elnc=eln/ga pn2=((elnc**2+dm2**2-amp**2)/(2.*elnc))**2-dm2**2 @@ -11369,7 +11727,16 @@ SUBROUTINE ddp2(SRT,ISEED,PX,PY,PZ,DM1,PNX, pny=pnT*sin(fain) SIG=1 IF(X.GT.0)SIG=-1 - pnz=SIG*SQRT(pn**2-PNT**2) + +clin-9/2012: check argument in sqrt(): + scheck=pn**2-PNT**2 + if(scheck.lt.0) then + write(99,*) 'scheck23: ', scheck + scheck=0. + endif + pnz=SIG*SQRT(scheck) +c pnz=SIG*SQRT(pn**2-PNT**2) + en=sqrt(dm2**2+pnx**2+pny**2+pnz**2) * (2) the momentum for the pion ppx=-pnx @@ -11429,7 +11796,16 @@ SUBROUTINE ddrho(SRT,ISEED,PX,PY,PZ,DM1,PNX, * (2.1) estimate the maximum transverse momentum PTMAX2=(SRT**2-(DM1+DM2+AMP)**2)* 1 (SRT**2-(DM1-AMP-DM2)**2)/4./SRT**2 - PTMAX=SQRT(PTMAX2)*1./3. + +clin-9/2012: check argument in sqrt(): + scheck=PTMAX2 + if(scheck.lt.0) then + write(99,*) 'scheck24: ', scheck + scheck=0. + endif + PTMAX=SQRT(scheck)*1./3. +c PTMAX=SQRT(PTMAX2)*1./3. + 7 PT=PTR(PTMAX,ISEED) * (3) GENGRATE THE LONGITUDINAL MOMENTUM FOR DM1 * USING THE GIVEN DISTRIBUTION @@ -11478,7 +11854,16 @@ SUBROUTINE ddrho(SRT,ISEED,PX,PY,PZ,DM1,PNX, bx=-Px/eln by=-Py/eln bz=-Pz/eln - ga=1./sqrt(1.-bx**2-by**2-bz**2) + +clin-9/2012: check argument in sqrt(): + scheck=1.-bx**2-by**2-bz**2 + if(scheck.le.0) then + write(99,*) 'scheck25: ', scheck + stop + endif + ga=1./sqrt(scheck) +c ga=1./sqrt(1.-bx**2-by**2-bz**2) + elnc=eln/ga pn2=((elnc**2+dm2**2-amp**2)/(2.*elnc))**2-dm2**2 if(pn2.le.0)then @@ -11498,7 +11883,16 @@ SUBROUTINE ddrho(SRT,ISEED,PX,PY,PZ,DM1,PNX, pny=pnT*sin(fain) SIG=1 IF(X.GT.0)SIG=-1 - pnz=SIG*SQRT(pn**2-PNT**2) + +clin-9/2012: check argument in sqrt(): + scheck=pn**2-PNT**2 + if(scheck.lt.0) then + write(99,*) 'scheck26: ', scheck + scheck=0. + endif + pnz=SIG*SQRT(scheck) +c pnz=SIG*SQRT(pn**2-PNT**2) + en=sqrt(dm2**2+pnx**2+pny**2+pnz**2) * (2) the momentum for the rho ppx=-pnx @@ -11558,7 +11952,16 @@ SUBROUTINE pprho(SRT,ISEED,PX,PY,PZ,DM1,PNX, * (2.1) estimate the maximum transverse momentum PTMAX2=(SRT**2-(DM1+DM2+AMP)**2)* 1 (SRT**2-(DM1-AMP-DM2)**2)/4./SRT**2 - PTMAX=SQRT(PTMAX2)*1./3. + +clin-9/2012: check argument in sqrt(): + scheck=PTMAX2 + if(scheck.lt.0) then + write(99,*) 'scheck27: ', scheck + scheck=0. + endif + PTMAX=SQRT(scheck)*1./3. +c PTMAX=SQRT(PTMAX2)*1./3. + 7 PT=PTR(PTMAX,ISEED) * (3) GENGRATE THE LONGITUDINAL MOMENTUM FOR DM1 * USING THE GIVEN DISTRIBUTION @@ -11608,7 +12011,16 @@ SUBROUTINE pprho(SRT,ISEED,PX,PY,PZ,DM1,PNX, bx=-Px/eln by=-Py/eln bz=-Pz/eln - ga=1./sqrt(1.-bx**2-by**2-bz**2) + +clin-9/2012: check argument in sqrt(): + scheck=1.-bx**2-by**2-bz**2 + if(scheck.le.0) then + write(99,*) 'scheck28: ', scheck + stop + endif + ga=1./sqrt(scheck) +c ga=1./sqrt(1.-bx**2-by**2-bz**2) + elnc=eln/ga pn2=((elnc**2+dm2**2-amp**2)/(2.*elnc))**2-dm2**2 if(pn2.le.0)then @@ -11628,7 +12040,16 @@ SUBROUTINE pprho(SRT,ISEED,PX,PY,PZ,DM1,PNX, pny=pnT*sin(fain) SIG=1 IF(X.GT.0)SIG=-1 - pnz=SIG*SQRT(pn**2-PNT**2) + +clin-9/2012: check argument in sqrt(): + scheck=pn**2-PNT**2 + if(scheck.lt.0) then + write(99,*) 'scheck29: ', scheck + scheck=0. + endif + pnz=SIG*SQRT(scheck) +c pnz=SIG*SQRT(pn**2-PNT**2) + en=sqrt(dm2**2+pnx**2+pny**2+pnz**2) * (2) the momentum for the rho ppx=-pnx @@ -11682,7 +12103,16 @@ SUBROUTINE ppomga(SRT,ISEED,PX,PY,PZ,DM1,PNX, * (2.1) estimate the maximum transverse momentum PTMAX2=(SRT**2-(DM1+DM2+AMP)**2)* 1 (SRT**2-(DM1-AMP-DM2)**2)/4./SRT**2 - PTMAX=SQRT(PTMAX2)*1./3. + +clin-9/2012: check argument in sqrt(): + scheck=PTMAX2 + if(scheck.lt.0) then + write(99,*) 'scheck30: ', scheck + scheck=0. + endif + PTMAX=SQRT(scheck)*1./3. +c PTMAX=SQRT(PTMAX2)*1./3. + 7 PT=PTR(PTMAX,ISEED) * (3) GENGRATE THE LONGITUDINAL MOMENTUM FOR DM1 * USING THE GIVEN DISTRIBUTION @@ -11731,7 +12161,16 @@ SUBROUTINE ppomga(SRT,ISEED,PX,PY,PZ,DM1,PNX, bx=-Px/eln by=-Py/eln bz=-Pz/eln - ga=1./sqrt(1.-bx**2-by**2-bz**2) + +clin-9/2012: check argument in sqrt(): + scheck=1.-bx**2-by**2-bz**2 + if(scheck.le.0) then + write(99,*) 'scheck31: ', scheck + stop + endif + ga=1./sqrt(scheck) +c ga=1./sqrt(1.-bx**2-by**2-bz**2) + elnc=eln/ga pn2=((elnc**2+dm2**2-amp**2)/(2.*elnc))**2-dm2**2 if(pn2.le.0)then @@ -11751,7 +12190,16 @@ SUBROUTINE ppomga(SRT,ISEED,PX,PY,PZ,DM1,PNX, pny=pnT*sin(fain) SIG=1 IF(X.GT.0)SIG=-1 - pnz=SIG*SQRT(pn**2-PNT**2) + +clin-9/2012: check argument in sqrt(): + scheck=pn**2-PNT**2 + if(scheck.lt.0) then + write(99,*) 'scheck32: ', scheck + scheck=0. + endif + pnz=SIG*SQRT(scheck) +c pnz=SIG*SQRT(pn**2-PNT**2) + en=sqrt(dm2**2+pnx**2+pny**2+pnz**2) * (2) the momentum for the rho ppx=-pnx @@ -11780,7 +12228,6 @@ REAL FUNCTION RMASS(DMAX,ISEED) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED * THE MINIMUM MASS FOR DELTA DMIN = 1.078 * Delta(1232) production @@ -11822,7 +12269,6 @@ REAL FUNCTION RHOMAS(DMAX,ISEED) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED * THE MINIMUM MASS FOR DELTA DMIN = 0.28 * RHO(770) production @@ -12109,7 +12555,16 @@ real function pp2(srt) * 1.Calculate p(lab) from srt [GeV] * Formula used: DSQRT(s) = 2 m DSQRT(E_kin/(2m) + 1) c ekin = 2.*pmass*((srt/(2.*pmass))**2 - 1.) - plab=sqrt(((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2) + +clin-9/2012: check argument in sqrt(): + scheck=((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2 + if(scheck.lt.0) then + write(99,*) 'scheck33: ', scheck + scheck=0. + endif + plab=sqrt(scheck) +c plab=sqrt(((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2) + pmin=2. pmax=2050 if(plab.gt.pmax)then @@ -12146,7 +12601,16 @@ real function ppt(srt) * 1.Calculate p(lab) from srt [GeV] * Formula used: DSQRT(s) = 2 m DSQRT(E_kin/(2m) + 1) c ekin = 2.*pmass*((srt/(2.*pmass))**2 - 1.) - plab=sqrt(((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2) + +clin-9/2012: check argument in sqrt(): + scheck=((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2 + if(scheck.lt.0) then + write(99,*) 'scheck34: ', scheck + scheck=0. + endif + plab=sqrt(scheck) +c plab=sqrt(((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2) + pmin=3. pmax=2100 if ((plab .lt. pmin).or.(plab.gt.pmax)) then @@ -12329,7 +12793,16 @@ real function pplpk(srt) * find the center of mass energy corresponding to the given pm as * if Lambda+N+K are produced pplpk=0. - plab=sqrt(((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2) + +clin-9/2012: check argument in sqrt(): + scheck=((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2 + if(scheck.lt.0) then + write(99,*) 'scheck35: ', scheck + scheck=0. + endif + plab=sqrt(scheck) +c plab=sqrt(((srt**2-2.*pmass**2)/(2.*pmass))**2-pmass**2) + pmin=2.82 pmax=25.0 if(plab.gt.pmax)then @@ -13314,7 +13787,15 @@ SUBROUTINE CRPN(PX,PY,PZ,SRT,I1,I2, cc1=ptr(xptr,iseed) clin-10/25/02-end - c1=sqrt(pr**2-cc1**2)/pr +clin-9/2012: check argument in sqrt(): + scheck=pr**2-cc1**2 + if(scheck.lt.0) then + write(99,*) 'scheck36: ', scheck + scheck=0. + endif + c1=sqrt(scheck)/pr +c c1=sqrt(pr**2-cc1**2)/pr + * C1 = 1.0 - 2.0 * RANART(NSEED) T1 = 2.0 * PI * RANART(NSEED) S1 = SQRT( 1.0 - C1**2 ) @@ -13482,7 +13963,15 @@ SUBROUTINE Crdir(PX,PY,PZ,SRT,I1,I2,IBLOCK) cc1=ptr(xptr,iseed) clin-10/25/02-end - c1=sqrt(pr**2-cc1**2)/pr +clin-9/2012: check argument in sqrt(): + scheck=pr**2-cc1**2 + if(scheck.lt.0) then + write(99,*) 'scheck37: ', scheck + scheck=0. + endif + c1=sqrt(scheck)/pr +c c1=sqrt(pr**2-cc1**2)/pr + T1 = 2.0 * PI * RANART(NSEED) S1 = SQRT( 1.0 - C1**2 ) CT1 = COS(T1) @@ -14146,7 +14635,15 @@ SUBROUTINE CRPD(PX,PY,PZ,SRT,I1,I2, cc1=ptr(xptr,iseed) clin-10/25/02-end - c1=sqrt(pr**2-cc1**2)/pr +clin-9/2012: check argument in sqrt(): + scheck=pr**2-cc1**2 + if(scheck.lt.0) then + write(99,*) 'scheck38: ', scheck + scheck=0. + endif + c1=sqrt(scheck)/pr +c c1=sqrt(pr**2-cc1**2)/pr + c C1 = 1.0 - 2.0 * RANART(NSEED) T1 = 2.0 * PI * RANART(NSEED) S1 = SQRT( 1.0 - C1**2 ) @@ -14998,7 +15495,15 @@ SUBROUTINE CRRD(PX,PY,PZ,SRT,I1,I2, cc1=ptr(xptr,iseed) clin-10/25/02-end - c1=sqrt(pr**2-cc1**2)/pr +clin-9/2012: check argument in sqrt(): + scheck=pr**2-cc1**2 + if(scheck.lt.0) then + write(99,*) 'scheck39: ', scheck + scheck=0. + endif + c1=sqrt(scheck)/pr +c c1=sqrt(pr**2-cc1**2)/pr + T1 = 2.0 * PI * RANART(NSEED) S1 = SQRT( 1.0 - C1**2 ) CT1 = COS(T1) @@ -15045,7 +15550,6 @@ SUBROUTINE Crlaba(PX,PY,PZ,SRT,brel,brsgm, COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - NT=NT c PX0=PX PY0=PY @@ -15361,9 +15865,8 @@ SUBROUTINE crkkpi(I1,I2,XSK1, XSK2, XSK3, XSK4, COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - - XSK11=XSK11 - IBLOCK=1907 + + IBLOCK=1907 X1 = RANART(NSEED) * SIGK XSK2 = XSK1 + XSK2 XSK3 = XSK2 + XSK3 @@ -15492,7 +15995,6 @@ SUBROUTINE Crkhyp(PX,PY,PZ,SRT,I1,I2, cc SAVE /RNDF77/ SAVE - XKY17=XKY17 PX0=PX PY0=PY PZ0=PZ @@ -15811,9 +16313,25 @@ SUBROUTINE Crkpla(PX,PY,PZ,EC,SRT,spika, sig2 = 372.378 if(srt .gt. aphi+aka)then pff = sqrt((srt**2-(aphi+aka)**2)*(srt**2-(aphi-aka)**2)) + +clin-9/2012: check argument in sqrt(): + scheck=pdd + if(scheck.le.0) then + write(99,*) 'scheck40: ', scheck + stop + endif + XKP9 = sig1*pff/sqrt(pdd)*1./32./pi/srt**2 if(srt .gt. aphi+aks)then pff = sqrt((srt**2-(aphi+aks)**2)*(srt**2-(aphi-aks)**2)) + +clin-9/2012: check argument in sqrt(): + scheck=pdd + if(scheck.le.0) then + write(99,*) 'scheck41: ', scheck + stop + endif + XKP10 = sig2*pff/sqrt(pdd)*3./32./pi/srt**2 endif endif @@ -16282,8 +16800,25 @@ SUBROUTINE Crksph(PX,PY,PZ,EC,SRT, c sig11=sig1*(6./dnr)*(srt**2-(aphi+aka)**2)* c & (srt**2-(aphi-aka)**2)/(srt**2-(e(i1)+e(i2))**2)/ c & (srt**2-(e(i1)-e(i2))**2) - pii = sqrt((srt**2-(e(i1)+e(i2))**2)*(srt**2-(e(i1)-e(i2))**2)) - pff = sqrt((srt**2-(aphi+aka)**2)*(srt**2-(aphi-aka)**2)) + +clin-9/2012: check argument in sqrt(): + scheck=(srt**2-(e(i1)+e(i2))**2)*(srt**2-(e(i1)-e(i2))**2) + if(scheck.le.0) then + write(99,*) 'scheck42: ', scheck + stop + endif + pii=sqrt(scheck) +c pii = sqrt((srt**2-(e(i1)+e(i2))**2)*(srt**2-(e(i1)-e(i2))**2)) + +clin-9/2012: check argument in sqrt(): + scheck=(srt**2-(aphi+aka)**2)*(srt**2-(aphi-aka)**2) + if(scheck.lt.0) then + write(99,*) 'scheck43: ', scheck + scheck=0. + endif + pff = sqrt(scheck) +c pff = sqrt((srt**2-(aphi+aka)**2)*(srt**2-(aphi-aka)**2)) + sig11 = sig1*pff/pii*6./dnr/32./pi/srt**2 c if(srt .gt. aphi+aks)then @@ -16440,7 +16975,16 @@ SUBROUTINE bbkaon(ic,SRT,PX,PY,PZ,ana,PlX, bx=-pkx/eln by=-pky/eln bz=-pkz/eln - ga=1./sqrt(1.-bx**2-by**2-bz**2) + +clin-9/2012: check argument in sqrt(): + scheck=1.-bx**2-by**2-bz**2 + if(scheck.le.0) then + write(99,*) 'scheck44: ', scheck + stop + endif + ga=1./sqrt(scheck) +c ga=1./sqrt(1.-bx**2-by**2-bz**2) + elnc=eln/ga pn2=((elnc**2+ana**2-ala**2)/(2.*elnc))**2-ana**2 if(pn2.le.0.)pn2=1.e-09 @@ -16956,7 +17500,16 @@ SUBROUTINE rotate(PX0,PY0,PZ0,px,py,pz) ELSE T2=ATAN2(PY0,PX0) END IF - S2 = SQRT( 1.0 - C2**2 ) + +clin-9/2012: check argument in sqrt(): + scheck=1.0 - C2**2 + if(scheck.lt.0) then + write(99,*) 'scheck45: ', scheck + scheck=0. + endif + S2=sqrt(scheck) +c S2 = SQRT( 1.0 - C2**2 ) + CT2 = COS(T2) ST2 = SIN(T2) * the momentum, polar and azimuthal angles of the momentum to be rotated @@ -16968,7 +17521,16 @@ SUBROUTINE rotate(PX0,PY0,PZ0,px,py,pz) ELSE T1=ATAN2(PY,PX) ENDIF - S1 = SQRT( 1.0 - C1**2 ) + +clin-9/2012: check argument in sqrt(): + scheck=1.0 - C1**2 + if(scheck.lt.0) then + write(99,*) 'scheck46: ', scheck + scheck=0. + endif + S1=sqrt(scheck) +c S1 = SQRT( 1.0 - C1**2 ) + CT1 = COS(T1) ST1 = SIN(T1) SS = C2 * S1 * CT1 + S2 * C1 @@ -16998,7 +17560,6 @@ real function Xpp(srt) &31.6,25.9,24.0,23.1, &24.0,28.3,33.6,41.5,47/ - xpp=0. pmass=0.9383 * 1.Calculate E_kin(lab) [MeV] from srt [GeV] * Formula used: DSQRT(s) = 2 m DSQRT(E_kin/(2m) + 1) @@ -17053,7 +17614,6 @@ real function Xnp(srt) data xarray / 410.,270.,214.5,130.,78.,53.5, &41.6,35.9,34.2,34.3,34.9/ - xnp=0. pmass=0.9383 * 1.Calculate E_kin(lab) [MeV] from srt [GeV] * Formula used: DSQRT(s) = 2 m DSQRT(E_kin/(2m) + 1) @@ -17100,7 +17660,6 @@ function ptr(ptmax,iseed) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED ptr=0. if(ptmax.le.1.e-02)then ptr=ptmax @@ -17687,7 +18246,6 @@ real function dirct1(srt) &6.420560,6.431045,6.441367,6.451529,6.461533,6.471386,6.481091, &6.490650,6.476413,6.297259,6.097826/ - dirct1=0 if (srt .lt. earray(1)) then dirct1 = 0.00001 return @@ -17770,7 +18328,6 @@ real function dirct2(srt) &5.552088,5.593359,5.634520,5.675570,5.716515,5.757356,5.798093, &5.838732,5.879272,5.919717,5.960068,5.980941/ - dirct2=0. if (srt .lt. earray(1)) then dirct2 = 0.00001 return @@ -18116,7 +18673,6 @@ subroutine Rmasdd(srt,am10,am20, COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED amn=0.94 amp=0.14 * the maximum mass for resonance 1 @@ -18241,7 +18797,6 @@ subroutine flow(nt) cc SAVE /input1/ SAVE *----------------------------------------------------------------------* - NT=NT ycut1=-2.6 ycut2=2.6 DY=0.2 @@ -18399,7 +18954,6 @@ subroutine pbarfs(srt,npion,iseed) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED C the factorial coefficients in the pion no. distribution * from n=2 to 6 calculated use the formula in the reference factor(2)=1. @@ -18590,7 +19144,16 @@ SUBROUTINE XKKANN(SRT, XSK1, XSK2, XSK3, XSK4, XSK5, c* K+ + K- --> phi fwdp = 1.68*(aphi**2-4.*aka**2)**1.5/6./aphi/aphi - pkaon=0.5*sqrt(srt**2-4.0*aka**2) + +clin-9/2012: check argument in sqrt(): + scheck=srt**2-4.0*aka**2 + if(scheck.le.0) then + write(99,*) 'scheck47: ', scheck + stop + endif + pkaon=0.5*sqrt(scheck) +c pkaon=0.5*sqrt(srt**2-4.0*aka**2) + XSK11 = 30.*3.14159*0.1973**2*(aphi*fwdp)**2/ & ((srt**2-aphi**2)**2+(aphi*fwdp)**2)/pkaon**2 c @@ -18964,7 +19527,15 @@ SUBROUTINE PHIMES(I1, I2, SRT, XSK1, XSK2, XSK3, XSK4, XSK5, c !! mb, elastic XSK1 = 5.0 - pii = sqrt((S-(em1+em2)**2)*(S-(em1-em2)**2)) +clin-9/2012: check argument in sqrt(): + scheck=(S-(em1+em2)**2)*(S-(em1-em2)**2) + if(scheck.le.0) then + write(99,*) 'scheck48: ', scheck + stop + endif + pii=sqrt(scheck) +c pii = sqrt((S-(em1+em2)**2)*(S-(em1-em2)**2)) + * phi + K(-bar) channel if( lb1.eq.23.or.lb2.eq.23 .or. lb1.eq.21.or.lb2.eq.21 )then if(srt .gt. (ap1+akap))then @@ -19676,7 +20247,7 @@ SUBROUTINE bbarfs(lbb1,lbb2,ei1,ei2,iblock,iseed) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED + c determine which final BbarB channel occurs: rd=RANART(NSEED) wsum=0. @@ -20021,7 +20592,6 @@ real function ptor(srt) * for rho rho -> pi pi, assumed a constant cross section (in mb) real function rtop(srt) ***************************************** - srt=srt rtop=5. return END @@ -20040,7 +20610,7 @@ SUBROUTINE pi2ro2(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - iseed=iseed + if((lb(i1).ge.3.and.lb(i1).le.5) 1 .and.(lb(i2).ge.3.and.lb(i2).le.5)) then iblock=1850 @@ -20104,7 +20674,7 @@ real function ptoe(srt) * for eta eta -> pi pi, assumed a constant cross section (in mb) real function etop(srt) ***************************************** - srt=srt + c eta equilibration: c most important channel is found to be pi pi <-> pi eta, then c rho pi <-> rho eta. @@ -20127,7 +20697,6 @@ SUBROUTINE pi2et2(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed) cc SAVE /RNDF77/ SAVE - iseed=iseed if((lb(i1).ge.3.and.lb(i1).le.5) 1 .and.(lb(i2).ge.3.and.lb(i2).le.5)) then iblock=1860 @@ -20194,7 +20763,7 @@ real function pptope(srt) * for pi eta -> pi pi, assumed a constant cross section (in mb) real function petopp(srt) ***************************************** - srt=srt + c eta equilibration: petopp=5. return @@ -20215,7 +20784,6 @@ SUBROUTINE pi3eta(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed) cc SAVE /RNDF77/ SAVE - ISEED=ISEED if((lb(i1).ge.3.and.lb(i1).le.5) 1 .and.(lb(i2).ge.3.and.lb(i2).le.5)) then iblock=1870 @@ -20286,7 +20854,7 @@ real function rptore(srt) * for rho eta -> rho pi, assumed a constant cross section (in mb) real function retorp(srt) ***************************************** - srt=srt + c eta equilibration: retorp=5. return @@ -20306,7 +20874,7 @@ SUBROUTINE rpiret(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ SAVE - ISEED=ISEED + if((lb(i1).ge.25.and.lb(i1).le.27 1 .and.lb(i2).ge.3.and.lb(i2).le.5).or. 2 (lb(i1).ge.3.and.lb(i1).le.5 @@ -20376,7 +20944,7 @@ real function xop2oe(srt) * for omega eta -> omega pi, assumed a constant cross section (in mb) real function xoe2op(srt) ***************************************** - srt=srt + c eta equilibration: xoe2op=5. return @@ -20397,7 +20965,6 @@ SUBROUTINE opioet(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed) cc SAVE /RNDF77/ SAVE - iseed=iseed if((lb(i1).ge.3.and.lb(i1).le.5.and.lb(i2).eq.28).or. 1 (lb(i2).ge.3.and.lb(i2).le.5.and.lb(i1).eq.28)) then iblock=1890 @@ -20462,7 +21029,7 @@ real function eetorr(srt) * for rho rho -> eta eta, assumed a constant cross section (in mb) real function rrtoee(srt) ***************************************** - srt=srt + c eta equilibration: rrtoee=5. return @@ -20483,7 +21050,6 @@ SUBROUTINE ro2et2(i1,i2,lbb1,lbb2,ei1,ei2,iblock,iseed) cc SAVE /RNDF77/ SAVE - ISEED=ISEED if(lb(i1).ge.25.and.lb(i1).le.27.and. 1 lb(i2).ge.25.and.lb(i2).le.27) then iblock=1895 @@ -20640,6 +21206,8 @@ SUBROUTINE KSRESO(I1,I2) PARAMETER (MAXSTR=150001,MAXR=1, 1 AMN=0.939457,AMP=0.93828, 2 AP1=0.13496,AP2=0.13957,AM0=1.232,PI=3.1415926) +clin-9/2012: improve precision for argument in sqrt(): + double precision e10,e20,scheck,p1,p2,p3 COMMON /AA/ R(3,MAXSTR) cc SAVE /AA/ COMMON /BB/ P(3,MAXSTR) @@ -20661,8 +21229,18 @@ SUBROUTINE KSRESO(I1,I2) SAVE * 1. DETERMINE THE MOMENTUM COMPONENT OF THE K* IN THE CMS OF PI-K FRAME * WE LET I1 TO BE THE K* AND ABSORB I2 - E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) - E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) + +clin-9/2012: improve precision for argument in sqrt(): +c E10=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2) +c E20=SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2) + E10=dSQRT(dble(E(I1))**2+dble(P(1,I1))**2 + 1 +dble(P(2,I1))**2+dble(P(3,I1))**2) + E20=dSQRT(dble(E(I2))**2+dble(P(1,I2))**2 + 1 +dble(P(2,I2))**2+dble(P(3,I2))**2) + p1=dble(P(1,I1))+dble(P(1,I2)) + p2=dble(P(2,I1))+dble(P(2,I2)) + p3=dble(P(3,I1))+dble(P(3,I2)) + IF(LB(I2) .EQ. 21 .OR. LB(I2) .EQ. 23) THEN E(I1)=0. I=I2 @@ -20679,7 +21257,18 @@ SUBROUTINE KSRESO(I1,I2) P(2,I)=P(2,I1)+P(2,I2) P(3,I)=P(3,I1)+P(3,I2) * 2. DETERMINE THE MASS OF K* BY USING THE REACTION KINEMATICS - DM=SQRT((E10+E20)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2) + +clin-9/2012: check argument in sqrt(): + scheck=(E10+E20)**2-p1**2-p2**2-p3**2 + if(scheck.lt.0) then + write(99,*) 'scheck49: ',scheck + write(99,*) 'scheck49',scheck,E10,E20,P(1,I),P(2,I),P(3,I) + write(99,*) 'scheck49-1',E(I1),P(1,I1),P(2,I1),P(3,I1) + write(99,*) 'scheck49-2',E(I2),P(1,I2),P(2,I2),P(3,I2) + endif + DM=sqrt(sngl(scheck)) +c DM=SQRT((E10+E20)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2) + E(I)=DM RETURN END @@ -20749,14 +21338,11 @@ SUBROUTINE pertur(PX,PY,PZ,SRT,IRUN,I1,I2,nt,kp,icont) 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) SAVE - kp=kp - nt=nt px0 = px py0 = py pz0 = pz LB1 = LB(I1) - lb1abs = iabs(lb1) EM1 = E(I1) X1 = R(1,I1) Y1 = R(2,I1) @@ -20807,7 +21393,7 @@ SUBROUTINE pertur(PX,PY,PZ,SRT,IRUN,I1,I2,nt,kp,icont) c---------------------------------------------------- * for process: K-bar + L(S) --> Ca + pi * -60 if(lb1abs.ge.14 .and. lb1abs.le.17)then +60 if(iabs(lb1).ge.14 .and. iabs(lb1).le.17)then asap = e(i1) akap = e(i2) idp = i1 @@ -20892,7 +21478,7 @@ SUBROUTINE pertur(PX,PY,PZ,SRT,IRUN,I1,I2,nt,kp,icont) c---------------------------------------------------- * for process: Cas(bar) + K_bar(K) --> Om(bar) + pi !! eta * -70 if(lb1abs.eq.40 .or. lb1abs.eq.41)then +70 if(iabs(lb1).eq.40 .or. iabs(lb1).eq.41)then acap = e(i1) akap = e(i2) idp = i1 @@ -20967,7 +21553,7 @@ SUBROUTINE pertur(PX,PY,PZ,SRT,IRUN,I1,I2,nt,kp,icont) c----------------------------------------------------------- * for process: Ca + pi/eta --> K-bar + L(S) * -90 if(lb1abs.eq.40 .or. lb1abs.eq.41)then +90 if(iabs(lb1).eq.40 .or. iabs(lb1).eq.41)then acap = e(i1) app = e(i2) idp = i1 @@ -21383,7 +21969,6 @@ SUBROUTINE distc0(drmax,deltr0,DT, COMMON /BG/ BETAX,BETAY,BETAZ,GAMMA cc SAVE /BG/ SAVE - deltr0=deltr0 Ifirst=-1 E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2) *NOW PARTICLES ARE CLOSE ENOUGH TO EACH OTHER ! @@ -21458,7 +22043,16 @@ subroutine sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal) c endif c s=srt**2 - pinitial=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt + +clin-9/2012: check argument in sqrt(): + scheck=(s-(em1+em2)**2)*(s-(em1-em2)**2) + if(scheck.le.0) then + write(99,*) 'scheck50: ', scheck + stop + endif + pinitial=sqrt(scheck)/2./srt +c pinitial=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt + fs=fnndpi(s) c Determine isospin and spin factors for the ratio between c BB->Deuteron+Meson and Deuteron+Meson->BB cross sections: @@ -21538,23 +22132,23 @@ subroutine sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal) else c for other BB initial states (spin- and isospin averaged): if(idxsec.eq.1) then -c 1: assume the same |matrix element|**2 (after averaging over initial +c 1: assume the same |matrix element|**2/s (after averaging over initial c spins and isospins) for B+B -> deuteron+meson at the same sqrt(s); sbbdpi=fs*pfinal/pinitial*3./16. elseif(idxsec.eq.2.or.idxsec.eq.4) then threshold=amax1(xmd+xmm,em1+em2) snew=(srt-threshold+srt0)**2 if(idxsec.eq.2) then -c 2: assume the same |matrix element|**2 for B+B -> deuteron+meson +c 2: assume the same |matrix element|**2/s for B+B -> deuteron+meson c at the same sqrt(s)-threshold: sbbdpi=fnndpi(snew)*pfinal/pinitial*3./16. elseif(idxsec.eq.4) then -c 4: assume the same |matrix element|**2 for B+B <- deuteron+meson +c 4: assume the same |matrix element|**2/s for B+B <- deuteron+meson c at the same sqrt(s)-threshold: sbbdpi=fnndpi(snew)*pfinal/pinitial/6.*pifactor endif elseif(idxsec.eq.3) then -c 3: assume the same |matrix element|**2 for B+B <- deuteron+meson +c 3: assume the same |matrix element|**2/s for B+B <- deuteron+meson c at the same sqrt(s): sbbdpi=fs*pfinal/pinitial/6.*pifactor endif @@ -21764,11 +22358,11 @@ subroutine sdmbb(SRT,sdm,ianti) endif clin-9/2008 For elastic collisions: if(idxsec.eq.1.or.idxsec.eq.3) then -c 1/3: assume the same |matrix element|**2 (after averaging over initial +c 1/3: assume the same |matrix element|**2/s (after averaging over initial c spins and isospins) for d+Meson elastic at the same sqrt(s); sdmel=fdpiel(s) elseif(idxsec.eq.2.or.idxsec.eq.4) then -c 2/4: assume the same |matrix element|**2 (after averaging over initial +c 2/4: assume the same |matrix element|**2/s (after averaging over initial c spins and isospins) for d+Meson elastic at the same sqrt(s)-threshold: threshold=em1+em2 snew=(srt-threshold+srt0)**2 @@ -21803,23 +22397,23 @@ subroutine sdmbb(SRT,sdm,ianti) if(srt.gt.(xmnn1+xmnn2)) then pfinal=sqrt((s-(xmnn1+xmnn2)**2)*(s-(xmnn1-xmnn2)**2))/2./srt if(idxsec.eq.1) then -c 1: assume the same |matrix element|**2 (after averaging over initial +c 1: assume the same |matrix element|**2/s (after averaging over initial c spins and isospins) for B+B -> deuteron+meson at the same sqrt(s); sdmnn=fs*pfinal/pinitial*3./16.*xnnfactor elseif(idxsec.eq.2.or.idxsec.eq.4) then threshold=amax1(xmnn1+xmnn2,em1+em2) snew=(srt-threshold+srt0)**2 if(idxsec.eq.2) then -c 2: assume the same |matrix element|**2 for B+B -> deuteron+meson +c 2: assume the same |matrix element|**2/s for B+B -> deuteron+meson c at the same sqrt(s)-threshold: sdmnn=fnndpi(snew)*pfinal/pinitial*3./16.*xnnfactor elseif(idxsec.eq.4) then -c 4: assume the same |matrix element|**2 for B+B <- deuteron+meson +c 4: assume the same |matrix element|**2/s for B+B <- deuteron+meson c at the same sqrt(s)-threshold: sdmnn=fnndpi(snew)*pfinal/pinitial/6. endif elseif(idxsec.eq.3) then -c 3: assume the same |matrix element|**2 for B+B <- deuteron+meson +c 3: assume the same |matrix element|**2/s for B+B <- deuteron+meson c at the same sqrt(s): sdmnn=fs*pfinal/pinitial/6. endif @@ -22125,7 +22719,16 @@ SUBROUTINE crdmbb(PX,PY,PZ,SRT,I1,I2,IBLOCK, write(91,*) ' d+',lbm,' (pert dbar M elastic) @nt=',nt 1 ,' @prob=',dpertp(ideut) endif - pfinal=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt + +clin-9/2012: check argument in sqrt(): + scheck=(s-(em1+em2)**2)*(s-(em1-em2)**2) + if(scheck.lt.0) then + write(99,*) 'scheck51: ', scheck + scheck=0. + endif + pfinal=sqrt(scheck)/2./srt +c pfinal=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt + CALL dmelangle(pxn,pyn,pzn,pfinal) CALL ROTATE(PX,PY,PZ,Pxn,Pyn,Pzn) EdCM=SQRT(E(ideut)**2+Pxn**2+Pyn**2+Pzn**2) @@ -22234,8 +22837,16 @@ SUBROUTINE crdmbb(PX,PY,PZ,SRT,I1,I2,IBLOCK, E(I2)=xmb2 lb1=lb(i1) lb2=lb(i2) - pfinal=sqrt((s-(xmb1+xmb2)**2)*(s-(xmb1-xmb2)**2))/2./srt -c + +clin-9/2012: check argument in sqrt(): + scheck=(s-(xmb1+xmb2)**2)*(s-(xmb1-xmb2)**2) + if(scheck.lt.0) then + write(99,*) 'scheck52: ', scheck + scheck=0. + endif + pfinal=sqrt(scheck)/2./srt +c pfinal=sqrt((s-(xmb1+xmb2)**2)*(s-(xmb1-xmb2)**2))/2./srt + if(iblock.eq.502) then CALL dmangle(pxn,pyn,pzn,nt,ianti,pfinal,lbm) elseif(iblock.eq.504) then @@ -22353,11 +22964,11 @@ subroutine sdbelastic(SRT,sdb) s=srt**2 c For elastic collisions: if(idxsec.eq.1.or.idxsec.eq.3) then -c 1/3: assume the same |matrix element|**2 (after averaging over initial +c 1/3: assume the same |matrix element|**2/s (after averaging over initial c spins and isospins) for d+Baryon elastic at the same sqrt(s); sdbel=fdbel(s) elseif(idxsec.eq.2.or.idxsec.eq.4) then -c 2/4: assume the same |matrix element|**2 (after averaging over initial +c 2/4: assume the same |matrix element|**2/s (after averaging over initial c spins and isospins) for d+Baryon elastic at the same sqrt(s)-threshold: threshold=em1+em2 snew=(srt-threshold+srt0)**2 @@ -22413,7 +23024,16 @@ SUBROUTINE crdbel(PX,PY,PZ,SRT,I1,I2,IBLOCK, 1 ,' @prob=',dpertp(ideut),p(1,idb),p(2,idb) 2 ,p(1,ideut),p(2,ideut) endif - pfinal=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt + +clin-9/2012: check argument in sqrt(): + scheck=(s-(em1+em2)**2)*(s-(em1-em2)**2) + if(scheck.lt.0) then + write(99,*) 'scheck53: ', scheck + scheck=0. + endif + pfinal=sqrt(scheck)/2./srt +c pfinal=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt + CALL dbelangle(pxn,pyn,pzn,pfinal) CALL ROTATE(PX,PY,PZ,Pxn,Pyn,Pzn) EdCM=SQRT(E(ideut)**2+Pxn**2+Pyn**2+Pzn**2) @@ -22446,7 +23066,15 @@ SUBROUTINE crdbel(PX,PY,PZ,SRT,I1,I2,IBLOCK, write (91,*) ' d+',lbb,' (regular dbar Bbar elastic) @evt#', 1 iaevt,' @nt=',nt,' lb1,2=',lb1,lb2 endif - pfinal=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt +clin-9/2012: check argument in sqrt(): + scheck=(s-(em1+em2)**2)*(s-(em1-em2)**2) + if(scheck.lt.0) then + write(99,*) 'scheck54: ', scheck + scheck=0. + endif + pfinal=sqrt(scheck)/2./srt +c pfinal=sqrt((s-(em1+em2)**2)*(s-(em1-em2)**2))/2./srt + CALL dbelangle(pxn,pyn,pzn,pfinal) * ROTATE THE MOMENTA OF PARTICLES IN THE CMS OF P1+P2 c (This is not needed for isotropic distributions) diff --git a/GeneratorInterface/AMPTInterface/src/hijing1.383_ampt.f b/GeneratorInterface/AMPTInterface/src/hijing1.383_ampt.f index 9b467f0168acb..73a08ed1a3505 100755 --- a/GeneratorInterface/AMPTInterface/src/hijing1.383_ampt.f +++ b/GeneratorInterface/AMPTInterface/src/hijing1.383_ampt.f @@ -324,6 +324,9 @@ SUBROUTINE HIJING(FRAME,BMIN0,BMAX0) common /phiHJ/iphirp,phiRP SAVE +ctuos-Need this line for correct b range-06/2015 + DOUBLE PRECISION BMIN0,BMAX0 + BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35)) BMIN=MIN(BMIN0,BMAX) IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN @@ -2061,7 +2064,7 @@ SUBROUTINE HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT) c (to limit strangeness enhancement when string tension is strongly c increased due to using a very low value of parameter b in Lund c symmetric splitting function as done in arXiv:1403.6321): -c PARJ(2)=min(PARJ(2),0.4) + PARJ(2)=min(PARJ(2),0.4) C ******** set up for jetset IF(FRAME.EQ.'LAB') THEN diff --git a/GeneratorInterface/AMPTInterface/src/hipyset1.35.f b/GeneratorInterface/AMPTInterface/src/hipyset1.35.f index 016edf4eda61f..87c0e4afdd07c 100755 --- a/GeneratorInterface/AMPTInterface/src/hipyset1.35.f +++ b/GeneratorInterface/AMPTInterface/src/hipyset1.35.f @@ -15,9 +15,6 @@ C C********************************************************************* -cms -cms gsfs 8/2009 Renamed common block PYINT4 due to conflict with something in CMSSW -cms SUBROUTINE LU2ENT(IP,KF1,KF2,PECM) C...Purpose: to store two partons/particles in their CM frame, @@ -500,9 +497,6 @@ SUBROUTINE LUPREP(IP) SAVE /LUDAT3/ DIMENSION DPS(5),DPC(5),UE(3) - ic1=0 - ic2=0 - kci=0 C...Rearrange parton shower product listing along strings: begin loop. I1=N DO 130 MQGST=1,2 @@ -668,7 +662,9 @@ SUBROUTINE LUPREP(IP) IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) GOTO 260 C...Perform two-particle decay of jet system, if possible. - IF(PECM.GE.0.02d0*DPC(4)) THEN +clin-5/2012: +c IF(PECM.GE.0.02d0*DPC(4)) THEN + IF(dble(PECM).GE.0.02d0*DPC(4)) THEN PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2- & (P(N+2,5)-P(N+3,5))**2))/(2.*PECM) UE(3)=2.*RLU(0)-1. @@ -855,11 +851,6 @@ SUBROUTINE LUSTRF(IP) FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3) DFOUR(I,J)=DP(I,4)*DP(J,4)-DP(I,1)*DP(J,1)-DP(I,2)*DP(J,2)- &DP(I,3)*DP(J,3) - - ir=0 - in3=0 - jr=0 - prev=0 C...Reset counters. Identify parton system. MSTJ(91)=0 @@ -1190,7 +1181,9 @@ SUBROUTINE LUSTRF(IP) DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2) DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2) DHC12=DFOUR(1,2) - IF(DHC12.LE.1E-2) THEN +clin-5/2012: +c IF(DHC12.LE.1E-2) THEN + IF(DHC12.LE.1D-2) THEN P(IN(1)+2,4)=P(IN(1)+2,3) P(IN(1)+2,1)=0. IN(1)=IN(1)+4 @@ -1250,7 +1243,9 @@ SUBROUTINE LUSTRF(IP) C...Junction strings: solve (m2, Gamma) equation system for energies. DHS1=DHM(3)*DHG(4)-DHM(4)*DHG(3) - IF(ABS(DHS1).LT.1E-4) GOTO 270 +clin-5/2012: +c IF(ABS(DHS1).LT.1E-4) GOTO 270 + IF(DABS(DHS1).LT.1D-4) GOTO 270 DHS2=DHM(4)*(dble(GAM(3))-DHG(1))-DHM(2)*DHG(3)-DHG(4)* &(dble(P(I,5))**2-DHM(1))+DHG(2)*DHM(3) DHS3=DHM(2)*(dble(GAM(3))-DHG(1)) @@ -1593,7 +1588,9 @@ SUBROUTINE LUSTRF(IP) DP(1,4)=DSQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2) DP(2,4)=DSQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2) DHC12=DFOUR(1,2) - IF(DHC12.LE.1E-2) THEN +clin-5/2012: +c IF(DHC12.LE.1E-2) THEN + IF(DHC12.LE.1D-2) THEN P(IN(JT)+2,4)=P(IN(JT)+2,3) P(IN(JT)+2,JT)=0. IN(JT)=IN(JT)+4*JS @@ -1655,7 +1652,9 @@ SUBROUTINE LUSTRF(IP) C...Solve (m2, Gamma) equation system for energies taken. DHS1=DHM(JR+1)*DHG(4)-DHM(4)*DHG(JR+1) - IF(ABS(DHS1).LT.1E-4) GOTO 550 +clin-5/2012: +c IF(ABS(DHS1).LT.1E-4) GOTO 550 + IF(DABS(DHS1).LT.1D-4) GOTO 550 DHS2=DHM(4)*(dble(GAM(3))-DHG(1))-DHM(JT+1)*DHG(JR+1)-DHG(4)* &(dble(P(I,5))**2-DHM(1))+DHG(JT+1)*DHM(JR+1) DHS3=DHM(JT+1)*(dble(GAM(3))-DHG(1))-DHG(JT+1) @@ -1836,7 +1835,6 @@ SUBROUTINE LUINDF(IP) DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3), &KFLO(2),PXO(2),PYO(2),WO(2) - pw=0. C...Reset counters. Identify parton system and take copy. Check flavour. NSAV=N NJET=0 @@ -2110,8 +2108,7 @@ SUBROUTINE LUINDF(IP) IF(NREM.GT.0) GOTO 280 C...Compensate for missing momentum in global scheme (3 options). - 320 MTMP = MOD(MSTJ(3),5) - IF(MTMP.NE.0.AND.MTMP.NE.4) THEN + 320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN DO 330 J=1,3 PSI(J)=0. DO 330 I=NSAV+NJET+1,N @@ -2264,19 +2261,6 @@ SUBROUTINE LUDECY(IP) &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA) C...Initial values. - - idc=0 - pqt=0. - hatu=0. - hmp1=0. - im=0 - kfam=0 - wtmax=0. - pmes=0. - pmst=0. - wt=0. - pmr=0. - NTRY=0 NSAV=N KFA=IABS(K(IP,2)) @@ -2953,16 +2937,6 @@ SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF) COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) SAVE /LUDAT2/ - par3m=0. - par4m=0. - pardm=0. - pars0=0. - pars1=0. - pars2=0. - parsm=0. - kmul=0 - ktab3=0 - C...Default flavour values. Input consistency checks. KF1A=IABS(KFL1) KF2A=IABS(KFL2) @@ -3304,10 +3278,6 @@ SUBROUTINE LUZDIS(KFL1,KFL2,PR,Z) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) SAVE /LUDAT1/ - zdiv=0. - fint=0. - zdivc=0. - C...Check if heavy flavour fragmentation. KFLA=IABS(KFL1) KFLB=IABS(KFL2) @@ -3421,24 +3391,6 @@ SUBROUTINE LUSHOW(IP1,IP2,QMAX) SAVE /LUDAT2/ DIMENSION PMTH(5,40),PS(5),PMA(4),PMSD(4),IEP(4),IPA(4), &KFLA(4),KFLD(4),KFL(4),ITRY(4),ISI(4),ISL(4),DP(4),DPT(5,4) - - npa=0 - kflm=0 - pem=0. - pmed=0. - fbre=0. - pm2=0. - ped=0. - zm=0. - pa1s=0. - pa2s=0. - pa3s=0. - pts=0. - pzm=0. - pmls=0. - pt=0. - hazip=0. - hazic=0. C...Initialization of cutoff masses etc. IF(MSTJ(41).LE.0.OR.(MSTJ(41).EQ.1.AND.QMAX.LE.PARJ(82)).OR. @@ -4143,7 +4095,9 @@ SUBROUTINE LUSHOW(IP1,IP2,QMAX) 370 DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/DPMM DPT(4,4)=DSQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2) DPT(5,4)=DSQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2) - IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1*PARJ(82)) THEN +clin-5/2012: +c IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1*PARJ(82)) THEN + IF(sngl(MIN(DPT(4,4),DPT(5,4))).GT.(0.1*PARJ(82))) THEN CAD=sngl((DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+ & DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))) IF(MAZIP.NE.0) THEN @@ -4277,12 +4231,6 @@ SUBROUTINE LUBOEI(NSAV) DIMENSION DPS(4),KFBE(9),NBE(0:9),BEI(100) DATA KFBE/211,-211,111,321,-321,130,310,221,331/ - pmhq=0. - qdel=0. - nbin=0 - beex=0. - bert=0. - C...Boost event to overall CM frame. Calculate CM energy. IF((MSTJ(51).NE.1.AND.MSTJ(51).NE.2).OR.N-NSAV.LE.1) RETURN DO 100 J=1,4 @@ -4419,8 +4367,6 @@ FUNCTION ULMASS(KF) SAVE /LUDAT1/ COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) SAVE /LUDAT2/ - - pmspl=0. C...Reset variables. Compressed code. ULMASS=0. @@ -4929,7 +4875,6 @@ FUNCTION RLU(IDUM) &(RRLU98,RRLU(98)),(RRLU99,RRLU(99)),(RRLU00,RRLU(100)) C...Initialize generation from given seed. - IDUM=IDUM IF(MRLU2.EQ.0) THEN IJ=MOD(MRLU1/30082,31329) KL=MOD(MRLU1,30082) @@ -5025,7 +4970,9 @@ SUBROUTINE LUROBO(THE,PHI,BEX,BEY,BEZ) ENDIF C...Rotate, typically from z axis to direction (theta,phi). - IF(THE**2+PHI**2.GT.1E-20) THEN +clin-5/2012: +c IF(THE**2+PHI**2.GT.1E-20) THEN + IF((THE**2+PHI**2).GT.1E-20) THEN ROT(1,1)=COS(THE)*COS(PHI) ROT(1,2)=-SIN(PHI) ROT(1,3)=SIN(THE)*COS(PHI) @@ -5047,7 +4994,9 @@ SUBROUTINE LUROBO(THE,PHI,BEX,BEY,BEZ) ENDIF C...Boost, typically from rest to momentum/energy=beta. - IF(DBX**2+DBY**2+DBZ**2.GT.1E-20) THEN +clin-5/2012: +c IF(DBX**2+DBY**2+DBZ**2.GT.1E-20) THEN + IF((DBX**2+DBY**2+DBZ**2).GT.1D-20) THEN DB=SQRT(DBX**2+DBY**2+DBZ**2) IF(DB.GT.0.99999999D0) THEN C...Rescale boost vector if too close to unity. @@ -5094,7 +5043,7 @@ SUBROUTINE HIROBO(THE,PHI,BEX,BEY,BEZ) SAVE /LUJETS/ COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) SAVE /LUDAT1/ - DIMENSION ROT(3,3),PR(3),DP(4) + DIMENSION ROT(3,3),PR(3),VR(3),DP(4),DV(4) C...Find range of rotation/boost. Convert boost to double precision. IMIN=1 @@ -5112,7 +5061,9 @@ SUBROUTINE HIROBO(THE,PHI,BEX,BEY,BEZ) ENDIF C...Rotate, typically from z axis to direction (theta,phi). - IF(THE**2+PHI**2.GT.1E-20) THEN +clin-5/2012: +c IF(THE**2+PHI**2.GT.1E-20) THEN + IF((THE**2+PHI**2).GT.1E-20) THEN ROT(1,1)=COS(THE)*COS(PHI) ROT(1,2)=-SIN(PHI) ROT(1,3)=SIN(THE)*COS(PHI) @@ -5132,7 +5083,9 @@ SUBROUTINE HIROBO(THE,PHI,BEX,BEY,BEZ) ENDIF C...Boost, typically from rest to momentum/energy=beta. - IF(DBX**2+DBY**2+DBZ**2.GT.1E-20) THEN +clin-5/2012: +c IF(DBX**2+DBY**2+DBZ**2.GT.1E-20) THEN + IF((DBX**2+DBY**2+DBZ**2).GT.1D-20) THEN DB=SQRT(DBX**2+DBY**2+DBZ**2) IF(DB.GT.0.99999999D0) THEN C...Rescale boost vector if too close to unity. @@ -5382,7 +5335,6 @@ SUBROUTINE LULIST(MLIST) DATA CHMO/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep', &'Oct','Nov','Dec'/,CHDL/'(())',' ','()','!!','<>','==','(==)'/ - CHMO(1)=CHMO(1) C...Initialization printout: version number and date of last change. C IF(MLIST.EQ.0.OR.MSTU(12).EQ.1) THEN C WRITE(MSTU(11),1000) MSTU(181),MSTU(182),MSTU(185), @@ -6117,8 +6069,10 @@ SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN) CHARACTER CHFRAM*8,CHBEAM*8,CHTARG*8,CHMO(12)*3,CHLH(2)*6 DATA CHMO/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep', &'Oct','Nov','Dec'/, CHLH/'lepton','hadron'/ + +clin-12/2012 correct NN differential cross section in HIJING: + WRITE(MSTU(11),*) 'In PYINIT: BEAM,TARGET= ',BEAM,TARGET - CHMO(1)=CHMO(1) C...Write headers. C IF(MSTP(122).GE.1) WRITE(MSTU(11),1000) MSTP(181),MSTP(182), C &MSTP(185),CHMO(MSTP(184)),MSTP(183) @@ -6822,15 +6776,13 @@ SUBROUTINE PYINRE SAVE /PYINT1/ COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2) SAVE /PYINT2/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ COMMON/PYINT6/PROC(0:200) CHARACTER PROC*28 SAVE /PYINT6/ DIMENSION WDTP(0:40),WDTE(0:40,0:5) - kc=0 - C...Calculate full and effective widths of gauge bosons. AEM=PARU(101) XW=PARU(102) @@ -7096,8 +7048,8 @@ SUBROUTINE PYMAXI SAVE /PYINT2/ COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000) SAVE /PYINT3/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3) SAVE /PYINT5/ COMMON/PYINT6/PROC(0:200) @@ -7108,20 +7060,7 @@ SUBROUTINE PYMAXI &NAREL(6),WTREL(6),WTMAT(6,6),COEFU(6),IACCMX(4),SIGSMX(4), &SIGSSM(3) DATA CVAR/'tau ','tau''','y* ','cth '/ - - taur1=0. - gamr1=0. - taur2=0. - gamr2=0. - atau3=0. - atau4=0. - atau5=0. - atau6=0. - ioff=0 - vvar=0. - vdel=0. - vmar=0. - + C...Select subprocess to study: skip cases not applicable. VINT(143)=1. VINT(144)=1. @@ -7674,8 +7613,8 @@ SUBROUTINE PYRAND SAVE /PYINT2/ COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000) SAVE /PYINT3/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3) SAVE /PYINT5/ @@ -8057,14 +7996,11 @@ SUBROUTINE PYSCAT SAVE /PYINT2/ COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000) SAVE /PYINT3/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3) SAVE /PYINT5/ DIMENSION WDTP(0:40),WDTE(0:40,0:5),PMQ(2),Z(2),CTHE(2),PHI(2) - - kflq=0 - phir=0. C...Choice of subprocess, number of documentation lines. ISUB=MINT(1) @@ -9097,14 +9033,6 @@ SUBROUTINE PYSSPA(IPU1,IPU2) &XFS(2,-6:6),XFA(-6:6),XFB(-6:6),XFN(-6:6),WTAP(-6:6),WTSF(-6:6), &THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),DPB(4) - tevb=0. - kfla=0 - z=0. - the2t=0. - ipo=0 - dmsma=0. - dpt2=0. - C...Calculate maximum virtuality and check that evolution possible. IPUS1=IPU1 IPUS2=IPU2 @@ -9505,13 +9433,6 @@ SUBROUTINE PYMULT(MMUL) DIMENSION NMUL(20),SIGM(20),KSTR(500,2) SAVE XT2,XT2FAC,XC2,XTS,IRBIN,RBIN,NMUL,SIGM - xf=0. - yf=0. - deltab=0. - ist1=0 - ist2=0 - istm=0 - C...Initialization of multiple interaction treatment. IF(MMUL.EQ.1) THEN IF(MSTP(122).GE.1) WRITE(MSTU(11),1000) MSTP(82) @@ -9967,16 +9888,6 @@ SUBROUTINE PYREMN(IPU1,IPU2) SAVE /PYINT1/ DIMENSION KFLCH(2),KFLSP(2),CHI(2),PMS(6),IS(2),ROBO(5) - iq=0 - ipu=0 - shs=0. - jpt=0 - peh=0. - pzh=0. - pei=0. - pzi=0. - pz=0. - C...Special case for lepton-lepton interaction. IF(MINT(43).EQ.1) THEN DO 100 JT=1,2 @@ -10312,8 +10223,8 @@ SUBROUTINE PYRESD SAVE /PYINT1/ COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2) SAVE /PYINT2/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ DIMENSION IREF(10,6),KDCY(2),KFL1(2),KFL2(2),NSD(2),ILIN(6), &COUP(6,4),PK(6,4),PKK(6,6),CTHE(2),PHI(2),WDTP(0:40), &WDTE(0:40,0:5) @@ -10330,10 +10241,6 @@ SUBROUTINE PYRESD &2.d0*DT*DU*(DT*DU/(D34*D56)-2.d0*(1.d0/D34+1.d0/D56)*(DT+DU)+ &2.d0*(D34/D56+D56/D34)) - i12=0 - wt=0. - wtmax=0. - C...Define initial two objects, initialize loop. ISUB=MINT(1) SH=VINT(44) @@ -10815,8 +10722,6 @@ SUBROUTINE PYDIFF COMMON/PYINT1/MINT(400),VINT(400) SAVE /PYINT1/ - chi=0. - C...Reset K, P and V vectors. Store incoming particles. DO 100 JT=1,MSTP(126)+10 I=MINT(83)+JT @@ -10999,26 +10904,9 @@ SUBROUTINE PYWIDT(KFLR,RMAS,WDTP,WDTE) SAVE /PYPARS/ COMMON/PYINT1/MINT(400),VINT(400) SAVE /PYINT1/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ DIMENSION WDTP(0:40),WDTE(0:40,0:5) - - wid2=0. - ai=0. - ggi=0. - gzi=0. - zzi=0. - ggf=0. - gzf=0. - zzf=0. - ej=0. - vj=0. - gzpi=0. - zzpi=0. - zpzpi=0. - gzpf=0. - zzpf=0. - zpzpf=0. C...Some common constants. KFLA=IABS(KFLR) @@ -11656,12 +11544,6 @@ SUBROUTINE PYKLIM(ILIM) SAVE /PYINT1/ COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2) SAVE /PYINT2/ - - tau=0. - rm3=0. - rm4=0. - be34=0. - st2eff=0. C...Common kinematical expressions. ISUB=MINT(1) @@ -11716,9 +11598,8 @@ SUBROUTINE PYKLIM(ILIM) ETA4=LOG(MIN(1.E10,MAX(1.E-10,EXPET4))) ETALAR=MAX(ETA3,ETA4) ETASMA=MIN(ETA3,ETA4) - 100 TMP1=(1.+RM3-RM4)*COSH(YST) - CTS3=((1.+RM3-RM4)*SINH(YST)+BE34*COSH(YST)*CTH)/ - & SQRT((TMP1+BE34*SINH(YST)*CTH)**2-4.*RM3) + 100 CTS3=((1.+RM3-RM4)*SINH(YST)+BE34*COSH(YST)*CTH)/ + & SQRT(((1.+RM3-RM4)*COSH(YST)+BE34*SINH(YST)*CTH)**2-4.*RM3) CTS4=((1.-RM3+RM4)*SINH(YST)-BE34*COSH(YST)*CTH)/ & SQRT(((1.-RM3+RM4)*COSH(YST)-BE34*SINH(YST)*CTH)**2-4.*RM4) CTSLAR=MAX(CTS3,CTS4) @@ -11925,10 +11806,6 @@ SUBROUTINE PYKMAP(IVAR,MVAR,VVAR) COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2) SAVE /PYINT2/ - taure=0. - gamre=0. - cth=0. - C...Convert VVAR to tau variable. ISUB=MINT(1) IF(IVAR.EQ.1) THEN @@ -12110,36 +11987,12 @@ SUBROUTINE PYSIGH(NCHN,SIGS) SAVE /PYINT2/ COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000) SAVE /PYINT3/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3) SAVE /PYINT5/ DIMENSION X(2),XPQ(-6:6),KFAC(2,-40:40),WDTP(0:40),WDTE(0:40,0:5) - as=0. - faca=0. - min1=0. - max1=0. - min2=0. - max2=0. - mina=0. - maxa=0. - sqmz=0. - gmmz=0. - sqmw=0. - gmmw=0. - sqmh=0. - gmmh=0. - sqmzp=0. - gmmzp=0. - sqmhc=0. - gmmhc=0. - sqmr=0. - gmmr=0. - aem=0. - xw=0. - comfac=0. - C...Reset number of channels and cross-section. NCHN=0 SIGS=0. @@ -14401,9 +14254,6 @@ SUBROUTINE PYSTFU(KF,X,Q2,XPQ,JBT) clin-10/25/02 get rid of argument usage mismatch in PYGAMM(): c EULBT(X,Y)=PYGAMM(X)*PYGAMM(Y)/PYGAMM(X+Y) - vx=0. - bbr2=0. - C...Reset structure functions, check x and hadron flavour. ALAM=0. DO 100 KFL=-6,6 @@ -14815,11 +14665,10 @@ FUNCTION PYW1AU(EPS,IREIM) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) SAVE /LUDAT1/ - ASINH(X)=LOG(X+SQRT(X**2+1.)) +clin-8/2014: +c ASINH(X)=LOG(X+SQRT(X**2+1.)) ACOSH(X)=LOG(X+SQRT(X**2-1.)) - pyw1au=0. - IF(EPS.LT.0.) THEN W1RE=2.*SQRT(1.-EPS)*ASINH(SQRT(-1./EPS)) W1IM=0. @@ -14847,11 +14696,10 @@ FUNCTION PYW2AU(EPS,IREIM) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) SAVE /LUDAT1/ - ASINH(X)=LOG(X+SQRT(X**2+1.)) +clin-8/2014: +c ASINH(X)=LOG(X+SQRT(X**2+1.)) ACOSH(X)=LOG(X+SQRT(X**2-1.)) - pyw2au=0. - IF(EPS.LT.0.) THEN W2RE=4.*(ASINH(SQRT(-1./EPS)))**2 W2IM=0. @@ -14879,9 +14727,6 @@ FUNCTION PYI3AU(BE,EPS,IREIM) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) SAVE /LUDAT1/ - pyi3au=0. - ga=0. - IF(EPS.LT.1.) GA=0.5*(1.+SQRT(1.-EPS)) IF(EPS.LT.0.) THEN @@ -14934,8 +14779,6 @@ FUNCTION PYSPEN(XREIN,XIMIN,IREIM) & 0.000000E+00, 7.575757E-02, 0.000000E+00, &-2.531135E-01, 0.000000E+00, 1.166667E+00/ - pyspen=0. - XRE=XREIN XIM=XIMIN IF(ABS(1.-XRE).LT.1.E-6.AND.ABS(XIM).LT.1.E-6) THEN @@ -15020,8 +14863,8 @@ BLOCK DATA PYDATA SAVE /PYINT2/ COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000) SAVE /PYINT3/ - COMMON/PYINT4A/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) - SAVE /PYINT4A/ + COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3) + SAVE /PYINT4/ COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3) SAVE /PYINT5/ COMMON/PYINT6/PROC(0:200) @@ -15336,9 +15179,6 @@ SUBROUTINE PYSTFE(KF,X,Q2,XPQ) DATA HEADER/'Tung evolution package has been invoked'/ DATA INIT/0/ - KF=KF - HEADER=HEADER - CHDFLM(1)=CHDFLM(1) C...Proton structure functions from Diemoz, Ferroni, Longo, Martinelli. C...Allowed variable range 10 GeV2 < Q2 < 1E8 GeV2, 5E-5 < x < .95. IF(MSTP(51).GE.11.AND.MSTP(51).LE.13.AND.MSTP(52).LE.1) THEN diff --git a/GeneratorInterface/AMPTInterface/src/linana.f b/GeneratorInterface/AMPTInterface/src/linana.f index 6f97a4cd11a25..56e7f18e0be0e 100755 --- a/GeneratorInterface/AMPTInterface/src/linana.f +++ b/GeneratorInterface/AMPTInterface/src/linana.f @@ -3,9 +3,6 @@ c 10/26/01 update freezeout positions in case of interactions: clin-3/2009 Note: freezeout spacetime values cannot be trusted for K0S & K0L c as K0S/K0L are converted from K+/K- by hand at the end of hadron cascade. -cms -cms dlw & gsfs Commented out output file writing -cms subroutine hbtout(nnew,nt,ntmax) c PARAMETER (MAXSTR=150001,MAXR=1) @@ -32,7 +29,7 @@ subroutine hbtout(nnew,nt,ntmax) cc SAVE /tdecay/ COMMON /AREVT/ IAEVT, IARUN, MISS cc SAVE /AREVT/ - common/snn/efrm,npart1,npart2 + common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG cc SAVE /snn/ COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP cc SAVE /HJGLBR/ @@ -44,6 +41,8 @@ subroutine hbtout(nnew,nt,ntmax) COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) EXTERNAL IARFLV, INVFLV common /para8/ idpert,npertd,idxsec +clin-2/2012: + common /phiHJ/iphirp,phiRP SAVE c do 1001 i=1,max0(nlast,nnew) @@ -164,12 +163,10 @@ subroutine hbtout(nnew,nt,ntmax) nlast=nnew ctest off look inside each NT timestep (for debugging purpose): c do ip=1,nlast -c if(nt.eq.5000) then -c write(95,*) ' p ',nt,ip,INVFLV(lblast(ip)),plast(1,ip), -c 1 plast(2,ip),plast(3,ip),plast(4,ip),dplast(ip) -c write(95,*) ' x ',nt,ip,INVFLV(lblast(ip)),xlast(1,ip), -c 1 xlast(2,ip),xlast(3,ip),xlast(4,ip),dplast(ip) -c endif +c write(99,*) ' p ',nt,ip,lblast(ip),plast(1,ip), +c 1 plast(2,ip),plast(3,ip),plast(4,ip),dplast(ip) +c write(99,*) ' x ',nt,ip,lblast(ip),xlast(1,ip), +c 1 xlast(2,ip),xlast(3,ip),xlast(4,ip),dplast(ip) c enddo c if(nt.eq.ntmax) then @@ -184,8 +181,11 @@ subroutine hbtout(nnew,nt,ntmax) c c write(16,190) IAEVT,IARUN,nlast,bimp,npart1,npart2, c 1 NELP,NINP,NELT,NINTHJ - write(16,190) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2, - 1 NELP,NINP,NELT,NINTHJ +clin-2/2012: +c write(16,190) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2, +c 1 NELP,NINP,NELT,NINTHJ + write(16,191) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2, + 1 NELP,NINP,NELT,NINTHJ,phiRP clin-5/2008 write out perturbatively-produced particles (deuterons only): if(idpert.eq.1.or.idpert.eq.2) 1 write(90,190) IAEVT,IARUN,ndpert,bimp,npart1,npart2, @@ -197,9 +197,15 @@ subroutine hbtout(nnew,nt,ntmax) c that are limited by the time-resolution (DT) of the hadron cascade, c freezeout time of spectator projectile or target nucleons is written as c DT as they are read at the 1st timestep and then propagated to time DT: - if(plast(1,ip).eq.0.and.plast(2,ip).eq.0 - 1 .and.(sqrt(plast(3,ip)**2+plast(4,ip)**2)*2/HINT1(1)) - 2 .gt.0.99.and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then +c +clin-9/2011 determine spectator nucleons consistently +c if(plast(1,ip).eq.0.and.plast(2,ip).eq.0 +c 1 .and.(sqrt(plast(3,ip)**2+plast(4,ip)**2)*2/HINT1(1)) +c 2 .gt.0.99.and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then + if(abs(plast(1,ip)).le.epsiPt.and.abs(plast(2,ip)).le.epsiPt + 1 .and.(plast(3,ip).gt.amax1(0.,PZPROJ-epsiPz) + 2 .or.plast(3,ip).lt.(-PZTARG+epsiPz)) + 3 .and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then clin-5/2008 perturbatively-produced particles (currently only deuterons) c are written to ana/ampt_pert.dat (without the column for the mass); c ana/ampt.dat has regularly-produced particles (including deuterons); @@ -259,6 +265,7 @@ subroutine hbtout(nnew,nt,ntmax) if(ioscar.eq.1) call hoscar endif 190 format(3(i7),f10.4,5x,6(i4)) + 191 format(3(i7),f10.4,5x,6(i4),5x,f7.4) clin-3/2009 improve the output accuracy of Pz 200 format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,f8.2)) 201 format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,e8.2)) @@ -283,22 +290,21 @@ SUBROUTINE decomp(px0,py0,pz0,xm0,i,itq1) COMMON/RNDF77/NSEED cc SAVE /RNDF77/ COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 - common/embed/iembed,pxqembd,pyqembd,xembd,yembd + common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd, + 1 psembd,tmaxembd,phidecomp SAVE c - itq1=itq1 dcth=dble(RANART(NSEED))*2.d0-1.d0 dPHI=dble(RANART(NSEED)*HIPR1(40))*2.d0 clin-6/2009 Added if embedding a high-Pt quark pair after string melting: - if(iembed.eq.1) then + if(iembed.ge.1.and.iembed.le.4) then c Decompose the parent high-Pt pion to q and qbar with an internal momentum c parallel to the pion direction so that one parton has ~the same hight Pt c and the other parton has a very soft Pt: c Note: htop() decomposes a meson to q as it(1) followed by qbar as it(2): - if(i.eq.(natt-1).or.i.eq.natt) then + if(i.eq.(natt-2*nsembd).or.i.eq.(natt-2*nsembd-1)) then dcth=0.d0 - dphi=dble(acos(pxqembd/sqrt(pxqembd**2+pyqembd**2))) - if(pyqembd.lt.0) dphi=HIPR1(40)*2.d0-dphi + dphi=dble(phidecomp) endif endif c @@ -347,7 +353,6 @@ SUBROUTINE HTOP 1 GXSGS,GYSGS,GZSGS,FTSGS dimension it(4) COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4) -cwei DOUBLE PRECISION PATT cc SAVE /HMAIN2/ COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 cc SAVE /HMAIN1/ @@ -389,6 +394,7 @@ SUBROUTINE HTOP cc SAVE /precpa/ common /para7/ ioscar,nsmbbbar,nsmmeson COMMON /AREVT/ IAEVT, IARUN, MISS + common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG SAVE c npar=0 @@ -402,10 +408,15 @@ SUBROUTINE HTOP id=ITYPAR(i) idabs=iabs(id) i2=MOD(idabs/10,10) - if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i) - 1 .ge.(HINT1(1)/2*0.99).and. - 2 ((id.eq.2112).or.(id.eq.2212))) then -c proj or targ nucleons without interactions, do not enter ZPC: +clin-9/2011 determine spectator nucleons consistently +c if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i) +c 1 .ge.(HINT1(1)/2*0.99).and. +c 2 .and.(id.eq.2112.or.id.eq.2212)) then + if(abs(PXAR(i)).le.epsiPt.and.abs(PYAR(i)).le.epsiPt + 1 .and.(PZAR(i).gt.amax1(0.,PZPROJ-epsiPz) + 2 .or.PZAR(i).lt.(-PZTARG+epsiPz)) + 3 .and.(id.eq.2112.or.id.eq.2212)) then +c spectator proj or targ nucleons without interactions, do not enter ZPC: elseif(idabs.gt.1000.and.i2.ne.0) then c baryons to be converted to q/qbar: nsmbbbar=nsmbbbar+1 @@ -446,9 +457,14 @@ SUBROUTINE HTOP it(3)=0 it(4)=0 c - if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i) - 1 .ge.(HINT1(1)/2*0.99).and.((id.eq.2112).or.(id.eq.2212))) then -c proj or targ nucleons without interactions, do not enter ZPC: +clin-9/2011 determine spectator nucleons consistently +c if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i) +c 1 .ge.(HINT1(1)/2*0.99).and.((id.eq.2112).or.(id.eq.2212))) then + if(abs(PXAR(i)).le.epsiPt.and.abs(PYAR(i)).le.epsiPt + 1 .and.(PZAR(i).gt.amax1(0.,PZPROJ-epsiPz) + 2 .or.PZAR(i).lt.(-PZTARG+epsiPz)) + 3 .and.(id.eq.2112.or.id.eq.2212)) then +c spectator proj or targ nucleons without interactions, do not enter ZPC: inozpc=1 elseif(idabs.gt.1000.and.i2.ne.0) then c baryons: @@ -583,8 +599,8 @@ SUBROUTINE HTOP vyp0(npar)=dble(patt(i,2)/patt(i,4)) vzp0(npar)=dble(patt(i,3)/patt(i,4)) 1001 continue -cyy 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2)) -cyy 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2)) + 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2)) + 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2)) c if((isoft.eq.4.or.isoft.eq.5) 1 .and.iabs(it(2)).gt.1000) then @@ -649,6 +665,8 @@ SUBROUTINE PTOH 1 GXSGS,GYSGS,GZSGS,FTSGS DOUBLE PRECISION xmdiag,px1,py1,pz1,e1,px2,py2,pz2,e2, 1 px3,py3,pz3,e3,xmpair,etot +clin-9/2012: improve precision for argument in sqrt(): + DOUBLE PRECISION p1,p2,p3 common /loclco/gxp(3),gyp(3),gzp(3),ftp(3), 1 pxp(3),pyp(3),pzp(3),pep(3),pmp(3) cc SAVE /loclco/ @@ -689,9 +707,10 @@ SUBROUTINE PTOH COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP COMMON /AREVT/ IAEVT, IARUN, MISS common /para7/ ioscar,nsmbbbar,nsmmeson +clin-5/2011 + common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT c dimension xmdiag(MAXSTR),indx(MAXSTR),ndiag(MAXSTR) -cwei DOUBLE PRECISION PATT SAVE c call coales @@ -873,6 +892,17 @@ SUBROUTINE PTOH nrhoch=nrhoch+1 endif endif +clin-7/2011-check charm hadron flavors: +c if(k1abs.eq.4.or.k2abs.eq.4) then +c if(k3.eq.0) then +c write(99,*) iaevt,k1,k2,kf,xmpair, +c 1 ULMASS(iabs(kf)),ULMASS(iabs(kf)+2),isg +c else +c write(99,*) iaevt,k1,k2,k3,kf,xmpair, +c 1 ULMASS(iabs(kf)),ULMASS(iabs(kf)+2),isg +c endif +c endif +clin-7/2011-end 1001 CONTINUE c assume Npi0=(Npi+ + Npi-)/2, Nrho0=(Nrho+ + Nrho-)/2 on the average: if(nuudd.ne.0) then @@ -945,6 +975,11 @@ SUBROUTINE PTOH PATT(inatt,2)=PYAR(inatt) PATT(inatt,3)=PZAR(inatt) etot=e1+e2 +clin-9/2012: improve precision for argument in sqrt(): + p1=px1+px2 + p2=py1+py2 + p3=pz1+pz2 +c elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3) 1 then PX3=PXSGS(ISG,3) @@ -958,8 +993,22 @@ SUBROUTINE PTOH PATT(inatt,2)=PYAR(inatt) PATT(inatt,3)=PZAR(inatt) etot=e1+e2+e3 +clin-9/2012: improve precision for argument in sqrt(): + p1=px1+px2+px3 + p2=py1+py2+py3 + p3=pz1+pz2+pz3 +c endif XMAR(inatt)=ULMASS(ITYPAR(inatt)) +clin-5/2011-add finite width to resonances (rho,omega,eta,K*,phi,Delta) after formation: + kf=KATT(inatt,1) + if(kf.eq.113.or.abs(kf).eq.213.or.kf.eq.221.or.kf.eq.223 + 1 .or.abs(kf).eq.313.or.abs(kf).eq.323.or.kf.eq.333 + 2 .or.abs(kf).eq.1114.or.abs(kf).eq.2114 + 3 .or.abs(kf).eq.2214.or.abs(kf).eq.2224) then + XMAR(inatt)=resmass(kf) + endif +c PEAR(inatt)=sqrt(PXAR(inatt)**2+PYAR(inatt)**2 1 +PZAR(inatt)**2+XMAR(inatt)**2) PATT(inatt,4)=PEAR(inatt) @@ -988,9 +1037,14 @@ SUBROUTINE PTOH gyavg0=gyavg0+gyp0(i)/ipartn gzavg0=gzavg0+gzp0(i)/ipartn 1005 CONTINUE - bex=dble(PXAR(inatt))/etot - bey=dble(PYAR(inatt))/etot - bez=dble(PZAR(inatt))/etot +clin-9/2012: improve precision for argument in sqrt(): +c bex=dble(PXAR(inatt))/etot +c bey=dble(PYAR(inatt))/etot +c bez=dble(PZAR(inatt))/etot + bex=p1/etot + bey=p2/etot + bez=p3/etot +c beta2 = bex ** 2 + bey ** 2 + bez ** 2 gam = 1.d0 / dsqrt(1.d0 - beta2) if(beta2.ge.0.9999999999999d0) then @@ -1004,14 +1058,16 @@ SUBROUTINE PTOH FTAR(inatt)=sngl(enenew) clin 4/19/2006 write out parton info after coalescence: if(ioscar.eq.3) then - WRITE (85, 312) K2SGS(ISG,1),px1,py1,pz1,PMSGS(ISG,1), - 1 inatt,katt(inatt,1) + WRITE (85, 313) K2SGS(ISG,1),px1,py1,pz1,PMSGS(ISG,1), + 1 inatt,katt(inatt,1),xmar(inatt) WRITE (85, 312) K2SGS(ISG,2),px2,py2,pz2,PMSGS(ISG,2), 1 inatt,katt(inatt,1) if(NJSGS(ISG).eq.3) WRITE (85, 312) K2SGS(ISG,3), 1 px3,py3,pz3,PMSGS(ISG,3),inatt,katt(inatt,1) endif 312 FORMAT(I6,4(1X,F10.3),1X,I6,1X,I6) +clin-5/02/2011 + 313 FORMAT(I6,4(1X,F10.3),1X,I6,1X,I6,1X,F10.3) c endif 1006 CONTINUE @@ -1022,6 +1078,56 @@ SUBROUTINE PTOH RETURN END +c======================================================================= +clin-5/2011-add finite width to resonances (rho,omega,eta,K*,phi,Delta) after formation: + FUNCTION resmass(kf) + + PARAMETER (arho=0.775,aomega=0.783,aeta=0.548,aks=0.894, + 1 aphi=1.019,adelta=1.232) + PARAMETER (wrho=0.149,womega=0.00849,weta=1.30E-6,wks=0.0498, + 1 wphi=0.00426,wdelta=0.118) + common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT + COMMON/RNDF77/NSEED + SAVE + + if(kf.eq.113.or.abs(kf).eq.213) then + amass=arho + wid=wrho + elseif(kf.eq.221) then + amass=aeta + wid=weta + elseif(kf.eq.223) then + amass=aomega + wid=womega + elseif(abs(kf).eq.313.or.abs(kf).eq.323) then + amass=aks + wid=wks + elseif(kf.eq.333) then + amass=aphi + wid=wphi + elseif(abs(kf).eq.1114.or.abs(kf).eq.2114 + 1 .or.abs(kf).eq.2214.or.abs(kf).eq.2224) then + amass=adelta + wid=wdelta + endif + dmin=amass-2*wid + dmax=amass+2*wid +c Delta mass needs to be big enough to decay to N+pi: + if(amass.eq.adelta) dmin=1.078 +c + FM=1. + NTRY1=0 + 10 DM = RANART(NSEED) * (DMAX-DMIN) + DMIN + NTRY1=NTRY1+1 + fmass=(amass*wid)**2/((DM**2-amass**2)**2+(amass*wid)**2) +check write (99,*) ntry1,kf,amass,wid,fmass,DM + IF((RANART(NSEED) .GT. FMASS/FM).AND. (NTRY1.LE.10)) GOTO 10 +c + resmass=DM + + RETURN + END + c======================================================================= SUBROUTINE coales @@ -1458,7 +1564,7 @@ subroutine hoscar parameter (MAXSTR=150001,AMN=0.939457,AMP=0.93828) character*8 code, reffra, FRAME character*25 amptvn - common/snn/efrm,npart1,npart2 + common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG cc SAVE /snn/ common /lastt/itimeh,bimp cc SAVE /lastt/ @@ -1530,9 +1636,8 @@ subroutine getnp cc SAVE /HMAIN2/ COMMON /HPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50) cc SAVE /HPARNT/ - common/snn/efrm,npart1,npart2 + common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG cc SAVE /snn/ -cwei DOUBLE PRECISION PATT SAVE if(NATT.eq.0) then @@ -1543,20 +1648,27 @@ subroutine getnp c PZPROJ=SQRT(HINT1(6)**2-HINT1(8)**2) PZTARG=SQRT(HINT1(7)**2-HINT1(9)**2) - epsiln=0.01 + epsiPz=0.01 +clin-9/2011-add Pt tolerance in determining spectator nucleons +c (affect string melting runs when LAB frame is used): + epsiPt=1e-6 c nspec1=0 nspec2=0 DO 1000 I = 1, NATT +clin-9/2011 determine spectator nucleons consistently +c if((KATT(I,1).eq.2112.or.KATT(I,1).eq.2212) +c 1 .and.PATT(I, 1).eq.0.and.PATT(I, 2).eq.0) then if((KATT(I,1).eq.2112.or.KATT(I,1).eq.2212) - 1 .and.PATT(I, 1).eq.0.and.PATT(I, 2).eq.0) then - if(PATT(I, 3).gt.amax1(0.,PZPROJ-epsiln)) then + 1 .and.abs(PATT(I, 1)).le.epsiPt + 2 .and.abs(PATT(I, 2)).le.epsiPt) then + if(PATT(I, 3).gt.amax1(0.,PZPROJ-epsiPz)) then nspec1=nspec1+1 - elseif(PATT(I, 3).lt.(-PZTARG+epsiln)) then + elseif(PATT(I, 3).lt.(-PZTARG+epsiPz)) then nspec2=nspec2+1 endif endif - 1000 CONTINUE + 1000 CONTINUE npart1=IHNT2(1)-nspec1 npart2=IHNT2(3)-nspec2 @@ -1565,7 +1677,9 @@ subroutine getnp c======================================================================= c 2/18/03 use PYTHIA to decay eta,rho,omega,k*,phi and Delta - subroutine resdec(i1,nt,nnn,wid,idecay) +c 4/2012 added pi0 decay flag: +c ipion=0: resonance or pi0 in lb(i1); >0: pi0 in lpion(ipion). + subroutine resdec(i1,nt,nnn,wid,idecay,ipion) PARAMETER (hbarc=0.19733) PARAMETER (AK0=0.498,APICH=0.140,API0=0.135,AN=0.940,ADDM=0.02) @@ -1608,9 +1722,15 @@ subroutine resdec(i1,nt,nnn,wid,idecay) 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) cc SAVE /RNDF77/ + common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy SAVE irun=idecay - if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27 +clin-4/2012 for option of pi0 decay: + if(nt.eq.ntmax.and.ipi0dcy.eq.1 + & .and.((lb1.eq.4.and.ipion.eq.0).or.ipion.ge.1)) then + kf=111 +c if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27 + elseif(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27 & .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30 & .or.lb1.eq.24.or.(iabs(lb1).ge.6.and.iabs(lb1).le.9) & .or.iabs(lb1).eq.16) then @@ -1628,10 +1748,13 @@ subroutine resdec(i1,nt,nnn,wid,idecay) K(IP,5)=0 c K(IP,2)=kf +clin-4/2012 for option of pi0 decay: + if(ipion.eq.0) then +c P(IP,1)=px1 P(IP,2)=py1 P(IP,3)=pz1 - em1a=em1 +c em1a=em1 c eta or omega in ART may be below or too close to (pi+pi-pi0) mass, c causing LUDECY error,thus increase their mass ADDM above this thresh, c noting that rho (m=0.281) too close to 2pi thrshold fails to decay: @@ -1655,6 +1778,21 @@ subroutine resdec(i1,nt,nnn,wid,idecay) P(IP,5)=em1 clin-5/2008: dpdecp=dpertp(i1) +clin-4/2012 for option of pi0 decay: + elseif(nt.eq.ntmax.and.ipi0dcy.eq.1.and.ipion.ge.1) then + P(IP,1)=PPION(1,ipion,IRUN) + P(IP,2)=PPION(2,ipion,IRUN) + P(IP,3)=PPION(3,ipion,IRUN) + P(IP,5)=EPION(ipion,IRUN) + P(IP,4)=SQRT(P(IP,5)**2+P(IP,1)**2+P(IP,2)**2+P(IP,3)**2) + dpdecp=dppion(ipion,IRUN) +ctest off +c write(99,*) P(IP,4), P(IP,5), dpdecp, ipion, wid + else + print *, 'stopped in resdec() a' + stop + endif +c call ludecy(IP) c add decay time to daughter's formation time at the last timestep: if(nt.eq.ntmax) then @@ -1663,19 +1801,32 @@ subroutine resdec(i1,nt,nnn,wid,idecay) ndaut=n-nsav if(ndaut.le.1) then write(10,*) 'note: ndaut(<1)=',ndaut - write(89,*) 'note: ndaut(<1)=',ndaut call lulist(2) stop endif c lorentz boost: - taudcy=taudcy*e1/em1 - tfnl=tfnl+taudcy - xfnl=xfnl+px1/e1*taudcy - yfnl=yfnl+py1/e1*taudcy - zfnl=zfnl+pz1/e1*taudcy -c at the last timestep, assign rho, K0S or eta (decay daughter) +clin-4/2012 for option of pi0 decay: + if(ipion.eq.0) then + taudcy=taudcy*e1/em1 + tfnl=tfnl+taudcy + xfnl=xfnl+px1/e1*taudcy + yfnl=yfnl+py1/e1*taudcy + zfnl=zfnl+pz1/e1*taudcy + elseif(ipion.ge.1) then + taudcy=taudcy*P(IP,4)/P(IP,5) + tfnl=tfdpi(ipion,IRUN)+taudcy + xfnl=RPION(1,ipion,IRUN)+P(IP,1)/P(IP,4)*taudcy + yfnl=RPION(2,ipion,IRUN)+P(IP,2)/P(IP,4)*taudcy + zfnl=RPION(3,ipion,IRUN)+P(IP,3)/P(IP,4)*taudcy + else + print *, 'stopped in resdec() b',ipion,wid,P(ip,4) + stop + endif +c at the last timestep, assign rho, K0S or eta (decay daughter) c to lb(i1) only (not to lpion) in order to decay them again: - if(n.ge.(nsav+2)) then +clin-4/2012 for option of pi0 decay: +c if(n.ge.(nsav+2)) then + if(n.ge.(nsav+2).and.ipion.eq.0) then do 1001 idau=nsav+2,n kdaut=K(idau,2) if(kdaut.eq.221.or.kdaut.eq.113 @@ -1716,10 +1867,6 @@ subroutine resdec(i1,nt,nnn,wid,idecay) c 1 write(93,*) 'resdec(): nt=',nt,enet-e1,lb1 endif -cyy 200 format(a20,3(1x,i6)) -cyy 210 format(i6,5(1x,f8.3)) -cyy 220 format(a2,i5,5(1x,f8.3)) - do 1003 idau=nsav+1,n kdaut=K(idau,2) lbdaut=IARFLV(kdaut) @@ -1742,13 +1889,28 @@ subroutine resdec(i1,nt,nnn,wid,idecay) endif c if(idau.eq.(nsav+1)) then - LB(i1)=lbdaut - E(i1)=p(idau,5) - px1n=p(idau,1) - py1n=p(idau,2) - pz1n=p(idau,3) +clin-4/2012 for option of pi0 decay: + if(ipion.eq.0) then + LB(i1)=lbdaut + E(i1)=p(idau,5) + px1n=p(idau,1) + py1n=p(idau,2) + pz1n=p(idau,3) clin-5/2008: - dp1n=dpdecp + dp1n=dpdecp + elseif(ipion.ge.1) then + LPION(ipion,IRUN)=lbdaut + EPION(ipion,IRUN)=p(idau,5) + PPION(1,ipion,IRUN)=p(idau,1) + PPION(2,ipion,IRUN)=p(idau,2) + PPION(3,ipion,IRUN)=p(idau,3) + RPION(1,ipion,IRUN)=xfnl + RPION(2,ipion,IRUN)=yfnl + RPION(3,ipion,IRUN)=zfnl + tfdpi(ipion,IRUN)=tfnl + dppion(ipion,IRUN)=dpdecp + endif +c else nnn=nnn+1 LPION(NNN,IRUN)=lbdaut @@ -1764,7 +1926,6 @@ subroutine resdec(i1,nt,nnn,wid,idecay) dppion(NNN,IRUN)=dpdecp endif 1003 continue -cyy 230 format(a2,i5,5(1x,e8.2)) return end @@ -1923,7 +2084,7 @@ subroutine inifrz return end -clin-5/2009 ctest on v2 analysis +clin-5/2009 v2 analysis c======================================================================= c idd=0,1,2,3 specifies different subroutines for partonic flow analysis. subroutine flowp(idd) @@ -1990,72 +2151,72 @@ subroutine flowp(idd) nfile(1)=60 nfile(2)=64 nfile(3)=20 -cms OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(1)+1, -cms 1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(1)+2, -cms 1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(1)+3, -cms 1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(2)+1, -cms 1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(2)+2, -cms 1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(2)+3, -cms 1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(3)+1, -cms 1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(3)+2, -cms 1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN') -cms OPEN (nfile(3)+3, -cms 1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN') -cms OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN') -cms write(49, *) ' ievt, v2p, v2p_y2, v2p_y1' -c -cms OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN') -cms OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN') -cms OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN') -cms OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN') -cms write(88, *) ' ievt, v2h, v2h_y2, v2h_y1' + OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN') + OPEN (nfile(1)+1, + 1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN') + OPEN (nfile(1)+2, + 1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN') + OPEN (nfile(1)+3, + 1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN') + OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN') + OPEN (nfile(2)+1, + 1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN') + OPEN (nfile(2)+2, + 1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN') + OPEN (nfile(2)+3, + 1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN') + OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN') + OPEN (nfile(3)+1, + 1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN') + OPEN (nfile(3)+2, + 1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN') + OPEN (nfile(3)+3, + 1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN') + OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN') + write(49, *) ' ievt, v2p, v2p_y2, v2p_y1' +c + OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN') + OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN') + OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN') + OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN') + write(88, *) ' ievt, v2h, v2h_y2, v2h_y1' csp07/05 nlfile(1)=70 nlfile(2)=72 nlfile(3)=74 -cms OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN') -cms OPEN (nlfile(1)+1, -cms 1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN') -cms OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN') -cms OPEN (nlfile(2)+1, -cms 1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN') -cms OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN') -cms OPEN (nlfile(3)+1, -cms 1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN') + OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN') + OPEN (nlfile(1)+1, + 1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN') + OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN') + OPEN (nlfile(2)+1, + 1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN') + OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN') + OPEN (nlfile(3)+1, + 1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN') nsfile(1)=76 nsfile(2)=78 nsfile(3)=80 -cms OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN') -cms OPEN (nsfile(1)+1, -cms 1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN') -cms OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN') -cms OPEN (nsfile(2)+1, -cms 1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN') -cms OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN') -cms OPEN (nsfile(3)+1, -cms 1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN') + OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN') + OPEN (nsfile(1)+1, + 1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN') + OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN') + OPEN (nsfile(2)+1, + 1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN') + OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN') + OPEN (nsfile(3)+1, + 1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN') nmfile(1)=82 nmfile(2)=83 nmfile(3)=84 -cms OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN') -cms OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN') -cms OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN') + OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN') + OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN') + OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN') clin-11/27/00 for animation: if(nevent.eq.1) then -cms OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN') -cms write(91,*) ntmax, dt -cms OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN') -cms OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN') +c OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN') +c write(91,*) ntmax, dt +c OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN') + OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN') endif c itimep=0 @@ -2182,15 +2343,21 @@ subroutine flowp(idd) endif endif do 50 iy=1,iloop - if(pt2.gt.0.) then +clin-5/2012: +c if(pt2.gt.0.) then + if(pt2.gt.0d0) then v2prtn=(PX5(I)**2-PY5(I)**2)/pt2 - if(abs(v2prtn).gt.1.) +clin-5/2012: +c if(abs(v2prtn).gt.1.) + if(dabs(v2prtn).gt.1d0) 1 write(nfile(iy),*) 'v2prtn>1',v2prtn v2p(ianp,iy)=v2p(ianp,iy)+v2prtn v2p2(ianp,iy)=v2p2(ianp,iy)+v2prtn**2 endif xperp2=GX5(I)**2+GY5(I)**2 - if(xperp2.gt.0.) +clin-5/2012: +c if(xperp2.gt.0.) + if(xperp2.gt.0d0) 1 s2p(ianp,iy)=s2p(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2 xnpart(ianp,iy)=xnpart(ianp,iy)+1d0 etp(ianp,iy)=etp(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) @@ -2201,7 +2368,9 @@ subroutine flowp(idd) if(FT5(I).le.t2time) then v2pf(ianp,iy)=v2pf(ianp,iy)+v2prtn v2pf2(ianp,iy)=v2pf2(ianp,iy)+v2prtn**2 - if(xperp2.gt.0.) +clin-5/2012: +c if(xperp2.gt.0.) + if(xperp2.gt.0d0) 1 s2pf(ianp,iy)=s2pf(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2 xnpf(ianp,iy)=xnpf(ianp,iy)+1d0 etpf(ianp,iy)=etpf(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) @@ -2289,13 +2458,17 @@ subroutine flowp(idd) if(FT5(I).lt.tscatt(ianp+1) 1 .and.FT5(I).ge.tscatt(ianp)) then do 1011 iy=1,iloop - if(pt2.gt.0.) then +clin-5/2012: +c if(pt2.gt.0.) then + if(pt2.gt.0d0) then v2prtn=(PX5(I)**2-PY5(I)**2)/pt2 v2pfrz(ianp,iy)=v2pfrz(ianp,iy)+v2prtn v2p2fz(ianp,iy)=v2p2fz(ianp,iy)+v2prtn**2 endif xperp2=GX5(I)**2+GY5(I)**2 - if(xperp2.gt.0.) s2pfrz(ianp,iy)= +clin-5/2012: +c if(xperp2.gt.0.) s2pfrz(ianp,iy)= + if(xperp2.gt.0d0) s2pfrz(ianp,iy)= 1 s2pfrz(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2 etpfrz(ianp,iy)=etpfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) xnpfrz(ianp,iy)=xnpfrz(ianp,iy)+1d0 @@ -2321,26 +2494,28 @@ subroutine flowp(idd) c idd=2: calculate average partonic elliptic flow, called from artdri.f, else if(idd.eq.2) then do 1012 i=1,3 -cms write(nfile(i),*) ' tsp, v2p, v2p2, '// -cms 1 ' s2p, etp, xmult, nevt, xnctot' -cms write ((nfile(i)+1),*) ' tsp, v2pf, v2pf2, '// -cms 1 ' s2pf, etpf, xnform, xcrate' -cms write ((nfile(i)+2),*) ' tsp, v2pa, v2pa2, '// -cms 1 ' s2pa, etpa, xmult_ap, xnform, nevt' -cms write ((nfile(i)+3),*) ' tsph, v2ph, v2ph2, '// -cms 1 ' s2ph, etph, xmult_(ap/2+h),xmult_ap/2,nevt' + write(nfile(i),*) ' tsp, v2p, v2p2, '// + 1 ' s2p, etp, xmult, nevt, xnctot' + write ((nfile(i)+1),*) ' tsp, v2pf, v2pf2, '// + 1 ' s2pf, etpf, xnform, xcrate' + write ((nfile(i)+2),*) ' tsp, v2pa, v2pa2, '// + 1 ' s2pa, etpa, xmult_ap, xnform, nevt' + write ((nfile(i)+3),*) ' tsph, v2ph, v2ph2, '// + 1 ' s2ph, etph, xmult_(ap/2+h),xmult_ap/2,nevt' csp -cms write(nlfile(i),*) ' tsp, v2, s2, etp, xmul' -cms write(nsfile(i),*) ' tsp, v2, s2, etp, xmul' -cms write(nlfile(i)+1,*) ' tsp, v2, s2, etp, xmul' -cms write(nsfile(i)+1,*) ' tsp, v2, s2, etp, xmul' + write(nlfile(i),*) ' tsp, v2, s2, etp, xmul' + write(nsfile(i),*) ' tsp, v2, s2, etp, xmul' + write(nlfile(i)+1,*) ' tsp, v2, s2, etp, xmul' + write(nsfile(i)+1,*) ' tsp, v2, s2, etp, xmul' c 1012 continue clin-3/30/00 ensemble average & variance of v2 (over particles in all events): do 150 ii=1, 30 if(nevt(ii).eq.0) goto 150 do 1014 iy=1,3 - if(xnpart(ii,iy).gt.1) then +clin-5/2012: +c if(xnpart(ii,iy).gt.1) then + if(xnpart(ii,iy).gt.1d0) then v2p(ii,iy)=v2p(ii,iy)/xnpart(ii,iy) v2p2(ii,iy)=dsqrt((v2p2(ii,iy)/xnpart(ii,iy) 1 -v2p(ii,iy)**2)/(xnpart(ii,iy)-1)) @@ -2357,13 +2532,15 @@ subroutine flowp(idd) xnctot=xnctot+xncoll(inum) 1013 continue if(nevt(1).ne.0) xnctot=xnctot/nevt(1) -cms write (nfile(iy),200) tsp(ii),v2p(ii,iy), -cms 1 v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot + write (nfile(iy),200) tsp(ii),v2p(ii,iy), + 1 v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot endif if(nevt(ii).ne.0) 1 xcrate=xncoll(ii)/(tsp(ii+1)-tsp(ii))/nevt(ii) c - if(xnpf(ii,iy).gt.1) then +clin-5/2012: +c if(xnpf(ii,iy).gt.1) then + if(xnpf(ii,iy).gt.1d0) then v2pf(ii,iy)=v2pf(ii,iy)/xnpf(ii,iy) v2pf2(ii,iy)=dsqrt((v2pf2(ii,iy)/xnpf(ii,iy) 1 -v2pf(ii,iy)**2)/(xnpf(ii,iy)-1)) @@ -2374,41 +2551,49 @@ subroutine flowp(idd) etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii)) etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii)) c -cms write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy), -cms 1 v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate + write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy), + 1 v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate endif csp - if(xnpl(ii,iy).gt.1) then +clin-5/2012: +c if(xnpl(ii,iy).gt.1) then + if(xnpl(ii,iy).gt.1d0) then v2pl(ii,iy)=v2pl(ii,iy)/xnpl(ii,iy) s2pl(ii,iy)=s2pl(ii,iy)/xnpl(ii,iy) xmult=dble(xnpl(ii,iy)/dble(nevt(ii))) etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii)) -cms write (nlfile(iy),201) tsp(ii),v2pl(ii,iy), -cms 1 s2pl(ii,iy),etpl(ii,iy),xmult + write (nlfile(iy),201) tsp(ii),v2pl(ii,iy), + 1 s2pl(ii,iy),etpl(ii,iy),xmult endif - if(xnps(ii,iy).gt.1) then +clin-5/2012: +c if(xnps(ii,iy).gt.1) then + if(xnps(ii,iy).gt.1d0) then v2ps(ii,iy)=v2ps(ii,iy)/xnps(ii,iy) s2ps(ii,iy)=s2ps(ii,iy)/xnps(ii,iy) xmult=dble(xnps(ii,iy)/dble(nevt(ii))) etps(ii,iy)=etps(ii,iy)/dble(nevt(ii)) -cms write (nsfile(iy),201) tsp(ii),v2ps(ii,iy), -cms 1 s2ps(ii,iy),etps(ii,iy),xmult + write (nsfile(iy),201) tsp(ii),v2ps(ii,iy), + 1 s2ps(ii,iy),etps(ii,iy),xmult endif - if(xnplf(ii,iy).gt.1) then +clin-5/2012: +c if(xnplf(ii,iy).gt.1) then + if(xnplf(ii,iy).gt.1d0) then v2plf(ii,iy)=v2plf(ii,iy)/xnplf(ii,iy) s2plf(ii,iy)=s2plf(ii,iy)/xnplf(ii,iy) xmult=dble(xnplf(ii,iy)/dble(nevt(ii))) etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii)) -cms write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy), -cms 1 s2plf(ii,iy),etplf(ii,iy),xmult + write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy), + 1 s2plf(ii,iy),etplf(ii,iy),xmult endif - if(xnpsf(ii,iy).gt.1) then +clin-5/2012: +c if(xnpsf(ii,iy).gt.1) then + if(xnpsf(ii,iy).gt.1d0) then v2psf(ii,iy)=v2psf(ii,iy)/xnpsf(ii,iy) s2psf(ii,iy)=s2psf(ii,iy)/xnpsf(ii,iy) xmult=dble(xnpsf(ii,iy)/dble(nevt(ii))) etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii)) -cms write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy), -cms 1 s2psf(ii,iy),etpsf(ii,iy),xmult + write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy), + 1 s2psf(ii,iy),etpsf(ii,iy),xmult endif csp-end 1014 continue @@ -2422,9 +2607,9 @@ subroutine flowp(idd) yra = 1d0 if(iy .eq. 2)yra = 2d0 do 1015 i=1,50 -cms WRITE(nmfile(iy),251) BMT*dble(I - 0.5), -cms & SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra, -cms & SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra + WRITE(nmfile(iy),251) BMT*dble(I - 0.5), + & SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra, + & SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra 1015 continue 1016 continue csp07/05 end @@ -2433,7 +2618,9 @@ subroutine flowp(idd) do 1017 iy=1,3 v2pavg(iy)=v2psum(iy)/nevt(30) v2var0=v2p2sm(iy)/nevt(30)-v2pavg(iy)**2 - if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0) +clin-5/2012: +c if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0) + if(v2var0.gt.0d0) varv2p(iy)=dsqrt(v2var0) 1017 continue write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): avg=', v2pavg write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): var=', varv2p @@ -2462,7 +2649,9 @@ subroutine flowp(idd) s2pact=0d0 etpact=0d0 xnacti=0d0 - if(xnpf(ianp,iy).gt.1) then +clin-5/2012: +c if(xnpf(ianp,iy).gt.1) then + if(xnpf(ianp,iy).gt.1d0) then c reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons: v2pact=v2pf(ianp,iy)*xnpf(ianp,iy) v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1) @@ -2485,16 +2674,18 @@ subroutine flowp(idd) etph=etpact xnp2=xnpact/2d0 c - if(xnpact.gt.1.and.nevt(ianp).ne.0) then +clin-5/2012: +c if(xnpact.gt.1.and.nevt(ianp).ne.0) then + if(xnpact.gt.1d0.and.nevt(ianp).ne.0) then v2pact=v2pact/xnpact v2p2ac=dsqrt((v2p2ac/xnpact 1 -v2pact**2)/(xnpact-1)) s2pact=s2pact/xnpact xnacti=dble(xnpact/dble(nevt(ianp))) etpact=etpact/dble(nevt(ianp)) -cms write (nfile(iy)+2, 250) tsp(ianp),v2pact, -cms 1 v2p2ac,s2pact,etpact,xnacti, -cms 2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp) + write (nfile(iy)+2, 250) tsp(ianp),v2pact, + 1 v2p2ac,s2pact,etpact,xnacti, + 2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp) endif endif c To calculate combined v2 for active partons plus formed hadrons, @@ -2509,16 +2700,18 @@ subroutine flowp(idd) etph=etph+eth(ianh,iy)*dble(nevent)*shadr xnph=xnpact+xnhadr(ianh,iy)*shadr xnp2h=xnp2+xnhadr(ianh,iy)*shadr - if(xnph.gt.1) then +clin-5/2012: +c if(xnph.gt.1) then + if(xnph.gt.1d0) then v2ph=v2ph/xnph v2ph2=dsqrt((v2ph2/xnph-v2ph**2)/(xnph-1)) s2ph=s2ph/xnph etph=etph/dble(nevt(ianp)) xnp2=xnp2/dble(nevt(ianp)) xnp2h=xnp2h/dble(nevent) -cms if(tsp(ianp).le.(ntmax*dt)) -cms 1 write (nfile(iy)+3, 250) tsp(ianp),v2ph, -cms 2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp) + if(tsp(ianp).le.(ntmax*dt)) + 1 write (nfile(iy)+3, 250) tsp(ianp),v2ph, + 2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp) endif c 400 continue @@ -2539,24 +2732,24 @@ subroutine flowp(idd) xnpact=xnpf(ianp,iy) 500 continue 550 continue -cms close (620) -cms close (630) + close (620) + close (630) do 1021 nf=1,3 do 1020 ifile=0,3 -cms close(nfile(nf)+ifile) + close(nfile(nf)+ifile) 1020 continue 1021 continue do 1022 nf=1,3 -cms close(740+nf) + close(740+nf) 1022 continue endif -cyy 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2) -cyy 210 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2)) + 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2) + 210 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2)) 240 format(a30,3(2x,f9.5)) -cyy 250 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6) + 250 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6) csp -cyy 201 format(2x,f5.2,4(2x,f9.2)) -cyy 251 format(5e15.5) + 201 format(2x,f5.2,4(2x,f9.2)) + 251 format(5e15.5) c return end @@ -2897,8 +3090,8 @@ subroutine frztm(NEVNT,idd) dnfrz(ii)=0d0 detfrz(ii)=0d0 1001 continue -cms OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN') -cms OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN') + OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN') + OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN') else if(idd.eq.1) then DO 100 ip = 1, MUL do 1002 ii=1,30 @@ -2927,10 +3120,10 @@ subroutine frztm(NEVNT,idd) 1002 continue 100 continue else if(idd.eq.2) then -cms write (86,*) ' t, np, dnp/dt, etp '// -cms 1 ' detp/dt' -cms write (87,*) ' t, nf, dnf/dt, etf '// -cms 1 ' detf/dt' + write (86,*) ' t, np, dnp/dt, etp '// + 1 ' detp/dt' + write (87,*) ' t, nf, dnf/dt, etf '// + 1 ' detf/dt' do 1003 ii=1,30 xnp=xnprod(ii)/dble(NEVNT) xnf=xnfrz(ii)/dble(NEVNT) @@ -2940,13 +3133,13 @@ subroutine frztm(NEVNT,idd) dxnf=dnfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii)) detp=detpro(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii)) detf=detfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii)) -cms write (86, 200) -cms 1 tsf(ii+1),xnp,dxnp,etp,detp -cms write (87, 200) -cms 1 tsf(ii+1),xnf,dxnf,etf,detf + write (86, 200) + 1 tsf(ii+1),xnp,dxnp,etp,detp + write (87, 200) + 1 tsf(ii+1),xnf,dxnf,etf,detf 1003 continue endif -cyy 200 format(2x,f9.2,4(2x,f10.2)) + 200 format(2x,f9.2,4(2x,f10.2)) c return end @@ -2954,7 +3147,9 @@ subroutine frztm(NEVNT,idd) c======================================================================= clin-6/2009 write out initial minijet information c before propagating to its formation time: - subroutine minijet_out(BB) +clin-2/2012: +c subroutine minijet_out(BB) + subroutine minijet_out(BB,phiRP) PARAMETER (MAXSTR=150001) COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) COMMON/hjcrdn/YP(3,300),YT(3,300) @@ -2968,7 +3163,7 @@ subroutine minijet_out(BB) & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) COMMON /AREVT/ IAEVT, IARUN, MISS common /para7/ ioscar,nsmbbbar,nsmmeson - common/phidcy/iphidcy,pttrig,ntrig,maxmiss + common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy SAVE ntrig=0 do I = 1, IHNT2(1) @@ -2999,10 +3194,13 @@ subroutine minijet_out(BB) ityp=KFPJ(I,J) c write out not only gluons: c if(ityp.ne.21) goto 1007 - gx=YP(1,I)+0.5*BB - gy=YP(2,I) - gz=0d0 - ft=0d0 +clin-2/2012: +c gx=YP(1,I)+0.5*BB +c gy=YP(2,I) + gx=YP(1,I)+0.5*BB*cos(phiRP) + gy=YP(2,I)+0.5*BB*sin(phiRP) + gz=0. + ft=0. px=PJPX(I,J) py=PJPY(I,J) pz=PJPZ(I,J) @@ -3021,10 +3219,13 @@ subroutine minijet_out(BB) DO 1009 J = 1, NTJ(I) ityp=KFTJ(I,J) c if(ityp.ne.21) goto 1009 - gx=YT(1,I)-0.5*BB - gy=YT(2,I) - gz=0d0 - ft=0d0 +clin-2/2012: +c gx=YT(1,I)-0.5*BB +c gy=YT(2,I) + gx=YT(1,I)-0.5*BB*cos(phiRP) + gy=YT(2,I)-0.5*BB*sin(phiRP) + gz=0. + ft=0. px=PJTX(I,J) py=PJTY(I,J) pz=PJTZ(I,J) @@ -3045,8 +3246,8 @@ subroutine minijet_out(BB) c if(ityp.ne.21) goto 1011 gx=0.5*(YP(1,IASG(I,1))+YT(1,IASG(I,2))) gy=0.5*(YP(2,IASG(I,1))+YT(2,IASG(I,2))) - gz=0d0 - ft=0d0 + gz=0. + ft=0. px=PXSG(I,J) py=PYSG(I,J) pz=PZSG(I,J) @@ -3073,8 +3274,10 @@ subroutine minijet_out(BB) c by generating the internal quark and antiquark momentum parallel to c the pion momentum (in order to produce a high-Pt and a low Pt parton): subroutine embedHighPt - PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140) - common/embed/iembed,pxqembd,pyqembd,xembd,yembd + PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140,pi0mass=0.135, + 1 pi=3.1415926,nxymax=10001) + common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd, + 1 psembd,tmaxembd,phidecomp COMMON/RNDF77/NSEED COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4) @@ -3082,10 +3285,27 @@ subroutine embedHighPt & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR), & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR), & XMAR(MAXSTR) -cwei DOUBLE PRECISION PATT + common/anim/nevent,isoft,isflag,izpc + COMMON /AREVT/ IAEVT, IARUN, MISS + common/xyembed/nxyjet,xyjet(nxymax,2) SAVE c - if(iembed.ne.1) return + if(iembed.eq.1.or.iembed.eq.2) then + xjet=xembd + yjet=yembd + elseif(iembed.eq.3.or.iembed.eq.4) then + if(nevent.le.nxyjet) then + read(97,*) xjet,yjet + else + ixy=mod(IAEVT,nxyjet) + if(ixy.eq.0) ixy=nxyjet + xjet=xyjet(ixy,1) + yjet=xyjet(ixy,2) + endif + else + return + endif +c ptq=sqrt(pxqembd**2+pyqembd**2) if(ptq.lt.(pichmass/2.)) then print *, 'Embedded quark transverse momentum is too small' @@ -3096,10 +3316,10 @@ subroutine embedHighPt c Flavor content for the charged pion that contains the leading quark: if(idqembd.eq.1) then idqsoft=-2 - idpi=-211 + idpi1=-211 elseif(idqembd.eq.2) then idqsoft=-1 - idpi=211 + idpi1=211 else print *, 'Wrong quark flavor embedded' stop @@ -3111,19 +3331,32 @@ subroutine embedHighPt 1 -sqrt((xmq**2+ptq**2)*(pichmass**4 2 -2.*pichmass**2*(xmq**2+xmqsoft**2)+(xmq**2-xmqsoft**2)**2))) 3 /(2.*xmq**2) - pxpi=ptpi*pxqembd/ptq - pypi=ptpi*pyqembd/ptq + if(iembed.eq.1.or.iembed.eq.3) then + pxpi1=ptpi*pxqembd/ptq + pypi1=ptpi*pyqembd/ptq + phidecomp=acos(pxqembd/ptq) + if(pyqembd.lt.0) phidecomp=2.*pi-phidecomp + else + phidecomp=2.*pi*RANART(NSEED) + pxpi1=ptpi*cos(phidecomp) + pypi1=ptpi*sin(phidecomp) + endif c Embedded quark/antiquark are assumed to have pz=0: - pzpi=0. + pzpi1=0. c Insert the two parent charged pions, c ipion=1 for the pion containing the leading quark, c ipion=2 for the pion containing the leading antiquark of the same flavor: do ipion=1,2 - if(ipion.eq.2) then - idpi=-idpi - pxpi=-pxpi - pypi=-pypi - pzpi=-pzpi + if(ipion.eq.1) then + idpi=idpi1 + pxpi=pxpi1 + pypi=pypi1 + pzpi=pzpi1 + elseif(ipion.eq.2) then + idpi=-idpi1 + pxpi=-pxpi1 + pypi=-pypi1 + pzpi=-pzpi1 endif NATT=NATT+1 KATT(NATT,1)=idpi @@ -3134,8 +3367,8 @@ subroutine embedHighPt PATT(NATT,3)=pzpi PATT(NATT,4)=sqrt(pxpi**2+pypi**2+pzpi**2+pichmass**2) EATT=EATT+PATT(NATT,4) - GXAR(NATT)=xembd - GYAR(NATT)=yembd + GXAR(NATT)=xjet + GYAR(NATT)=yjet GZAR(NATT)=0. FTAR(NATT)=0. ITYPAR(NATT)=KATT(NATT,1) @@ -3145,6 +3378,62 @@ subroutine embedHighPt PEAR(NATT)=PATT(NATT,4) XMAR(NATT)=pichmass enddo +c +clin-8/2009 +c Randomly embed a number of soft pions around each high-Pt quark in pair: + if(nsembd.gt.0) then + do ipion=1,2 + do ispion=1,nsembd + idsart=3+int(3*RANART(NSEED)) + if(idsart.eq.3) then + pimass=pichmass + idpis=-211 + elseif(idsart.eq.4) then + pimass=pi0mass + idpis=111 + else + pimass=pichmass + idpis=211 + endif + NATT=NATT+1 + KATT(NATT,1)=idpis + KATT(NATT,2)=40 + KATT(NATT,3)=0 +c theta: relative angle between soft pion & associated high-Pt q or qbar, +c generate theta and phi uniformly: +c Note: it is not generated uniformly in solid angle because that gives +c a valley at theta=0, unlike the jet-like correlation (a peak at theta=0). + theta=tmaxembd*RANART(NSEED) + phi=2.*pi*RANART(NSEED) + pxspi=psembd*sin(theta)*cos(phi) + pyspi=psembd*sin(theta)*sin(phi) + pzspi=psembd*cos(theta) + if(ipion.eq.1) then + call rotate(pxpi1,pypi1,pzpi1,pxspi,pyspi,pzspi) + else + call rotate(-pxpi1,-pypi1,-pzpi1,pxspi,pyspi,pzspi) + endif +ctest off +c write(99,*) "2 ",pxspi,pyspi,pzspi + PATT(NATT,1)=pxspi + PATT(NATT,2)=pyspi + PATT(NATT,3)=pzspi + PATT(NATT,4)=sqrt(psembd**2+pimass**2) + EATT=EATT+PATT(NATT,4) + GXAR(NATT)=xjet + GYAR(NATT)=yjet + GZAR(NATT)=0. + FTAR(NATT)=0. + ITYPAR(NATT)=KATT(NATT,1) + PXAR(NATT)=PATT(NATT,1) + PYAR(NATT)=PATT(NATT,2) + PZAR(NATT)=PATT(NATT,3) + PEAR(NATT)=PATT(NATT,4) + XMAR(NATT)=pimass + enddo + enddo + endif +clin-8/2009-end c return end diff --git a/GeneratorInterface/AMPTInterface/src/main.f b/GeneratorInterface/AMPTInterface/src/main.f index 7ede4b9f567ee..5c18c380e17f2 100755 --- a/GeneratorInterface/AMPTInterface/src/main.f +++ b/GeneratorInterface/AMPTInterface/src/main.f @@ -1,7 +1,7 @@ c.....driver program for A Multi-Phase Transport model SUBROUTINE AMPT(FRAME0,BMIN,BMAX) c - double precision xmp, xmu, alpha, rscut2, cutof2 + double precision xmp, xmu, alpha, rscut2, cutof2, dshadow double precision smearp,smearh,dpcoal,drcoal,ecritl cgsfs added following line to match C++ call double precision BMIN, BMAX @@ -13,39 +13,42 @@ SUBROUTINE AMPT(FRAME0,BMIN,BMAX) cgsfs added to match specification in AMPTSET character*25 amptvn - COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 COMMON /HPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) - COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) + COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50),CARINT1(100), IAINT2(50) COMMON /AROUT/ IOUT COMMON /AREVT/ IAEVT, IARUN, MISS COMMON /smearz/smearp,smearh COMMON/RNDF77/NSEED common/anim/nevent,isoft,isflag,izpc c parton coalescence radii in case of string melting: - common /coal/dpcoal,drcoal,ecritl - common/snn/efrm,npart1,npart2 -c initialization value for parton cascade: - common /para2/ xmp, xmu, alpha, rscut2, cutof2 - common /para7/ ioscar,nsmbbbar,nsmmeson - common /para8/ idpert,npertd,idxsec - common /rndm3/ iseedp -c initialization value for hadron cascade: + common/coal/dpcoal,drcoal,ecritl + common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG +c initialization value for parton cascade: + common /para2/xmp,xmu,alpha,rscut2,cutof2 + common/para7/ioscar,nsmbbbar,nsmmeson + common/para8/idpert,npertd,idxsec + common/rndm3/iseedp +c initialization value for hadron cascade: COMMON /RUN/ NUM - common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT - COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, - & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL + common/input1/MASSPR,MASSTA,ISEED,IAVOID,DT + COMMON/INPUT2/ILAB,MANYB,NTMAX,ICOLL,INSYS,IPOT,MODE, + & IMOMEN,NFREQ,ICFLOW,ICRHO,ICOU,KPOTEN,KMUL common/oscar1/iap,izp,iat,izt common/oscar2/FRAME,amptvn common/resdcy/NSAV,iksdcy clin-6/2009: -c common/phidcy/iphidcy - common/phidcy/iphidcy,pttrig,ntrig,maxmiss - common/embed/iembed,pxqembd,pyqembd,xembd,yembd - - EXTERNAL HIDATA, PYDATA, LUDATA, ARDATA, PPBDAT, zpcbdt - SAVE +c common/phidcy/iphidcy + common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy + common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd, + 1 psembd,tmaxembd,phidecomp + common/cmsflag/dshadow,ishadow +clin-2/2012 allow random orientation of reaction plane: + common /phiHJ/iphirp,phiRP + + EXTERNAL HIDATA,PYDATA,LUDATA,ARDATA,PPBDAT,zpcbdt + SAVE c**************** FRAME=FRAME0 @@ -55,7 +58,7 @@ SUBROUTINE AMPT(FRAME0,BMIN,BMAX) K=1 100 CALL HIJING(FRAME, BMIN, BMAX) - IAINT2(1) = NATT + IAINT2(1) = NATT c evaluate Npart (from primary NN collisions) for both proj and targ: @@ -66,14 +69,14 @@ SUBROUTINE AMPT(FRAME0,BMIN,BMAX) c still repeat the event to get an interaction c (this may have an additional "trigger" effect): if(NATT.eq.0) then - imiss=imiss+1 - if(imiss.le.20) then + imiss=imiss+1 + if(imiss.le.20) then write(6,*) 'repeated event: natt=0,j,imiss=',j,imiss - goto 100 - else + goto 100 + else write(6,*) 'missed event: natt=0,j=',j goto 2000 - endif + endif endif c.....ART initialization and run CALL ARINI @@ -97,8 +100,7 @@ SUBROUTINE AMPT(FRAME0,BMIN,BMAX) END SUBROUTINE AMPTSET(EFRM0,FRAME0,PROJ0,TARG0,IAP0,IZP0,IAT0,IZT0) -c -cgsfs added following line to match C++ call +c cgsfs added following line to match C++ call double precision EFRM0 double precision xmp, xmu, alpha, rscut2, cutof2 double precision smearp,smearh,dpcoal,drcoal,ecritl @@ -106,39 +108,41 @@ SUBROUTINE AMPTSET(EFRM0,FRAME0,PROJ0,TARG0,IAP0,IZP0,IAT0,IZT0) CHARACTER FRAME*8,PROJ*8,TARG*8 character*25 amptvn COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 - COMMON /HPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50) + COMMON /HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) - COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) - COMMON /AROUT/ IOUT - COMMON /AREVT/ IAEVT, IARUN, MISS - COMMON /smearz/smearp,smearh + COMMON/ARPRNT/ARPAR1(100),IAPAR2(50),ARINT1(100),IAINT2(50) + COMMON/AROUT/IOUT + COMMON/AREVT/IAEVT,IARUN,MISS + COMMON/smearz/smearp,smearh COMMON/RNDF77/NSEED common/anim/nevent,isoft,isflag,izpc -c parton coalescence radii in case of string melting: - common /coal/dpcoal,drcoal,ecritl - common/snn/efrm,npart1,npart2 -c initialization value for parton cascade: - common /para2/ xmp, xmu, alpha, rscut2, cutof2 - common /para7/ ioscar,nsmbbbar,nsmmeson - common /para8/ idpert,npertd,idxsec - common /rndm3/ iseedp -c initialization value for hadron cascade: - COMMON /RUN/ NUM - common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT - COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, - & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL +c parton coalescence radii in case of string melting: + common/coal/dpcoal,drcoal,ecritl + common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG +c initialization value for parton cascade: + common/para2/xmp,xmu,alpha,rscut2,cutof2 + common/para7/ioscar,nsmbbbar,nsmmeson + common/para8/idpert,npertd,idxsec + common/rndm3/iseedp +c initialization value for hadron cascade: + COMMON/RUN/NUM + common/input1/MASSPR,MASSTA,ISEED,IAVOID,DT + COMMON/INPUT2/ILAB,MANYB,NTMAX,ICOLL,INSYS,IPOT,MODE, + & IMOMEN,NFREQ,ICFLOW,ICRHO,ICOU,KPOTEN,KMUL common/oscar1/iap,izp,iat,izt common/oscar2/FRAME,amptvn common/resdcy/NSAV,iksdcy clin-6/2009: c common/phidcy/iphidcy - common/phidcy/iphidcy,pttrig,ntrig,maxmiss - common/embed/iembed,pxqembd,pyqembd,xembd,yembd + common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy + common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd, + 1 psembd,tmaxembd,phidecomp common/popcorn/ipop - EXTERNAL HIDATA, PYDATA, LUDATA, ARDATA, PPBDAT, zpcbdt - SAVE + EXTERNAL HIDATA,PYDATA, LUDATA,ARDATA,PPBDAT,zpcbdt + SAVE c**************** + EFRM=EFRM0 FRAME=FRAME0 PROJ=PROJ0 @@ -155,26 +159,26 @@ SUBROUTINE AMPTSET(EFRM0,FRAME0,PROJ0,TARG0,IAP0,IZP0,IAT0,IZT0) c Trigger Pt of high-pt jets in HIJING: c HIPR1(10)=7. c - if(isoft.eq.1) then - amptvn = '1.25 (Default)' + amptvn = '1.26t5 (Default)' elseif(isoft.eq.4) then - amptvn = '2.25 (StringMelting)' + amptvn = '2.26t5 (StringMelting)' else amptvn = 'Test-Only' endif - - WRITE(*,50) amptvn + WRITE(6,50) amptvn + WRITE(12,50) amptvn 50 FORMAT(' '/ &11X,'##################################################'/1X, &10X,'# AMPT (A Multi-Phase Transport) model #'/1X, - &10X,'# Version ',a20, ' #'/1X, - &10X,'# 06/25/2009 #'/1X, + &10X,'# Version ',a25, ' #'/1X, + &10X,'# 4/06/2015 #'/1X, &10X,'##################################################'/1X, &10X,' ') -c an odd number is needed for the random number generator: - if(mod(NSEED,2).eq.0) NSEED=NSEED+1 +clin-5/2015 an odd number is needed for the random number generator: +c if(mod(NSEED,2).eq.0) NSEED=NSEED+1 + NSEED=2*NSEED+1 c 9/26/03 random number generator for f77 compiler: CALL SRAND(NSEED) c @@ -188,7 +192,8 @@ SUBROUTINE AMPTSET(EFRM0,FRAME0,PROJ0,TARG0,IAP0,IZP0,IAT0,IZT0) smearp=0d0 IAmax=max(iap,iat) smearh=1.2d0*IAmax**0.3333d0/(dble(EFRM)/2/0.938d0) -cgsfs Restored this call which was missing + nevent=NEVNT +c CALL HIJSET(EFRM, FRAME, PROJ, TARG, IAP, IZP, IAT, IZT) CALL ARTSET CALL INIZPC diff --git a/GeneratorInterface/AMPTInterface/src/zpc.f b/GeneratorInterface/AMPTInterface/src/zpc.f index f3450ebb036d9..4f8eb833c5a66 100755 --- a/GeneratorInterface/AMPTInterface/src/zpc.f +++ b/GeneratorInterface/AMPTInterface/src/zpc.f @@ -4,9 +4,6 @@ SUBROUTINE ZPCMN c Version: 1.0.1 c Author: Bin Zhang c (suggestions, problems -> bzhang@nt1.phys.columbia.edu) -cms -cms dlw & gsfs Comments our writing of output files -cms implicit double precision (a-h, o-z) clin-4/20/01 PARAMETER (NMAXGL = 16000) parameter (MAXPTN=400001) @@ -186,7 +183,6 @@ subroutine readpa cc SAVE /rndm3/ SAVE - str=str iseed=iseedp c this is the initialization file containing the initial values of c the parameters @@ -197,7 +193,7 @@ subroutine readpa c this is the final data file containing general info about the cascade cbz1/31/99 c open (6, file = 'zpc.res', status = 'unknown') -cms open (25, file = 'ana/zpc.res', status = 'unknown') +cms open (25, file = 'ana/zpc.res', status = 'unknown') cbz1/31/99end c this is the input file containing initial particle records @@ -209,8 +205,8 @@ subroutine readpa cbz1/31/99 c open (8, file = 'zpc.oscar', status = 'unknown') if(ioscar.eq.1) then -cms open (26, file = 'ana/parton.oscar', status = 'unknown') -cms open (19, file = 'ana/hadron.oscar', status = 'unknown') +cms open (26, file = 'ana/parton.oscar', status = 'unknown') +cms open (19, file = 'ana/hadron.oscar', status = 'unknown') endif cbz1/31/99end @@ -235,7 +231,7 @@ subroutine readpa ireflg=1 cbz1/31/99 IF (ireflg .EQ. 0) THEN -cms OPEN (27, FILE = 'zpc.inp', STATUS = 'UNKNOWN') +cms OPEN (27, FILE = 'zpc.inp', STATUS = 'UNKNOWN') END IF cbz1/31/99end c read (29, *) str, igeflg @@ -494,8 +490,8 @@ subroutine genei c check if it's necessary to adjust array size 'adarr' if (mul .ge. MAXPTN .or. mul .eq. 0) then - print *, 'event',ievt,'has',mul,'number of gluon', - & 'adjusting counting is necessary' +cms print *, 'event',ievt,'has',mul,'number of gluon', +cms & 'adjusting counting is necessary' stop 'adarr' end if @@ -525,7 +521,7 @@ subroutine posit2(x, y) implicit double precision (a-h, o-z) -c external ran1 +cc external ran1 common /ilist3/ size1, size2, size3, v1, v2, v3, size cc SAVE /ilist3/ @@ -895,8 +891,8 @@ subroutine inirec 1002 continue clin-6/2009 -cms if(isoft.eq.1.and.(ioscar.eq.2.or.ioscar.eq.3)) -cms 1 write(92,*) iaevt,miss,mul + if(isoft.eq.1.and.(ioscar.eq.2.or.ioscar.eq.3)) + 1 write(92,*) iaevt,miss,mul do 1003 i = 1, mul energy = e(i) @@ -935,15 +931,15 @@ subroutine inirec if(ioscar.eq.2.or.ioscar.eq.3) then if(dmax1(abs(gx(i)),abs(gy(i)), 1 abs(gz(i)),abs(ft(i))).lt.9999) then -cms write(92,200) ityp(i),px(i),py(i),pz(i),xmass(i), -cms 1 gx(i),gy(i),gz(i),ft(i) + write(92,200) ityp(i),px(i),py(i),pz(i),xmass(i), + 1 gx(i),gy(i),gz(i),ft(i) else -cms write(92,201) ityp(i),px(i),py(i),pz(i),xmass(i), -cms 1 gx(i),gy(i),gz(i),ft(i) + write(92,201) ityp(i),px(i),py(i),pz(i),xmass(i), + 1 gx(i),gy(i),gz(i),ft(i) endif endif -cyy 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2)) -cyy 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2)) + 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2)) + 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2)) c 1003 continue @@ -1216,24 +1212,24 @@ subroutine zpcrun(*) c if (mod(ictype, 2) .eq. 0) call scat(t, iscat, jscat) if (mod(ictype, 2) .eq. 0) then if(ioscar.eq.3) then -cms write(95,*) 'event,miss,iscat,jscat=',iaevt,miss,iscat,jscat + write(95,*) 'event,miss,iscat,jscat=',iaevt,miss,iscat,jscat if(dmax1(abs(gx(iscat)),abs(gy(iscat)), 1 abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)), 2 abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat))) 3 .lt.9999) then -cms write(95,200) ityp(iscat),px(iscat),py(iscat), -cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), -cms 2 gz(iscat),ft(iscat) -cms write(95,200) ityp(jscat),px(jscat),py(jscat), -cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), -cms 2 gz(jscat),ft(jscat) + write(95,200) ityp(iscat),px(iscat),py(iscat), + 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), + 2 gz(iscat),ft(iscat) + write(95,200) ityp(jscat),px(jscat),py(jscat), + 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), + 2 gz(jscat),ft(jscat) else -cms write(95,201) ityp(iscat),px(iscat),py(iscat), -cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), -cms 2 gz(iscat),ft(iscat) -cms write(95,201) ityp(jscat),px(jscat),py(jscat), -cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), -cms 2 gz(jscat),ft(jscat) + write(95,201) ityp(iscat),px(iscat),py(iscat), + 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), + 2 gz(iscat),ft(iscat) + write(95,201) ityp(jscat),px(jscat),py(jscat), + 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), + 2 gz(jscat),ft(jscat) endif endif c @@ -1244,24 +1240,24 @@ subroutine zpcrun(*) 1 abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)), 2 abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat))) 3 .lt.9999) then -cms write(95,200) ityp(iscat),px(iscat),py(iscat), -cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), -cms 2 gz(iscat),ft(iscat) -cms write(95,200) ityp(jscat),px(jscat),py(jscat), -cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), -cms 2 gz(jscat),ft(jscat) + write(95,200) ityp(iscat),px(iscat),py(iscat), + 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), + 2 gz(iscat),ft(iscat) + write(95,200) ityp(jscat),px(jscat),py(jscat), + 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), + 2 gz(jscat),ft(jscat) else -cms write(95,201) ityp(iscat),px(iscat),py(iscat), -cms 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), -cms 2 gz(iscat),ft(iscat) -cms write(95,201) ityp(jscat),px(jscat),py(jscat), -cms 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), -cms 2 gz(jscat),ft(jscat) + write(95,201) ityp(iscat),px(iscat),py(iscat), + 1 pz(iscat),xmass(iscat),gx(iscat),gy(iscat), + 2 gz(iscat),ft(iscat) + write(95,201) ityp(jscat),px(jscat),py(jscat), + 1 pz(jscat),xmass(jscat),gx(jscat),gy(jscat), + 2 gz(jscat),ft(jscat) endif endif endif -cyy format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2)) -cyy format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2)) + 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2)) + 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2)) end if @@ -1856,7 +1852,6 @@ subroutine newpos(t, i) cc SAVE /ilist5/ SAVE - t=t dt1 = ct(i) - ft(i) gx(i) = gx(i) + vx(i) * dt1 @@ -1930,7 +1925,6 @@ subroutine newmom(t) cc SAVE /frzprc/ SAVE - t=t clin-6/06/02 no momentum change for partons already frozen out, c however, spatial upgrade is needed to ensure overall system freezeout: if(isoft.eq.5) then @@ -2138,8 +2132,6 @@ subroutine getht(iscat, jscat, pp2, that) cc SAVE /rndm3/ SAVE - iscat=iscat - jscat=jscat iseed=iseedp xmu2 = (hbarc * xmu) ** 2 xmp2 = xmp ** 2 @@ -2474,10 +2466,6 @@ subroutine wallc2(i, i1, i2, i3, t, tmin) cc SAVE /ilist5/ SAVE - i1=i1 - i2=i2 - i3=i3 - t=t x1p = gx(i) x2p = gy(i) x3p = gz(i) @@ -3609,7 +3597,6 @@ subroutine fixtim(l, t, tmin1, tmin, nc) cc SAVE /ilist5/ SAVE - t=t k = nc if (tmin .lt. tmin1) then ot(l) = tmin @@ -5823,7 +5810,6 @@ subroutine chcell(il, i1, i2, i3, last0, t, tmin, nc) cc SAVE /ilist4/ SAVE - t=t if (iconfg .eq. 3 .or. iconfg .eq. 5) then jj = ichkpt do 1001 j = 1, jj @@ -6127,13 +6113,13 @@ subroutine zpca2 c & ',\n\t freezeout time=', t, c & ',\n\t ending at the ', number, 'th random number', c & ',\n\t ending collision iff=', iff -cms WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN -cms WRITE (25, *) ' number of operations = ', iopern -cms WRITE (25, *) ' number of collisions between particles = ', -cms & icolln -cms WRITE (25, *) ' freezeout time=', t -cms WRITE (25, *) ' ending at the ', number, 'th random number' -cms WRITE (25, *) ' ending collision iff=', iff + WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN + WRITE (25, *) ' number of operations = ', iopern + WRITE (25, *) ' number of collisions between particles = ', + & icolln + WRITE (25, *) ' freezeout time=', t + WRITE (25, *) ' ending at the ', number, 'th random number' + WRITE (25, *) ' ending collision iff=', iff return end @@ -6180,9 +6166,6 @@ subroutine zpca2a cc SAVE /ana2/ common /ana4/ fdetdy(24), fdndy(24), fdndpt(12) cc SAVE /ana4/ - - logical iwrite - data iwrite / .false. / SAVE do 1004 i = 1, ichkpt @@ -6229,14 +6212,13 @@ subroutine zpca2a 1004 continue do 1005 ian = 1, 12 - if ( iwrite ) then if (dn(ian) .eq. 0d0 .or. dn1(ian) .eq. 0d0 .or. & dn2(ian) .eq. 0d0) then - print *, 'event=', ievt - print *, 'dn(', ian, ')=', dn(ian), 'dn1(', ian, - & ')=', dn1(ian), 'dn2(', ian, ')=', dn2(ian) +clin-9/2012 suppress output: +c print *, 'event=', ievt +c print *, 'dn(', ian, ')=', dn(ian), 'dn1(', ian, +c & ')=', dn1(ian), 'dn2(', ian, ')=', dn2(ian) end if - endif detdy(ian) = detdy(ian) + det(ian) if (dn(ian) .ne. 0) then detdn(ian) = detdn(ian) + det(ian) / dn(ian) @@ -6314,8 +6296,8 @@ subroutine zpca2c c file header if (nff .eq. 0) then -cms write (26, 101) 'OSCAR1997A' -cms write (26, 101) 'final_id_p_x' + write (26, 101) 'OSCAR1997A' + write (26, 101) 'final_id_p_x' code = 'ZPC' versn = '1.0.1' aproj = -1 @@ -6325,8 +6307,8 @@ subroutine zpca2c reffra = 'cm' ebeam = 0d0 ntestp = 1 -cms write (26, 102) code, versn, aproj, zproj, atarg, ztarg, -cms & reffra, ebeam, ntestp + write (26, 102) code, versn, aproj, zproj, atarg, ztarg, + & reffra, ebeam, ntestp nff = 1 event = 1 bimp = 0d0 @@ -6336,22 +6318,22 @@ subroutine zpca2c c comment c event header -cms write (26, 103) event, mul, bimp, phi + write (26, 103) event, mul, bimp, phi c particles do 99 i = 1, mul -cms write (26, 104) i, ityp(i), -cms & px(i), py(i), pz(i), e(i), xmass(i), -cms & gx(i), gy(i), gz(i), ft(i) + write (26, 104) i, ityp(i), + & px(i), py(i), pz(i), e(i), xmass(i), + & gx(i), gy(i), gz(i), ft(i) 99 continue event = event + 1 -cyy 101 format (a12) -cyy 102 format (2(a8, 2x), '(',i3, ',',i6, ')+(',i3, ',', i6, ')', -cyy & 2x, a4, 2x, e10.4, 2x, i8) -cyy 103 format (i10, 2x, i10, 2x, f8.3, 2x, f8.3) -cyy 104 format (i10, 2x, i10, 2x, 9(e12.6, 2x)) + 101 format (a12) + 102 format (2(a8, 2x), '(',i3, ',',i6, ')+(',i3, ',', i6, ')', + & 2x, a4, 2x, e10.4, 2x, i8) + 103 format (i10, 2x, i10, 2x, f8.3, 2x, f8.3) + 104 format (i10, 2x, i10, 2x, 9(e12.6, 2x)) return end @@ -6415,17 +6397,17 @@ subroutine zpcou2 cc SAVE /ana3/ SAVE c -cms open (28, file = 'ana4/em.dat', status = 'unknown') + open (28, file = 'ana4/em.dat', status = 'unknown') vol = 1000.d0 * size1 * size2 * size3 ntotal = nevnt * nsbrun do 1002 ian = 1, 12 -cms write (28, *) '*** for time ', ts(ian), 'fm(s)' + write (28, *) '*** for time ', ts(ian), 'fm(s)' do 1001 i = 1, 4 -cms write (28, *) em(i, 1, ian) / vol / ntotal, -cms & em(i, 2, ian) / vol / ntotal, -cms & em(i, 3, ian) / vol / ntotal, -cms & em(i, 4, ian) / vol / ntotal + write (28, *) em(i, 1, ian) / vol / ntotal, + & em(i, 2, ian) / vol / ntotal, + & em(i, 3, ian) / vol / ntotal, + & em(i, 4, ian) / vol / ntotal 1001 continue 1002 continue @@ -6621,25 +6603,25 @@ subroutine zprota(xn1, xn2, xn3, theta, v1, v2, v3) end c double precision function ran1(idum) - +c c* return a uniform random deviate between 0.0 and 1.0. set idum to c* any negative value to initialize or reinitialize the sequence. - +c c implicit double precision (a-h, o-z) - +c c dimension r(97) - +c c common /rndm1/ number -cc SAVE /rndm1/ +ccc SAVE /rndm1/ c parameter (m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 1d0 / m1) c parameter (m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 1d0 / m2) c parameter (m3 = 243000, ia3 = 4561, ic3 = 51349) -clin-6/23/00 save ix1-3: -clin-10/30/02 r unsaved, causing wrong values for ran1 when compiled with f77: -cc SAVE ix1,ix2,ix3,r +cclin-6/23/00 save ix1-3: +cclin-10/30/02 r unsaved, causing wrong values for ran1 when compiled with f77: +ccc SAVE ix1,ix2,ix3,r c SAVE c data iff/0/ - +c c if (idum .lt. 0 .or. iff .eq. 0) then c iff = 1 c ix1 = mod(ic1 - idum, m1) @@ -6657,17 +6639,17 @@ subroutine zprota(xn1, xn2, xn3, theta, v1, v2, v3) c ix1 = mod(ia1 * ix1 + ic1, m1) c ix2 = mod(ia2 * ix2 + ic2, m2) c ix3 = mod(ia3 * ix3 + ic3, m3) -clin-7/01/02 j = 1 + (97 * i x 3) / m3 +cclin-7/01/02 j = 1 + (97 * i x 3) / m3 c j=1+(97*ix3)/m3 -clin-4/2008: -c if (j .gt. 97 .or. j .lt. 1) pause +cclin-4/2008: +cc if (j .gt. 97 .or. j .lt. 1) pause c if (j .gt. 97 .or. j .lt. 1) print *, 'In zpc ran1, j<1 or j>97',j c ran1 = r(j) c r(j) = (dble(ix1) + dble(ix2) * rm2) * rm1 - -clin-6/23/00 check random number generator: +c +cclin-6/23/00 check random number generator: c number = number + 1 -c if(number.le.100000) write(99,*) 'number, ran1=', number,ran1 - +cc if(number.le.100000) write(99,*) 'number, ran1=', number,ran1 +c c return c end