Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
medium-simple.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file medium-simple.f
1 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 C++ Copyright (C) 2017 Korinna C. Zapp [Korinna.Zapp@cern.ch] ++
3 C++ ++
4 C++ This file is part of JEWEL 2.2.0 ++
5 C++ ++
6 C++ The JEWEL homepage is jewel.hepforge.org ++
7 C++ ++
8 C++ The medium model was partly implemented by Jochen Klein. ++
9 C++ Raghav Kunnawalkam Elayavalli helped with the implementation ++
10 C++ of the V+jet processes. ++
11 C++ ++
12 C++ Please follow the MCnet GUIDELINES and cite Eur.Phys.J. C74 ++
13 C++ (2014) no.2, 2762 [arXiv:1311.0048] for the code and ++
14 C++ JHEP 1303 (2013) 080 [arXiv:1212.1599] and ++
15 C++ optionally EPJC 60 (2009) 617 [arXiv:0804.3568] for the ++
16 C++ physics. The reference for V+jet processes is EPJC 76 (2016) ++
17 C++ no.12 695 [arXiv:1608.03099] and for recoil effects it is ++
18 C++ arXiv:1707.01539.
19 C++ ++
20 C++ JEWEL relies heavily on PYTHIA 6 for the event generation. The ++
21 C++ modified version of PYTHIA 6.4.25 that is distributed with ++
22 C++ JEWEL is, however, not an official PYTHIA release and must not ++
23 C++ be used for anything else. Please refer to results as ++
24 C++ "JEWEL+PYTHIA". ++
25 C++ ++
26 C++ JEWEL also uses code provided by S. Zhang and J. M. Jing ++
27 C++ (Computation of Special Functions, John Wiley & Sons, New York, ++
28 C++ 1996 and http://jin.ece.illinois.edu) for computing the ++
29 C++ exponential integral Ei(x). ++
30 C++ ++
31 C++ ++
32 C++ JEWEL is free software; you can redistribute it and/or ++
33 C++ modify it under the terms of the GNU General Public License ++
34 C++ as published by the Free Software Foundation; either version 2 ++
35 C++ of the License, or (at your option) any later version. ++
36 C++ ++
37 C++ JEWEL is distributed in the hope that it will be useful, ++
38 C++ but WITHOUT ANY WARRANTY; without even the implied warranty of ++
39 C++ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ++
40 C++ GNU General Public License for more details. ++
41 C++ ++
42 C++ You should have received a copy of the GNU General Public ++
43 C++ License along with this program; if not, write to the Free ++
44 C++ Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, ++
45 C++ MA 02110-1301 USA ++
46 C++ ++
47 C++ Linking JEWEL statically or dynamically with other modules is ++
48 C++ making a combined work based on JEWEL. Thus, the terms and ++
49 C++ conditions of the GNU General Public License cover the whole ++
50 C++ combination. ++
51 C++ ++
52 C++ In addition, as a special exception, I give you permission to ++
53 C++ combine JEWEL with the code for the computation of special ++
54 C++ functions provided by S. Zhang and J. M. Jing. You may copy and ++
55 C++ distribute such a system following the terms of the GNU GPL for ++
56 C++ JEWEL and the licenses of the other code concerned, provided ++
57 C++ that you include the source code of that other code when and as ++
58 C++ the GNU GPL requires distribution of source code. ++
59 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
60 
61  SUBROUTINE medinit(FILE,id,etam,mass)
62  IMPLICIT NONE
63 C--medium parameters
64  common/medparam/centrmin,centrmax,breal,centr,rau,nf
65  INTEGER nf
66  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
67  common/medparamint/taui,ti,tc,d3,zeta3,d,
68  &n0,sigmann,a,woodssaxon
69  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
70  &sigmann
71  INTEGER a
72  LOGICAL woodssaxon
73 C--max rapidity
74  common/rapmax2/etamax2
75  double precision etamax2
76 C--longitudinal boost of momentum distribution
77  common/boostmed/boost
78  logical boost
79 C--factor to vary Debye mass
80  common/mdfac/mdfactor,mdscalefac
81  DOUBLE PRECISION mdfactor,mdscalefac
82 C--nuclear thickness function
83  COMMON /thickfnc/ rmax,ta(100,2)
84  DOUBLE PRECISION rmax,ta
85 C--geometrical cross section
86  COMMON /crosssec/ impmax,cross(200,3)
87  DOUBLE PRECISION impmax,cross
88 C--identifier of log file
89  common/logfile/logfid
90  integer logfid
91 
92  DATA rau/10./
93  DATA d3/0.9d0/
94  DATA zeta3/1.2d0/
95 C--local variables
96  INTEGER i,lun,pos,ios,id,mass
97  double precision etam
98  CHARACTER*100 buffer,label,tempbuf
99  CHARACTER*80 file
100  character firstchar
101  logical fileexist
102 
103  etamax2 = etam
104  logfid = id
105 
106  ios=0
107  lun=77
108 
109 C--default settings
110  taui=0.6d0
111  ti=0.36d0
112  tc=0.17d0
113  woodssaxon=.true.
114  centrmin=0.d0
115  centrmax=10.d0
116  nf=3
117  a=mass
118  n0=0.17d0
119  d=0.54d0
120  sigmann=6.2
121  mdfactor=0.45d0
122  mdscalefac=0.9d0
123  boost = .true.
124 
125 C--read settings from file
126  write(logfid,*)
127  inquire(file=file,exist=fileexist)
128  if(fileexist)then
129  write(logfid,*)'Reading medium parameters from ',file
130  OPEN(unit=lun,file=file,status='old',err=10)
131  do 20 i=1,1000
132  READ(lun, '(A)', iostat=ios) buffer
133  if (ios.ne.0) goto 30
134  firstchar = buffer(1:1)
135  if (firstchar.eq.'#') goto 20
136  pos=scan(buffer,' ')
137  label=buffer(1:pos)
138  buffer=buffer(pos+1:)
139  IF (label=="TAUI")THEN
140  READ(buffer,*,iostat=ios) taui
141  ELSE IF (label=="TI") THEN
142  READ(buffer,*,iostat=ios) ti
143  ELSE IF (label=="TC") THEN
144  READ(buffer,*,iostat=ios) tc
145  ELSE IF (label=="WOODSSAXON") THEN
146  READ(buffer,*,iostat=ios) woodssaxon
147  ELSE IF (label=="CENTRMIN") THEN
148  READ(buffer,*,iostat=ios) centrmin
149  ELSE IF (label=="CENTRMAX") THEN
150  READ(buffer,*,iostat=ios) centrmax
151  ELSE IF (label=="NF") THEN
152  READ(buffer,*,iostat=ios) nf
153  ELSE IF (label=="N0") THEN
154  READ(buffer,*,iostat=ios) n0
155  ELSE IF (label=="D") THEN
156  READ(buffer,*,iostat=ios) d
157  ELSE IF (label=="SIGMANN") THEN
158  READ(buffer,*,iostat=ios) sigmann
159  ELSE IF (label=="MDFACTOR") THEN
160  READ(buffer,*,iostat=ios) mdfactor
161  ELSE IF (label=="MDSCALEFAC") THEN
162  READ(buffer,*,iostat=ios) mdscalefac
163  else
164  write(logfid,*)'unknown label ',label
165  endif
166  20 continue
167 
168  30 close(lun,status='keep')
169  write(logfid,*)'...done'
170  goto 40
171 
172  10 write(logfid,*)'Could not open medium parameter file, '//
173  & 'will run with default settings.'
174 
175  else
176  write(logfid,*)'No medium parameter file found, '//
177  & 'will run with default settings.'
178  endif
179 
180  40 write(logfid,*)'using parameters:'
181  write(logfid,*)'TAUI =',taui
182  write(logfid,*)'TI =',ti
183  write(logfid,*)'TC =',tc
184  write(logfid,*)'WOODSSAXON =',woodssaxon
185  write(logfid,*)'CENTRMIN =',centrmin
186  write(logfid,*)'CENTRMAX =',centrmax
187  write(logfid,*)'NF =',nf
188  write(logfid,*)'A =',a
189  write(logfid,*)'N0 =',n0
190  write(logfid,*)'D =',d
191  write(logfid,*)'SIGMANN =',sigmann
192  write(logfid,*)'MDFACTOR =',mdfactor
193  write(logfid,*)'MDSCALEFAC =',mdscalefac
194  write(logfid,*)
195  write(logfid,*)
196  write(logfid,*)
197 
198 C--calculate T_A(x,y)
199  CALL calcta
200 C--calculate geometrical cross section
201  CALL calcxsection
202 
203  END
204 
205 
206 
207  SUBROUTINE mednextevt
208  IMPLICIT NONE
209 C--medium parameters
210  common/medparam/centrmin,centrmax,breal,centr,rau,nf
211  INTEGER nf
212  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
213  common/medparamint/taui,ti,tc,d3,zeta3,d,
214  &n0,sigmann,a,woodssaxon
215  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
216  &sigmann
217  INTEGER a
218  LOGICAL woodssaxon
219 C--geometrical cross section
220  COMMON /crosssec/ impmax,cross(200,3)
221  DOUBLE PRECISION impmax,cross
222 C--local variables
223  integer i,j
224  DOUBLE PRECISION pyr,r,b1,b2,gettemp
225 
226 C--pick an impact parameter
227  r=(pyr(0)*(centrmax-centrmin)+centrmin)/100.
228  i=0
229  do 130 j=1,200
230  if ((r-cross(j,3)/cross(200,3)).ge.0.) then
231  i=i+1
232  else
233  goto 132
234  endif
235  130 continue
236  132 continue
237  b1 = (i-1)*0.1d0
238  b2 = i*0.1d0
239  breal = (b2*(cross(i,3)/cross(200,3)-r)
240  & +b1*(r-cross(i+1,3)/cross(200,3)))/
241  & (cross(i,3)/cross(200,3)-cross(i+1,3)/cross(200,3))
242  centr = r;
243  END
244 
245  double precision function getcentrality()
246  implicit none
247  common/medparam/centrmin,centrmax,breal,centr,rau,nf
248  INTEGER nf
249  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
250  getcentrality=centr
251  end
252 
253 
254 
255  SUBROUTINE pickvtx(X,Y)
256  IMPLICIT NONE
257  DOUBLE PRECISION x,y
258 C--medium parameters
259  common/medparam/centrmin,centrmax,breal,centr,rau,nf
260  INTEGER nf
261  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
262 C--local variables
263  DOUBLE PRECISION x1,x2,y1,y2,z,xval,yval,zval,nthick,pyr
264 
265  x1=breal/2.-rau
266  x2=rau-breal/2.
267  y1=-sqrt(4*rau**2-breal**2)/2.
268  y2=sqrt(4*rau**2-breal**2)/2.
269  131 xval=pyr(0)*(x2-x1)+x1
270  yval=pyr(0)*(y2-y1)+y1
271  IF((nthick(xval-breal/2.,yval).EQ.0.d0).OR.
272  & nthick(xval+breal/2.,yval).EQ.0.d0) goto 131
273  zval=pyr(0)*nthick(-breal/2.,0d0)*nthick(breal/2.,0d0)
274  z=nthick(xval-breal/2.,yval)*nthick(xval+breal/2.,yval)
275  IF(zval.GT.z) goto 131
276  x=xval
277  y=yval
278  END
279 
280  SUBROUTINE setb(BVAL)
281  IMPLICIT NONE
282 C--medium parameters
283  common/medparam/centrmin,centrmax,breal,centr,rau,nf
284  INTEGER nf
285  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
286  DOUBLE PRECISION bval
287  breal=bval
288  END
289 
290 
291 
292  SUBROUTINE getscatterer(X,Y,Z,T,TYPE,PX,PY,PZ,E,MS)
293  IMPLICIT NONE
294 C--medium parameters
295  common/medparam/centrmin,centrmax,breal,centr,rau,nf
296  INTEGER nf
297  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
298 C--internal medium parameters
299  common/medparamint/taui,ti,tc,d3,zeta3,d,
300  &n0,sigmann,a,woodssaxon
301  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
302  &sigmann
303  INTEGER a
304  LOGICAL woodssaxon
305 C--longitudinal boost of momentum distribution
306  common/boostmed/boost
307  logical boost
308 C--function calls
309  DOUBLE PRECISION gettemp,getmd,getmom,getms
310 C--identifier of log file
311  common/logfile/logfid
312  integer logfid
313 C--local variables
314  DOUBLE PRECISION x,y,z,t,ms,px,py,pz,e,md,temp
315  INTEGER type
316  DOUBLE PRECISION r,pyr,pmax,wt,tau,theta,phi,pi,p,ys,pz2,e2
317  DATA pi/3.141592653589793d0/
318 
319  r=pyr(0)
320  IF(r.LT.(2.*12.*nf*d3/3.)/(2.*12.*nf*d3/3.+3.*16.*zeta3/2.))THEN
321  type=2
322  ELSE
323  type=21
324  ENDIF
325  ms=getms(x,y,z,t)
326  md=getmd(x,y,z,t)
327  temp=gettemp(x,y,z,t)
328  tau=sqrt(t**2-z**2)
329  if (boost) then
330  ys = 0.5*log((t+z)/(t-z))
331  else
332  ys = 0.d0
333  endif
334  pmax = 10.*temp
335 
336  IF(temp.LT.1.d-2)THEN
337  write(logfid,*)'asking for a scattering centre without medium:'
338  write(logfid,*)'at (x,y,z,t)=',x,y,z,t
339  write(logfid,*)'making one up to continue but '//
340  & 'something is wrong!'
341  type=21
342  px=0.d0
343  py=0.d0
344  pz=0.d0
345  ms=getms(0.d0,0.d0,0.d0,0.d0)
346  md=getmd(0.d0,0.d0,0.d0,0.d0)
347  e=sqrt(px**2+py**2+pz**2+ms**2)
348  RETURN
349  ENDIF
350 
351  10 p = pyr(0)**0.3333333*pmax
352  e2 = sqrt(p**2+ms**2)
353  if (type.eq.2) then
354  wt = (exp(ms/temp)-1.)/(exp(e2/temp)-1.)
355  else
356  wt = (exp(ms/temp)+1.)/(exp(e2/temp)+1.)
357  endif
358  if (wt.gt.1.) write(logfid,*)'Error in getscatterer: weight = ',wt
359  if (wt.lt.0.) write(logfid,*)'Error in getscatterer: weight = ',wt
360  if (pyr(0).gt.wt) goto 10
361  phi = pyr(0)*2.*pi
362  theta = -acos(2.*pyr(0)-1.)+pi
363  px = p*sin(theta)*cos(phi)
364  py = p*sin(theta)*sin(phi)
365  pz2 = p*cos(theta)
366  e = cosh(ys)*e2 + sinh(ys)*pz2
367  pz = sinh(ys)*e2 + cosh(ys)*pz2
368  END
369 
370 
371  SUBROUTINE avscatcen(X,Y,Z,T,PX,PY,PZ,E,m)
372  IMPLICIT NONE
373 C--longitudinal boost of momentum distribution
374  common/boostmed/boost
375  logical boost
376 C--max rapidity
377  common/rapmax2/etamax2
378  double precision etamax2
379 C--local variables
380  double precision x,y,z,t,px,py,pz,e,getms,m,ys
381  if (boost) then
382  ys = 0.5*log((t+z)/(t-z))
383  if ((z.eq.0.d0).and.(t.eq.0.d0)) ys =0.d0
384  if (ys.gt.etamax2) ys=etamax2
385  if (ys.lt.-etamax2) ys=-etamax2
386  else
387  ys = 0.d0
388  endif
389  m = getms(x,y,z,t)
390  e = m*cosh(ys)
391  px = 0.d0
392  py = 0.d0
393  pz = m*sinh(ys)
394  end
395 
396 
397  SUBROUTINE maxscatcen(PX,PY,PZ,E,m)
398  IMPLICIT NONE
399 C--longitudinal boost of momentum distribution
400  common/boostmed/boost
401  logical boost
402 C--max rapidity
403  common/rapmax2/etamax2
404  double precision etamax2
405 C--local variables
406  double precision px,py,pz,e,getmsmax,m,ys
407  if (boost) then
408  ys = etamax2
409  else
410  ys = 0.d0
411  endif
412  m = getmsmax()
413  e = m*cosh(ys)
414  px = 0.d0
415  py = 0.d0
416  pz = m*sinh(ys)
417  end
418 
419 
420 
421  DOUBLE PRECISION FUNCTION getmd(X1,Y1,Z1,T1)
422  IMPLICIT NONE
423 C--factor to vary Debye mass
424  common/mdfac/mdfactor,mdscalefac
425  DOUBLE PRECISION mdfactor,mdscalefac
426  DOUBLE PRECISION x1,y1,z1,t1,gettemp
427  getmd=mdscalefac*3.*gettemp(x1,y1,z1,t1)
428  getmd=max(getmd,mdfactor)
429  END
430 
431 
432 
433  DOUBLE PRECISION FUNCTION getms(X2,Y2,Z2,T2)
434  IMPLICIT NONE
435  DOUBLE PRECISION x2,y2,z2,t2,getmd
436  getms=getmd(x2,y2,z2,t2)/sqrt(2.)
437  END
438 
439 
440 
441  DOUBLE PRECISION FUNCTION getneff(X3,Y3,Z3,T3)
442  IMPLICIT NONE
443  common/medparam/centrmin,centrmax,breal,centr,rau,nf
444  INTEGER nf
445  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
446  common/medparamint/taui,ti,tc,d3,zeta3,d,
447  &n0,sigmann,a,woodssaxon
448  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
449  &sigmann
450  INTEGER a
451  LOGICAL woodssaxon
452 C-- local variables
453  DOUBLE PRECISION x3,y3,z3,t3,pi,gettemp,tau,cosheta
454  DATA pi/3.141592653589793d0/
455  tau = sqrt(t3**2-z3**2)
456  cosheta = t3/tau
457  getneff=(2.*6.*nf*d3*2./3. + 16.*zeta3*3./2.)
458  & *gettemp(x3,y3,z3,t3)**3/pi**2
459  getneff = getneff/cosheta
460  END
461 
462 
463 
464  DOUBLE PRECISION FUNCTION gettemp(X4,Y4,Z4,T4)
465  IMPLICIT NONE
466 C--medium parameters
467  common/medparam/centrmin,centrmax,breal,centr,rau,nf
468  INTEGER nf
469  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
470  common/medparamint/taui,ti,tc,d3,zeta3,d,
471  &n0,sigmann,a,woodssaxon
472  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
473  &sigmann
474  INTEGER a
475  LOGICAL woodssaxon
476 C--max rapidity
477  common/rapmax2/etamax2
478  double precision etamax2
479 C--local variables
480  DOUBLE PRECISION x4,y4,z4,t4,tau,npart,eps0,epsin,tempin,pi,
481  &nthick,ys
482  DATA pi/3.141592653589793d0/
483 
484  gettemp=0.d0
485 
486  IF(abs(z4).GT.t4)RETURN
487 
488  tau=sqrt(t4**2-z4**2)
489 C--check for overlap region
490  IF((nthick(x4-breal/2.,y4).EQ.0.d0).OR.
491  &nthick(x4+breal/2.,y4).EQ.0.d0) RETURN
492 
493  ys = 0.5*log((t4+z4)/(t4-z4))
494  if (abs(ys).gt.etamax2) return
495 C--determine initial temperature at transverse position
496  IF(woodssaxon)THEN
497  eps0=(16.*8.+7.*2.*6.*nf)*pi**2*ti**4/240.
498  epsin=eps0*npart(x4-breal/2.,y4,x4+breal/2.,y4)
499  & *pi*rau**2/(2.*a)
500  tempin=(epsin*240./(pi**2*(16.*8.+7.*2.*6.*nf)))**0.25
501  ELSE
502  tempin=ti
503  ENDIF
504 C--calculate temperature if before initial time
505  IF(tau.LE.taui)THEN
506  gettemp=tempin*tau/taui
507  ELSE
508 C--evolve temperature
509  gettemp=tempin*(taui/tau)**0.3333
510  ENDIF
511  IF(gettemp.LT.tc) gettemp=0.d0
512  END
513 
514 
515 
516  DOUBLE PRECISION FUNCTION gettempmax()
517  IMPLICIT NONE
518 C--medium parameters
519  common/medparam/centrmin,centrmax,breal,centr,rau,nf
520  INTEGER nf
521  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
522  common/medparamint/taui,ti,tc,d3,zeta3,d,
523  &n0,sigmann,a,woodssaxon
524  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
525  &sigmann
526  INTEGER a
527  LOGICAL woodssaxon
528 C--function call
529  DOUBLE PRECISION gettemp
530  gettempmax=gettemp(0.d0,0.d0,0.d0,taui)
531  END
532 
533 
534 
535  DOUBLE PRECISION FUNCTION getmdmax()
536  IMPLICIT NONE
537 C--factor to vary Debye mass
538  common/mdfac/mdfactor,mdscalefac
539  DOUBLE PRECISION mdfactor,mdscalefac
540  DOUBLE PRECISION gettempmax
541  getmdmax=mdscalefac*3.*gettempmax()
542  getmdmax=max(getmdmax,mdfactor)
543  END
544 
545 
546 
547  DOUBLE PRECISION FUNCTION getmdmin()
548  IMPLICIT NONE
549 C--medium parameters
550  common/medparam/centrmin,centrmax,breal,centr,rau,nf
551  INTEGER nf
552  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
553  common/medparamint/taui,ti,tc,d3,zeta3,d,
554  &n0,sigmann,a,woodssaxon
555  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
556  &sigmann
557  INTEGER a
558  LOGICAL woodssaxon
559 C--factor to vary Debye mass
560  common/mdfac/mdfactor,mdscalefac
561  DOUBLE PRECISION mdfactor,mdscalefac
562  DOUBLE PRECISION gettempmax
563  getmdmin=mdscalefac*3.*tc
564  getmdmin=max(getmdmin,mdfactor)
565  END
566 
567 
568 
569  DOUBLE PRECISION FUNCTION getmsmax()
570  IMPLICIT NONE
571  DOUBLE PRECISION getmdmax,sqrt
572  getmsmax=getmdmax()/sqrt(2.d0)
573  END
574 
575 
576 
577  DOUBLE PRECISION FUNCTION getnatmdmin()
578  IMPLICIT NONE
579 C--medium parameters
580  common/medparam/centrmin,centrmax,breal,centr,rau,nf
581  INTEGER nf
582  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
583  common/medparamint/taui,ti,tc,d3,zeta3,d,
584  &n0,sigmann,a,woodssaxon
585  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
586  &sigmann
587  INTEGER a
588  LOGICAL woodssaxon
589 C--max rapidity
590  common/rapmax2/etamax2
591  double precision etamax2
592 C--factor to vary Debye mass
593  common/mdfac/mdfactor,mdscalefac
594  DOUBLE PRECISION mdfactor,mdscalefac,pi
595  DATA pi/3.141592653589793d0/
596 C--local variables
597  DOUBLE PRECISION t,getmdmin
598  t=getmdmin()/(mdscalefac*3.)
599  getnatmdmin=(2.*6.*nf*d3*2./3. + 16.*zeta3*3./2.)
600  & *t**3/pi**2
601  END
602 
603 
604 
605  DOUBLE PRECISION FUNCTION getltimemax()
606  IMPLICIT NONE
607 C--medium parameters
608  common/medparam/centrmin,centrmax,breal,centr,rau,nf
609  INTEGER nf
610  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
611  common/medparamint/taui,ti,tc,d3,zeta3,d,
612  &n0,sigmann,a,woodssaxon
613  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
614  &sigmann
615  INTEGER a
616  LOGICAL woodssaxon
617 C--max rapidity
618  common/rapmax2/etamax2
619  double precision etamax2
620 C--function call
621  DOUBLE PRECISION gettempmax
622  getltimemax=taui*(gettempmax()/tc)**3*cosh(etamax2)
623  END
624 
625 
626 
627  DOUBLE PRECISION FUNCTION getneffmax()
628  IMPLICIT NONE
629  common/medparam/centrmin,centrmax,breal,centr,rau,nf
630  INTEGER nf
631  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
632  common/medparamint/taui,ti,tc,d3,zeta3,d,
633  &n0,sigmann,a,woodssaxon
634  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
635  &sigmann
636  INTEGER a
637  LOGICAL woodssaxon
638 C--max rapidity
639  common/rapmax2/etamax2
640  double precision etamax2
641 C-- local variables
642  DOUBLE PRECISION pi,gettempmax
643  DATA pi/3.141592653589793d0/
644  getneffmax=(2.*6.*nf*d3*2./3. + 16.*zeta3*3./2.)
645  & *gettempmax()**3/pi**2
646  END
647 
648 
649 
650  DOUBLE PRECISION FUNCTION npart(XX1,YY1,XX2,YY2)
651  IMPLICIT NONE
652  common/medparamint/taui,ti,tc,d3,zeta3,d,
653  &n0,sigmann,a,woodssaxon
654  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
655  &sigmann
656  INTEGER a
657  LOGICAL woodssaxon
658 C--local variables
659  DOUBLE PRECISION xx1,yy1,xx2,yy2,nthick
660 
661  npart = nthick(xx1,yy1)*(1.-exp(-sigmann*nthick(xx2,yy2))) +
662  & nthick(xx2,yy2)*(1.-exp(-sigmann*nthick(xx1,yy1)))
663  END
664 
665 
666 
667  DOUBLE PRECISION FUNCTION nthick(X1,Y1)
668  IMPLICIT NONE
669 C--medium parameters
670  common/medparam/centrmin,centrmax,breal,centr,rau,nf
671  INTEGER nf
672  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
673  common/medparamint/taui,ti,tc,d3,zeta3,d,
674  &n0,sigmann,a,woodssaxon
675  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
676  &sigmann
677  INTEGER a
678  LOGICAL woodssaxon
679 C--identifier of log file
680  common/logfile/logfid
681  integer logfid
682 C--nuclear thickness function
683  COMMON /thickfnc/ rmax,ta(100,2)
684  DOUBLE PRECISION rmax,ta
685  INTEGER line,lmin,lmax,i
686  DOUBLE PRECISION x1,y1,xa(4),ya(4),y,dy,r,c,b,delta
687 
688  r=sqrt(x1**2+y1**2)
689  IF(r.GT.ta(100,1))THEN
690  nthick=0.
691  ELSE
692  line=int(r*99.d0/ta(100,1)+1)
693  lmin=max(line,1)
694  lmin=min(lmin,99)
695  IF((r.LT.ta(lmin,1)).OR.(r.GT.ta(lmin+1,1)))
696  & write(logfid,*)line,lmin,r,ta(lmin,1),ta(lmin+1,1)
697  xa(1)=ta(lmin,1)
698  xa(2)=ta(lmin+1,1)
699  ya(1)=ta(lmin,2)
700  ya(2)=ta(lmin+1,2)
701  c=(ya(2)-ya(1))/(xa(2)-xa(1))
702  b=ya(1)-c*xa(1)
703  nthick=c*r+b
704  ENDIF
705  END
706 
707 
708 
709  SUBROUTINE calcta()
710  IMPLICIT NONE
711 C--medium parameters
712  common/medparam/centrmin,centrmax,breal,centr,rau,nf
713  INTEGER nf
714  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
715  common/medparamint/taui,ti,tc,d3,zeta3,d,
716  &n0,sigmann,a,woodssaxon
717  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
718  &sigmann
719  INTEGER a
720  LOGICAL woodssaxon
721 C-- nuclear thickness function
722  COMMON /thickfnc/ rmax,ta(100,2)
723  DOUBLE PRECISION rmax,ta
724 C--variables for integration
725  common/integ/b,r
726  DOUBLE PRECISION b,r
727 C--local variables
728  INTEGER nsteps,i
729  DOUBLE PRECISION eps,hfirst,y
730 
731  nsteps=100
732  eps=1.e-4
733  hfirst=0.1d0
734 
735  r=1.12*a**(0.33333)-0.86*a**(-0.33333)
736  rmax=2.*r
737 
738  DO 10 i=1,nsteps
739 C--set transverse position
740  b=(i-1)*2.d0*r/nsteps
741  y=0.d0
742 C--integrate along longitudinal line
743  CALL odeint(y,-2*r,2*r,eps,hfirst,0.d0,101)
744  ta(i,1)=b
745  ta(i,2)=y
746  10 CONTINUE
747  END
748 
749 
750 
751  SUBROUTINE calcxsection()
752  IMPLICIT NONE
753 C--medium parameters
754  common/medparam/centrmin,centrmax,breal,centr,rau,nf
755  INTEGER nf
756  DOUBLE PRECISION centrmin,centrmax,breal,centr,rau
757  common/medparamint/taui,ti,tc,d3,zeta3,d,
758  &n0,sigmann,a,woodssaxon
759  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
760  &sigmann
761  INTEGER a
762  LOGICAL woodssaxon
763 C-- geometrical cross section
764  COMMON /crosssec/ impmax,cross(200,3)
765  DOUBLE PRECISION impmax,cross
766 C--local variables
767  INTEGER ix,iy,ib
768  DOUBLE PRECISION b,p,prod,x,y,nthick,npart,pprev
769 
770  pprev=0.
771  DO 30 ib=1,200
772  b=0.1d0*ib
773  prod=1.d0
774  DO 10 ix=1,100
775  DO 20 iy=1,100
776  x=-20.d0+ix*0.4d0
777  y=-20.d0+iy*0.4d0
778  prod=prod*
779  &exp(-nthick(x+b/2.d0,y)*sigmann)**(0.16d0*nthick(x-b/2.d0,y))
780  20 CONTINUE
781  10 CONTINUE
782  p=(1.d0-prod)*8.8d0/14.d0*b
783  cross(ib,1)=b
784  cross(ib,2)=p
785  if (ib.eq.1) then
786  cross(ib,3)=0.
787  else
788  cross(ib,3)=cross(ib-1,3)+(p+pprev)/2.*0.1
789  endif
790  pprev=p
791  30 CONTINUE
792  impmax=19.95
793  END
794 
795 
796 
797  DOUBLE PRECISION FUNCTION medderiv(XVAL,W)
798  IMPLICIT NONE
799  DOUBLE PRECISION xval
800  INTEGER w
801 C--medium parameters
802  common/medparamint/taui,ti,tc,d3,zeta3,d,
803  &n0,sigmann,a,woodssaxon
804  DOUBLE PRECISION taui,ti,tc,alpha,beta,gamma,d3,zeta3,d,n0,
805  &sigmann
806  INTEGER a
807  LOGICAL woodssaxon
808 C--variables for integration
809  common/integ/b,r
810  DOUBLE PRECISION b,r
811 
812  IF (w.EQ.1) THEN
813 C--XVAL corresponds to z-coordinate
814  medderiv=n0/(1+exp((sqrt(b**2+xval**2)-r)/d))
815  ELSE
816  medderiv=0.d0
817  ENDIF
818  END