IAP GITLAB

Commit dc345649 authored by Ralf Ulrich's avatar Ralf Ulrich

sibyll2.3c03

git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@7392 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 81e47526
......@@ -7,7 +7,7 @@ C SSSSSS IIIIIII BBBBB YY LLLLLLL LLLLLLL
C=======================================================================
C Code for SIBYLL: hadronic interaction Monte Carlo event generator
C=======================================================================
C Version 2.3c01 (Jun-01-2017, modified Sept-05-2017)
C Version 2.3c03 (Jun-01-2017, modified Aug-22-2019)
C
C with CHARM production
C
......@@ -31,13 +31,13 @@ C paolo.lipari@roma1.infn.it
C friehn@lip.pt
C stanev@bartol.udel.edu
C
C differences to Sibyll 2.3 include:
C last changes relative to Sibyll 2.3c:
C * added cross section tables for hadron-nitrogen and hadron-oxygen
C (changed S_CCSIG common)
C * no remnant in high mass diff. events (pi0-had scattering)
C * repaired had-nuc. cross section routine for kaon beams
C routine remains inactive in ordinary calls.
C
C *extend eta' and phi decay to include prompt muons
C *allow initialization from file
C *explicit hyperon production in 2-particle fireballs
C *forbid baryon pair production immediately next to leading baryon
C *adjusted string tension to help with central prod. in CMS
C=======================================================================
SUBROUTINE SIBYLL (K_beam, IATARG, Ecm)
......@@ -427,7 +427,7 @@ C. important information for the generation of events
C.
C PARAMETER (NS_max = 20, NH_max = 80)
C COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
C & SSIGN(61,3),SSIGNSD(61,3) ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
C & SSIGN(61,3,3),SSIGNSD(61,3,3) ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
C.
C. NSQS = number of energy points (61 is current version)
C. ASQSMIN = log_10 [sqrt(s) GeV] minimum value
......@@ -438,8 +438,10 @@ C.
C. SSIG(J,1) inelastic cross section for pp interaction
C. at energy: sqrt(s)(GeV) = 10**[ASQSMIN+DASQS*(J-1)]
C. SSIG(J,2) inelastic cross section for pi-p interaction
C. SSIGN(J,1) inelastic cross section for p-Air interaction
C. SSIGN(J,2) inelastic cross section for pi-Air interaction
C. SSIGN(J,1,1) inelastic cross section for p-Air interaction
C. SSIGN(J,2,1) inelastic cross section for pi-Air interaction
C. SSIGN(J,1,2) inelastic cross section for p-Nitrogen interaction
C. SSIGN(J,2,2) inelastic cross section for pi-Nitrogen interaction
C.
C. PJETC(n_s,n_j,J,1) Cumulative probability distribution
C. for the production of n_s soft interactions and
......@@ -469,7 +471,7 @@ C-----------------------------------------------------------------------
* /,' ','| F. RIEHN et al., Proc. 35th Int. Cosmic Ray Conf.|',
* /,' ','| Bexco, Busan, Korea, cont. 301 (2017) |',
* /,' ','| |',
* /,' ','| last modifications: F. Riehn (09/18/2017) |',
* /,' ','| last modifications: F. Riehn (08/22/2019) |',
* /,' ','====================================================',
* /)
......@@ -1434,26 +1436,33 @@ C----------------------------------------------------------------
DOUBLE PRECISION AM,AM2
COMMON /S_MASS1/ AM(99), AM2(99)
INTEGER ICHP,ISTR,IBAR
COMMON /S_CHP/ ICHP(99), ISTR(99), IBAR(99)
INTEGER ICHM
COMMON /S_CHM/ ICHM(99)
CHARACTER*6 NAMP
COMMON /S_CNAM/ NAMP (0:99)
SAVE
WRITE(LUN,50)
50 FORMAT(/,2X,16X,'SIBYLL PARTICLE TABLE:',/,2x,56('-'))
50 FORMAT(/,2X,16X,'SIBYLL PARTICLE TABLE:',/,2x,80('-'))
WRITE(LUN,100)
100 FORMAT(2X,'Particle',4X,'SIB PID',6x,'SIB2PDG',6x,'SIB2PDG^-1',
& 4x,'MASS'/, 2X,56('-'))
& 4x,'MASS',4x,'STRG',4x,'CHRM',4x,'BRYN'/, 2X,80('-'))
DO J=1,99
IA = ISIB_PID2PDG( j )
IF(IA.ne.0)THEN
ISIBPDG2PIDIA=ISIB_PDG2PID( IA )
ELSE
WRITE(LUN,*) 'PDG conversion not found!'
ENDIF
WRITE (LUN,120) NAMP(J), J, IA, ISIBPDG2PIDIA, AM(J)
WRITE(LUN,'(1X,A,I2)') 'PDG conversion not found! pid=', j
ENDIF
WRITE (LUN,120) NAMP(J), J, IA, ISIBPDG2PIDIA, AM(J), ISTR(J),
& ICHM(J), IBAR(J)
ENDDO
120 FORMAT(4X,A6,4X,I4,7X,I7,8X,I4,5X,F9.3)
120 FORMAT(4X,A6,4X,I4,7X,I7,8X,I4,5X,F9.3,3(6X,I2))
END
......@@ -4446,8 +4455,11 @@ C--------------------------------------------------------------------
DOUBLE PRECISION FACN
DIMENSION FACN(3:10)
COMMON /SIB_FAC/ FACN
DOUBLE PRECISION EPS100
SAVE
DATA EPS100 /1.D-100/
DO J=0,NH_max
DO I=0,NS_max
P_int(I,J) = 0.D0
......@@ -4528,16 +4540,20 @@ C--------------------------------------------------------------------
chi2_hard_12 = (1.D0-al1+ga1)*(1.D0-al2-ga2)*chi2_hard
chi2_hard_21 = (1.D0-al1-ga1)*(1.D0-al2+ga2)*chi2_hard
ef_11 = max(-0.5D0*(chi2_soft_11+chi2_hard_11),log(EPS100))
ef_22 = max(-0.5D0*(chi2_soft_22+chi2_hard_22),log(EPS100))
ef_12 = max(-0.5D0*(chi2_soft_12+chi2_hard_12),log(EPS100))
ef_21 = max(-0.5D0*(chi2_soft_21+chi2_hard_21),log(EPS100))
ef_11 = dexp(ef_11)
ef_22 = dexp(ef_22)
ef_12 = dexp(ef_12)
ef_21 = dexp(ef_21)
ef_11 = dexp(-0.5D0*(chi2_soft_11+chi2_hard_11))
ef_22 = dexp(-0.5D0*(chi2_soft_22+chi2_hard_22))
ef_12 = dexp(-0.5D0*(chi2_soft_12+chi2_hard_12))
ef_21 = dexp(-0.5D0*(chi2_soft_21+chi2_hard_21))
esf_11 = ef_11**2
esf_22 = ef_22**2
esf_12 = ef_12**2
esf_21 = ef_21**2
esf_11 = max(ef_11,EPS100)**2
esf_22 = max(ef_22,EPS100)**2
esf_12 = max(ef_12,EPS100)**2
esf_21 = max(ef_21,EPS100)**2
F_ine = B*(1.D0 - fe_11*esf_11 - fe_12*esf_12
& - fe_21*esf_21 - fe_22*esf_22)
......@@ -4571,24 +4587,24 @@ C--------------------------------------------------------------------
soft_rec_22 = 1.D0/chi2_soft_22
soft_rec_12 = 1.D0/chi2_soft_12
soft_rec_21 = 1.D0/chi2_soft_21
chi2_hard_11 = max(chi2_hard_11,EPS10)
chi2_hard_22 = max(chi2_hard_22,EPS10)
chi2_hard_12 = max(chi2_hard_12,EPS10)
chi2_hard_21 = max(chi2_hard_21,EPS10)
chi2_hard_11 = max(chi2_hard_11,EPS100)
chi2_hard_22 = max(chi2_hard_22,EPS100)
chi2_hard_12 = max(chi2_hard_12,EPS100)
chi2_hard_21 = max(chi2_hard_21,EPS100)
DO I=0,I0MAX
soft_rec_11 = soft_rec_11*chi2_soft_11
soft_rec_22 = soft_rec_22*chi2_soft_22
soft_rec_12 = soft_rec_12*chi2_soft_12
soft_rec_21 = soft_rec_21*chi2_soft_21
hard_rec_11 = 1.D0/chi2_hard_11
hard_rec_22 = 1.D0/chi2_hard_22
hard_rec_12 = 1.D0/chi2_hard_12
hard_rec_21 = 1.D0/chi2_hard_21
soft_rec_11 = max(soft_rec_11*chi2_soft_11,EPS100)
soft_rec_22 = max(soft_rec_22*chi2_soft_22,EPS100)
soft_rec_12 = max(soft_rec_12*chi2_soft_12,EPS100)
soft_rec_21 = max(soft_rec_21*chi2_soft_21,EPS100)
hard_rec_11 = max(1.D0/chi2_hard_11,EPS100)
hard_rec_22 = max(1.D0/chi2_hard_22,EPS100)
hard_rec_12 = max(1.D0/chi2_hard_12,EPS100)
hard_rec_21 = max(1.D0/chi2_hard_21,EPS100)
DO J=0,J0MAX
hard_rec_11 = hard_rec_11*chi2_hard_11
hard_rec_22 = hard_rec_22*chi2_hard_22
hard_rec_12 = hard_rec_12*chi2_hard_12
hard_rec_21 = hard_rec_21*chi2_hard_21
hard_rec_11 = max(hard_rec_11*chi2_hard_11,EPS100)
hard_rec_22 = max(hard_rec_22*chi2_hard_22,EPS100)
hard_rec_12 = max(hard_rec_12*chi2_hard_12,EPS100)
hard_rec_21 = max(hard_rec_21*chi2_hard_21,EPS100)
P_int(I,J) = P_int(I,J)
& + fe_11*soft_rec_11*hard_rec_11*fac_11
& + fe_22*soft_rec_22*hard_rec_22*fac_22
......@@ -5511,11 +5527,14 @@ c external type declarations
c local type declarations
DOUBLE PRECISION SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO,
& SIGPROD,SIGBDIF,S_RNDM,S,PF,PB,PD,P0,P1,P2,R
& SIGPROD,SIGBDIF,SIGELA,S_RNDM,S,PF,PB,PD,P0,P1,P2,R
DIMENSION SIGDIF(3)
INTEGER K
SAVE
IF(NDEBUG.gt.0)
&WRITE(LUN,*)'SIB_START_EV:', SQS, L, IA, IAFLG, NW, JDIF
C...sample number of wounded nucleons
c read hadron-nucleon cross section from table
CALL SIB_SIGMA_HP(L,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
......@@ -5525,8 +5544,8 @@ c read hadron-nucleon cross section from table
IF(IPAR(12).eq.3)THEN
c distinguish between nuclear cross sections..
IF(IAFLG.eq.0)THEN
c if target is nucleus calc. hadron-nucleus cross section (slow)
CALL SIB_SIGMA_HNUC(L,IA,SQS,SIGprod,SIGbdif)
c if target is nucleus calc. hadron-nucleus cross section
CALL SIB_SIGMA_HNUC(L,IA,SQS,SIGprod,SIGbdif,SIGela)
ELSE
c if target is air read hadron-air cross section from table
CALL SIB_SIGMA_HAIR(L,SQS,SIGprod,SIGbdif)
......@@ -10878,12 +10897,12 @@ c COMMON /S_DEBUG/ Ncall, Ndebug, Lun
INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION STR_mass_val, STR_mass_val_hyp, STR_mass_sea
COMMON /S_CUTOFF/ STR_mass_val, STR_mass_val_hyp, STR_mass_sea
......@@ -10985,12 +11004,12 @@ C-----------------------------------------------------------------------
INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B,
& SSIG_RHO
COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
......@@ -11179,12 +11198,12 @@ C-----------------------------------------------------------------------
INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DIMENSION PJ(2),PS(2),PW(2)
SAVE
......@@ -11207,7 +11226,7 @@ C-----------------------------------------------------------------------
DO K=1,2
PW(K) = ATARGET*SSIG(J,K)/SSIGN(J,K)
PW(K) = ATARGET*SSIG(J,K)/SSIGN(J,K,1)
PJ(K) = 0.D0
PS(K) = 0.D0
......@@ -11227,8 +11246,8 @@ C-----------------------------------------------------------------------
ENDDO
WRITE(LUN,20) SQS,SSIG(J,1),SSIGN(J,1),PS(1),PJ(1),PW(1)
& ,SSIG(J,2),SSIGN(J,2),PS(2),PJ(2),PW(2)
WRITE(LUN,20) SQS,SSIG(J,1),SSIGN(J,1,1),PS(1),PJ(1),PW(1)
& ,SSIG(J,2),SSIGN(J,2,1),PS(2),PJ(2),PW(2)
ENDDO
......@@ -11241,10 +11260,11 @@ C-----------------------------------------------------------------------
C=======================================================================
SUBROUTINE SIG_AIR_INI
SUBROUTINE SIG_AIR_INI
C-----------------------------------------------------------------------
C...Initialize the cross section and interaction lengths on air
C... Initialize the cross section and interaction lengths on air,
C. nitrogen and oxygen
C. (this version initializes p-air, pi-air, and K-air cross sections)
C.
C. also calculates the low mass beam diffraction cross section in hAir \FR
......@@ -11256,12 +11276,12 @@ C-----------------------------------------------------------------------
INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
COMMON /GLAUB_SCR/ XI_MAX , ALAM(61)
INTEGER NCALL, NDEBUG, LUN
......@@ -11272,15 +11292,15 @@ C-----------------------------------------------------------------------
INTEGER IPAR
COMMON /S_CFLAFR/ PAR(NPAR_max), IPAR(NIPAR_max)
DIMENSION SIGDIF(3)
DIMENSION ITARGC(3)
CHARACTER*3 TARGN
DIMENSION TARGN(3)
SAVE
DATA AVOG /6.0221367D-04/
DATA ATARGET /14.514D0/
DATA ITARGC /0,14,16/
DATA TARGN /'air','nit','oxy'/
c PRINT *,'inel. screening in hadron - nucleus interactions'
c ALAM = 0.5
c PRINT *,'const. coupling: ', ALAM
IF ( IPAR(12).GT.0 ) THEN
if (ndebug.gt.0) then
WRITE(LUN,*) ' SIG_AIR_INI:'
......@@ -11290,60 +11310,79 @@ c PRINT *,'const. coupling: ', ALAM
if (ndebug.gt.0)WRITE(LUN,*)' low mass Xi_max: ' , XI_MAX
ENDIF
C...target loop (air, N, O)
DO IK=1,3
IAT = ITARGC(IK)
WRITE(6,*) 'SIG_AIR_INI: initializing target: (i,A)',
& ik, IAT, TARGN(IK) , '..'
C...particle loop (p, pi, K)
DO K=1,3
DO K=1,3
if (NDEBUG .gt. 0 ) then
WRITE(LUN,'(/,1X,A,A)')
& 'Table: J, sqs, SIGtot, SIGprod, SIG_SD,',
& ' Lambda '
WRITE(LUN,*)
& '-------------------------------------------------',
& '-------------'
endif
DO J=1,NSQS
if (NDEBUG .gt. 0 ) then
WRITE(6,'(/,1X,A,A)')
& 'Table: J, IK, sqs, SIGtot, SIGprod, SIG_SD,',
&' Lambda '
WRITE(6,*)
& '-------------------------------------------------',
& '-------------'
endif
DO J=1,NSQS
ASQS = ASQSMIN + DASQS*DBLE(J-1)
SQS = 10.D0**ASQS
ASQS = ASQSMIN + DASQS*DBLE(J-1)
SQS = 10.D0**ASQS
IF (K.EQ.1) THEN
IF (K.EQ.1) THEN
c Goulianos param. from GAP-2012-056, Mx**2s = 0.02
c against PDG elastic cross section
CALL SIB_HADCS1(K,SQS,SIGT1,SIGEL1,SIGINEL1,SLOPE1,RHO1)
SIGEFF = 0.68D0*(1.D0+36.D0/SQS**2)*
& dlog(0.6D0+XI_MAX/1.5D0*SQS**2)
ALAM(J) = dSQRT(SIGEFF/SIGEL1)
ENDIF
c CALL SIB_HADCSL(k,SQS,
c & SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
CALL SIB_SIGMA_HP(K,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
CALL SIG_H_AIR
& (SIGT, SLOPE, RHO, ALAM(J),
& SSIGT, SSIGEL, SSIGQE, SSIGSD, SSIGQSD)
if (ndebug .gt. 0 ) WRITE(LUN,'(1X,I2,1P,5E12.3)')
& K,SQS,SSIGT,SSIGT-SSIGQE,SSIGQSD,ALAM(J)
C particle production cross section
SSIGN(J,K) = SSIGT-SSIGQE
SSIGNSD(J,K) = SSIGQSD
ALINT(J,K) = 1.D0/(AVOG*SSIGn(j,K)/ATARGET)
ENDDO
CALL SIB_HADCS1
& (K,SQS,SIGT1,SIGEL1,SIGINEL1,SLOPE1,RHO1)
SIGEFF = 0.68D0*(1.D0+36.D0/SQS**2)*
& dlog(0.6D0+XI_MAX/1.5D0*SQS**2)
ALAM(J) = dSQRT(SIGEFF/SIGEL1)
ENDIF
CALL SIB_SIGMA_HP(K,SQS,
& SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
IF(IK.eq.1)THEN
c fixed O-N mixture
CALL SIG_H_AIR
& (SIGT, SLOPE, RHO, ALAM(J),
& SSIGT, SSIGEL, SSIGQE, SSIGSD, SSIGQSD)
ELSE
CALL SIG_H_NUC
& (IAT, SIGT, SLOPE, RHO, ALAM(J),
& SSIGT, SSIGEL, SSIGQE, SSIGSD, SSIGQSD)
ENDIF
if (ndebug .gt. 0 ) WRITE(6,'(1X,I2,1P,5E12.3)')
& K,SQS,SSIGT,SSIGT-SSIGQE,SSIGQSD,ALAM(J)
C particle production cross section
SSIGN(J,K,IK) = SSIGT-SSIGQE
c diffractive cross section
SSIGNSD(J,K,IK) = SSIGQSD
c elastic cross section
SSIGNEL(J,K,IK) = SSIGEL
c interaction length
IF(IK.eq.1)then
ALINT(J,K,IK) = 1.D0/(AVOG*SSIGn(j,K,IK)/ATARGET)
else
ALINT(J,K,IK) = 1.D0/(AVOG*SSIGn(j,K,IK)/IAT)
endif
ENDDO
ENDDO
if (ndebug .gt. 0 ) then
WRITE(6,'(/,1X,A)')
& ' SIG_AIR_INI: NUCLIB interaction lengths [g/cm**2]'
WRITE(6,*) 'target:', TARGN(IK)
WRITE(6,'(1X,A)')
& ' sqs, p-targ, pi-targ, K-targ'
DO J=1,NSQS
ASQS = ASQSMIN + DASQS*DBLE(J-1)
SQS = 10.D0**ASQS
WRITE(6,'(1X,1P,4E12.3)')
& SQS,ALINT(J,1,IK),ALINT(J,2,IK),ALINT(J,3,IK)
ENDDO
endif
ENDDO
if (ndebug .gt. 0 ) then
WRITE(LUN,'(/,1X,A)')
& ' SIG_AIR_INI: NUCLIB interaction lengths [g/cm**2]'
WRITE(LUN,'(1X,A)')
& ' sqs, p-air, pi-air, K-air'
DO J=1,NSQS
ASQS = ASQSMIN + DASQS*DBLE(J-1)
SQS = 10.D0**ASQS
WRITE(LUN,'(1X,1P,4E12.3)')
& SQS,ALINT(J,1),ALINT(J,2),ALINT(J,3)
ENDDO
endif
END
C=======================================================================
......@@ -11768,15 +11807,19 @@ c empirical parametrization
ALAM = dsqrt(SIGEFF/SSIGEL)
SSIGSD = 2.D0 * SIGEFF
elseif( IPARM.EQ.2 ) then
C parametrization by Goulianos
c lambda derived from proton interactions
CALL SIB_HADCS1(1,ECM,SIGT1,SSIGEL1,SIGINEL1,SLOPE1,RHO1)
C parametrization by Goulianos for diff. interaction
SIGEFF = 0.68D0*(1.D0+36.D0/Ecm**2)
& *dLOG(0.6D0+0.02D0/1.5D0*Ecm**2)
& *LOG(0.6D0+0.02D0/1.5D0*Ecm**2)
SIGEFF = MAX(0.D0,SIGEFF)
ALAM = dsqrt(SIGEFF/SSIGEL)
ALAM = sqrt(SIGEFF/SSIGEL1)
SSIGSD = 2.D0 * SIGEFF
elseif( IPARM.eq.3)then
C data from Paolo Lipari's note
......@@ -11835,6 +11878,7 @@ C. INPUT : SSIG (mbarn) total pp cross section
C. SLOPE (GeV**-2) elastic scattering slope for pp
C. ALPHA real/imaginary part of the forward pp elastic
C. scattering amplitude
C. ALAM coupling to inel. intermediat states
C. OUTPUT : SIGT = Total cross section
C. SIGEL = Elastic cross section
C. SIGQEL = Elastic + Quasi elastic cross section
......@@ -11864,6 +11908,55 @@ C......................................................................
C=======================================================================
SUBROUTINE SIG_H_NUC
+ (IAT,SSIG,SLOPE,ALPHA,ALAM,SIGT,SIGEL,SIGQE,SIGSD,SIGQSD)
C-----------------------------------------------------------------------
C**********************************************************************
C...Subroutine to compute hadron-nucleus cross sections
C. according to:
C. R.J. Glauber and G.Matthiae Nucl.Phys. B21, 135, (1970)
C.
C. INPUT : IAT nucleon number in target nucleus
C. SSIG (mbarn) total pp cross section
C. SLOPE (GeV**-2) elastic scattering slope for pp
C. ALPHA real/imaginary part of the forward pp elastic
C. scattering amplitude
C. OUTPUT : SIGT = Total cross section
C. SIGEL = Elastic cross section
C. SIGQEL = Elastic + Quasi elastic cross section
C. SIGSD = single diff. cross section (beam)
C. SIGQSD = Elastic + Quasi elastic SD cross section (beam)
C.
C......................................................................
IMPLICIT NONE
INTEGER IAT
DOUBLE PRECISION SSIG,SLOPE,ALPHA,ALAM
DOUBLE PRECISION SIG1,SIGEL1,SIGQE1,SIGSD1,SIGQSD1
DOUBLE PRECISION SIGT,SIGEL,SIGQE,SIGSD,SIGQSD
INTEGER NCALL, NDEBUG, LUN
COMMON /S_DEBUG/ NCALL, NDEBUG, LUN
SAVE
IF(IAT.eq.0.or.IAT.gt.18) THEN
WRITE(LUN,'(//,1X,A)')
& ' SIG_H_NUC: number of target nucleons too large!',
& ' (1<=IAT<=18)'
SIGT = -1.D0
STOP
ENDIF
CALL GLAUBER2
+ (IAT,SSIG,SLOPE,ALPHA,ALAM,SIG1,SIGEL1,SIGQE1,SIGSD1,SIGQSD1)
SIGT = SIG1
SIGEL = SIGEL1
SIGQE = SIGQE1
SIGSD = SIGSD1
SIGQSD = SIGQSD1
RETURN
END
C=======================================================================
SUBROUTINE GLAUBER2
+ (JA,SSIG,SLOPE,ALPHA,ALAM,SIGT,SIGEL,SIGQEL,SIGSD,SIGQSD)
......@@ -13773,12 +13866,12 @@ c external types
INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B,
& SSIG_RHO
COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
......@@ -13880,12 +13973,12 @@ c external types
INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
DOUBLE PRECISION SSIG_TOT,SSIG_SD1,SSIG_SD2,SSIG_DD,SSIG_B,
& SSIG_RHO
COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
......@@ -13958,12 +14051,12 @@ Cf2py double precision, intent(out) :: SIGprod,SIGbdif
INTEGER NS_max, NH_max
PARAMETER (NS_max = 20, NH_max = 80)
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
c external
DOUBLE PRECISION SQS,SIGPROD,SIGBDIF
......@@ -13992,13 +14085,12 @@ c internal
J1 = max(J1,1)
endif
T = (AL-1.D0)*10.D0 - DBLE(J1-1)
SIGprod = SSIGN(J1,L)*(1.D0-T) + SSIGN(J1+1,L)*T
SIGbdif = SSIGNSD(J1,L)*(1.D0-T) + SSIGNSD(J1+1,L)*T
SIGprod = SSIGN(J1,L,1)*(1.D0-T) + SSIGN(J1+1,L,1)*T
SIGbdif = SSIGNSD(J1,L,1)*(1.D0-T) + SSIGNSD(J1+1,L,1)*T
END
C=======================================================================
SUBROUTINE SIB_SIGMA_HNUC (L,IAT,SQS,SIGprod,SIGbdif)
SUBROUTINE SIB_SIGMA_HNUC (L,IAT,SQS,SIGprod,SIGbdif,SIGela)
C-----------------------------------------------------------------------
C calculate Hadron-nucleus cross sections
......@@ -14012,6 +14104,7 @@ C SQS sqrt(s)
C
C output: SIGprod particle production cross section (mb)
C SIGbdif q.ela and ela beam diff. cross section
C SIGela elastic cross section
C-----------------------------------------------------------------------
Cf2py integer, intent(in) :: L,IAT
Cf2py double precision, intent(in) :: SQS
......@@ -14023,12 +14116,12 @@ Cf2py double precision, intent(out) :: SIGprod,SIGbdif
INTEGER NCALL, NDEBUG, LUN
COMMON /S_DEBUG/ NCALL, NDEBUG, LUN
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,ALINT,ASQSMIN,ASQSMAX,
& DASQS
DOUBLE PRECISION SSIG,PJETC,SSIGN,SSIGNSD,SSIGNEL,ALINT,ASQSMIN,
& ASQSMAX,DASQS
INTEGER NSQS
COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
& SSIGN(61,3), SSIGNSD(61,3), ALINT(61,3),
& ASQSMIN, ASQSMAX, DASQS, NSQS
& SSIGN(61,3,3), SSIGNSD(61,3,3), SSIGNEL(61,3,3),
& ALINT(61,3,3), ASQSMIN, ASQSMAX, DASQS, NSQS
C--------------------------------------------------------------------
C SIBYLL utility common blocks containing constants \FR'14
......@@ -14049,21 +14142,48 @@ C--------------------------------------------------------------------
+ SIGQSD,SIGPPT,SIGPPEL,SIGPPSD,ITG
c external
DOUBLE PRECISION SQS,SIGPROD,SIGBDIF
DOUBLE PRECISION SQS,SIGPROD,SIGBDIF,SIGELA
INTEGER L,IAT
c internal
DOUBLE PRECISION ALAM
INTEGER IPARM,ICSMOD
DOUBLE PRECISION ALAM,T,AL
INTEGER IPARM,ICSMOD,IK,J1
SAVE
IF(NSQS.LE.0) THEN
WRITE(LUN,'(//,1X,A)')
& ' SIB_SIGMA_HNUC: interpolation table not initialized.'
STOP
ENDIF
IF(IAT.ge.0.and.IAT.lt.19)THEN
IF(IAT.eq.0.or.IAT.eq.14.or.IAT.eq.16)THEN
c read cross section from table
IF(IAT.eq.0) THEN
IK=1
ELSEIF(IAT.eq.14)THEN
IK=2
ELSE
IK=3
ENDIF
AL = LOG10(SQS)
J1 = INT((AL - 1.D0)*10.D0 + 1)
if((j1.lt.1).or.(j1.gt.NSQS)) then
if (ndebug .gt. 0)
& write (LUN,'(1x,a,i3,1p,e12.3)')
& ' SIB_SIGMA_HNUC: energy out of range ',L,sqs
endif
if((j1.lt.1).or.(j1.ge.NSQS)) then
J1 = min(J1,NSQS-1)
J1 = max(J1,1)
endif
T = (AL-1.D0)*10.D0 - DBLE(J1-1)
SIGprod = SSIGN(J1,L,IK)*(1.D0-T) + SSIGN(J1+1,L,IK)*T
SIGbdif = SSIGNSD(J1,L,IK)*(1.D0-T) + SSIGNSD(J1+1,L,IK)*T
SIGela = SSIGNEL(J1,L,IK)*(1.D0-T) + SSIGNEL(J1+1,L,IK)*T
ELSEIF(IAT.lt.19)THEN
c calculate cross section
IF(ndebug.gt.0)THEN
WRITE(LUN,'(1X,A,2I4,F8.2)')
& 'SIB_SIGMA_HNUC: L,IAT,SQS:',L,IAT,SQS
ENDIF
c calculate hadron - nucleus cross section
c dummy arg, coupling derived from dif xsctn
ALAM = 1.D0
......@@ -14076,6 +14196,13 @@ C particle production cross section