generators/JEWEL/medium-simple.f

0001 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0002 C++ Copyright (C) 2017 Korinna C. Zapp [Korinna.Zapp@cern.ch]       ++
0003 C++                                                                 ++
0004 C++ This file is part of JEWEL 2.2.0                                ++
0005 C++                                                                 ++
0006 C++ The JEWEL homepage is jewel.hepforge.org                        ++
0007 C++                                                                 ++
0008 C++ The medium model was partly implemented by Jochen Klein.        ++
0009 C++ Raghav Kunnawalkam Elayavalli helped with the implementation    ++
0010 C++ of the V+jet processes.                                         ++
0011 C++                                                                 ++
0012 C++ Please follow the MCnet GUIDELINES and cite Eur.Phys.J. C74     ++
0013 C++ (2014) no.2, 2762 [arXiv:1311.0048] for the code and            ++
0014 C++ JHEP 1303 (2013) 080 [arXiv:1212.1599] and                      ++
0015 C++ optionally EPJC 60 (2009) 617 [arXiv:0804.3568] for the         ++
0016 C++ physics. The reference for V+jet processes is EPJC 76 (2016)    ++
0017 C++ no.12 695 [arXiv:1608.03099] and for recoil effects it is       ++
0018 C++ arXiv:1707.01539.
0019 C++                                                                 ++
0020 C++ JEWEL relies heavily on PYTHIA 6 for the event generation. The  ++
0021 C++ modified version of PYTHIA 6.4.25 that is distributed with      ++
0022 C++ JEWEL is, however, not an official PYTHIA release and must not  ++
0023 C++ be used for anything else. Please refer to results as           ++
0024 C++ "JEWEL+PYTHIA".                                                 ++
0025 C++                                                                 ++
0026 C++ JEWEL also uses code provided by S. Zhang and J. M. Jing        ++
0027 C++ (Computation of Special Functions, John Wiley & Sons, New York, ++
0028 C++ 1996 and http://jin.ece.illinois.edu) for computing the         ++
0029 C++ exponential integral Ei(x).                                     ++
0030 C++                                                                 ++
0031 C++                                                                 ++
0032 C++ JEWEL  is free software; you can redistribute it and/or         ++
0033 C++ modify it under the terms of the GNU General Public License     ++
0034 C++ as published by the Free Software Foundation; either version 2  ++
0035 C++ of the License, or (at your option) any later version.          ++
0036 C++                                                                 ++
0037 C++ JEWEL is distributed in the hope that it will be useful,        ++
0038 C++ but WITHOUT ANY WARRANTY; without even the implied warranty of  ++
0039 C++ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the    ++
0040 C++ GNU General Public License for more details.                    ++
0041 C++                                                                 ++
0042 C++ You should have received a copy of the GNU General Public       ++
0043 C++ License along with this program; if not, write to the Free      ++
0044 C++ Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, ++
0045 C++ MA 02110-1301 USA                                               ++
0046 C++                                                                 ++
0047 C++ Linking JEWEL statically or dynamically with other modules is   ++
0048 C++ making a combined work based on JEWEL. Thus, the terms and      ++
0049 C++ conditions of the GNU General Public License cover the whole    ++
0050 C++ combination.                                                    ++
0051 C++                                                                 ++
0052 C++ In addition, as a special exception, I give you permission to   ++
0053 C++ combine JEWEL with the code for the computation of special      ++
0054 C++ functions provided by S. Zhang and J. M. Jing. You may copy and ++
0055 C++ distribute such a system following the terms of the GNU GPL for ++
0056 C++ JEWEL and the licenses of the other code concerned, provided    ++
0057 C++ that you include the source code of that other code when and as ++
0058 C++ the GNU GPL requires distribution of source code.               ++
0059 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0060
0061       SUBROUTINE MEDINIT(FILE,id,etam,mass)
0062       IMPLICIT NONE
0063 C--medium parameters
0064       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0065       INTEGER NF
0066       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0067       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0068      &N0,SIGMANN,A,WOODSSAXON
0069       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0070      &SIGMANN
0071       INTEGER A
0072       LOGICAL WOODSSAXON
0073 C--max rapidity
0074         common/rapmax2/etamax2
0075         double precision etamax2
0076 C--longitudinal boost of momentum distribution
0077         common/boostmed/boost
0078         logical boost
0079 C--factor to vary Debye mass
0080         COMMON/MDFAC/MDFACTOR,MDSCALEFAC
0081         DOUBLE PRECISION MDFACTOR,MDSCALEFAC
0082 C--nuclear thickness function
0083       COMMON /THICKFNC/ RMAX,TA(100,2)
0084       DOUBLE PRECISION RMAX,TA
0085 C--geometrical cross section
0086       COMMON /CROSSSEC/ IMPMAX,CROSS(200,3)
0087       DOUBLE PRECISION IMPMAX,CROSS
0088 C--identifier of log file
0089         common/logfile/logfid
0090         integer logfid
0091
0092       DATA RAU/10./
0093       DATA D3/0.9d0/
0094       DATA ZETA3/1.2d0/
0095 C--local variables
0096       INTEGER I,LUN,POS,IOS,id,mass
0097         double precision etam
0098       CHARACTER*100 BUFFER,LABEL,tempbuf
0099         CHARACTER*80 FILE
0100         character firstchar
0101         logical fileexist
0102
0103         etamax2 = etam
0104         logfid = id
0105
0106       IOS=0
0107       LUN=77
0108
0109 C--default settings
0110       TAUI=0.6d0
0111       TI=0.36d0
0112       TC=0.17d0
0113       WOODSSAXON=.TRUE.
0114       CENTRMIN=0.d0
0115       CENTRMAX=10.d0
0116       NF=3
0117       A=mass
0118       N0=0.17d0
0119       D=0.54d0
0120       SIGMANN=6.2
0121         MDFACTOR=0.45d0
0122         MDSCALEFAC=0.9d0
0123         boost = .true.
0124
0125 C--read settings from file
0126         write(logfid,*)
0127         inquire(file=FILE,exist=fileexist)
0128         if(fileexist)then
0129         write(logfid,*)'Reading medium parameters from ',FILE
0130         OPEN(unit=LUN,file=FILE,status='old',err=10)
0131           do 20 i=1,1000
0132           READ(LUN, '(A)', iostat=ios) BUFFER
0133             if (ios.ne.0) goto 30
0134             firstchar = buffer(1:1)
0135             if (firstchar.eq.'#') goto 20
0136           POS=SCAN(BUFFER,' ')
0137           LABEL=BUFFER(1:POS)
0138           BUFFER=BUFFER(POS+1:)
0139           IF (LABEL=="TAUI")THEN
0140             READ(BUFFER,*,IOSTAT=IOS) TAUI
0141           ELSE IF (LABEL=="TI") THEN
0142             READ(BUFFER,*,IOSTAT=IOS) TI
0143           ELSE IF (LABEL=="TC") THEN
0144             READ(BUFFER,*,IOSTAT=IOS) TC
0145           ELSE IF (LABEL=="WOODSSAXON") THEN
0146             READ(BUFFER,*,IOSTAT=IOS) WOODSSAXON
0147           ELSE IF (LABEL=="CENTRMIN") THEN
0148             READ(BUFFER,*,IOSTAT=IOS) CENTRMIN
0149           ELSE IF (LABEL=="CENTRMAX") THEN
0150             READ(BUFFER,*,IOSTAT=IOS) CENTRMAX
0151           ELSE IF (LABEL=="NF") THEN
0152             READ(BUFFER,*,IOSTAT=IOS) NF
0153           ELSE IF (LABEL=="N0") THEN
0154             READ(BUFFER,*,IOSTAT=IOS) N0
0155           ELSE IF (LABEL=="D") THEN
0156             READ(BUFFER,*,IOSTAT=IOS) D
0157           ELSE IF (LABEL=="SIGMANN") THEN
0158             READ(BUFFER,*,IOSTAT=IOS) SIGMANN
0159           ELSE IF (LABEL=="MDFACTOR") THEN
0160             READ(BUFFER,*,IOSTAT=IOS) MDFACTOR
0161           ELSE IF (LABEL=="MDSCALEFAC") THEN
0162             READ(BUFFER,*,IOSTAT=IOS) MDSCALEFAC
0163             else
0164               write(logfid,*)'unknown label ',label
0165             endif
0166  20       continue
0167
0168  30       close(LUN,status='keep')
0169           write(logfid,*)'...done'
0170           goto 40
0171
0172  10     write(logfid,*)'Could not open medium parameter file, '//
0173      &  'will run with default settings.'
0174
0175         else
0176           write(logfid,*)'No medium parameter file found, '//
0177      &  'will run with default settings.'
0178         endif
0179
0180  40   write(logfid,*)'using parameters:'
0181       write(logfid,*)'TAUI       =',TAUI
0182       write(logfid,*)'TI         =',TI
0183       write(logfid,*)'TC         =',TC
0184       write(logfid,*)'WOODSSAXON =',WOODSSAXON
0185       write(logfid,*)'CENTRMIN   =',CENTRMIN
0186       write(logfid,*)'CENTRMAX   =',CENTRMAX
0187       write(logfid,*)'NF         =',NF
0188       write(logfid,*)'A          =',A
0189       write(logfid,*)'N0         =',N0
0190       write(logfid,*)'D          =',D
0191       write(logfid,*)'SIGMANN    =',SIGMANN
0192       write(logfid,*)'MDFACTOR   =',MDFACTOR
0193       write(logfid,*)'MDSCALEFAC =',MDSCALEFAC
0194         write(logfid,*)
0195         write(logfid,*)
0196         write(logfid,*)
0197
0198 C--calculate T_A(x,y)
0199       CALL CALCTA
0200 C--calculate geometrical cross section
0201       CALL CALCXSECTION
0202
0203       END
0204
0205
0206
0207       SUBROUTINE MEDNEXTEVT
0208       IMPLICIT NONE
0209 C--medium parameters
0210       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0211       INTEGER NF
0212       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0213       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0214      &N0,SIGMANN,A,WOODSSAXON
0215       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0216      &SIGMANN
0217       INTEGER A
0218       LOGICAL WOODSSAXON
0219 C--geometrical cross section
0220       COMMON /CROSSSEC/ IMPMAX,CROSS(200,3)
0221       DOUBLE PRECISION IMPMAX,CROSS
0222 C--local variables
0223       integer i,j
0224       DOUBLE PRECISION PYR,R,b1,b2,gettemp
0225
0226 C--pick an impact parameter
0227       r=(pyr(0)*(centrmax-centrmin)+centrmin)/100.
0228       i=0
0229       do 130 j=1,200
0230        if ((r-cross(j,3)/cross(200,3)).ge.0.) then
0231         i=i+1
0232        else
0233         goto 132
0234        endif
0235  130  continue
0236  132  continue
0237       b1 = (i-1)*0.1d0
0238       b2 = i*0.1d0
0239       breal = (b2*(cross(i,3)/cross(200,3)-r)
0240      &      +b1*(r-cross(i+1,3)/cross(200,3)))/
0241      &  (cross(i,3)/cross(200,3)-cross(i+1,3)/cross(200,3))
0242       centr = r;
0243       END
0244
0245       double precision function getcentrality()
0246       implicit none
0247       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0248       INTEGER NF
0249       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0250       getcentrality=centr
0251       end
0252
0253
0254
0255       SUBROUTINE PICKVTX(X,Y)
0256       IMPLICIT NONE
0257       DOUBLE PRECISION X,Y
0258 C--medium parameters
0259       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0260       INTEGER NF
0261       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0262 C--local variables
0263       DOUBLE PRECISION X1,X2,Y1,Y2,Z,XVAL,YVAL,ZVAL,NTHICK,PYR
0264
0265       X1=BREAL/2.-RAU
0266       X2=RAU-BREAL/2.
0267       Y1=-SQRT(4*RAU**2-BREAL**2)/2.
0268       Y2=SQRT(4*RAU**2-BREAL**2)/2.
0269  131  XVAL=PYR(0)*(X2-X1)+X1
0270       YVAL=PYR(0)*(Y2-Y1)+Y1
0271       IF((NTHICK(XVAL-BREAL/2.,YVAL).EQ.0.d0).OR.
0272      &     NTHICK(XVAL+BREAL/2.,YVAL).EQ.0.d0) GOTO 131
0273       ZVAL=PYR(0)*NTHICK(-BREAL/2.,0d0)*NTHICK(BREAL/2.,0d0)
0274       Z=NTHICK(XVAL-BREAL/2.,YVAL)*NTHICK(XVAL+BREAL/2.,YVAL)
0275       IF(ZVAL.GT.Z) GOTO 131
0276       X=XVAL
0277       Y=YVAL
0278       END
0279
0280         SUBROUTINE SETB(BVAL)
0281         IMPLICIT NONE
0282 C--medium parameters
0283       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0284       INTEGER NF
0285       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0286         DOUBLE PRECISION BVAL
0287         BREAL=BVAL
0288         END
0289
0290
0291
0292       SUBROUTINE GETSCATTERER(X,Y,Z,T,TYPE,PX,PY,PZ,E,MS)
0293       IMPLICIT NONE
0294 C--medium parameters
0295       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0296       INTEGER NF
0297       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0298 C--internal medium parameters
0299       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0300      &N0,SIGMANN,A,WOODSSAXON
0301       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0302      &SIGMANN
0303       INTEGER A
0304       LOGICAL WOODSSAXON
0305 C--longitudinal boost of momentum distribution
0306         common/boostmed/boost
0307         logical boost
0308 C--function calls
0309       DOUBLE PRECISION GETTEMP,GETMD,GETMOM,GETMS
0310 C--identifier of log file
0311         common/logfile/logfid
0312         integer logfid
0313 C--local variables
0314       DOUBLE PRECISION X,Y,Z,T,MS,PX,PY,PZ,E,MD,TEMP
0315       INTEGER TYPE
0316       DOUBLE PRECISION R,PYR,pmax,wt,tau,theta,phi,pi,p,ys,pz2,e2
0317       DATA PI/3.141592653589793d0/
0318
0319       R=PYR(0)
0320       IF(R.LT.(2.*12.*NF*D3/3.)/(2.*12.*NF*D3/3.+3.*16.*ZETA3/2.))THEN
0321          TYPE=2
0322       ELSE
0323          TYPE=21
0324       ENDIF
0325       MS=GETMS(X,Y,Z,T)
0326       MD=GETMD(X,Y,Z,T)
0327       TEMP=GETTEMP(X,Y,Z,T)
0328         tau=sqrt(t**2-z**2)
0329         if (boost) then
0330           ys = 0.5*log((t+z)/(t-z))
0331         else
0332           ys = 0.d0
0333         endif
0334         pmax = 10.*temp
0335
0336       IF(TEMP.LT.1.D-2)THEN
0337        write(logfid,*)'asking for a scattering centre without medium:'
0338        write(logfid,*)'at (x,y,z,t)=',X,Y,Z,T
0339        write(logfid,*)'making one up to continue but '//
0340      &  'something is wrong!'
0341        TYPE=21
0342        PX=0.d0
0343        PY=0.d0
0344        PZ=0.d0
0345        MS=GETMS(0.d0,0.d0,0.d0,0.d0)
0346        MD=GETMD(0.d0,0.d0,0.d0,0.d0)
0347        E=SQRT(PX**2+PY**2+PZ**2+MS**2)
0348        RETURN
0349       ENDIF
0350
0351  10     p = pyr(0)**0.3333333*pmax
0352         E2 = sqrt(p**2+ms**2)
0353         if (type.eq.2) then
0354           wt = (exp(ms/temp)-1.)/(exp(E2/temp)-1.)
0355         else
0356           wt = (exp(ms/temp)+1.)/(exp(E2/temp)+1.)
0357         endif
0358         if (wt.gt.1.) write(logfid,*)'Error in getscatterer: weight = ',wt
0359         if (wt.lt.0.) write(logfid,*)'Error in getscatterer: weight = ',wt
0360         if (pyr(0).gt.wt) goto 10
0361         phi = pyr(0)*2.*pi
0362         theta = -acos(2.*pyr(0)-1.)+pi
0363         px  = p*sin(theta)*cos(phi)
0364         py  = p*sin(theta)*sin(phi)
0365         pz2 = p*cos(theta)
0366         E   = cosh(ys)*E2 + sinh(ys)*pz2
0367         pz  = sinh(ys)*E2 + cosh(ys)*pz2
0368       END
0369
0370
0371       SUBROUTINE AVSCATCEN(X,Y,Z,T,PX,PY,PZ,E,m)
0372       IMPLICIT NONE
0373 C--longitudinal boost of momentum distribution
0374         common/boostmed/boost
0375         logical boost
0376 C--max rapidity
0377         common/rapmax2/etamax2
0378         double precision etamax2
0379 C--local variables
0380         double precision x,y,z,t,px,py,pz,e,getms,m,ys
0381         if (boost) then
0382           ys = 0.5*log((t+z)/(t-z))
0383           if ((z.eq.0.d0).and.(t.eq.0.d0)) ys =0.d0
0384           if (ys.gt.etamax2) ys=etamax2
0385           if (ys.lt.-etamax2) ys=-etamax2
0386         else
0387           ys = 0.d0
0388         endif
0389         m  = getms(x,y,z,t)
0390         e  = m*cosh(ys)
0391         px = 0.d0
0392         py = 0.d0
0393         pz = m*sinh(ys)
0394         end
0395
0396
0397       SUBROUTINE maxscatcen(PX,PY,PZ,E,m)
0398       IMPLICIT NONE
0399 C--longitudinal boost of momentum distribution
0400         common/boostmed/boost
0401         logical boost
0402 C--max rapidity
0403         common/rapmax2/etamax2
0404         double precision etamax2
0405 C--local variables
0406         double precision px,py,pz,e,getmsmax,m,ys
0407         if (boost) then
0408           ys = etamax2
0409         else
0410           ys = 0.d0
0411         endif
0412         m  = getmsmax()
0413         e  = m*cosh(ys)
0414         px = 0.d0
0415         py = 0.d0
0416         pz = m*sinh(ys)
0417         end
0418
0419
0420
0421       DOUBLE PRECISION FUNCTION GETMD(X1,Y1,Z1,T1)
0422       IMPLICIT NONE
0423 C--factor to vary Debye mass
0424         COMMON/MDFAC/MDFACTOR,MDSCALEFAC
0425         DOUBLE PRECISION MDFACTOR,MDSCALEFAC
0426       DOUBLE PRECISION X1,Y1,Z1,T1,GETTEMP
0427       GETMD=MDSCALEFAC*3.*GETTEMP(X1,Y1,Z1,T1)
0428       GETMD=MAX(GETMD,MDFACTOR)
0429       END
0430
0431
0432
0433       DOUBLE PRECISION FUNCTION GETMS(X2,Y2,Z2,T2)
0434       IMPLICIT NONE
0435       DOUBLE PRECISION X2,Y2,Z2,T2,GETMD
0436       GETMS=GETMD(X2,Y2,Z2,T2)/SQRT(2.)
0437       END
0438
0439
0440
0441       DOUBLE PRECISION FUNCTION GETNEFF(X3,Y3,Z3,T3)
0442       IMPLICIT NONE
0443       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0444       INTEGER NF
0445       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0446       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0447      &N0,SIGMANN,A,WOODSSAXON
0448       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0449      &SIGMANN
0450       INTEGER A
0451       LOGICAL WOODSSAXON
0452 C--   local variables
0453       DOUBLE PRECISION X3,Y3,Z3,T3,PI,GETTEMP,tau,cosheta
0454       DATA PI/3.141592653589793d0/
0455         tau = sqrt(t3**2-z3**2)
0456         cosheta = t3/tau
0457       GETNEFF=(2.*6.*NF*D3*2./3. + 16.*ZETA3*3./2.)
0458      &     *GETTEMP(X3,Y3,Z3,T3)**3/PI**2
0459         getneff = getneff/cosheta
0460       END
0461
0462
0463
0464       DOUBLE PRECISION FUNCTION GETTEMP(X4,Y4,Z4,T4)
0465       IMPLICIT NONE
0466 C--medium parameters
0467       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0468       INTEGER NF
0469       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0470       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0471      &N0,SIGMANN,A,WOODSSAXON
0472       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0473      &SIGMANN
0474       INTEGER A
0475       LOGICAL WOODSSAXON
0476 C--max rapidity
0477         common/rapmax2/etamax2
0478         double precision etamax2
0479 C--local variables
0480       DOUBLE PRECISION X4,Y4,Z4,T4,TAU,NPART,EPS0,EPSIN,TEMPIN,PI,
0481      &NTHICK,ys
0482       DATA PI/3.141592653589793d0/
0483
0484       GETTEMP=0.D0
0485
0486       IF(ABS(Z4).GT.T4)RETURN
0487
0488       TAU=SQRT(T4**2-Z4**2)
0489 C--check for overlap region
0490       IF((NTHICK(X4-BREAL/2.,Y4).EQ.0.d0).OR.
0491      &NTHICK(X4+BREAL/2.,Y4).EQ.0.d0) RETURN
0492
0493         ys = 0.5*log((t4+z4)/(t4-z4))
0494         if (abs(ys).gt.etamax2) return
0495 C--determine initial temperature at transverse position
0496       IF(WOODSSAXON)THEN
0497          EPS0=(16.*8.+7.*2.*6.*NF)*PI**2*TI**4/240.
0498          EPSIN=EPS0*NPART(X4-BREAL/2.,Y4,X4+BREAL/2.,Y4)
0499      &        *PI*RAU**2/(2.*A)
0500          TEMPIN=(EPSIN*240./(PI**2*(16.*8.+7.*2.*6.*NF)))**0.25
0501       ELSE
0502          TEMPIN=TI
0503       ENDIF
0504 C--calculate temperature if before initial time
0505       IF(TAU.LE.TAUI)THEN
0506          GETTEMP=TEMPIN*TAU/TAUI
0507       ELSE
0508 C--evolve temperature
0509        GETTEMP=TEMPIN*(TAUI/TAU)**0.3333
0510       ENDIF
0511       IF(GETTEMP.LT.TC) GETTEMP=0.d0
0512       END
0513
0514
0515
0516       DOUBLE PRECISION FUNCTION GETTEMPMAX()
0517       IMPLICIT NONE
0518 C--medium parameters
0519       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0520       INTEGER NF
0521       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0522       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0523      &N0,SIGMANN,A,WOODSSAXON
0524       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0525      &SIGMANN
0526       INTEGER A
0527       LOGICAL WOODSSAXON
0528 C--function call
0529       DOUBLE PRECISION GETTEMP
0530       GETTEMPMAX=GETTEMP(0.D0,0.D0,0.D0,TAUI)
0531       END
0532
0533
0534
0535       DOUBLE PRECISION FUNCTION GETMDMAX()
0536       IMPLICIT NONE
0537 C--factor to vary Debye mass
0538         COMMON/MDFAC/MDFACTOR,MDSCALEFAC
0539         DOUBLE PRECISION MDFACTOR,MDSCALEFAC
0540       DOUBLE PRECISION GETTEMPMAX
0541       GETMDMAX=MDSCALEFAC*3.*GETTEMPMAX()
0542       GETMDMAX=MAX(GETMDMAX,MDFACTOR)
0543       END
0544
0545
0546
0547       DOUBLE PRECISION FUNCTION GETMDMIN()
0548       IMPLICIT NONE
0549 C--medium parameters
0550       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0551       INTEGER NF
0552       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0553       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0554      &N0,SIGMANN,A,WOODSSAXON
0555       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0556      &SIGMANN
0557       INTEGER A
0558       LOGICAL WOODSSAXON
0559 C--factor to vary Debye mass
0560         COMMON/MDFAC/MDFACTOR,MDSCALEFAC
0561         DOUBLE PRECISION MDFACTOR,MDSCALEFAC
0562       DOUBLE PRECISION GETTEMPMAX
0563         GETMDMIN=MDSCALEFAC*3.*TC
0564       GETMDMIN=MAX(GETMDMIN,MDFACTOR)
0565       END
0566
0567
0568
0569       DOUBLE PRECISION FUNCTION GETMSMAX()
0570       IMPLICIT NONE
0571       DOUBLE PRECISION GETMDMAX,SQRT
0572       GETMSMAX=GETMDMAX()/SQRT(2.D0)
0573       END
0574
0575
0576
0577         DOUBLE PRECISION FUNCTION GETNATMDMIN()
0578         IMPLICIT NONE
0579 C--medium parameters
0580       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0581       INTEGER NF
0582       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0583       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0584      &N0,SIGMANN,A,WOODSSAXON
0585       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0586      &SIGMANN
0587       INTEGER A
0588       LOGICAL WOODSSAXON
0589 C--max rapidity
0590         common/rapmax2/etamax2
0591         double precision etamax2
0592 C--factor to vary Debye mass
0593         COMMON/MDFAC/MDFACTOR,MDSCALEFAC
0594         DOUBLE PRECISION MDFACTOR,MDSCALEFAC,PI
0595       DATA PI/3.141592653589793d0/
0596 C--local variables
0597         DOUBLE PRECISION T,GETMDMIN
0598         T=GETMDMIN()/(MDSCALEFAC*3.)
0599       GETNATMDMIN=(2.*6.*NF*D3*2./3. + 16.*ZETA3*3./2.)
0600      &     *T**3/PI**2
0601         END
0602
0603
0604
0605         DOUBLE PRECISION FUNCTION GETLTIMEMAX()
0606         IMPLICIT NONE
0607 C--medium parameters
0608       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0609       INTEGER NF
0610       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0611       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0612      &N0,SIGMANN,A,WOODSSAXON
0613       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0614      &SIGMANN
0615       INTEGER A
0616       LOGICAL WOODSSAXON
0617 C--max rapidity
0618         common/rapmax2/etamax2
0619         double precision etamax2
0620 C--function call
0621       DOUBLE PRECISION GETTEMPMAX
0622         GETLTIMEMAX=TAUI*(GETTEMPMAX()/TC)**3*cosh(etamax2)
0623         END
0624
0625
0626
0627       DOUBLE PRECISION FUNCTION GETNEFFMAX()
0628       IMPLICIT NONE
0629       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0630       INTEGER NF
0631       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0632       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0633      &N0,SIGMANN,A,WOODSSAXON
0634       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0635      &SIGMANN
0636       INTEGER A
0637       LOGICAL WOODSSAXON
0638 C--max rapidity
0639         common/rapmax2/etamax2
0640         double precision etamax2
0641 C--   local variables
0642       DOUBLE PRECISION PI,GETTEMPMAX
0643       DATA PI/3.141592653589793d0/
0644       GETNEFFMAX=(2.*6.*NF*D3*2./3. + 16.*ZETA3*3./2.)
0645      &     *GETTEMPMAX()**3/PI**2
0646       END
0647
0648
0649
0650       DOUBLE PRECISION FUNCTION NPART(XX1,YY1,XX2,YY2)
0651       IMPLICIT NONE
0652       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0653      &N0,SIGMANN,A,WOODSSAXON
0654       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0655      &SIGMANN
0656       INTEGER A
0657       LOGICAL WOODSSAXON
0658 C--local variables
0659       DOUBLE PRECISION XX1,YY1,XX2,YY2,NTHICK
0660
0661       NPART = NTHICK(XX1,YY1)*(1.-EXP(-SIGMANN*NTHICK(XX2,YY2))) +
0662      &        NTHICK(XX2,YY2)*(1.-EXP(-SIGMANN*NTHICK(XX1,YY1)))
0663       END
0664
0665
0666
0667       DOUBLE PRECISION FUNCTION NTHICK(X1,Y1)
0668       IMPLICIT NONE
0669 C--medium parameters
0670       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0671       INTEGER NF
0672       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0673       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0674      &N0,SIGMANN,A,WOODSSAXON
0675       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0676      &SIGMANN
0677       INTEGER A
0678       LOGICAL WOODSSAXON
0679 C--identifier of log file
0680         common/logfile/logfid
0681         integer logfid
0682 C--nuclear thickness function
0683       COMMON /THICKFNC/ RMAX,TA(100,2)
0684       DOUBLE PRECISION RMAX,TA
0685       INTEGER LINE,LMIN,LMAX,I
0686       DOUBLE PRECISION X1,Y1,XA(4),YA(4),Y,DY,R,C,B,DELTA
0687
0688       R=SQRT(X1**2+Y1**2)
0689       IF(R.GT.TA(100,1))THEN
0690          NTHICK=0.
0691       ELSE
0692          LINE=INT(R*99.d0/TA(100,1)+1)
0693          LMIN=MAX(LINE,1)
0694          LMIN=MIN(LMIN,99)
0695          IF((R.LT.TA(LMIN,1)).OR.(R.GT.TA(LMIN+1,1)))
0696      &  write(logfid,*)LINE,LMIN,R,TA(LMIN,1),TA(LMIN+1,1)
0697          XA(1)=TA(LMIN,1)
0698          XA(2)=TA(LMIN+1,1)
0699          YA(1)=TA(LMIN,2)
0700          YA(2)=TA(LMIN+1,2)
0701          C=(YA(2)-YA(1))/(XA(2)-XA(1))
0702          B=YA(1)-C*XA(1)
0703          NTHICK=C*R+B
0704       ENDIF
0705       END
0706
0707
0708
0709       SUBROUTINE CALCTA()
0710       IMPLICIT NONE
0711 C--medium parameters
0712       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0713       INTEGER NF
0714       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0715       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0716      &N0,SIGMANN,A,WOODSSAXON
0717       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0718      &SIGMANN
0719       INTEGER A
0720       LOGICAL WOODSSAXON
0721 C--   nuclear thickness function
0722       COMMON /THICKFNC/ RMAX,TA(100,2)
0723       DOUBLE PRECISION RMAX,TA
0724 C--variables for integration
0725       COMMON/INTEG/B,R
0726       DOUBLE PRECISION B,R
0727 C--local variables
0728       INTEGER NSTEPS,I
0729       DOUBLE PRECISION EPS,HFIRST,Y
0730
0731       NSTEPS=100
0732       EPS=1.E-4
0733       HFIRST=0.1D0
0734
0735       R=1.12*A**(0.33333)-0.86*A**(-0.33333)
0736       RMAX=2.*R
0737
0738       DO 10 I=1,NSTEPS
0739 C--set transverse position
0740        B=(I-1)*2.D0*R/NSTEPS
0741        Y=0.D0
0742 C--integrate along longitudinal line
0743        CALL ODEINT(Y,-2*R,2*R,EPS,HFIRST,0.d0,101)
0744        TA(I,1)=B
0745        TA(I,2)=Y
0746  10   CONTINUE
0747       END
0748
0749
0750
0751       SUBROUTINE CALCXSECTION()
0752       IMPLICIT NONE
0753 C--medium parameters
0754       COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
0755       INTEGER NF
0756       DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
0757       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0758      &N0,SIGMANN,A,WOODSSAXON
0759       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0760      &SIGMANN
0761       INTEGER A
0762       LOGICAL WOODSSAXON
0763 C--   geometrical cross section
0764       COMMON /CROSSSEC/ IMPMAX,CROSS(200,3)
0765       DOUBLE PRECISION IMPMAX,CROSS
0766 C--local variables
0767       INTEGER IX,IY,IB
0768       DOUBLE PRECISION B,P,PROD,X,Y,NTHICK,NPART,pprev
0769
0770       pprev=0.
0771       DO 30 IB=1,200
0772        B=0.1d0*IB
0773        PROD=1.d0
0774        DO 10 IX=1,100
0775         DO 20 IY=1,100
0776          X=-20.d0+IX*0.4d0
0777          Y=-20.d0+IY*0.4d0
0778          PROD=PROD*
0779      &EXP(-NTHICK(X+B/2.D0,Y)*SIGMANN)**(0.16d0*NTHICK(X-B/2.D0,Y))
0780  20     CONTINUE
0781  10    CONTINUE
0782        P=(1.D0-PROD)*8.8D0/14.D0*B
0783        CROSS(IB,1)=B
0784        CROSS(IB,2)=P
0785        if (ib.eq.1) then
0786         cross(ib,3)=0.
0787        else
0788         cross(ib,3)=cross(ib-1,3)+(p+pprev)/2.*0.1
0789        endif
0790        pprev=p
0791  30   CONTINUE
0792       IMPMAX=19.95
0793       END
0794
0795
0796
0797       DOUBLE PRECISION FUNCTION MEDDERIV(XVAL,W)
0798       IMPLICIT NONE
0799       DOUBLE PRECISION XVAL
0800       INTEGER W
0801 C--medium parameters
0802       COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
0803      &N0,SIGMANN,A,WOODSSAXON
0804       DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
0805      &SIGMANN
0806       INTEGER A
0807       LOGICAL WOODSSAXON
0808 C--variables for integration
0809       COMMON/INTEG/B,R
0810       DOUBLE PRECISION B,R
0811
0812       IF (W.EQ.1) THEN
0813 C--XVAL corresponds to z-coordinate
0814        MEDDERIV=N0/(1+EXP((SQRT(B**2+XVAL**2)-R)/D))
0815       ELSE
0816        MEDDERIV=0.D0
0817       ENDIF
0818       END