IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/cxroot
1 result
Show changes
Commits on Source (6)
Hide whitespace changes
Inline Side-by-side
Makefile
View file @ cfcd301f
......@@ -34,8 +34,8 @@ CFGDIR = $(CONEX_PREFIX)/cfg
# general compiling options
MC3D =
#MC3D = -D__MC3D__ -D__COAST__
LOWMEMORY = -D__SORT_FOR_ENERGY__
#LOWMEMORY= -D__SAVEMEMO__ -D__SORT_FOR_ENERGY__ -D__CXDEBUG__
#LOWMEMORY = -D__SORT_FOR_ENERGY__
LOWMEMORY= -D__SAVEMEMO__ -D__SORT_FOR_ENERGY__ -D__CXDEBUG__
CPPFLAGS = -D__URQMD__ -D__CXSUB__ $(LOWMEMORY) $(CONEX_EXTENSIONS) $(MC3D) -fPIC
# compiler selection
COMP := `which gfortran 2>&1 | grep -v which`
......
src/conex_mod.F
View file @ cfcd301f
......@@ -871,7 +871,7 @@ C-----------------------------------------------------------------------
END
#endif
 
#ifndef __CORSIKA8__
#if !__CXCORSIKA__ && !__CORSIKA8__
#ifdef __NEXUS__
c-----------------------------------------------------------------------
......@@ -22022,9 +22022,9 @@ c-----------------------------------------------------------------------
INCLUDE '(PHNCCM)' !
INCLUDE '(CTITLE)' !
 
C FLUKA 2011.2 BLOCK DATA PROGRAMS
C FLUKA 2020.4 BLOCK DATA PROGRAMS
EXTERNAL BDINPT,BDTRNS,BDHDR1,BDHDR2,BDHDR3,BDPART,BDPRDC,
& BDNOPT,BDEVAP,BDPREE
& BDNOPT,BDEVAP,BDPREE,BDESJE
 
DIMENSION WHAT (6)
CHARACTER SDUM*10
......@@ -22197,7 +22197,6 @@ c-----------------------------------------------------------------------
if(isx.ge.4)call cxalist('Determine FLUKA Production&',0,0,0)
#endif
nptlxs=0
np=0
 
 
C CONVERT PARTICLE TYPE TO FLUKA
......@@ -22224,6 +22223,10 @@ C USE THE FLUKA-INTERNALLY USED DIRECTION COSINES:
IJ = KPROJ
C CALCULATE THE MOMENTUM WITH FLUKA MASSES
POO = PPROJ
C RESET STACK INDEX
NP=0
NP0=0
NPHEAV=0
C NOW INTERACTION IS PERFORMED
CALL EVENTV( IJ,POO,EKIN1,TXX,TYY,TZZ,WEE,MMMAT )
c write(ifck,*)'target used',ibtar,ichtar,mmmat,ibres
......@@ -48864,7 +48867,6 @@ c------------------------------------------------------------------------------
end
 
#if !__CXCORSIKA__ && !__CORSIKA8__
c--------------------------------------------------------------------
double precision function qgran(b10)
c--------------------------------------------------------------------
......@@ -48876,14 +48878,9 @@ c--------------------------------------------------------------------
return
end
 
integer function size(array)
double precision array(*)
size=int(array(1))
end
#endif
 
#endif
......@@ -49172,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
......@@ -49208,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
......@@ -49233,17 +49230,17 @@ 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
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
......@@ -49253,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
 
......@@ -49315,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
......@@ -66672,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
......@@ -66854,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=======================================================================
......@@ -69157,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
......@@ -69250,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)
......@@ -69410,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-----------------------------------------------------------------------
......@@ -69426,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)
......@@ -69489,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) /
......@@ -69672,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
......@@ -69765,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)
......@@ -69937,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
......@@ -69982,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
......@@ -70115,17 +70088,34 @@ C...pi-p inelastic cross sections (mbarn)
SIGT = 13.63D0*S**DELDL+(36.02D0+27.56D0)/2.D0*S**EPSDL
 
