]>
Commit | Line | Data |
---|---|---|
9ef1c2d9 | 1 | c----------------------------------------------------------------------- |
2 | subroutine conaa(iret) | |
3 | c----------------------------------------------------------------------- | |
4 | c determines interaction configuration | |
5 | c----------------------------------------------------------------------- | |
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 | ||
22 | c initialisations | |
23 | c --------------- | |
24 | ||
25 | vel=tanh(ypjtl-yhaha)+tanh(yhaha) | |
26 | ||
27 | c determine phi, bimp, coll, iproj, itarg, x/y/zproj, x/y/ztarg | |
28 | c --------------------------------------------------------------- | |
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 | |
140 | c 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 | ||
180 | 11 continue | |
181 | 12 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 | ||
199 | c determine coord | |
200 | c --------------- | |
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 | ||
213 | c exit | |
214 | c ---- | |
215 | 1000 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 | ||
230 | 1001 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 | ||
242 | c----------------------------------------------------------------------- | |
243 | subroutine xGeometry(iii) | |
244 | c----------------------------------------------------------------------- | |
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 | ||
337 | c----------------------------------------------------------------------- | |
338 | function conbmx() | |
339 | c----------------------------------------------------------------------- | |
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 | ||
353 | 10 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 | ||
381 | c----------------------------------------------------------------------- | |
382 | function conbmxndif() | |
383 | c----------------------------------------------------------------------- | |
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 | ||
401 | 10 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 | ||
430 | c----------------------------------------------------------------------- | |
431 | function conbmxdif() | |
432 | c----------------------------------------------------------------------- | |
433 | c find b to have (1-exp(-om))pmax=pdiff | |
434 | c----------------------------------------------------------------------- | |
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 | ||
457 | c----------------------------------------------------------------------- | |
458 | subroutine conre | |
459 | c----------------------------------------------------------------------- | |
460 | c initializes remnants | |
461 | c----------------------------------------------------------------------- | |
462 | include "epos.incems" | |
463 | include 'epos.inc' | |
464 | ||
465 | call utpri('conre ',ish,ishini,6) | |
466 | ||
467 | c proj | |
468 | c ---- | |
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 | ||
493 | c targ | |
494 | c ---- | |
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 | ||
523 | c----------------------------------------------------------------------- | |
524 | subroutine conrl | |
525 | c----------------------------------------------------------------------- | |
526 | c initializes target remnant in case of appl lepton | |
527 | c----------------------------------------------------------------------- | |
528 | include "epos.incems" | |
529 | common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx | |
530 | common/hadr2/iomodl,idproj,idtarg,wexcit | |
531 | ||
532 | c targ | |
533 | c ---- | |
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 | ||
556 | c----------------------------------------------------------------------- | |
557 | subroutine conwr | |
558 | c----------------------------------------------------------------------- | |
559 | c writes /cptl/ | |
560 | c----------------------------------------------------------------------- | |
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 | ||
579 | c write /cptl/ | |
580 | c ------------ | |
581 | nptl=0 | |
582 | ||
583 | if(iokoll.eq.1)then ! precisely matarg collisions | |
584 | ||
585 | nptl=nptl+1 | |
586 | do 3 i=1,4 | |
587 | 3 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 | |
596 | 4 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 | |
602 | 1 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 | |
632 | c for visualisation uncomment | |
633 | c-c tivptl(1,nptl)=-100 | |
634 | tivptl(2,nptl)= ainfin | |
635 | 6 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 | |
662 | c for visualisation uncomment | |
663 | c-c tivptl(1,nptl)=-100 | |
664 | tivptl(2,nptl)= ainfin | |
665 | 7 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 | ||
729 | c exit | |
730 | c ---- | |
731 | ||
732 | call utprix('conwr ',ish,ishini,6) | |
733 | return | |
734 | end | |
735 | ||
736 | c------------------------------------------------------------------------ | |
737 | subroutine conxyz(ch,n,x,y,z,ynuc) | |
738 | c----------------------------------------------------------------------- | |
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 | ||
834 | c...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 | ||
845 | c...lorentz trafo | |
846 | ||
847 | do i=1,massnr | |
848 | z(i)=z(i)/cosh(ynuc) | |
849 | enddo | |
850 | ||
851 | return | |
852 | end | |
853 | ||
854 | c----------------------------------------------------------------------- | |
855 | subroutine conini | |
856 | c----------------------------------------------------------------------- | |
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 | ||
882 | c----------------------------------------------------------------------- | |
883 | subroutine xConNuclDens(iii) | |
884 | c----------------------------------------------------------------------- | |
885 | c plots distribution of nucleons in nuclei | |
886 | c iii = 1 (proj) or 2 (targ) | |
887 | c----------------------------------------------------------------------- | |
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 | ||
937 | c----------------------------------------------------------------------- | |
938 | subroutine xConThick(iii) | |
939 | c----------------------------------------------------------------------- | |
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 | |
982 | c----------------------------------------------------------------------- | |
983 | function conrho(iii,r) | |
984 | c----------------------------------------------------------------------- | |
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 |