v3.0.02 update (#50)

Author: ArcsaberHep

PACIAE.sh

@@ -2,7 +2,7 @@
      Written by Markdown language.
                  By Anke at CCNU on 17/10/2022 
                                                
-                    Last updated on 13/03/2024 
+                    Last updated on 15/10/2024 
  -->
 
 # The parton and hadron cascade model PACIAE
@@ -122,7 +122,20 @@ The released code are available on https://github.com/ArcsaberHep/PACIAE/release
 ## Update Notes:
 
 <!----------------------------------------------------------------------------->
-...(waiting for updating)
+<!-- ...(waiting for updating) -->
+
+<!----------------------------------------------------------------------------->
+### <font color=red> 04-10/2024: </font> In version PACIAE 3.0.02
+
+- In "main_30.f" and "parini_30.f", the *NN* total cross section "para1_1" has been recalculated in the program by subroutine "PAXTOT" via Donnachie-Landshoff parameterization instead of the manual input.
+- In "main_30.f" and "parini_30.f", the adjustment of energy conservation stemming from diquark breaking-up and recovering has been cancelled. It may cause errors in string fragmentation and abnormal hardening of intermediate and high $p_T$. The cancellation may lead to 0.1-1% energy violation.
+- In "parini_30.f", a new "D-framework" has been introduced to take the leading-proton effect into account, mainly for ***p*A** collisions.
+- In "parini_30.f", for **AA** and ***p*A** collisions, now the initial positions of projectile and target nucleon/nuclei are setting at (b/2, 0, -20) and (-b/2, 0, 20), respectively.
+- In "parcas_30.f", the analytical form of parton-parton cross sections have been introduced replacing the old numerical method of integral cross section calculation, which may cause abnormal sharp peaks at forward and backward rapidities. The mass correction of cross section for heavy c and b quarks has been taken into account now.
+- An experimental feature, string fragmentation after the gluon-splitting and the deexcitation of excited quarks, has been introduced in order to enhance the flow effect in string fragmentation model.
+- Some statements optimized in program.
+- Bug fixed.
+- "PACIAE.sh" updated.
 
 <!----------------------------------------------------------------------------->
 ### <font color=red> 01-03/2024: </font> In version PACIAE 3.0.01

README.md

@@ -30,8 +30,11 @@ subroutine coales(ijk,neve,nnout,nap,nat,nzp,nzt,i_coal)
         common/sa6_p/ithroq_p,ithrob_p,ich_p,non6_p,throe_p(4)   ! 201104 300623 Lei
         common/sa18/i_deex,n_deex_step,i_pT,i_pT_max,a_FF,aPS_c,aPS_b   ! 280823 Lei
         common/sa24/adj1(40),nnstop,non24,zstop
+        common/sa28/nstr,nstra(kszj),nstrv(kszj),nstr0,
+     &   nstr1,nstr1a(kszj),nstr1v(kszj)   ! 030620
         common/sa36/nglu,nongu,kglu(kszj,5),pglu(kszj,5),vglu(kszj,5) ! 220822
         common/sa37/nth,npadth,kth(kszj,5),pth(kszj,5),vth(kszj,5)   ! 150922
+        common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
         common/sbh/nbh,nonh,kbh(kszj,5),pbh(kszj,5),vbh(kszj,5)
         common/syspar/ipden,itden,suppm,suptm,suppc,suptc,r0p,r0t,
      c   napp,natt,nzpp,nztt,pio
@@ -58,19 +61,21 @@ subroutine coales(ijk,neve,nnout,nap,nat,nzp,nzt,i_coal)
 
 c-------------------------------------------------------------------------------
 c----------------------------   Junctions removing   ---------------------------
+        if( INT(adj1(40)).eq.3 )then   ! 070223
 c220822     Remove junctions.
-        jb = 0
-2010    do i1=jb+1,N,1  ! i1 loop
-            kf   = K(i1,2)
-            kfab = ABS(kf)
-            if(kfab.ne.88)then
-                jb = jb + 1
-                goto 2020
-            endif
-            call updad_pyj(N,i1+1,1)   ! 090922 'updad_pyj' in sfm_30.f
-            N = N - 1
-            goto 2010
-2020    enddo   ! i1 loop
+            jb = 0
+2010        do i1=jb+1,N,1  ! i1 loop
+                kf   = K(i1,2)
+                kfab = ABS(kf)
+                if(kfab.ne.88)then
+                    jb = jb + 1
+                    goto 2020
+                endif
+                call updad_pyj(N,i1+1,1)   ! 090922 'updad_pyj' in sfm_30.f
+                N = N - 1
+                goto 2010
+2020        enddo   ! i1 loop
+        endif   ! 070223
 c----------------------------   Junctions removing   ---------------------------
 c-------------------------------------------------------------------------------
 
@@ -83,6 +88,7 @@ subroutine coales(ijk,neve,nnout,nap,nat,nzp,nzt,i_coal)
 c-------------------------------------------------------------------------------
 c-----------------------------   Gluon splitting   -----------------------------
 c220122
+        n00 = N   ! Original total entries in PYJETS
 c       Move gluons from 'pyjest' to 'sa36'.
         call remo_glu
 c       Break-up gluon (with E_g>2E_u in 'sa36') -> qqbar string 
@@ -100,20 +106,15 @@ subroutine coales(ijk,neve,nnout,nap,nat,nzp,nzt,i_coal)
 c-------------------------------------------------------------------------------
 
 
-c250823 Debug mode.   ! 250823 Lei
-c       Do g-splitting only, without q-deexcitation.
-        if( INT(adj1(12)).eq.3 ) goto 1000
-
-
 c-------------------------------------------------------------------------------
 c---------------------------   Quark deexcitation   ----------------------------
 c280822 energetic q (qbar) de-excitation
-        n00   = N   ! Original total entries in PYJETS
         i_call_deex  = 0
         i_daught_gen = 1   ! the #-th newly produced daughter qqbar
         n_deex = 0   ! the number of successful deexcitation
         jb = 0
         n0 = N   ! Current total entries in PYJETS
+        if( i_deex_gen.eq.0 ) goto 900   ! 300324 Lei
 700     continue
         do i1=jb+1,n0,1
             kf0   = K(i1,2)
@@ -147,12 +148,37 @@ subroutine coales(ijk,neve,nnout,nap,nat,nzp,nzt,i_coal)
 c300623 Shares the 4-momentum in 'throe_p' among partons.   ! 300623 Lei
         call share_p_PYJETS   ! 300623 Lei
 c220122
+
+c300324 Lei
+        if( INT(adj1(12)).eq.3 )then
+c       Records the location (line numbers) of new strings (qqbar).
+            do i1=n00+1,N,2
+                nstr1 = nstr1 + 1
+                nstr1a( nstr1 ) = i1
+                nstr1v( nstr1 ) = i1 + 1
+            end do
+            nstr0 = nstr1
+c       Updates "sbe".
+            do ii=n00+1,N,1
+                nbe = nbe + 1
+                do jj=1,5,1
+                    kbe( nbe, jj ) = K(ii,jj)
+                    pbe( nbe, jj ) = P(ii,jj)
+                    vbe( nbe, jj ) = V(ii,jj)
+                end do
+            end do
+        end if
+c300324 Lei
+
 c---------------------------   Quark deexcitation   ----------------------------
 c-------------------------------------------------------------------------------
 
 
 c       Just do the g-splitting and quark deexcitation, without real coalescence
-1000    if( i_coal.eq.0 ) return   ! 300623 Lei For adj12 = 2
+c250324
+1000    continue
+        if( i_coal.eq.0 ) return   ! 300623 Lei For adj12 = 2
+c250324
 
 
 c-------------------------------------------------------------------------------
@@ -506,7 +532,7 @@ subroutine hadpro(rrp,iphas)   ! 080324
         common/sa6_p/ithroq_p,ithrob_p,ich_p,non6_p,throe_p(4)   ! 201104 300623 Lei
         common/sa24/adj1(40),nnstop,non24,zstop
         common/sa37/nth,npadth,kth(kszj,5),pth(kszj,5),vth(kszj,5)   ! 150922
-        common/coal1/bmrat   ! ratio of baryon to meson
+        common/coal1/bmrat,i_mm   ! ratio of baryon to meson
         dimension pc(4),rc(3),iar(3),rcp(3)
         dimension psu(3),peo(5),pnn(kszj,5)
         dimension numb(3)   ! 110905
@@ -575,7 +601,7 @@ subroutine hadpro(rrp,iphas)   ! 080324
 c110324 Lei
             KF_in_1 = kf1
             KF_in_2 = kf2
-c           Exchanges KFSs of qbar and q to ensure the first one is q.
+c       Exchanges KFs of qbar and q to ensure the first one is q.
             if( kf1.lt.0 )then
                 KF_in_1 = kf2
                 KF_in_2 = kf1
@@ -584,13 +610,13 @@ subroutine hadpro(rrp,iphas)   ! 080324
 c110324 Lei
             if(isucc.eq.0) goto 500
 
-c           Phase space adjudgment.
+c       Phase space adjudgment.
             if( iphas.ne.0 )then 
             call phas(i1,i2,0,isucc,2,iphas)
             if( isucc.eq.0 ) goto 500   ! fail
             endif
 
-c           Proceed for success
+c       Proceed for success
             imes = imes + 1
             nme  = nme  + 1
 
@@ -659,6 +685,12 @@ subroutine hadpro(rrp,iphas)   ! 080324
                 call findb(kf1,kf2,kf3,cm,kfii,amasi,isucc,1)
                 if(isucc.eq.0) goto 600
 
+c       Phase space adjudgment.
+                if( iphas.ne.0 )then 
+                call phas(i1,i2,i3,isucc,3,iphas)
+                if( isucc.eq.0 ) goto 600   ! fail
+                endif
+
 c       Proceed for success.
                 ibarp = ibarp + 1
                 nba   = nba   + 1
@@ -735,6 +767,12 @@ subroutine hadpro(rrp,iphas)   ! 080324
             call findb(-kf1,-kf2,-kf3,cm,kfii,amasi,isucc,-1)
             if(isucc.eq.0) goto 700   ! 110324 Lei
 
+c       Phase space adjudgment.
+            if( iphas.ne.0 )then 
+            call phas(i1,i2,i3,isucc,3,iphas)
+            if( isucc.eq.0 ) goto 700   ! fail
+            endif
+
             ibarm = ibarm + 1
             nba = nba + 1
 
@@ -2424,13 +2462,21 @@ subroutine tdgaus(v,pmax,np,pp)
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         subroutine remo_glu   ! 160822
 c       moves gluons from 'pyjets' to 'sa36'
+c290324 Removes gluons from 'sbe', too.   ! 290324 Lei
         IMPLICIT DOUBLE PRECISION(A-H, O-Z)
         IMPLICIT INTEGER(I-N)
         INTEGER PYK,PYCHGE,PYCOMP
         PARAMETER (KSZJ=80000)
         COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
-        common/sa36/nglu,nongu,kglu(kszj,5),pglu(kszj,5),vglu(kszj,5)
         common/sa1/kjp21,non1,bp,iii,neve,nout,nosc
+        common/sa24/adj1(40),nnstop,non24,zstop   ! 170205
+        common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio   ! 220110
+        common/sa28/nstr,nstra(kszj),nstrv(kszj),nstr0,
+     &   nstr1,nstr1a(kszj),nstr1v(kszj)   ! 030620
+        common/sa36/nglu,nongu,kglu(kszj,5),pglu(kszj,5),vglu(kszj,5)
+        common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
+
+
         nglu=0
         if(iii.eq.1)then   ! 300623 Lei
         kglu=0
@@ -2466,9 +2512,62 @@ subroutine remo_glu   ! 160822
         enddo
         enddo
         N=N-1
+c290324 Lei
+        if( INT(adj1(12)).eq.3 )then
+c       Adjusts line numbers of components of broken diquarks.
+        do i_diq = 1, idi, 1
+            if( i1.lt.ifcom(i_diq) ) ifcom(i_diq) = ifcom(i_diq) - 1
+            if( i1.lt.  npt(i_diq) )   npt(i_diq) =   npt(i_diq) - 1
+        end do
+        do j = i1+1, N+1, 1
+            ndiq( j-1 ) = ndiq(j)
+        end do
+c       Adjusts line numbers of the string-locating.
+        do i_string = 1, nstr1, 1
+            if( i1.le.nstr1a(i_string) ) 
+     &          nstr1a(i_string) = nstr1a(i_string) - 1
+            if( i1.le.nstr1v(i_string) ) 
+     &          nstr1v(i_string) = nstr1v(i_string) - 1
+        end do
+        end if
+c290324 Lei
         goto 201
 202     enddo   ! do loop
 203     continue
+
+
+c290324 Lei
+        if( INT(adj1(12)).eq.3 )then
+c       Removes gluons from 'sbe'.
+        jb=0
+301     do i1=jb+1,nbe
+        kf=kbe(i1,2)
+        kfab=iabs(kf)
+        eng=pbe(i1,4)
+        if(kfab.ne.21)then   ! stay
+        jb=jb+1
+        goto 302
+        endif
+        if(i1.eq.nbe)then
+        nbe=nbe-1
+        goto 303
+        endif
+c       move particle list one step downward from i1+1 to nbe
+        do jj=1,5
+        do j=i1+1,nbe,1
+        kbe(j-1,jj)=kbe(j,jj)
+        pbe(j-1,jj)=pbe(j,jj)
+        vbe(j-1,jj)=vbe(j,jj)
+        enddo
+        enddo
+        nbe=nbe-1
+        goto 301
+302     enddo   ! do loop
+303     continue
+        end if
+c290324 Lei
+
+
         return
         end
 

coales_30.f

@@ -145,9 +145,7 @@ subroutine hadcas(ijk,neve,nout,time_had,ijkk)
 c        neutron, ...
 c       initial particle list is compsed of the arraies in common block 
 c        'sa1_h', 'tau' and 'ishp' in 'sa8_h', and 'numb' in 'sa9_h'
-c       call prt_sa1_h(nsa) ! sa 
         call filt_h
-c       call prt_sa1_h(nsa) ! sa
 c280910 time=0.
 
 c       calculate position of center of mass of the system. distance of a 
@@ -1107,7 +1105,7 @@ subroutine coinel(l,l1,ss,b,pi,pj,icp,pii,pjj,lc,tc,tw,winel,
 c       Total 112 channels: 385-496
 c--------------------------------------------
 c       D + pion/rho
-c       Single production channel.
+c       Single production channels.
 c
 c       Charge ++, --
 c       D+- + pion+-/rho+-
@@ -1132,7 +1130,7 @@ subroutine coinel(l,l1,ss,b,pi,pj,icp,pii,pjj,lc,tc,tw,winel,
 c       400. D*bar0 + rho+  -->  Dbar0 + pi+
 c--------------------------------------------
 c       D + pion/rho
-c       Dual production chennels.
+c       Dual production channels.
 c
 c       Charge 00 
 c       D0 + pi0/rho0
@@ -1248,7 +1246,7 @@ subroutine coinel(l,l1,ss,b,pi,pj,icp,pii,pjj,lc,tc,tw,winel,
 c============================================
 c031223 Lei D meson induced reactions.
 
-c       489-592: null.
+c       497-592: null.
 
 c       593. lambda- + p to K*+ + Omega
 c       594. lambda- + n to K*0 + Omega
@@ -10896,7 +10894,7 @@ subroutine DMeson_coll( iA,iB,KF_A,KF_B,E_AB,icp,lc,tc,tw,ioo )
         iCharge_D = PYCHGE( KF_D )
         iCharge_M = PYCHGE( KF_M )
 
-c       Single production chennel: 16 chennels, 385-400.
+c       Single production channels: 16 channels, 385-400.
         if( (iCharge_D.gt.0 .AND. iCharge_M.gt.0) .OR.   ! ++
      &      (iCharge_D.lt.0 .AND. iCharge_M.lt.0) .OR.   ! --
      &      (iCharge_D.eq.0 .AND. iCharge_M.gt.0 .AND.   ! 0+, KF-+
@@ -10905,7 +10903,7 @@ subroutine DMeson_coll( iA,iB,KF_A,KF_B,E_AB,icp,lc,tc,tw,ioo )
      &            KF_D.gt.0 .AND.      KF_M.lt.0)        ! 0-, KF+-
      &    )then
             call Dmeson1(iA,iB,KF_D,KF_M,E_AB,icp,lc,tc,tw,ioo,i_order)
-c       Dual production chennels: 64 channels, 401-464,
+c       Dual production channels: 64 channels, 401-464,
 c        else if( (iCharge_D.eq.0 .AND. iCharge_M.eq.0) .OR.   ! 00
 c     &           (iCharge_D.gt.0 .AND. iCharge_M.lt.0) .OR.   ! +-
 c     &           (iCharge_D.lt.0 .AND. iCharge_M.gt.0) .OR.   ! -+
@@ -11936,8 +11934,12 @@ subroutine upddep(ij,lc,tc,tw,time,iii)
         i2=i1
         goto 3001
 3000    enddo
+c200624 Lei
+3001    continue
         if(i2.gt.0)then   ! 280224 Lei
-3001    do i1=i2,kfmax
+c3001   do i1=i2,kfmax
+        do i1=i2,kfmax
+c200624 Lei
         numb(i1)=numb(i1)-1
         enddo
         end if   ! 280224 Lei