diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index d6162179ca3440e788a40a9a6ab0c18edb7d44f6..2835a7043cc80101ef22d1301af981b5735de6bb 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -9,6 +9,7 @@ add_subdirectory (TrackWriter) add_subdirectory (TrackingLine) add_subdirectory (HadronicElasticModel) add_subdirectory (EnergyLoss) +add_subdirectory (UrQMD) #add_custom_target(CORSIKAprocesses) add_library (CORSIKAprocesses INTERFACE) @@ -20,3 +21,4 @@ endif (PYTHIA8_FOUND) add_dependencies(CORSIKAprocesses ProcessStackInspector) add_dependencies(CORSIKAprocesses ProcessTrackingLine) add_dependencies(CORSIKAprocesses ProcessEnergyLoss) +add_dependencies(CORSIKAprocesses ProcessUrQMD) diff --git a/Processes/UrQMD/CMakeLists.txt b/Processes/UrQMD/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..3ba71c209f3d3e746458c2a80e4c1623f4684981 --- /dev/null +++ b/Processes/UrQMD/CMakeLists.txt @@ -0,0 +1,94 @@ +set ( + MODEL_SOURCES + UrQMD.cc + urqmdInterface.F + addpart.f + angdis.f + anndec.f + blockres.f + boxprg.f + cascinit.f + coload.f + dectim.f + delpart.f + detbal.f + dwidth.f + error.f + getmass.f + getspin.f + init.f + iso.f + ityp2pdg.f + jdecay2.f + make22.f + numrec.f + output.f + paulibl.f + proppot.f + saveinfo.f + scatter.f + siglookup.f + string.f + tabinit.f + urqmd.f + whichres.f +) + +set ( + MODEL_HEADERS + UrQMD.h + ) + +set ( + MODEL_NAMESPACE + corsika/process/urqmd + ) + +add_library (ProcessUrQMD STATIC ${MODEL_SOURCES}) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessUrQMD ${MODEL_NAMESPACE} ${MODEL_HEADERS}) + +set_target_properties ( + ProcessUrQMD + PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION 1 +# PUBLIC_HEADER "${MODEL_HEADERS}" + ) + +# target dependencies on other libraries (also the header onlys) +target_link_libraries ( + ProcessUrQMD + CORSIKAparticles + CORSIKAunits + CORSIKAgeometry + CORSIKArandom + ) + +target_include_directories ( + ProcessUrQMD + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) + +install ( + TARGETS ProcessUrQMD + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib +# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE} + ) + + +# -------------------- +# code unit testing +add_executable (testUrQMD + testUrQMD.cc + ${MODEL_HEADERS} + ) + +target_link_libraries ( + testUrQMD + ProcessUrQMD + CORSIKAthirdparty # for catch2 + ) +CORSIKA_ADD_TEST(testUrQMD) diff --git a/Processes/UrQMD/README b/Processes/UrQMD/README new file mode 100644 index 0000000000000000000000000000000000000000..3fa900fee63555075a318e21c5ce38b2d1a3c27c --- /dev/null +++ b/Processes/UrQMD/README @@ -0,0 +1,49 @@ +additional files needed for building (from UrQMD v1.3cr): + + Copyright + README + UrQMD-1.3.1-xs.dat + addpart.f + angdis.f + anndec.f + blockres.f + boxinc.f + boxprg.f + cascinit.f + colltab.f + coload.f + comnorm.f + comres.f + coms.f + comstr.f + comwid.f + conex.h + conex.incnex + dectim.f + delpart.f + detbal.f + dwidth.f + error.f + freezeout.f + getmass.f + getspin.f + init.f + inputs.f + iso.f + ityp2pdg.f + jdecay2.f + make22.f + newpart.f + numrec.f + options.f + outcom.f + output.f + paulibl.f + proppot.f + saveinfo.f + scatter.f + siglookup.f + string.f + tabinit.f + urqmd.f + whichres.f diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc new file mode 100644 index 0000000000000000000000000000000000000000..92e5d53e4fba0e6a18479e0479dcfd346eb67cdf --- /dev/null +++ b/Processes/UrQMD/UrQMD.cc @@ -0,0 +1,15 @@ +#include <corsika/process/urqmd/UrQMD.h> +#include <corsika/random/RNGManager.h> +#include <random> + +corsika::process::UrQMD::UrQMD::UrQMD() { + iniurqmd_(); +} + +double ranf_(int*) { + static corsika::random::RNG& rng = + corsika::random::RNGManager::GetInstance().GetRandomStream("ranf"); + + std::uniform_real_distribution<double> dist; + return dist(rng); +} diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h new file mode 100644 index 0000000000000000000000000000000000000000..83d0d5c42416fa766bcde00f5f431701796306a0 --- /dev/null +++ b/Processes/UrQMD/UrQMD.h @@ -0,0 +1,16 @@ +#ifndef _Processes_UrQMD_UrQMD_h +#define _Processes_UrQMD_UrQMD_h + +extern "C" { + void iniurqmd_(); + double ranf_(int*); +} + +namespace corsika::process::UrQMD { + class UrQMD { + public: + UrQMD(); + }; +} + +#endif diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc new file mode 100644 index 0000000000000000000000000000000000000000..d50b88d41dbc292969dd6884e92744d718b9bf13 --- /dev/null +++ b/Processes/UrQMD/testUrQMD.cc @@ -0,0 +1,21 @@ +/* + * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <corsika/process/urqmd/UrQMD.h> + +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one + // cpp file +#include <catch2/catch.hpp> + +using namespace corsika::process::UrQMD; + +TEST_CASE("UrQMD") { + UrQMD proc; +} diff --git a/Processes/UrQMD/urqmdInterface.F b/Processes/UrQMD/urqmdInterface.F new file mode 100644 index 0000000000000000000000000000000000000000..f22e5d94d7d83f490dabf3f6433a98e918fe30db --- /dev/null +++ b/Processes/UrQMD/urqmdInterface.F @@ -0,0 +1,433 @@ +c 18.11.2011 Link routines between UrQMD 1.3 and CONEX. +c author T. Pierog based on CORSIKA and EPOS link to UrQMD + +c adapted by M. Reininghaus for linking UrQMD to CORSIKA 8 (Apr 2019) + +#ifdef __STD__ +#define __GHEISHA__ +#define __QGSJET__ +#define __ANALYSIS__ +#endif + + +c----------------------------------------------------------------------- + subroutine IniUrQMD +c----------------------------------------------------------------------- +c Primary initialization for UrQMD 1.31 +c----------------------------------------------------------------------- + implicit none +c CONEX includes +#include "conex.h" +#include "conex.incnex" +#ifndef __CXCORSIKA__ + character*500 furqdat + integer ifurqdat, nfurqdat + common/urqfname/ furqdat, ifurqdat, nfurqdat + + include 'boxinc.f' + include 'inputs.f' + include 'options.f' + +c commons from coms.f + integer Ap, At, Zp, Zt, npart, nbar, nmes, ctag + integer nsteps,ranseed,event,eos,dectag,uid_cnt + integer NHardRes,NSoftRes,NDecRes,NElColl,NBlColl + common /sys/ npart, nbar, nmes, ctag,nsteps,uid_cnt, + + ranseed,event,Ap,At,Zp,Zt,eos,dectag, + + NHardRes,NSoftRes,NDecRes,NElColl,NBlColl + +c local + INTEGER i,io,ia,ie,id + CHARACTER CTPStrg(numctp)*60, CTOStrng(numcto)*60 + integer mxie,mxid,mxia + parameter (mxie=41,mxid=10,mxia=3) + character adum + double precision sig_u1,ekdummy + integer iamaxu,idmaxu,iemaxu + common /cxs_u1/ sig_u1(mxie,mxid,mxia),iamaxu,idmaxu,iemaxu + double precision xs(3),bim(3) + common /cxs_u2/ xs + integer iudebug + data bim/6.d0,6.d0,7.d0/ + integer init + data init/0/ + SAVE + + if(init.ge.1)return + init=init+1 +#ifdef __CXDEBUG__ + call utisx1('iniurqmd ',4) + write(*,'(a)')'initialize URQMD ...' +#endif + +C----------------------------------------------------------------------- + + IF ( isx.ge.2 ) THEN + IUDEBUG = isx-1 + ELSE + IUDEBUG = 0 + ENDIF + + WRITE (*,*) + $ '############################################################' + WRITE (*,*) + $ '## ##' + WRITE (*,*) + $ '## UrQMD 1.3.1 University of Frankfurt ##' + WRITE (*,*) + $ '## urqmd@th.physik.uni-frankfurt.de ##' + WRITE (*,*) + $ '## ##' + WRITE (*,*) + $ '############################################################' + WRITE (*,*) + $ '## ##' + WRITE (*,*) + $ '## please cite when using this model: ##' + WRITE (*,*) + $ '## S.A.Bass et al. Prog.Part.Nucl.Phys. 41 (1998) 225 ##' + WRITE (*,*) + $ '## M.Bleicher et al. J.Phys. G25 (1999) 1859 ##' + WRITE (*,*) + $ '## ##' + WRITE (*,*) + $ '############################################################' + +C SET THE 'LARGE' CROSS-SECTIONS FOR ALL 3 TARGET ELEMENTS + DO I = 1, 3 + XS(I) = 10.D0 * PI * BIM(I)**2 + ENDDO + +C SET NMAX TO DEFAULT VALUE + call set0 + call params + +C THIS IS THE SUBSTITUE FOR THE URQMD INPUT ROUTINE +C INITIALIZE COUNTERS + boxflag = 0 + mbflag = 0 + edens = 0.d0 + para = 0 + solid = 0 + mbox = 0 + io = 0 + +C THE FOLLOWING FLAGS CHECK, WHETHER ALL NECESSARY INPUT IS GIVEN +C PROJECTILE + prspflg = 0 +C TARGET + trspflg = 0 +C + srtflag = 0 + firstev = 0 +C EXCITATION FUNCTION + nsrt = 1 + npb = 1 + efuncflag = 0 +C DEFAULT NUMBER OF EVENTS + nevents = 1 +C DEFAULT NUMBER OF TIMESTEPS + nsteps = 1000 + +C SKIP CONDITIONS ON UNIT 13, 14, 15, 16 & 18 +C SUPPRESS ALL OUTPUT + bf13 = .true. + bf14 = .true. + bf15 = .true. + bf16 = .true. + bf18 = .true. + bf19 = .true. + bf20 = .true. +C SET DEBUG OUTPUT DEPENDING ON CHOSEN DEBUG LEVEL +C SET THE OUTPUT OF UNITS 13, 14, 15 TO THE DEBUG OUTPUT UNIT + IF ( IUDEBUG .EQ. 1 ) THEN + bf13 = .true. + bf14 = .false. + call uounit(14,IFCK) + bf15 = .true. + ELSEIF ( IUDEBUG .EQ. 2 ) THEN + bf13 = .false. + call uounit(13,IFCK) + bf14 = .true. + bf15 = .true. + ELSEIF ( IUDEBUG .GT. 2 ) THEN + bf13 = .true. + bf14 = .true. + bf15 = .false. + call uounit(15,IFCK) + ENDIF + do i = 1, numcto + CTOdc(i) = ' ' + enddo + do i = 1, numctp + CTPdc(i) = ' ' + enddo + do i = 1, maxstables + stabvec(i) = 0 + enddo + nstable = 0 + +C DEFAULT SETTINGS FOR CTParam AND CTOption +C DEFAULT SETTINGS FOR CTParam + CTParam(1)=1.d0 + CTPStrg(1)='scaling factor for decay-width' + CTParam(2)=0.52d0 + CTPStrg(2)='used for minimal stringmass & el/inel cut in makestr' + CTParam(3)=2.d0 + CTPStrg(3)='velocity exponent for modified AQM' + CTParam(4)=0.3d0 + CTPStrg(4)='transverse pion mass, used in make22 & strexct' + CTParam(5)=0.d0 + CTPStrg(5)='probabil. for quark rearrangement in cluster' + CTParam(6)=0.37d0 + CTPstrg(6)='strangeness probability' + CTParam(7)=0.d0 + CTPStrg(7)='charm probability (not yet implemented in UQMD)' + CTParam(8)=0.093d0 + CTPStrg(8)='probability to create a diquark' + CTParam(9)=0.35d0 + CTPStrg(9)='kinetic energy cut off for last string break' + CTParam(10)=0.25d0 + CTPStrg(10)='min. kinetic energy for hadron in string' + CTParam(11)=0.d0 + CTPStrg(11)='fraction of non groundstate resonances' + CTParam(12)=.5d0 + CTPStrg(12)='probability for rho 770 in String' + CTParam(13)=.27d0 + CTPStrg(13)='probability for rho 1450 (rest->rho1700)' + CTParam(14)=.49d0 + CTPStrg(14)='probability for omega 782' + CTParam(15)=.27d0 + CTPStrg(15)='probability for omega 1420(rest->om1600)' + CTParam(16)=1.0d0 + CTPStrg(16)='mass cut betw. rho770 and rho 1450' + CTParam(17)=1.6d0 + CTPSTRG(17)='mass cut betw. rho1450 and rho1700' + CTParam(18)=.85d0 + CTPStrg(18)='mass cut betw. om 782 and om1420' + CTParam(19)=1.55d0 + CTPStrg(19)='mass cut betw. om1420 and om1600' + CTParam(20)=0.0d0 + CTPStrg(20)=' distance for second projectile' + CTParam(21)=0.0d0 + CTPStrg(21)=' deformation parameter' + + CTParam(25)=.9d0 + CTPStrg(25)=' probability for diquark not to break' + CTParam(26)=50.d0 + CTPStrg(26)=' maximum trials to get string masses' + CTParam(27)=1.d0 + CTPStrg(27)=' scaling factor for xmin in string excitation' + CTParam(28)=1.d0 + CTPStrg(28)=' scaling factor for transverse fermi motion' + CTParam(29)=0.4d0 + CTPStrg(29)=' single strange di-quark suppression factor ' + CTParam(30)=1.5d0 + CTPStrg(30)=' radius offset for initialization ' + CTParam(31)=1.6d0 + CTPStrg(31)=' sigma of gaussian for tranverse momentum tranfer ' + CTParam(32)=0.d0 + CTPStrg(32)=' alpha-1 for valence quark distribution ' + CTParam(33)=2.5d0 + CTPStrg(33)=' betav for valence quark distribution (DPM)' + CTParam(34)=0.1d0 + CTPStrg(34)=' minimal x multiplied with ecm ' + CTParam(35)=3.0d0 + CTPStrg(35)=' offset for cut for the FSM ' + CTParam(36)=0.275d0 + CTPStrg(36)=' fragmentation function parameter a ' + CTParam(37)=0.42d0 + CTPStrg(37)=' fragmentation function parameter b ' + CTParam(38)=1.08d0 + CTPStrg(38)=' diquark pt scaling factor ' + CTParam(39)=0.8d0 + CTPStrg(39)=' strange quark pt scaling factor ' + CTParam(40)=0.5d0 + CTPStrg(40)=' betas-1 for valence quark distribution (LEM)' + CTParam(41)=0.d0 + CTPStrg(41)=' distance of initialization' + CTParam(42)=0.55d0 + CTPStrg(42)=' width of gaussian -> pt in string-fragmentation ' + CTParam(43)=5.d0 + CTPStrg(43)=' maximum kinetic energy in mesonic clustr ' + CTParam(44)=0.8d0 + CTPStrg(44)=' prob. of double vs. single excitation for AQM inel.' + CTParam(45)=0.5d0 + CTPStrg(45)=' offset for minimal mass generation of strings' + CTParam(46)=800000.d0 + CTPStrg(46)=' maximal number of rejections for initialization' + CTParam(47)=1.0d0 + CTPStrg(47)=' field feynman fragmentation funct. param. a' + CTParam(48)=2.0d0 + CTPStrg(48)=' field feynman fragmentation funct. param. b' + + CTParam(50)=1.d0 + CTPStrg(50)=' enhancement factor for 0- mesons' + CTParam(51)=1.d0 + CTPStrg(51)=' enhancement factor for 1- mesons' + CTParam(52)=1.d0 + CTPStrg(52)=' enhancement factor for 0+ mesons' + CTParam(53)=1.d0 + CTPStrg(53)=' enhancement factor for 1+ mesons' + CTParam(54)=1.d0 + CTPStrg(54)=' enhancement factor for 2+ mesons' + CTParam(55)=1.d0 + CTPStrg(55)=' enhancement factor for 1+-mesons' + CTParam(56)=1.d0 + CTPStrg(56)=' enhancement factor for 1-*mesons' + CTParam(57)=1.d0 + CTPStrg(57)=' enhancement factor for 1-*mesons' + CTParam(58)=1.d0 + CTPStrg(58)=' scaling factor for DP time-delay' + +C DEFAULT SETTINGS FOR CTOption + CTOption(1)=1 ! hjd1 + CTOStrng(1)=' resonance widths are mass dependent ' + CTOption(2)=0 + CTOStrng(2)=' conservation of scattering plane' + CTOption(3)=0 + CTOStrng(3)=' use modified detailed balance' + CTOption(4)=0 + CTOStrng(4)=' no initial conf. output ' + CTOption(5)=0 + CTOStrng(5)=' fixed random impact parameter' + CTOption(6)=0 + CTOStrng(6)=' no first collisions inside proj/target' + CTOption(7)=0 + CTOStrng(7)=' elastic cross-section enabled (<>0:total=inelast)' + CTOption(8)=0 + CTOStrng(8)=' extrapolate branching ratios ' + CTOption(9)=0 + CTOStrng(9)=' use tabulated pp cross-sections ' + CTOption(10)=0 + CTOStrng(10)=' enable Pauli Blocker' + CTOption(11)=0 + CTOStrng(11)=' mass reduction for cascade initialization' + CTOption(12)=0 + CTOStrng(12)=' string condition =0 (.ne.0 no strings)' + CTOption(13)=0 + CTOStrng(13)=' enhanced file16 output ' + CTOption(14)=0 + CTOStrng(14)=' cos(the) is distributet between -1..1 ' + CTOption(15)=0 + CTOStrng(15)=' allow mm&mb-scattering' + CTOption(16)=0 + CTOStrng(16)=' propagate without collisions' + CTOption(17)=0 + CTOStrng(17)=' colload after every timestep ' + CTOption(18)=0 + CTOStrng(18)=' final decay of unstable particles' + CTOption(19)=0 + CTOStrng(19)=' allow bbar annihilaion' + CTOption(20)=0 + CTOStrng(20)=' dont generate e+e- instead of bbar' + CTOption(21)=0 + CTOStrng(21)=' use field feynman frgm. function' + CTOption(22)=1 + CTOStrng(22)=' use lund excitation function' + CTOption(23)=0 + CTOStrng(23)=' lorentz contraction of projectile & targed' + CTOption(24)=2 ! 1 is default 2 means fast method + CTOStrng(24)=' Wood-Saxon initialization' + CTOption(25)=0 + CTOStrng(25)=' phase space corrections for resonance mass' + CTOption(26)=0 + CTOStrng(26)=' use z -> 1-z for diquark-pairs' + CTOption(27)=1 ! hjd1 + CTOStrng(27)=' reference frame (1=target, 2=projectile, else=cms)' + CTOption(28)=0 + CTOStrng(28)=' propagate spectators also ' + CTOption(29)=2 + CTOStrng(29)=' no transverse momentum in clustr ' + CTOption(30)=1 + CTOStrng(30)=' frozen fermi motion ' + CTOption(31)=0 + CTOStrng(31)=' reduced mass spectrum in string' + CTOption(32)=0 + CTOStrng(32)=' masses are distributed acc. to m-dep. widths' + CTOption(33)=0 + CTOStrng(33)=' use tables & m-dep. for pmean in fprwdt & fwidth' + CTOption(34)=1 + CTOStrng(34)=' lifetme according to m-dep. width' + CTOption(35)=1 + CTOStrng(35)=' generate high precision tables' + CTOption(36)=0 + CTOStrng(36)=' normalize Breit-Wigners with m.dep. widths ' + CTOption(37)=0 + CTOStrng(37)=' heavy quarks form di-quark clusters' + CTOption(38)=0 + CTOStrng(38)=' scale p-pbar to b-bbar with equal p_lab ' + CTOption(39)=0 + CTOStrng(39)=' dont call pauliblocker' + CTOption(40)=0 + CTOStrng(40)=' read old fort.14 file ' + CTOption(41)=0 + CTOStrng(41)=' generate extended output for cto40' + CTOption(42)=0 + CTOStrng(42)=' hadrons now have color fluctuations' + CTOption(43)=0 + CTOStrng(43)=' dont generate dimuon intead of dielectron output' + CTOption(44)=0 + CTOStrng(44)=' not used at the moment' + CTOption(45)=0 + CTOStrng(45)=' not used at the moment' + +C INITIALIZE ARRAYS FOR SPECIAL PRO/TAR COMBINATIONS + do i = 1, 2 + spityp(i) = 0 + spiso3(i) = 0 + enddo + +C INITIALIZE ARRAYS FOR SPECIAL PARTICLES + EoS = 0 + +C READ CROSS-SECTION FILES +Cdh CALL URQREC() + +C INITIALIZES SOME ARRAYS + call strini ! initialize mixing angles for meson-multipletts + call loginit + + IF ( CTOption(33) .EQ. 0 .OR. CTOption(9) .EQ. 0 ) THEN + call loadwtab(io) + IF ( IUDEBUG .GT. 0 ) WRITE(IFCK,*) 'URQINI: AFTER LOADWTAB' + ENDIF + +C READ URQMD TOTAL CROSS SECTION TABLE +c +c ie=1..41 E=10.0**(float(ie)/10-1.0-0.05) (bin-middle) +c id=1..9 p,ap,n,an,pi+,pi-,K+,K-,KS +c ia=1..3 N,O,Ar +c + if(ifurqdat.eq.1)then + OPEN(UNIT=76,FILE=furqdat(1:nfurqdat),STATUS='OLD') + else + OPEN(UNIT=76,FILE='UrQMD-1.3.1-xs.dat',STATUS='OLD') + endif + read(76,*) adum,iamaxu,idmaxu,iemaxu + do ia=1,iamaxu + do id=1,idmaxu + do ie=1,iemaxu + read(76,*) ekdummy,sig_u1(ie,id,ia) + enddo + read(76,*) + read(76,*) + enddo + enddo + close(76) + +C IN CASE OF CASCADE MODE, THE POTENTIALS NEED NOT BE CALCULATED + +C CALCULATE NORMALIZATION OF RESONANCES DISTRIBUTION... + call norm_init +#endif + + + xsegymin=0.25d0 + +#ifdef __CXDEBUG__ + call utisx2 +#endif + + end +