SIGEL = R*SIGT
SIGINEL = SIGT-SIGEL
SLOPE = SIGT**2/(16.D0*PI*SIGEL)/CMBARN
RHO = 0.D0
 
c ICSPA=4 reserved for CONEX_EXTENSION
c ELSE IF(ICSPA.EQ.4) THEN
* cross section from 2014 Review of Particle Physics
ELSE IF(ICSPA.EQ.5) THEN
c elastic slope not included in fit
c taking slope parameterization from sigma_pp Donnie.-Landshoff
ALPHAP = 0.25D0
SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
CALL SIG_RPP2014(2,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(3,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
CALL SIB_SIGMA_HP(2,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
factMod=1.d0
nkaon=4
call modifiercx(factMod, e0*1.d3, nkaon)
npion=2
call modifiercx(factMod, e0*1.d3,npion)
sigt=factMod*sigt
sigdif(1)=factMod*sigdif(1)
sigdif(2)=factMod*sigdif(2)
......@@ -70137,25 +70127,9 @@ 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
c ICSPA=4 reserved for CONEX_EXTENSION
c ELSE IF(ICSPA.EQ.4) THEN
#endif
c CONEX_EXTENSIONS
 
* cross section from 2014 Review of Particle Physics
ELSE IF(ICSPA.EQ.5) THEN
c elastic slope not included in fit
c taking slope parameterization from sigma_pp Donnie.-Landshoff
ALPHAP = 0.25D0
SLOPE = 8.5D0+4.D0*ALPHAP*LOG(SQS)
CALL SIG_RPP2014(2,1,SQS,SLOPE,SIGT,SIGEL,SIGINEL,RHO)
ENDIF
 
 
......@@ -70311,6 +70285,29 @@ c taking slope parameterization from sigma_pp Donnie.-Landshoff
CALL SIG_RPP2014(3,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(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
ENDIF
 
RETURN
......@@ -157633,6 +157630,9 @@ c-----------------------------------------------------------------------
 
c Conex common
 
double precision EgyHiLoLim
integer MCModel,MCleModel,ilowegy
common/cxmodel/EgyHiLoLim,MCModel,MCleModel,ilowegy
double precision airz,aira,airw,airavz,airava,airi
common/cxair/airz(3),aira(3),airw(3),airavz,airava,airi(3)
common/xsfiles/ixsfcx
......@@ -157700,7 +157700,24 @@ c air
enddo
airavanxs=sngl(airava)
airavznxs=sngl(airavz)
#ifdef __BUGCORS__
airanxs(1)=14.
airanxs(2)=16.
airanxs(3)=40.
airznxs(1)=7.
airznxs(2)=8.
airznxs(3)=18.
airwnxs(1)=0.78479D0
airwnxs(2)=0.21052D0
airwnxs(3)=0.00469D0
airavznxs=0.
airavanxs=14.543
do i=1,3
airavznxs=airavznxs+airwnxs(i)*airznxs(i)
enddo
#endif
iappl = iapplxs
nevent = neventxs
iframe = iframexs
src/conex_sub.F
View file @ cfcd301f
......@@ -164,10 +164,10 @@ c To give priority to CE for EM (same results in hybrid or MC mode), if elow > e
endif
thetas=theta
phisho=phi
XmaxP=Xmaximum
#if __CXCORSIKA__ || __CORSIKA8__
XminP=hstart !temporary use of XminP in meter (starting altitude to define minimum height in Initialize2)
#endif
XmaxP=Xmaximum
altitude=dimpact
call ranfgt(seed) !get seed before shower
#ifdef __CXCORSIKA__
......@@ -203,9 +203,6 @@ c CORSIKA 8 treat the hadronic shower, so only cascade equations are compu
call InitialParticleSho(id) !called only to setup stacks
lxfirst=.true. !first interaction in CORSIKA 8
lxfirstIn=.false. !do not fix first interaction
c ethin=thin*eprima !egs4
c isxegs=isx
c ZBOUND=zshmax
return
 
entry ConexCascade() !call after source function is filled
......@@ -3155,10 +3152,10 @@ c-----------------------------------------------------------------------
IQIN=id
if(abs(IQIN).gt.1)then
#else
if(id.eq. 12)then !electron
if(id.eq. 12)then !electron
if(mode.eq.3.and.ekin.le.eecut)mc2ce=.true.
IQIN=-1
elseif(id.eq.-12)then !positron
elseif(id.eq.-12)then !positron
if(mode.eq.3.and.ekin.le.eecut)mc2ce=.true.
IQIN= 1
elseif(abs(id).eq. 10)then !photon
......@@ -3166,7 +3163,7 @@ c-----------------------------------------------------------------------
IQIN= 0
else !neutrino : lost energy
#endif
if(iwrt.ge.2)call Profana(dptl(13)-0.000001d0*dzHa,zshmax,
if(iwrt.ge.2)call Profana(dptl(13)-0.000001d0*dzHa,zshmax,
& dptl(4),dptl(4),dptl(11),999,2)
#ifdef __CXDEBUG__
etotsource=etotsource-dptl(4)*dptl(11)
......@@ -6413,7 +6410,6 @@ c write(*,*)'ROOT output file : ', DSN(1:nfnrt),'.root'
 
end
 
 
c$$$c-----------------------------------------------------------------------
c$$$ subroutine AddDataPath(file)
......@@ -10514,7 +10510,7 @@ c low energy model oriented input
if(ilowegy.eq.1)then
 
iehlim1=iehlim-1
if(ifwle.gt.0)then
write(6,'(a,a)')'read LE model table from: ',fnwle(1:nfnwle)
luseCompress = (index(fnwle(1:nfnwle), ".bz2") .eq. nfnwle-3)
......@@ -10758,12 +10754,12 @@ c high energy model oriented input
if(n1maxi.gt.n1max+2.or.n2maxi.gt.n2max+1)
& stop'***** n1/n2 mismatch (high) ***** '
if(.not.go)then
iemin=max(iemin,iemin1) !set minimum to minimum in table
iehlim=iemin !only HE model
iemin=max(iemin,iemin1) !set minimum to minimum in table
iehlim=iemin !only HE model
endif
if((iemin1.gt.iemin.and..not.go).or.iehlim.lt.iemin1)
& stop'***** iemin too small (high) ***** '
iemax=iemax1 !set maximum to maximum in table
& stop'***** iemin too small (high) ***** '
iemax=iemax1 !set maximum to maximum in table
if(iehlim.lt.iemin1)stop'***** iehlim too small (high) ***** '
if (luseCompress.ne.0) then
do i1 = 1, n1max
......@@ -10975,10 +10971,10 @@ c Pions and pions 0 from charged Kaons
close(ifdkz)
endif
else
write(6,*)'Table dkz is not defined for hadron cascade !'
write(6,*)'file: ', fndkz(1:nfndkz)
write(6,*)'lzma=', luseCompress
stop
write(6,*)'Table dkz is not defined for hadron cascade !'
write(6,*)'file: ', fndkz(1:nfndkz)
write(6,*)'lzma=', luseCompress
stop
endif
 
c Pions and pions 0 from Kaon Long
......@@ -11012,8 +11008,8 @@ c Pions and pions 0 from Kaon short
else
open(ifdks,file=fndks(1:nfndks),status='old')
read(ifdks,*) aks,aks0
close(ifdks)
endif
close(ifdks)
else
write(6,*)'Table dks is not defined for hadron cascade !'
stop
......@@ -11036,7 +11032,7 @@ c Muons from Charge Kaon, Kaon Long and Charged Pions
endif
else
write(6,*)'Table dkm is not defined for hadron cascade !'
stop
stop
endif
 
 
......@@ -11523,19 +11519,9 @@ c Muons and hadrons from gamma by muon pair production or photonuclear effect
 
if(mode.ge.7)then
if(ifdkg.gt.0)then
luseCompress = (index(fndkg(1:nfndkg), ".bz2") .eq. nfndkg-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndkg(1:nfndkg))
call CorDataFillArray(agpr, size(agpr))
call CorDataFillArray(agne, size(agne))
call CorDataFillArray(agpi, size(agpi))
call CorDataFillArray(agmu, size(agmu))
call CorDataCloseFile()
else
open(ifdkg,file=fndkg(1:nfndkg),status='old')
read(ifdkg,*) agpr,agne,agpi,agmu
close(ifdkg)
endif
open(ifdkg,file=fndkg(1:nfndkg),status='old')
read(ifdkg,*) agpr,agne,agpi,agmu
close(ifdkg)
else
write(6,*)'Table dkg is not defined for hadron cascade !'
endif
......@@ -11551,20 +11537,8 @@ c low energy model oriented input
 
if(ifwle.gt.0)then
write(6,'(a,a)')'read LE model table from ',fnwle(1:nfnwle)
luseCompress = (index(fnwle(1:nfnwle), ".bz2") .eq. nfnwle-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fnwle(1:nfnwle))
call CorDataNextText(ppjver)
exmin0 = CorDataNextNumber()
iemin0 = CorDataNextNumber()
iemax0 = CorDataNextNumber()
n1max0 = CorDataNextNumber()
n2max0 = CorDataNextNumber()
nde = CorDataNextNumber()
else
open(ifwle,file=fnwle(1:nfnwle),status='old')
read(ifwle,*) ppjver,exmin0,iemin0,iemax0,n1max0,n2max0,nde
else
open(ifwle,file=fnwle(1:nfnwle),status='old')
read(ifwle,*) ppjver,exmin0,iemin0,iemax0,n1max0,n2max0,nde
if(nde.ne.nint(dnHa))stop'***** nde mismatch (low) ***** '
if(iemin0.gt.iemin)stop'***** iemin mismatch (low) ***** '
if(abs(dble(exmin0)-exmin).gt.1.d-4)
......@@ -11572,31 +11546,18 @@ c low energy model oriented input
if(n1maxi.gt.n1max0+2.or.n2maxi.gt.n2max0+1)
& stop'***** n1/n2 mismatch (low) ***** '
if(iehlim1.ge.iemin0.and.iehlim1.le.iemax0)then
if (luseCompress.ne.0) then
do i1 = 1, n1max0
do i2 = 1, n2max0
do i = iemin0, iemax0
do j = 1, i
ppj(j,i,i2,i1) = CorDataNextNumber()
enddo
enddo
enddo
enddo
call CorDataCloseFile()
else
read(ifwle,*)
$ ((((ppj(j,i,i2,i1)
$ ,j=1,i)
$ ,i=iemin0,iemax0)
$ ,i2=1,n2max0)
$ ,i1=1,n1max0)
close(ifwle)
endif
read(ifwle,*)
$ ((((ppj(j,i,i2,i1)
$ ,j=1,i)
$ ,i=iemin0,iemax0)
$ ,i2=1,n2max0)
$ ,i1=1,n1max0)
elseif(iehlim1.lt.iemin0)then
stop'***** iemin0 too big *****'
elseif(iehlim1.gt.iemax0)then
stop'***** iemax0 too small *****'
endif
close(ifwle)
go=.true.
else
write(6,*)'Table wle is not defined for hadron cascade !',ppjver
......@@ -11947,17 +11908,9 @@ c Read the decay tables for kaons
c Pions and pions 0 from charged Kaons
 
if(ifdkz.gt.0)then
luseCompress = (index(fndkz(1:nfndkz), ".bz2") .eq. nfndkz-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndkz(1:nfndkz))
call CorDataFillArray(akz, size(akz))
call CorDataFillArray(akz0, size(akz0))
call CorDataCloseFile()
else
open(ifdkz,file=fndkz(1:nfndkz),status='old')
read(ifdkz,*) akz,akz0
close(ifdkz)
endif
open(ifdkz,file=fndkz(1:nfndkz),status='old')
read(ifdkz,*) akz,akz0
close(ifdkz)
else
write(6,*)'Table dkz is not defined for hadron cascade !'
endif
......@@ -11965,17 +11918,9 @@ c Pions and pions 0 from charged Kaons
c Pions and pions 0 from Kaon Long
 
if(ifdkl.gt.0)then
luseCompress = (index(fndkl(1:nfndkl), ".bz2") .eq. nfndkl-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndkl(1:nfndkl))
call CorDataFillArray(akl, size(akl))
call CorDataFillArray(akl0, size(akl0))
call CorDataCloseFile()
else
open(ifdkl,file=fndkl(1:nfndkl),status='old')
read(ifdkl,*) akl,akl0
close(ifdkl)
endif
open(ifdkl,file=fndkl(1:nfndkl),status='old')
read(ifdkl,*) akl,akl0
close(ifdkl)
else
write(6,*)'Table dkl is not defined for hadron cascade !'
endif
......@@ -11983,17 +11928,9 @@ c Pions and pions 0 from Kaon Long
c Pions and pions 0 from Kaon short
 
if(ifdks.gt.0)then
luseCompress = (index(fndks(1:nfndks), ".bz2") .eq. nfndks-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndks(1:nfndks))
call CorDataFillArray(aks, size(aks))
call CorDataFillArray(aks0, size(aks0))
call CorDataCloseFile()
else
open(ifdks,file=fndks(1:nfndks),status='old')
read(ifdks,*) aks,aks0
close(ifdks)
endif
open(ifdks,file=fndks(1:nfndks),status='old')
read(ifdks,*) aks,aks0
close(ifdks)
else
write(6,*)'Table dks is not defined for hadron cascade !'
endif
......@@ -12001,18 +11938,9 @@ c Pions and pions 0 from Kaon short
c Muons from Charge Kaon, Kaon Long and Charged Pions
 
if(ifdkm.gt.0)then
luseCompress = (index(fndkm(1:nfndkm), ".bz2") .eq. nfndkm-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndkm(1:nfndkm))
call CorDataFillArray(akzm, size(akzm))
call CorDataFillArray(aklm, size(aklm))
call CorDataFillArray(apim, size(apim))
call CorDataCloseFile()
else
open(ifdkm,file=fndkm(1:nfndkm),status='old')
read(ifdkm,*) akzm,aklm,apim
close(ifdkm)
endif
open(ifdkm,file=fndkm(1:nfndkm),status='old')
read(ifdkm,*) akzm,aklm,apim
close(ifdkm)
else
write(6,*)'Table dkm is not defined for hadron cascade !'
endif
......@@ -12020,19 +11948,9 @@ c Muons from Charge Kaon, Kaon Long and Charged Pions
c Neutrinos from Charge Kaon, Kaon Long and Charged Pions
 
if(ifdkn.gt.0)then
luseCompress = (index(fndkn(1:nfndkn), ".bz2") .eq. nfndkn-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndkn(1:nfndkn))
call CorDataFillArray(akzn, size(akzn))
call CorDataFillArray(akln, size(akln))
call CorDataFillArray(apin, size(apin))
call CorDataFillArray(amun, size(amun))
call CorDataCloseFile()
else
open(ifdkn,file=fndkn(1:nfndkn),status='old')
read(ifdkn,*) akzn,akln,apin,amun
close(ifdkn)
endif
open(ifdkn,file=fndkn(1:nfndkn),status='old')
read(ifdkn,*) akzn,akln,apin,amun
close(ifdkn)
else
write(6,*)'Table dkn is not defined for hadron cascade !'
endif
......@@ -12192,20 +12110,9 @@ c Neutrino for Edepo
c Electrons from Charge Kaon or Kaon Long and Gammas from Neutral Pions
 
if(ifdke.gt.0)then
luseCompress = (index(fndke(1:nfndke), ".bz2") .eq. nfndke-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndke(1:nfndke))
call CorDataFillArray(akze, size(akze))
call CorDataFillArray(akle, size(akle))
call CorDataFillArray(ap0g, size(ap0g))
call CorDataFillArray(ap0e, size(ap0e))
call CorDataFillArray(amue, size(amue))
call CorDataCloseFile()
else
open(ifdke,file=fndke(1:nfndke),status='old')
read(ifdke,*) akze,akle,ap0g,ap0e,amue
close(ifdke)
endif
open(ifdke,file=fndke(1:nfndke),status='old')
read(ifdke,*) akze,akle,ap0g,ap0e,amue
close(ifdke)
else
write(6,*)'Table dke is not defined for hadron cascade !'
endif
......@@ -13320,7 +13227,7 @@ c if(ebal.lt.0.d0)write(*,*)'ha',k,zha(k),ebal,enpart(1),enpart(2) !there
edep=max(0d0,ebal-edepmn)
fact=0.44d0
if(xpo.gt.0d0.or.(edep/edepmn.gt.xxx1.and.edep.gt.xxx2
. .and.edep.gt.fact*edepmn))then
. .and.edep.gt.fact*edepmn))then
xpo=1d0
c edep=max(edep,sumEloss)
else
......@@ -14487,6 +14394,49 @@ c inelastic interaction path
rlamold=rlamold/(1.d0-w)
return
end
#if __CXCORSIKA__ || !__CXSUB__
c Dummy functions when we don't want to use compressed files (old CORSIKA and CONEX)
subroutine CorDataOpenFile(name)
character*256 name
write(*,*)'WARNING with',name
write(*,*)'This version can not read compressed data file !!!'
stop
end
subroutine CorDataCloseFile()
end
subroutine CorDataCanDeCompress(idum)
integer idum,idum2
idum2=idum
end
subroutine CorDataFillArray(dum,idum)
double precision dum,dum2
dimension dum(*)
integer idum,idum2
dum2=dum(1)
idum2=idum
end
subroutine CorDataNextText(dum)
character*20 dum,dum2
dum2=dum
end
double precision function CorDataNextNumber()
end
integer function size(array)
double precision array(*)
size=int(array(1))
end
#endif
c Original File : conex-lat.F
......@@ -23373,6 +23323,7 @@ C-----------------------------------------------------------------------
 
end
 
#endif
 
C-----------------------------------------------------------------------
......@@ -31633,6 +31584,7 @@ c
c RETURN
c END
c
c Original File : conex-xan.F
......@@ -36831,16 +36783,32 @@ c Energy distribution of gammas from pi0 decay
 
c Pi0 from Charge Kaon or Kaon Long
if(ifdkz.gt.0)then
open(ifdkz,file=fndkz(1:nfndkz),status='old')
read(ifdkz,*) akz,akz0
close(ifdkz)
luseCompress = (index(fndkz(1:nfndkz), ".bz2") .eq. nfndkz-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndkz(1:nfndkz))
call CorDataFillArray(akz, size(akz))
call CorDataFillArray(akz0, size(akz0))
call CorDataCloseFile()
else
open(ifdkz,file=fndkz(1:nfndkz),status='old')
read(ifdkz,*) akz,akz0
close(ifdkz)
endif
else
write(6,*)'Table dkz is not defined for hadron cascade !'
endif
if(ifdkl.gt.0)then
open(ifdkl,file=fndkl(1:nfndkl),status='old')
read(ifdkl,*) akl,akl0
close(ifdkl)
luseCompress = (index(fndkl(1:nfndkl), ".bz2") .eq. nfndkl-3)
if (luseCompress.ne.0) then
call CorDataOpenFile(fndkl(1:nfndkl))
call CorDataFillArray(akl, size(akl))
call CorDataFillArray(akl0, size(akl0))
call CorDataCloseFile()
else
open(ifdkl,file=fndkl(1:nfndkl),status='old')
read(ifdkl,*) akl,akl0
close(ifdkl)
endif
else
write(6,*)'Table dkl is not defined for hadron cascade !'
endif
......@@ -36852,9 +36820,9 @@ c Pi0 from Charge Kaon or Kaon Long
call CorDataFillArray(aks0, size(aks0))
call CorDataCloseFile()
else
open(ifdks,file=fndks(1:nfndks),status='old')
read(ifdks,*) aks,aks0
close(ifdks)
open(ifdks,file=fndks(1:nfndks),status='old')
read(ifdks,*) aks,aks0
close(ifdks)
endif
else
write(6,*)'Table dks is not defined for hadron cascade !'

data privacy - impressum