diff --git a/Makefile b/Makefile
index 05855b43d2a771b2f813767af090af9ecb8497a2..c20a1cee6f48e944e4b28df20d69c07febf5f149 100644
--- a/Makefile
+++ b/Makefile
@@ -4,7 +4,7 @@ VER := 7.60
# activate this, if you want to manipulate cross-sections and
# multiplicities etc., or you want to write/read particle lists
# during the simulation
-#CONEX_EXTENSIONS+=-DCONEX_EXTENSIONS #
+CONEX_EXTENSIONS+=-DCONEX_EXTENSIONS #
# system type and directory layout
SYSTEM := $(shell uname -s)
diff --git a/src/conex_mod.F b/src/conex_mod.F
index eac4f5aa75da6bc30defbc45affa903e3a4ccbb6..5a4ea78f7e89d5d4fb48d1dfe406f6541c303ace 100644
--- a/src/conex_mod.F
+++ b/src/conex_mod.F
@@ -49169,13 +49169,13 @@ C. SQS (GeV) is the center of mass energy of each
C. nucleon - nucleon cross section
C-----------------------------------------------------------------------
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+ IMPLICIT INTEGER(I-N)
#ifdef CONEX_EXTENSIONS
c Not needed anymore and create a conflict with common SIB_CST for "pi" defined in conex.h too
c#include "conex.h"
c include 'conex.incnex'
c CONEX_EXTENSIONS
#endif
- IMPLICIT INTEGER(I-N)
C The final particle output is contained in COMMON /S_PLIST/
C NP : number of final particles
C P(1:NP, 1:5) : 4-momenta + masses of the final particles
@@ -49205,12 +49205,12 @@ C--------------------------------------------------------------------
COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
+ ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
+ ,JJAEL(IAMAX), JJBEL(IAMAX)
+ COMMON /FRAGMENTS/ PPP(3,60)
#ifdef CONEX_EXTENSIONS
C additional common block for interaction history information
COMMON /S_NUCLST/ LBEGIN(IAMAX), LEND(IAMAX), LUSED
c CONEX_EXTENSIONS
#endif
- COMMON /FRAGMENTS/ PPP(3,60)
DIMENSION SIGDIF(3)
DIMENSION IAF(60)
DOUBLE PRECISION FOX
@@ -49230,18 +49230,18 @@ c select target IATARGET from air composition
C...Single nucleon (proton) case
IF (IAB .EQ. 1) THEN
+ NPA = 0
#ifdef CONEX_EXTENSIONS
LBEGIN(1) = NPA+1
#endif
- NPA = 0
CALL SIBYLL (13,IATARGET, ECM)
+ CALL DECSIB
+ DO J=1,NP
c#ifdef CONEX_EXTENSIONS
c if (doResampling) then
c call resampleSIBYLL(NP, xselab, idprojxs, SQS, xsamproj)
c endif
c#endif
- CALL DECSIB
- DO J=1,NP
LA = IABS(LLIST(J))
IF (LA .LT. 10000) THEN
NPA = NPA + 1
@@ -49250,11 +49250,11 @@ c#endif
PA(K,NPA) = P(J,K)
ENDDO
ENDIF
+ ENDDO
#ifdef CONEX_EXTENSIONS
LEND(1) = NPA
LUSED = 1
#endif
- ENDDO
RETURN
ENDIF
@@ -49312,13 +49312,13 @@ C...Elastically scattered fragments
C...Superimpose NB nucleon interactions
DO JJ=1,NB
CALL SIBYLL (13,IATARGET, ECM)
+ CALL DECSIB
#ifdef CONEX_EXTENSIONS
LBEGIN(JJ) = NPA+1
c if (doResampling) then
c call resampleSIBYLL(NP, xselab, idprojxs, SQS, xsamproj)
c endif
#endif
- CALL DECSIB
DO J=1,NP
LA = IABS(LLIST(J))
IF (LA .LT. 10000) THEN
@@ -66669,21 +66669,11 @@ c setup quarks: u,d,s
ELSEIF(IABS(KF).eq.78)THEN ! D*+
KP1 = 4
KP2 = -2
-#ifdef CONEX_EXTENSIONS
- COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC ! RU 04/2009:
- LOGICAL INITNUC ! See SIGNUC_INI_INI for details
-#endif
ELSEIF(IABS(KF).eq.79)THEN ! D*-
KP1 = 2
KP2 = -4
-#ifdef CONEX_EXTENSIONS
- if (INITNUC) then
- call SIGNUC_INI_INI ! VERY SLOW (~hour)
- endif
-#endif
-
ELSEIF(IABS(KF).eq.80)THEN ! D*0
KP1 = 4
KP2 = -1
@@ -66851,11 +66841,21 @@ C-----------------------------------------------------------------------
C...Initialization for the generation of nucleus-nucleus interactions
C. INPUT : E0 (TeV) Energy per nucleon of the beam nucleus
C........................................................................
+#ifdef CONEX_EXTENSIONS
+ COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC ! RU 04/2009:
+ LOGICAL INITNUC ! See SIGNUC_INI_INI for details
+#endif
SAVE
CALL NUC_GEOM_INI ! nucleus profiles
CALL SIGMA_INI ! initialize pp cross sections
+#ifdef CONEX_EXTENSIONS
+ if (INITNUC) then
+ call SIGNUC_INI_INI ! VERY SLOW (~hour)
+ endif
+#endif
+
RETURN
END
C=======================================================================
@@ -69154,83 +69154,6 @@ C..........................................................................
M1AEL(NAEL) = M1AEL(NAEL)+1
M1BEL(NBEL) = M1BEL(NBEL)+1
ELSE
-#ifdef CONEX_EXTENSIONS
-
- SUBROUTINE DO_SIGNUC_INI_INI(useInitNuc)
- COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC
- LOGICAL INITNUC
- LOGICAL useInitNuc
- SAVE
- INITNUC=useInitNuc
- END
-
- SUBROUTINE SIGNUC_INI_INI
-c.....This routine is added 04/2009 by R. Ulrich (ralf.ulrich@kit.edu) to
-c. compute nucleus-air interaction cross sections from first principles
-c. of pp-scattering.
-c. The tabulations provided in SIGNUC_INI are thus not used !
- COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC
- LOGICAL INITNUC
- SAVE
- double precision factMod, factModPP
- DIMENSION AA(6)
- DATA NE /6/
- DATA AA /1.,2.,3.,4.,5.,6./
-c ru ... output sibyll p-p parameters ...
-c write(*,*) 'E0/D:SQS/D:SIGT/D:SIGEL/D:SIGINEL/D:SLOPE/D:RHO/D'
-c write(*,*) '# sigma_pp'
-c do i = 1, 150
-c E0 = 10.**(-3. + 0.09 * i) ! TeV
-c ESQS = SQRT(2000.*0.938*E0)
-c call SIGMA_PP(E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
-c write(*,*) E0,ESQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO
-c enddo
-c ru ...
- write (*,*), 'SIGNUC_INI_INI: ',
- + 'Initialize SIBYLL nucleus-air cross sections (up to 1 hour)'
- do I=1,NE ! do energy bins
- ASQS = AA(I)
- E0 = 10**(ASQS*2.0) / 1.876E+03
- call SIGMA_PP(E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
- factMod=1.d0
- nproton=1
- call modifiercx(factMod, e0*1.d3, nproton)
- factModPP=1.d0
- call sibyllf19(factModPP, factMod)
- SIGT=factModPP*SIGT
- SIGINEL=factModPP*SIGINEL
-c rescaling variant A
-c sigel=factModPP**2*sigel
-c rescaling variant B
- SIGEL=factModPP*SIGEL
- SLOPE=factModPP*SLOPE
-
- do J=2,56 ! do nuclei from A=2..56
- NINT=300000
-c NINT=1000
- 10 call SIGMA_AIR(J,SIGINEL,SIGEL,NINT,
- + SIGMA,DSIGMA,SIGQE,DSIGQE)
- if (DSIGMA/SIGMA.gt.1.e-2) then
-c if (DSIGMA/SIGMA.gt.1.e2) then
- write (*,*) 'SIGNUC_INI_INI: precision too low: ',
- + DSIGMA/SIGMA
- write (*,*) 'SIGNUC_INI_INI: inrease precision to: ',
- + NINT
- NINT=NINT*3 ! increase precision
- goto 10
- endif
- SIGMATAB(I,J)=SIGMA
- SIGQETAB(I,J)=SIGQE
- write (*,*) 'E0,SQS,A,SIG,dSIG,SIGEQ,dSIGQE', E0,AA(I),J,
- + SIGMA,DSIGMA,SIGQE,DSIGQE
- enddo
- enddo
- end
-
-#endif
-c CONEX_EXTENSIONS
-
-
M2AEL(NAEL) = M2AEL(NAEL)+1
M2BEL(NBEL) = M2BEL(NBEL)+1
ENDIF
@@ -69247,10 +69170,6 @@ c CONEX_EXTENSIONS
DO J=1,IB
PROBB(J) = DBLE(MMB(J))/DBLE(M)
DPROBB(J) = SQRT(DBLE(MMB(J)))/DBLE(M)
-#ifdef CONEX_EXTENSIONS
- COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC ! RU 04/2009:
- LOGICAL INITNUC ! See SIGNUC_INI_INI for details
-#endif
ENDDO
DO J=1,IA*IB
PROBI(J) = DBLE(MMI(J))/DBLE(M)
@@ -69407,6 +69326,83 @@ C*=============================================================
C. Nucleus-nucleus cross sections
C=======================================================================
+#ifdef CONEX_EXTENSIONS
+
+ SUBROUTINE DO_SIGNUC_INI_INI(useInitNuc)
+ COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC
+ LOGICAL INITNUC
+ LOGICAL useInitNuc
+ SAVE
+ INITNUC=useInitNuc
+ END
+
+ SUBROUTINE SIGNUC_INI_INI
+c.....This routine is added 04/2009 by R. Ulrich (ralf.ulrich@kit.edu) to
+c. compute nucleus-air interaction cross sections from first principles
+c. of pp-scattering.
+c. The tabulations provided in SIGNUC_INI are thus not used !
+ COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC
+ LOGICAL INITNUC
+ SAVE
+ double precision factMod, factModPP
+ DIMENSION AA(6)
+ DATA NE /6/
+ DATA AA /1.,2.,3.,4.,5.,6./
+c ru ... output sibyll p-p parameters ...
+c write(*,*) 'E0/D:SQS/D:SIGT/D:SIGEL/D:SIGINEL/D:SLOPE/D:RHO/D'
+c write(*,*) '# sigma_pp'
+c do i = 1, 150
+c E0 = 10.**(-3. + 0.09 * i) ! TeV
+c ESQS = SQRT(2000.*0.938*E0)
+c call SIGMA_PP(E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
+c write(*,*) E0,ESQS,SIGT,SIGEL,SIGINEL,SLOPE,RHO
+c enddo
+c ru ...
+ write (*,*), 'SIGNUC_INI_INI: ',
+ + 'Initialize SIBYLL nucleus-air cross sections (up to 1 hour)'
+ do I=1,NE ! do energy bins
+ ASQS = AA(I)
+ E0 = 10**(ASQS*2.0) / 1.876E+03
+ call SIGMA_PP(E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
+ factMod=1.d0
+ nproton=1
+ call modifiercx(factMod, e0*1.d3, nproton)
+ factModPP=1.d0
+ call sibyllf19(factModPP, factMod)
+ SIGT=factModPP*SIGT
+ SIGINEL=factModPP*SIGINEL
+c rescaling variant A
+c sigel=factModPP**2*sigel
+c rescaling variant B
+ SIGEL=factModPP*SIGEL
+ SLOPE=factModPP*SLOPE
+
+ do J=2,56 ! do nuclei from A=2..56
+ NINT=300000
+c NINT=1000
+ 10 call SIGMA_AIR(J,SIGINEL,SIGEL,NINT,
+ + SIGMA,DSIGMA,SIGQE,DSIGQE)
+ if (DSIGMA/SIGMA.gt.1.e-2) then
+c if (DSIGMA/SIGMA.gt.1.e2) then
+ write (*,*) 'SIGNUC_INI_INI: precision too low: ',
+ + DSIGMA/SIGMA
+ write (*,*) 'SIGNUC_INI_INI: inrease precision to: ',
+ + NINT
+ NINT=NINT*3 ! increase precision
+ goto 10
+ endif
+ SIGMATAB(I,J)=SIGMA
+ SIGQETAB(I,J)=SIGQE
+ write (*,*) 'E0,SQS,A,SIG,dSIG,SIGEQ,dSIGQE', E0,AA(I),J,
+ + SIGMA,DSIGMA,SIGQE,DSIGQE
+ enddo
+ enddo
+ end
+
+#endif
+c CONEX_EXTENSIONS
+
+
SUBROUTINE SIGNUC_INI (IA,E0)
C-----------------------------------------------------------------------
@@ -69423,6 +69419,10 @@ C.
C........................................................
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT INTEGER(I-N)
+#ifdef CONEX_EXTENSIONS
+ COMMON /CLENNNTAB/ SIGMATAB(6,56),SIGQETAB(6,56),INITNUC ! RU 04/2009:
+ LOGICAL INITNUC ! See SIGNUC_INI_INI for details
+#endif
COMMON /CLENNN/ SSIGNUC(60), ALNUC(60)
DIMENSION SIGMA(6,56), SIGQE(6,56)
DIMENSION AA(6)
@@ -69486,21 +69486,10 @@ C...Data on `inelastic-production' nucleus-air cross section
&1.228D+03,1.304D+03,1.438D+03,1.546D+03,1.714D+03,1.836D+03/
DATA (SIGMA(J, 28),J=1,6) /
&1.289D+03,1.370D+03,1.451D+03,1.597D+03,1.754D+03,1.913D+03/
-#ifdef CONEX_EXTENSIONS
- if (INITNUC) then
- S1 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
- + SIGMATAB(JE,JA),SIGMATAB(JE+1,JA),SIGMATAB(JE+2,JA))
- S2 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
- + SIGQETAB(JE,JA),SIGQETAB(JE+1,JA),SIGQETAB(JE+2,JA))
- else
-#endif
DATA (SIGMA(J, 29),J=1,6) /
&1.411D+03,1.469D+03,1.613D+03,1.777D+03,1.910D+03,2.075D+03/
DATA (SIGMA(J, 30),J=1,6) /
&1.347D+03,1.401D+03,1.498D+03,1.642D+03,1.816D+03,1.975D+03/
-#ifdef CONEX_EXTENSIONS
- endif
-#endif
DATA (SIGMA(J, 31),J=1,6) /
&1.359D+03,1.448D+03,1.551D+03,1.694D+03,1.858D+03,2.007D+03/
DATA (SIGMA(J, 32),J=1,6) /
@@ -69669,10 +69658,21 @@ C...Data on `inelastic-production' nucleus-air cross section
JE = MIN(INT((ASQS-AMIN)/DA)+1,NE-2)
DO JA=2,IA
ABEAM = DBLE(JA)
+#ifdef CONEX_EXTENSIONS
+ if (INITNUC) then
+ S1 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ + SIGMATAB(JE,JA),SIGMATAB(JE+1,JA),SIGMATAB(JE+2,JA))
+ S2 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ + SIGQETAB(JE,JA),SIGQETAB(JE+1,JA),SIGQETAB(JE+2,JA))
+ else
+#endif
S1 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ SIGMA(JE,JA),SIGMA(JE+1,JA),SIGMA(JE+2,JA))
S2 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ SIGQE(JE,JA),SIGQE(JE+1,JA),SIGQE(JE+2,JA))
+#ifdef CONEX_EXTENSIONS
+ endif
+#endif
SSIGNUC(JA) = S1 + S2
ALNUC(JA) = ATARGET/(AVOG*SSIGNUC(JA))
ENDDO
@@ -69762,27 +69762,6 @@ C..............................................................
if((yy(k) .gt. y) .eqv. (yy(n) .gt. yy(1))) goto 100
enddo
100 y2 = yy(k)
-#ifdef CONEX_EXTENSIONS
-
-* cross section as calculated in SIBYLL, but re-scaled (RU 02/2009)
-
- ELSE IF(ICSPA.EQ.4) THEN
- CALL SIB_SIGMA_HP(1,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
- factMod=1.d0
- nproton=1
- call modifiercx(factMod, e0*1.d3,nproton)
- sigt=factMod*sigt
- sigdif(1)=factMod*sigdif(1)
- sigdif(2)=factMod*sigdif(2)
- sigdif(3)=factMod*sigdif(3)
- siginel=factMod*siginel
-c rescaling variant A
-c sigel=factMod**2*sigel
-c rescaling variant B
- sigel=factMod*sigel
- slope=factMod*slope
-#endif
-
y1 = yy(k-1)
k0 = k-1
x1 = xmin + dx*float(k-2)
@@ -69934,30 +69913,6 @@ C...p-p inelastic cross sections (mbarn)
DELDL = 0.0808D0
EPSDL = -0.4525D0
S = SQS*SQS
-#ifdef CONEX_EXTENSIONS
-
-* cross section as calculated in SIBYLL, but re-scaled (RU 02/2009)
-
- ELSE IF(ICSPA.EQ.4) THEN
-
- CALL SIB_SIGMA_HP(2,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
- factMod=1.d0
- npion=2
- call modifiercx(factMod, e0*1.d3,npion)
- sigt=factMod*sigt
- sigdif(1)=factMod*sigdif(1)
- sigdif(2)=factMod*sigdif(2)
- sigdif(3)=factMod*sigdif(3)
- siginel=factMod*siginel
-c rescaling variant A
-c sigel=factMod**2*sigel
-c rescaling variant B
- sigel=factMod*sigel
- slope=factMod*slope
-
-#endif
-c CONEX_EXTENSIONS
-
SIGT = 21.7D0*S**DELDL+56.08D0*S**EPSDL
SIGEL = R*SIGT
@@ -69979,6 +69934,27 @@ c taking slope parameterization from sigma_pp Donnie.-Landshoff
CALL SIG_RPP2014(1,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)
+#ifdef CONEX_EXTENSIONS
+
+* cross section as calculated in SIBYLL, but re-scaled (RU 02/2009)
+
+ ELSE IF(ICSPA.EQ.4) THEN
+ CALL SIB_SIGMA_HP(1,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
+ factMod=1.d0
+ nproton=1
+ call modifiercx(factMod, e0*1.d3,nproton)
+ sigt=factMod*sigt
+ sigdif(1)=factMod*sigdif(1)
+ sigdif(2)=factMod*sigdif(2)
+ sigdif(3)=factMod*sigdif(3)
+ siginel=factMod*siginel
+c rescaling variant A
+c sigel=factMod**2*sigel
+c rescaling variant B
+ sigel=factMod*sigel
+ slope=factMod*slope
+#endif
+
ENDIF
RETURN
@@ -70112,29 +70088,6 @@ C...pi-p inelastic cross sections (mbarn)
SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL
SIGEL = R*SIGT
-
-#ifdef CONEX_EXTENSIONS
-
-* cross section as calculated in SIBYLL, but re-scaled (RU 02/2009)
-
- ELSE IF(ICSPA.EQ.4) THEN
-
- CALL SIB_SIGMA_HP(3,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
- factMod=1.d0
- nkaon=4
- call modifiercx(factMod, e0*1.d3, nkaon)
- sigt=factMod*sigt
- sigdif(1)=factMod*sigdif(1)
- sigdif(2)=factMod*sigdif(2)
- sigdif(3)=factMod*sigdif(3)
- siginel=factMod*siginel
-c rescaling variant A
-c sigel=factMod**2*sigel
-c rescaling variant B
- sigel=factMod*sigel
- slope=factMod*slope
-
-#endif
SIGINEL = SIGT-SIGEL
SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
RHO = 0.D0