]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-con-161.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes needed by gfortran, it doesn't perform lazy evaluation (Piotr)
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-con-161.f
CommitLineData
9ef1c2d9 1c-----------------------------------------------------------------------
2 subroutine conaa(iret)
3c-----------------------------------------------------------------------
4c determines interaction configuration
5c-----------------------------------------------------------------------
6
7 include 'epos.inc'
8 include 'epos.incems'
9 include 'epos.incsem'
10
11 common/geom/rmproj,rmtarg,bmax,bkmx
12 common/nucl3/phi,bimp
13 common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
14 * ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
15 common/cncl3/iactpr(mamx),iacttg(mamx)
16 common/cfacmss/facmss
17
18 call utpri('conaa ',ish,ishini,4)
19
20 iret=0
21
22c initialisations
23c ---------------
24
25 vel=tanh(ypjtl-yhaha)+tanh(yhaha)
26
27c determine phi, bimp, coll, iproj, itarg, x/y/zproj, x/y/ztarg
28c ---------------------------------------------------------------
29
30 if(iokoll.eq.1)then
31
32 koll=matarg
33 do k=1,koll
34 do n=1,4
35 coord(n,k)=0.
36 enddo
37 enddo
38 bimp=0
39 phi=0
40 do k=1,koll
41 bij=bkmx*sqrt(rangen())
42 bk(k)=bij
43 iproj(k)=1
44 itarg(k)=k
45 enddo
46 xproj(1)=0
47 yproj(1)=0
48 zproj(1)=0
49 iactpr(1)=1
50 do k=1,koll
51 phi=2.*pi*rangen()
52 xtarg(k)=bij*cos(phi)
53 ytarg(k)=bij*sin(phi)
54 ztarg(k)=0
55 iacttg(k)=1
56 enddo
57 lproj(1)=koll
58 lproj3(1)=0
59 do k=1,koll
60 ltarg(k)=1
61 kproj(1,k)=k
62 ktarg(k,1)=k
63 ltarg3(k)=0
64 if(iscreen.ne.0)then
65 lproj3(1)=lproj3(1)+1
66 ltarg3(k)=1
67 kproj3(1,lproj3(1))=k
68 ktarg3(k,1)=k
69 endif
70 enddo
71
72 elseif(maproj.eq.1.and.matarg.eq.1)then
73
74 b1=bminim
75 b2=amin1(bkmx,bmaxim)
76 if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
77 bimp=sqrt(b1*b1+(b2*b2-b1*b1)*rangen())
78 koll=1
79 do n=1,4
80 coord(n,1)=0.
81 enddo
82 bk(1)=bimp
83 iproj(1)=1
84 itarg(1)=1
85 phi=2.*pi*rangen()
86 xproj(1)=bimp*cos(phi)
87 yproj(1)=bimp*sin(phi)
88 zproj(1)=0
89 xtarg(1)=0
90 ytarg(1)=0
91 ztarg(1)=0
92 lproj(1)=1
93 ltarg(1)=1
94 if(iscreen.ne.0)then
95 lproj3(1)=1
96 ltarg3(1)=1
97 else
98 lproj3(1)=0
99 ltarg3(1)=0
100 endif
101 kproj(1,1)=1
102 ktarg(1,1)=1
103 kproj3(1,1)=1
104 ktarg3(1,1)=1
105 iactpr(1)=1 !int(rangen()+1./chad(iclpro)) !active proj.
106 iacttg(1)=1 !int(rangen()+1./chad(icltar)) !active targ.
107 if(iactpr(1).eq.0.and.iacttg(1).eq.0)iret=-1
108
109 else
110
111 call conxyz('p',mamx,xproj,yproj,zproj,ypjtl-yhaha)
112 call conxyz('t',mamx,xtarg,ytarg,ztarg,yhaha)
113
114 bx=0
115 by=0
116 if(maproj.gt.0)then
117 if(bimevt.lt.0)then
118 b1=bminim
119 b2=amin1(rmproj+rmtarg,bmaxim)
120 if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
121 bimp=sqrt(b1**2+(b2**2-b1**2)*rangen())
122 if(nbarray.gt.0)bimp=barray(mod(nrevt,nbarray)+1)
123 if(jpsi.gt.0)then
124 bimp=b1+(b2-b1)*(float(mod(nrevt,12))+rangen())/12.
125 bimevt=bimp
126 endif
127 phi=phimin+rangen()*(phimax-phimin)
128 else
129 phi=phievt
130 bimp=bimevt
131 endif
132 bx=cos(phi)*bimp
133 by=sin(phi)*bimp
134 endif
135 if(jpsi.lt.0)then !modify b
136 bx=xtarg(1)
137 by=ytarg(1)
138 endif
139 if(maproj.eq.0)goto1000
140c wactv=1./chad(2) !probability to be "active"
141 koll=0
142 do i=1,maproj
143 lproj(i)=0
144 lproj3(i)=0
145 iactpr(i)=1 !int(rangen()+wactv) !active proj. nucleons
146 enddo
147 do j=1,matarg
148 ltarg(j)=0
149 ltarg3(j)=0
150 iacttg(j)=1 !int(rangen()+wactv) !active targ. nucleons
151 enddo
152 do 12 i=1,maproj
153 if(iactpr(i).eq.0)goto 12
154 do 11 j=1,matarg
155 if(iacttg(j).eq.0)goto 11
156 if(jpsi.lt.0.and.ztarg(j).le.ztarg(1))goto11
157 bij=sqrt((xproj(i)+bx-xtarg(j))**2+(yproj(i)+by-ytarg(j))**2)
158 if(ish.ge.7)write(ifch,*)'i_p:',i,' i_t:',j,' b_ij:',bij
159 if(bij.gt.bkmx)goto 11
160
161 koll=koll+1
162 if(koll.gt.kollmx)call utstop('conaa: kollmx too small&')
163 bk(koll)=bij
164 bkx(koll)=xproj(i)+bx-xtarg(j)
165 bky(koll)=yproj(i)+by-ytarg(j)
166 iproj(koll)=i
167 itarg(koll)=j
168 lproj(i)=lproj(i)+1
169 ltarg(j)=ltarg(j)+1
170 kproj(i,lproj(i))=koll
171 ktarg(j,ltarg(j))=koll
172 if(iscreen.ne.0)then
173 if(bij.gt.bkmxndif)goto 11
174 lproj3(i)=lproj3(i)+1
175 ltarg3(j)=ltarg3(j)+1
176 kproj3(i,lproj3(i))=koll
177 ktarg3(j,ltarg3(j))=koll
178 endif
179
18011 continue
18112 continue
182
183 do k=1,koll
184 do n=1,4
185 coord(n,k)=0.
186 enddo
187 enddo
188
189 if(((iactpr(1).eq.0.and.maproj.eq.1)
190 &.or.(iacttg(1).eq.0.and.matarg.eq.1))
191 &.and.iret.eq.0)iret=-1
192
193 endif
194
195 if(ish.ge.3)write(ifch,*)'koll=',koll
196 if(koll.eq.0)goto 1001
197
198
199c determine coord
200c ---------------
201 do kl=1,koll
202 i=iproj(kl)
203 j=itarg(kl)
204 dist=ztarg(j)-zproj(i)
205 coord(1,kl)=(xproj(i)+xtarg(j))*0.5
206 coord(2,kl)=(yproj(i)+ytarg(j))*0.5
207 coord(3,kl)=(zproj(i)+ztarg(j))*0.5
208 coord(4,kl)=dist/vel
209 enddo
210
211
212
213c exit
214c ----
2151000 continue
216 if(ish.ge.5)then
217 write(ifch,*)'ia,x/y/zproj:'
218 do mm=1,maproj
219 write(ifch,*)mm,iactpr(mm),xproj(mm),yproj(mm),zproj(mm)
220 enddo
221 write(ifch,*)'ia,x/y/ztarg:'
222 do mm=1,matarg
223 write(ifch,*)mm,iacttg(mm),xtarg(mm),ytarg(mm),ztarg(mm)
224 enddo
225 write(ifch,*)'iret',iret
226 endif
227 call utprix('conaa ',ish,ishini,4)
228 return
229
2301001 continue !iret=1 causes redo of whole collision
231 iret=1
232 if(ish.ge.3)then
233 write(ifch,*)
234 write(ifch,*)'***** subroutine conaa:'
235 write(ifch,*)'***** no nucleon pair found --> no interaction'
236 write(ifch,*)
237 endif
238 goto 1000
239
240 end
241
242c-----------------------------------------------------------------------
243 subroutine xGeometry(iii)
244c-----------------------------------------------------------------------
245 include 'epos.inc'
246 include 'epos.incems'
247 include 'epos.incsem'
248 common/xgeom/nnn,naa(kollmx),nbb(kollmx)
249 character*5 fmt1,fmt2
250
251 if(iii.eq.0)then
252
253 do k=1,kollmx
254 naa(k)=0
255 enddo
256 nnn=0
257
258
259 elseif(iii.eq.1)then
260
261 ngl=0
262 do k=1,koll
263 r=bk(k)
264 if(r.le.sqrt(sigine/10./pi))ngl=ngl+1
265 enddo
266 if(ngl.ne.0)then
267 nnn=nnn+1
268 naa(ngl)=naa(ngl)+1
269 endif
270
271 elseif(iii.eq.2)then
272
273 if(xpar1.eq.0..and.xpar2.eq.0.)then
274 print*,'---------- xGeometry -----------------------'
275 return
276 endif
277 x1=1-0.01*xpar2
278 x2=1-0.01*xpar1
279 kmx=0
280 nbb(1)=naa(1)
281 do k=2,kollmx
282 if(naa(k).ne.0.)kmx=k
283 nbb(k)=nbb(k-1)+naa(k)
284 enddo
285 k1=0
286 k2=0
287 do k=1,kmx
288 x=nbb(k)/float(nnn)
289 if(x.lt.x1)k1=k
290 if(x.lt.x2)k2=k
291 enddo
292 k1=k1+1
293 k2=k2+1
294 x=0
295 av=0
296 su=0
297 do k=1,kmx
298 xb=x
299 x=nbb(k)/float(nnn)
300 y=naa(k)/float(nnn)
301 dx=x-xb
302 p=0.
303 if(k.eq.k1)then
304 p=(x-x1)/dx
305 p1=p
306 elseif(k.eq.k2)then
307 p=(x2-xb)/dx
308 p2=p
309 elseif(k.gt.k1.and.k.lt.k2)then
310 p=1
311 endif
312 av=av+y*p*k
313 su=su+y*p
314 enddo
315 av=av/su
316 n1=nint(100*p1)
317 n2=nint(100*p2)
318 if(n1.eq.0)then
319 k1=k1+1
320 n1=100
321 endif
322 if(k1.le.9)fmt1='i1,4x'
323 if(k1.gt.9.and.k1.le.99)fmt1='i2,3x'
324 if(k1.gt.99.and.k1.le.999)fmt1='i3,2x'
325 if(k1.gt.999.and.k1.le.9999)fmt1='i4,1x'
326 if(k2.le.9)fmt2='i1,4x'
327 if(k2.gt.9.and.k2.le.99)fmt2='i2,3x'
328 if(k2.gt.99.and.k2.le.999)fmt2='i3,2x'
329 if(k2.gt.999.and.k2.le.9999)fmt2='i4,1x'
330 write(6,'(i4,a,i5,a,'//fmt1//',i6,a,i5,a,'//fmt2//',5x,a,f8.2)')
331 & nint(xpar2),':MIN',n1,'%',k1
332 & ,nint(xpar1),':MAX',n2,'%',k2 ,'av:',av
333 endif
334
335 end
336
337c-----------------------------------------------------------------------
338 function conbmx()
339c-----------------------------------------------------------------------
340 double precision om1intbc,p,eps
341 include 'epos.inc'
342 include 'epos.incsem'
343
344 conbmx=0.
345 b1=0.
346 b2=7.
347 eps=5.0d-3
348 p=1.d0-dexp(-om1intbc(b2))
349 if(p.gt.2.d0*eps)return
350
351 ntry=0
352
35310 ntry=ntry+1
354 b=b1+.5*(b2-b1)
355
356 p=(1.d0-dexp(-om1intbc(b)))
357
358 if(p.gt.eps)then
359 if(p.gt.2.d0*eps)then
360 b1=b
361 else
362 conbmx=b
363 return
364 endif
365 else
366 if(p.lt.eps/5.d0)then
367 b2=b
368 else
369 conbmx=b
370 return
371 endif
372 endif
373
374 if(ntry.le.1000)goto 10
375 write(ifmt,*)'Too much try in conbmx ... bmax=',b
376 conbmx=b
377 return
378
379 end
380
381c-----------------------------------------------------------------------
382 function conbmxndif()
383c-----------------------------------------------------------------------
384 double precision om1intbc,p,eps
385 include 'epos.inc'
386 include 'epos.incsem'
387
388
389 iomegasave=iomega
390 iomega=2
391 conbmxndif=0.
392 b1=0.
393 b2=7.
394 bkmxndif=b2
395 eps=zbcut
396 p=1.d0-dexp(-om1intbc(b2))
397 if(p.gt.2.d0*eps)goto 100
398
399 ntry=0
400
40110 ntry=ntry+1
402 b=b1+.5*(b2-b1)
403
404 p=(1.d0-dexp(-om1intbc(b)))
405
406 if(p.gt.eps)then
407 if(p.gt.2.d0*eps)then
408 b1=b
409 else
410 conbmxndif=b
411 goto 100
412 endif
413 else
414 if(p.lt.eps/5.d0)then
415 b2=b
416 else
417 conbmxndif=b
418 goto 100
419 endif
420 endif
421
422 if(ntry.le.1000)goto 10
423 write(ifmt,*)'Too much try in conbmxndif ... bkmxndif=',b
424 conbmxndif=b
425 100 iomega=iomegasave
426 return
427
428 end
429
430c-----------------------------------------------------------------------
431 function conbmxdif()
432c-----------------------------------------------------------------------
433c find b to have (1-exp(-om))pmax=pdiff
434c-----------------------------------------------------------------------
435 double precision om1intbc,pmax,drootom,pdiff
436 include 'epos.inc'
437 include 'epos.incsem'
438
439 conbmxdif=0.
440 b1=0.
441 bmax=7.
442 iomegasave=iomega
443 iomega=2
444
445 eps=1.e-5
446 pmax=1.d0-dexp(-om1intbc(b1))
447 pdiff=facdif
448 if(pmax.gt.eps)then
449 conbmxdif=drootom(pdiff,pmax,bmax,eps)
450 endif
451 iomega=iomegasave
452
453 return
454
455 end
456
457c-----------------------------------------------------------------------
458 subroutine conre
459c-----------------------------------------------------------------------
460c initializes remnants
461c-----------------------------------------------------------------------
462 include "epos.incems"
463 include 'epos.inc'
464
465 call utpri('conre ',ish,ishini,6)
466
467c proj
468c ----
469 la=laproj
470 ma=iabs(maproj)
471 las=0
472 mas=0
473 do l=1,ma
474 if(la.lt.0)then
475 ia=iabs(idproj/10)
476 is=idproj/iabs(idproj)
477 if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
478 if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
479 if(ia.eq.213)id=1230*is
480 if(iabs(idproj).lt.20)id=idproj
481 else
482 id=1220
483 if(rangen().le.(la-las)*1./(ma-mas))id=1120
484 if(id.eq.1120)las=las+1
485 mas=mas+1
486 endif
487 ic1=idtrai(1,id,1)
488 ic2=idtrai(2,id,1)
489 icproj(1,l)=ic1
490 icproj(2,l)=ic2
491 enddo
492
493c targ
494c ----
495 la=latarg
496 ma=iabs(matarg)
497 las=0
498 mas=0
499 do l=1,ma
500 if(la.lt.0)then
501 ia=iabs(idtarg/10)
502 is=idtarg/iabs(idtarg)
503 if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
504 if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
505 if(ia.eq.213)id=1230*is
506 if(iabs(idtarg).lt.20)id=idtarg
507 else
508 id=1220
509 if(rangen().le.(la-las)*1./(ma-mas))id=1120
510 if(id.eq.1120)las=las+1
511 mas=mas+1
512 endif
513 ic1=idtrai(1,id,1)
514 ic2=idtrai(2,id,1)
515 ictarg(1,l)=ic1
516 ictarg(2,l)=ic2
517 enddo
518
519 call utprix('conre ',ish,ishini,6)
520 return
521 end
522
523c-----------------------------------------------------------------------
524 subroutine conrl
525c-----------------------------------------------------------------------
526c initializes target remnant in case of appl lepton
527c-----------------------------------------------------------------------
528 include "epos.incems"
529 common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx
530 common/hadr2/iomodl,idproj,idtarg,wexcit
531
532c targ
533c ----
534 la=latarg
535 ma=iabs(matarg)
536 las=0
537 mas=0
538 do l=1,ma
539 if(la.lt.0)then
540 id=idtarg
541 else
542 id=1220
543 if(rangen().le.(la-las)*1./(ma-mas))id=1120
544 if(id.eq.1120)las=las+1
545 mas=mas+1
546 endif
547 ic1=idtrai(1,id,1)
548 ic2=idtrai(2,id,1)
549 ictarg(1,l)=ic1
550 ictarg(2,l)=ic2
551 enddo
552
553 return
554 end
555
556c-----------------------------------------------------------------------
557 subroutine conwr
558c-----------------------------------------------------------------------
559c writes /cptl/
560c-----------------------------------------------------------------------
561 include "epos.inc"
562 include "epos.incems"
563 double precision XA(64,3),XB(64,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
564 COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
565 common/nucl3/phi,bimp
566 common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
567 *,xtarg(mamx),ytarg(mamx),ztarg(mamx)
568 parameter(iapmax=207)
569 double precision bqgs2,bmaxqgs2,bmaxnex2,bminnex2,xan,xbn
570 common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
571 *,bqgs2,bmaxqgs2,bmaxnex2,bminnex2
572 integer ic(2)
573
574 call utpri('conwr ',ish,ishini,6)
575
576 bx=cos(phi)*bimp
577 by=sin(phi)*bimp
578
579c write /cptl/
580c ------------
581 nptl=0
582
583 if(iokoll.eq.1)then ! precisely matarg collisions
584
585 nptl=nptl+1
586 do 3 i=1,4
5873 xorptl(i,nptl)=0
588 tivptl(1,nptl)=-ainfin
589 tivptl(2,nptl)=0
590 istptl(nptl)=1
591 iorptl(nptl)=-1
592 jorptl(nptl)=0
593 do 1 k=1,koll
594 nptl=nptl+1
595 do 4 i=1,4
5964 xorptl(i,nptl)=0
597 tivptl(1,nptl)=-ainfin
598 tivptl(2,nptl)=0
599 istptl(nptl)=1
600 iorptl(nptl)=-1
601 jorptl(nptl)=0
6021 continue
603
604 elseif(iappl.ne.7)then
605
606 do 6 i=1,maproj
607 nptl=nptl+1
608 istptl(nptl)=0
609 iorptl(nptl)=0
610 jorptl(nptl)=0
611 if(model.eq.2)then !QGSJet
612 xproj(i)=XA(i,1)
613 yproj(i)=XA(i,2)
614 zproj(i)=XA(i,3)
615 istptl(nptl)=1
616 iorptl(nptl)=-1
617 elseif(model.eq.7)then !QGSJetII
618 xproj(i)=xan(i,1)
619 yproj(i)=xan(i,2)
620 zproj(i)=xan(i,3)
621 istptl(nptl)=1
622 iorptl(nptl)=-1
623 elseif(model.ge.3)then !Gheisha, ...
624 istptl(nptl)=1
625 iorptl(nptl)=-1
626 endif
627 xorptl(1,nptl)=xproj(i)+bx/2
628 xorptl(2,nptl)=yproj(i)+by/2
629 xorptl(3,nptl)=zproj(i)
630 xorptl(4,nptl)=0
631 tivptl(1,nptl)=-ainfin
632c for visualisation uncomment
633c-c tivptl(1,nptl)=-100
634 tivptl(2,nptl)= ainfin
6356 continue
636 do 7 i=1,matarg
637 nptl=nptl+1
638 istptl(nptl)=0
639 iorptl(nptl)=0
640 jorptl(nptl)=0
641 if(model.eq.2)then !QGSJet
642 xtarg(i)=XB(i,1)
643 ytarg(i)=XB(i,2)
644 ztarg(i)=XB(i,3)
645 istptl(nptl)=1
646 iorptl(nptl)=-1
647 elseif(model.eq.7)then !QGSJetII
648 xtarg(i)=xbn(i,1)
649 ytarg(i)=xbn(i,2)
650 ztarg(i)=xbn(i,3)
651 istptl(nptl)=1
652 iorptl(nptl)=-1
653 elseif(model.ge.3)then !Gheisha, ...
654 istptl(nptl)=1
655 iorptl(nptl)=-1
656 endif
657 xorptl(1,nptl)=xtarg(i)-bx/2
658 xorptl(2,nptl)=ytarg(i)-by/2
659 xorptl(3,nptl)=ztarg(i)
660 xorptl(4,nptl)=0
661 tivptl(1,nptl)=-ainfin
662c for visualisation uncomment
663c-c tivptl(1,nptl)=-100
664 tivptl(2,nptl)= ainfin
6657 continue
666
667 endif
668
669 nptl=0
670 if(iappl.le.2)then
671 do i=1,maproj
672 nptl=nptl+1
673 ic(1)=icproj(1,i)
674 ic(2)=icproj(2,i)
675 id=idtra(ic,0,0,3)
676 call idmass(id,ams)
677 idptl(nptl)=id
678 pptl(1,nptl)=0.
679 pptl(2,nptl)=0.
680 pptl(3,nptl)=pnullx
681 pptl(4,nptl)=sqrt(pnullx**2+ams**2)
682 pptl(5,nptl)=ams
683 ifrptl(1,nptl)=0
684 ifrptl(2,nptl)=0
685 ityptl(nptl)=0
686 enddo
687 endif
688 if(iappl.ne.7)then
689 do i=1,matarg
690 nptl=nptl+1
691 ic(1)=ictarg(1,i)
692 ic(2)=ictarg(2,i)
693 id=idtra(ic,0,0,3)
694 call idmass(id,ams)
695 idptl(nptl)=id
696 pptl(1,nptl)=0.
697 pptl(2,nptl)=0.
698 pptl(3,nptl)=-pnullx
699 pptl(4,nptl)=sqrt(pnullx**2+ams**2)
700 pptl(5,nptl)=ams
701 ifrptl(1,nptl)=0
702 ifrptl(2,nptl)=0
703 ityptl(nptl)=0
704 enddo
705
706 else
707
708 nptl=nptl+1
709 id=idproj
710 call idmass(id,ams)
711 idptl(nptl)=id
712 pptl(1,nptl)=0.
713 pptl(2,nptl)=0.
714 pptl(3,nptl)=pnullx
715 pptl(4,nptl)=sqrt(pnullx**2+ams**2)
716 pptl(5,nptl)=ams
717 ifrptl(1,nptl)=0
718 ifrptl(2,nptl)=0
719 ityptl(nptl)=0
720 iorptl(nptl)=-1
721 jorptl(nptl)=0
722 istptl(nptl)=0
723 do 5 i=1,4
724 5 xorptl(i,nptl)=0
725 tivptl(1,nptl)=0
726 tivptl(2,nptl)=0
727 endif
728
729c exit
730c ----
731
732 call utprix('conwr ',ish,ishini,6)
733 return
734 end
735
736c------------------------------------------------------------------------
737 subroutine conxyz(ch,n,x,y,z,ynuc)
738c-----------------------------------------------------------------------
739 include 'epos.inc'
740
741 real x(n),y(n),z(n)
742 character ch*1
743
744 if(ch.eq.'p')then
745 massnr=maproj
746 iii=1
747 elseif(ch.eq.'t')then
748 massnr=matarg
749 iii=2
750 else
751 call utstop('conxyz: nucleus neither proj nor targ&')
752 endif
753
754 if(massnr.eq.0)return
755 if(massnr.gt.n)call utstop('conxyz: massnr.gt.n&')
756 if(massnr.eq.1)then
757 x(1)=0
758 y(1)=0
759 z(1)=0
760 return
761 endif
762
763 rad=radnuc(massnr)
764
765 if(massnr.ge.10)then !---wood-saxon density---
766
767 rad=rad/difnuc(massnr)
768 cr1=1.+3./rad+6./rad**2+6./rad**3
769 cr2=3./rad
770 cr3=3./rad+6./rad**2
771 do i=1,massnr
772 1 zuk=rangen()*cr1-1.
773 if(zuk.le.0.)then
774 tt=rad*(rangen()**.3333-1.)
775 elseif(zuk.le.cr2 )then
776 tt=-log(rangen())
777 elseif(zuk.lt.cr3 )then
778 tt=-log(rangen())-log(rangen())
779 else
780 tt=-log(rangen())-log(rangen())-log(rangen())
781 endif
782 if(rangen().gt.1./(1.+exp(-abs(tt))))goto 1
783 rim=tt+rad
784 zz=rim*(2.*rangen()-1.)
785 rim=sqrt(rim*rim-zz*zz)
786 z(i)=zz*difnuc(massnr)
787 call pscs(c,s)
788 x(i)=rim*c*difnuc(massnr)
789 y(i)=rim*s*difnuc(massnr)
790 enddo
791
792 elseif(massnr.ge.3)then ! ---gaussian density---
793
794 rad=rad*sqrt(2.*massnr/(massnr-1.)) !van hove simulation
795 do l=1,3
796 summ=0.
797 do i=1,massnr-1
798 j=massnr-i
799 aks=rad *(rangen()+rangen()+rangen()-1.5)
800 k=j+1
801 if(l.eq.1)x(k)=summ-aks*sqrt(float(j)/k)
802 if(l.eq.2)y(k)=summ-aks*sqrt(float(j)/k)
803 if(l.eq.3)z(k)=summ-aks*sqrt(float(j)/k)
804 summ=summ+aks/sqrt(float(j*k))
805 enddo
806 if(l.eq.1)x(1)=summ
807 if(l.eq.2)y(1)=summ
808 if(l.eq.3)z(1)=summ
809 enddo
810
811 elseif(massnr.eq.2)then ! ---deuteron---
812
813 !.........r=t*difnuc(massnr), t~exp(-2*t)*(1-exp(-a*t))
814 a=radnuc(massnr)
815 2 t=-0.5*alog(rangen()) !~exp(-2*t)
816 if(rangen().gt.(1-exp(-a*t))**2)goto2
817 r=t*difnuc(massnr)
818 zz=r*(2.*rangen()-1.)
819 call pscs(c,s)
820 rxy=sqrt(r*r-zz*zz)
821 z(1)=0.5*zz
822 x(1)=0.5*rxy*c
823 y(1)=0.5*rxy*s
824 z(2)=-z(1)
825 x(2)=-x(1)
826 y(2)=-y(1)
827
828 else
829
830 stop'conxyz: wrong massnr. '
831
832 endif
833
834c...plot preparation
835
836 rmax=(radnuc(massnr)+3)
837 drnucl(iii)=rmax/mxnucl
838 nrnucl(iii)=nrnucl(iii)+1
839 do i=1,massnr
840 r=sqrt(x(i)**2+y(i)**2+z(i)**2)
841 k=1+int(r/drnucl(iii))
842 if(k.le.mxnucl)rnucl(k,iii)=rnucl(k,iii)+1
843 enddo
844
845c...lorentz trafo
846
847 do i=1,massnr
848 z(i)=z(i)/cosh(ynuc)
849 enddo
850
851 return
852 end
853
854c-----------------------------------------------------------------------
855 subroutine conini
856c-----------------------------------------------------------------------
857 include 'epos.inc'
858
859 imax=max(maproj,matarg)
860 if(idtargin.eq.0)imax=max(imax,40)
861 do massnr=1,mamxx
862 dif=0.54
863 if(massnr.gt.imax.or.massnr.eq.1)then
864 dif=0
865 rad=0
866 elseif(massnr.eq.197)then
867 dif=0.562
868 rad=6.5
869 elseif(massnr.ge.10)then
870 rad=1.12*massnr**0.33333-0.86*massnr**(-0.33333)
871 elseif(massnr.ge.3)then
872 rad=.9*float(massnr)**.3333
873 elseif(massnr.eq.2)then
874 dif=4.316
875 rad=4.68
876 endif
877 difnuc(massnr)=dif
878 radnuc(massnr)=rad
879 enddo
880 end
881
882c-----------------------------------------------------------------------
883 subroutine xConNuclDens(iii)
884c-----------------------------------------------------------------------
885c plots distribution of nucleons in nuclei
886c iii = 1 (proj) or 2 (targ)
887c-----------------------------------------------------------------------
888 include 'epos.inc'
889 if(model.ne.1)return
890 if(iii.eq.1)then
891 massnr=maproj
892 elseif(iii.eq.2)then
893 massnr=matarg
894 endif
895 a=1./4.316
896 b=4.68
897 write(ifhi,'(a)') '!-----------------------------------------'
898 write(ifhi,'(a)') '! nuclear density '
899 write(ifhi,'(a)') '!-----------------------------------------'
900 write(ifhi,'(a)') 'openhisto'
901 if(massnr.ge.10)write(ifhi,'(a)')'htyp lin xmod lin ymod lin'
902 if(massnr.lt.10)write(ifhi,'(a)')'htyp lin xmod lin ymod log'
903 write(ifhi,'(a)') 'text 0 0 "title nuclear density"'
904 write(ifhi,'(a)') 'text 0.99 -0.15 " r (fm) "'
905 write(ifhi,'(a)') 'text 0 0 "yaxis rho(r)"'
906 write(ifhi,'(a,2e11.3)')'xrange',0.,mxnucl*drnucl(iii)
907 write(ifhi,'(3a)')'yrange',' 0 ',' auto'
908 write(ifhi,'(a)') 'array 2'
909 do j=1,mxnucl
910 r=(j-0.5)*drnucl(iii)
911 d=0.5*drnucl(iii)
912 write(ifhi,'(2e12.4)') r,rnucl(j,iii)/nrnucl(iii)/
913 * (4./3.*pi*((r+d)**3-(r-d)**3))
914 enddo
915 write(ifhi,'(a)') ' endarray'
916 write(ifhi,'(a)') 'closehisto plot 0-'
917 write(ifhi,'(a)') 'openhisto'
918 write(ifhi,'(a)') 'htyp lbo '
919 write(ifhi,'(a)') 'array 2'
920 do j=1,mxnucl
921 r=(j-0.5)*drnucl(iii)
922 rr=2*r
923 rho=1
924 if(massnr.eq.2)then
925 rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
926 elseif(massnr.eq.197)then
927 rho=0.16/(1+exp((r-6.5)/0.562))
928 elseif(massnr.ge.10)then
929 rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
930 endif
931 write(ifhi,'(2e12.4)') r,rho
932 enddo
933 write(ifhi,'(a)') ' endarray'
934 write(ifhi,'(a)') 'closehisto plot 0'
935 end
936
937c-----------------------------------------------------------------------
938 subroutine xConThick(iii)
939c-----------------------------------------------------------------------
940 ! plots sigma_pp *T_A (b) (=average number of collisions)
941 ! T_A = thickness function
942 ! iii = 1 (proj) or 2 (targ)
943 !----------------------------------------------------------------
944 include 'epos.inc'
945 parameter(iconimax=20,iconkmax=100)
946 real thick(2,0:iconimax)
947 imax=iconimax
948 kmax=iconkmax
949 if(model.ne.1)return
950 ramx=mxnucl*drnucl(iii)
951 do i=0,imax
952 x=i/float(imax)*ramx
953 sum=0
954 rho0=conrho(iii,0.)
955 h=ramx/kmax
956 do k=1,kmax
957 z=k*h
958 r=sqrt(x*x+z*z)
959 rho2=conrho(iii,r)
960 z=(k-0.5)*h
961 r=sqrt(x*x+z*z)
962 rho1=conrho(iii,r)
963 sum=sum+h/6.*(rho0+4*rho1+rho2)
964 rho0=rho2
965 enddo
966 sum=sum*2 ! integral fro -infty to +infty
967 thick(1,i)=x
968 thick(2,i)=sum
969 enddo
970 write(ifhi,'(a)') 'openhisto'
971 write(ifhi,'(a)') 'htyp lbo '
972 write(ifhi,'(a)') 'txt "xaxis b (fm)" '
973 write(ifhi,'(a)') 'txt "yaxis [s]?pp! T?A! (b) " '
974 write(ifhi,'(a)') 'array 2'
975 do i=0,imax
976 write(ifhi,'(2e12.4)') thick(1,i),sigine/10.*thick(2,i)
977 enddo
978 write(ifhi,'(a)') ' endarray'
979 write(ifhi,'(a)') 'closehisto plot 0'
980
981 end
982c-----------------------------------------------------------------------
983 function conrho(iii,r)
984c-----------------------------------------------------------------------
985 ! nuclear density
986 ! iii = 1 (proj) or 2 (targ)
987 !----------------------------------------------------------------
988 include 'epos.inc'
989 conrho=1.
990 if(model.ne.1)return
991 if(iii.eq.1)then
992 massnr=maproj
993 elseif(iii.eq.2)then
994 massnr=matarg
995 endif
996 a=1./4.316
997 b=4.68
998 rr=2*r
999 rho=1
1000 if(massnr.eq.2)then
1001 rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
1002 elseif(massnr.eq.197)then
1003 rho=0.16/(1+exp((r-6.5)/0.562))
1004 elseif(massnr.ge.10)then
1005 rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
1006 endif
1007 conrho=rho
1008 end
1009
1010
1011
1012
1013
1014
1015