diff --git a/src/conex_mod.F b/src/conex_mod.F index 122d893d1268c9c20175df14fa74786affe31041..a61f4abbd11f1426cf2849a84c69cf3c78a5c27c 100644 --- a/src/conex_mod.F +++ b/src/conex_mod.F @@ -871,7 +871,7 @@ C----------------------------------------------------------------------- END #endif -#ifndef __CORSIKA8__ +#if !__CXCORSIKA__ && !__CORSIKA8__ #ifdef __NEXUS__ c----------------------------------------------------------------------- @@ -48864,7 +48864,6 @@ c------------------------------------------------------------------------------ end #if !__CXCORSIKA__ && !__CORSIKA8__ - c-------------------------------------------------------------------- double precision function qgran(b10) c-------------------------------------------------------------------- @@ -48876,14 +48875,9 @@ c-------------------------------------------------------------------- return end - - - integer function size(array) - double precision array(*) - size=int(array(1)) - end #endif + #endif @@ -157633,6 +157627,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 diff --git a/src/conex_sub.F b/src/conex_sub.F index ab3e025bb050f7b367fde5390e297f0a209283e1..418dfe96bd51ee7c72a255df73d7f381c76c0172 100644 --- a/src/conex_sub.F +++ b/src/conex_sub.F @@ -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) + double precision dum,dum2 + dum2=dum + end + + 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 !'