]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-uti-165.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
minor fixes
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-uti-165.f
CommitLineData
9ef1c2d9 1c-----------------------------------------------------------------------
2 subroutine utresc(iret)
3c-----------------------------------------------------------------------
60ed90b3 4c if irescl=1 rescaling is done, otherwise the purpose of going through
9ef1c2d9 5c this routine is to not change the seed in case of no rescaling
6c-----------------------------------------------------------------------
7 include 'epos.inc'
8 include 'epos.incems'
9 double precision p1,esoll,ppp,seedp
10 dimension p1(5),p0(5,mamx+mamx),pini(mxptl)
11 logical force,nolead(mxptl),lspec(mxptl)
12 data scalmean/0./
13 save scalmean
14 call utpri('utresc',ish,ishini,4)
15
16 errlim=0.001 !max(0.001,1./engy)
60ed90b3 17
9ef1c2d9 18 iret=0
19 scal=1.
20 nptlpt=iabs(maproj)+iabs(matarg)
21 call ranfgt(seedp) !not to change the seed ...
22 if(nptl.le.nptlpt) goto 9999
60ed90b3 23
9ef1c2d9 24 esoll=0.d0
25 p1(1)=0.d0
26 p1(2)=0.d0
27 p1(3)=0.d0
28 p1(4)=0.d0
29 ipmax=4
30 imin=nptlpt+1
31 if(iappl.eq.1)then
32 imin=1
33 ipmax=2
34 endif
35 do i=1,nptlpt
36 nolead(i)=.false.
37 do j=1,5
38 p0(j,i)=pptl(j,i)
39 enddo
40c if(mod(istptl(i),10).eq.1)then
41 do j=ipmax+1,4
42 p1(j)=p1(j)+dble(pptl(j,i))
43 enddo
44c endif
45 enddo
46 do i=nptlpt+1,nptl
47 if(mod(istptl(i),10).eq.0)then
48 do j=1,ipmax
49 p1(j)=p1(j)+dble(pptl(j,i))
50 enddo
51 lspec(i)=.false.
52 if((ityptl(i).eq.45.and.maproj.ge.100)
53 & .or.(ityptl(i).eq.55.and.matarg.ge.100))lspec(i)=.true.
54 if((abs(pptl(3,i)/pnullx).le.0.9
55 & .and.abs(pptl(3,i)).gt.pptl(5,i)).or.lspec(i))then
56 nolead(i)=.true.
57c write(ifch,*)'nolead',i
58 else
59 nolead(i)=.false.
60c write(ifch,*)'lead',i
61 endif
62 endif
63 enddo
64 ppp=(p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)*p1(2)-p1(1)*p1(1)
65 if(ppp.gt.0.d0)then
66 p1(5)=dsqrt(ppp)
67 else
68 iret=1
69 write(ifch,*)'p1',p1(1),p1(2),p1(3),p1(4),ppp
70 if(ish.ge.2)write (ifch,*) 'Problem in utresc (0): redo event'
71 write (ifmt,*) 'Problem in utresc (0): redo event'
72c call utstop('utresc&')
60ed90b3 73 goto 9999
9ef1c2d9 74 endif
75 esoll=p1(5)
76 if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector: ',p1
77
78c trafo
79c -----
80 pmax=0.d0
81 do i=imin,nptl
82 if(mod(istptl(i),10).le.1)then
83 call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
84 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
85 endif
86 if(i.gt.nptlpt.and.nolead(i))pmax=max(pmax,abs(pptl(3,i)))
87 enddo
88
89 if(ish.ge.6)then
90 call alistf('list after boost&')
60ed90b3 91 endif
9ef1c2d9 92
93 if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------'
94
95c rescale momenta in rest frame
96c -----------------------------
97 scal=1.
98 diff0=0.
99c ndif0=1
100 ferr=0.05
101 force=.false.
102 npart=nptl-imin+1
103 do ipass=1,300
104 sum=0.
105 sum3=0.
106 difft=diff0
107 ndif=0
108 nfirst=int(rangen()*float(npart)) !start list randomly
109 do i=0,npart-1
110 j=imin+i+nfirst
111 if(j.gt.nptl)j=imin+i+nfirst-npart
112 if(mod(istptl(j),10).eq.0)then
113c modified particles
114 if(nolead(j))then
115c if(j.gt.nptlpt)then
116c if(abs(pptl(3,j))/pnullx.lt.0.9)then !not spectator or diffraction
117 if(scal.eq.1..and.diff0.eq.0.)then
118 ndif=ndif+1
119 pini(j)=pptl(3,j)
120 else
121 pptl3new=0.
60ed90b3 122 if( force .or.(
9ef1c2d9 123 & ityptl(j)/10.eq.4.or.ityptl(j)/10.eq.5
124 & ))then !try just remnant first
125 ndif=ndif+1
126 diff=sign(min(0.3*abs(pini(j)),
60ed90b3 127 & rangen()*abs(difft)),difft)
9ef1c2d9 128 pptl3new=scal*(pptl(3,j)-diff)
129c write(ifch,*)'par',j,pptl3new,pptl(3,j),diff,difft
130c & ,ndif,pmax,scal
131 if(abs(pptl3new).lt.pmax)then
132 if(abs(pptl3new-pini(j)).lt.ferr*abs(pini(j))
133 & .or.(lspec(j).and.abs(pptl3new).lt.abs(0.8*pini(j))))then !particle should not be too fast or too modified
134c write(ifch,*)'used'
135 difft=difft-diff
136 pptl(3,j)=scal*(pptl(3,j)-diff)
137 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
138 * +pptl(3,j)**2+pptl(5,j)**2)
139 endif
140 endif
141 endif
142 endif
143 endif
144c sum over all particles
145 sum=sum+pptl(4,j)
146 sum3=sum3+pptl(3,j)
147 endif
148 enddo
149
150 diff=sum3
151 scal=real(esoll)/(sum-diff)
152 if(ish.ge.6)write(ifch,*)
153 $ 'ipass,scal,diff/pnullx,e,esoll,pz,pmax,ndif,f:'
154 $ ,ipass,scal,diff/pnullx,sum,esoll,sum3,pnullx,ndif,force
60ed90b3 155 if(abs(scal-1.).le.errlim.and.abs(diff/pnullx).lt.errlim)
9ef1c2d9 156 $ goto 300
157 if(ndif.gt.0.and.(force.or.ipass.lt.150))then
158c ndif0=ndif
159 diff0=diff
60ed90b3 160 elseif(abs(scal-1.).le.1e-2.and.abs(diff/pnullx).lt.5e-2)then
9ef1c2d9 161 goto 300
162 elseif(force)then
163 if(ish.ge.2)
164 $ write(ifmt,*)'Warning in utresc: no more allowed particle'
165 goto 302
166 else
167 force=.true.
168 ferr=0.1
169 diff=0.
170 endif
171 301 enddo
172 302 if(ish.ge.1)then
173 call utmsg('utrescl')
174 write(ifch,*)'***** scal=',scal,diff
175 call utmsgf
176 endif
177 if(abs(scal)+abs(diff/pnullx).gt.2.5)then
178 write(ifmt,*)'Warning in utresc !'
179 write(ifch,*)'redo EPOS event ...'
180 iret=1
181 goto 9999
182 endif
60ed90b3 183
9ef1c2d9 184c trafo
185c -----
186 300 continue
187 do i=1,nptl
188 if(i.le.nptlpt)then
189 do j=1,5
190 pptl(j,i)=p0(j,i)
191 enddo
192 else
193 if(mod(istptl(i),10).le.1)then
194 call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
195 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
196 endif
197 endif
198 enddo
60ed90b3 199
9ef1c2d9 200 if(ish.ge.4)call alist('list after rescaling&',1,nptl)
60ed90b3 201
9ef1c2d9 202 9999 continue
203 if(ish.ge.2)then
204 scalmean=scalmean+scal
205 write(ifch,*)' average rescaling factor: ',scalmean
206 & /float(nrevt+1)
60ed90b3 207 endif
208 call ranfst(seedp) ! ... after this subroutine
9ef1c2d9 209 call utprix('utresc',ish,ishini,4)
210
211 end
212
213c-----------------------------------------------------------------------
214 subroutine utghost(iret)
215c-----------------------------------------------------------------------
216c if irescl=1 make particle on-shell if not
217c-----------------------------------------------------------------------
218 include 'epos.inc'
219 include 'epos.incems'
220 double precision seedp
221 call utpri('ughost',ish,ishini,4)
222
223 iret=0
224 nptlpt=iabs(maproj)+iabs(matarg)
225 call ranfgt(seedp) !not to change the seed ...
226 if(nptl.le.nptlpt) goto 9999
60ed90b3 227
9ef1c2d9 228 if(ish.ge.5)write(ifch,'(a)')'---------mark ghosts---------'
229
230c mark ghosts
231c -----------
60ed90b3 232 do j=nptlpt+1,nptl
9ef1c2d9 233 if(mod(istptl(j),10).eq.0.and.pptl(4,j).gt.0.d0)then
234 amass=pptl(5,j)
235 call idmass(idptl(j),amass)
236 if(abs(pptl(5,j)-amass).gt.0.01*amass)then
237 if(ish.ge.5)write(ifch,*)'wrong particle mass',j,idptl(j)
238 & ,pptl(5,j),amass
239 amass=pptl(5,j)
240 call idres(idptl(j),amass,idr,iadj,0)
241 if(idr.ne.0)then
242 pptl(5,j)=amass
243 idptl(j)=idr
244 else
245 call idmass(idptl(j),amass)
246 pptl(5,j)=amass
247 endif
248 call idtau(idptl(j),pptl(4,j),pptl(5,j),taugm)
249 tivptl(2,j)=tivptl(1,j)+taugm*(-alog(rangen()))
250 endif
251 if(abs((pptl(4,j)+pptl(3,j))*(pptl(4,j)-pptl(3,j))
252 $ -pptl(2,j)**2-pptl(1,j)**2-pptl(5,j)**2).gt.0.3
253 $ .and.abs(1.-abs(pptl(3,j))/pptl(4,j)).gt.0.01)then
254 !print*,'ghost',ityptl(j),idptl(j)
255 if(ish.ge.1)write(ifmt,*)'ghost:',j,idptl(j),ityptl(j)
60ed90b3 256 if(ish.ge.5)then
9ef1c2d9 257 write(ifch,'(a,$)')'ghost:'
258 call alistc("&",j,j)
259 endif
260 ityptl(j)=100+ityptl(j)/10
60ed90b3 261 endif
9ef1c2d9 262 elseif(mod(istptl(j),10).eq.0)then
263 if(ish.ge.1)then
264 write(ifmt,*)'Lost particle (E=0)'
265 write(ifch,*)'Lost particle (E=0) :'
266 call alistc("utghost&",j,j)
267 endif
268 istptl(j)=istptl(j)+2
269 endif
270 enddo
271
272 if(ish.ge.5)write(ifch,'(a)')'---------treat ghosts---------'
60ed90b3 273
9ef1c2d9 274c treat ghosts
275c ------------
276 ifirst=1
277 scal=1.
278 pfif=0.
279 efif=0.
280 ntry=0
281 132 nfif=0
282 psum=0
283 esum=0.
284 ntry=ntry+1
60ed90b3 285 do j=nptlpt+1,nptl
9ef1c2d9 286 if(mod(istptl(j),10).eq.0)then
287 if(ityptl(j).ge.101.and.ityptl(j).le.105)then
288 nfif=nfif+1
60ed90b3 289 if(ifirst.eq.1)then
9ef1c2d9 290 pfif=pfif+pptl(3,j)
291 if(pptl(4,j).gt.0.)efif=efif+pptl(4,j)
292 endif
60ed90b3 293 if(irescl.ge.1) then
9ef1c2d9 294 if(ifirst.gt.1)then
295 if(pptl(4,j).gt.0.)then
296 Einv=1./pptl(4,j)
297 amt=1.-(pptl(5,j)*Einv)**2+(pptl(1,j)*Einv)**2
298 $ +(pptl(2,j)*Einv)**2
299 else
300 amt=-1.
301 endif
302 if(amt.gt.0.)then
303 pptl(3,j)=sign(pptl(4,j),pptl(3,j))*sqrt(amt)
304 else
305 y=(rangen()+rangen()+rangen()+rangen()-2.)/2.*yhaha
306 y=sign(abs(y),pptl(3,j))
307 pptl(3,j)
308 $ =sqrt(pptl(5,j)**2+pptl(1,j)**2
309 $ +pptl(2,j)**2)*sinh(y)
310 pptl(4,j)
311 $ =sqrt(pptl(5,j)**2+pptl(1,j)**2
312 $ +pptl(2,j)**2)*cosh(y)
313 efif=efif+pptl(4,j)
314 endif
315 ifirst=0
316 else
317c do k=1,3
318 do k=3,3
319 pptl(k,j)=pptl(k,j)*scal
320 enddo
321 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
322 * +pptl(5,j)**2)
323 endif
324 endif
325 psum=psum+pptl(3,j)
326 esum=esum+pptl(4,j)
327 if(ish.ge.5)
328 $ write (ifch,*) 'nrevt,psum,esum,pfif,efif,nfif,scal'
329 $ ,nrevt,psum,esum,pfif,efif,nfif,scal
60ed90b3 330 endif
9ef1c2d9 331 endif
332 enddo
333 if ( ish.ge.5 ) write (ifch,*) 'tot',nfif,efif,pfif,esum,psum
334
335
336 if(nfif.gt.5.or.(esum.gt.0.05*engy.and.nfif.ne.1))then
337 if(ifirst.eq.0)then
338 do j=nptlpt+1,nptl
339 if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
340 if((psum-pfif)*(1.-scal).ge.0)
341 & pptl(3,j)=pptl(3,j)-(psum-pfif)/nfif
342 endif
343 enddo
344 else
345 ifirst=2
346 goto 132
347 endif
348 scal=efif/esum
349 if ( ish.ge.5 ) write (ifch,*) 'scal',scal
350 if ( abs(scal-1.) .gt. 0.05 ) then
351 if(ntry.le.1000)then
352 goto 132
353 else
354 iret=1
355 if(ish.ge.2)write (ifch,*) 'Problem in utghost : redo event'
356 if(ish.ge.1)write (ifmt,*) 'Problem in utghost : redo event'
357 goto 9999
358 endif
359 endif
360 else
361 do j=nptlpt+1,nptl
362 if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
363 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
364 * +pptl(5,j)**2)
365 endif
366 enddo
367 endif
368
369 if(ish.ge.5)write(ifch,'(a)')'---------Check Ghost list---------'
370
371c Check Ghost list
372
60ed90b3 373 if(ish.ge.5)then
374 do j=nptlpt+1,nptl
9ef1c2d9 375 if(mod(istptl(j),10).eq.0)then
376 if(ityptl(j).le.105.and.ityptl(j).ge.101)then
377 write(ifch,'(a,$)')'ghost:'
378 call alistc("&",j,j)
379 endif
380 endif
381 enddo
382 endif
383
60ed90b3 384
9ef1c2d9 385 9999 continue
60ed90b3 386 call ranfst(seedp) ! ... after this subroutine
9ef1c2d9 387 call utprix('ughost',ish,ishini,4)
388
389 end
390
391c-----------------------------------------------------------------------
392 subroutine utrsph(iret)
393c-----------------------------------------------------------------------
60ed90b3 394c if irescl=1 and ispherio=1 rescaling is done for particle used by
9ef1c2d9 395c spherio as initial condition.
396c-----------------------------------------------------------------------
397 include 'epos.inc'
398 include 'epos.incems'
399 double precision p1,esoll
400 dimension p1(5),p0(5,mamx+mamx)
401 call utpri('utrsph',ish,ishini,4)
402
403 errlim=0.0001
60ed90b3 404
9ef1c2d9 405 iret=0
406 nptlpt=iabs(maproj)+iabs(matarg)
407 if(nptl.le.nptlpt) goto 9999
60ed90b3 408
9ef1c2d9 409 esoll=0.d0
410 p1(1)=0.d0
411 p1(2)=0.d0
412 p1(3)=0.d0
413 p1(4)=0.d0
414 do i=nptlpt+1,nptl
415 if((istptl(i).le.11
416 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
60ed90b3 417 $ .or.istptl(i).eq.20.or.istptl(i).eq.21)then
9ef1c2d9 418 do j=1,2
419 p1(j)=p1(j)+dble(pptl(j,i))
420 enddo
421 endif
422 enddo
423 do i=1,nptlpt
424 do j=1,5
425 p0(j,i)=pptl(j,i)
426 enddo
427 do j=3,4
428 p1(j)=p1(j)+dble(pptl(j,i))
429 enddo
430 enddo
431 p1(5)=dsqrt((p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)**2.d0-p1(1)**2.d0)
432 esoll=p1(5)
433 if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector',p1
434
435c trafo
436c -----
437 do i=1,nptl
438 if((istptl(i).le.11
439 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
440 $ .or.istptl(i).eq.20.or.istptl(i).eq.21
60ed90b3 441 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then
9ef1c2d9 442 call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
443 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
444 endif
445 enddo
446
447
448 if(ish.ge.5)write(ifch,'(a)')'------------------'
449
450c rescale momenta in rest frame
451c -----------------------------
452
453 scal=1.
454 diff=0.
455 do ipass=1,1000
456 sum=0.
457 sum3=0.
458 ndif=0
60ed90b3 459 do j=1,nptl
9ef1c2d9 460 if((istptl(j).le.11
461 $ .and.(iorptl(j).ge.1.and.istptl(iorptl(j)).eq.41))
462 $ .or.istptl(j).eq.20.or.istptl(j).eq.21
60ed90b3 463 $ .or.(istptl(j).eq.0.and.j.le.nptlpt))then
9ef1c2d9 464 if(j.gt.nptlpt)then
465 ndif=ndif+1
466 pptl(3,j)=scal*(pptl(3,j)-diff)
467 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
468 * +pptl(5,j)**2)
469 endif
470 sum=sum+pptl(4,j)
471 sum3=sum3+pptl(3,j)
472 endif
473 enddo
474
475 diff=sum3/real(ndif)
476 scal=real(esoll)/sum
477 if(ish.ge.6)write(ifch,*)'ipass,scal,diff,e,esoll,pz,ndif:'
478 $ ,ipass,scal,diff,sum,esoll,sum3,ndif
479 if(abs(scal-1.).le.errlim.and.abs(diff).lt.10.*errlim) goto300
480 301 enddo
481 if(ish.ge.1)then
482 call utmsg('hresph')
483 write(ifch,*)'***** scal=',scal,diff
484 call utmsgf
485 endif
486
60ed90b3 487
9ef1c2d9 488c trafo
489c -----
490 300 continue
491c do i=nptlpt+1,nptl
492 do i=1,nptl
493 if((istptl(i).le.11
494 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
495 $ .or.istptl(i).eq.20.or.istptl(i).eq.21
60ed90b3 496 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then
9ef1c2d9 497 call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
498 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
499 endif
500 if(i.le.nptlpt)then
501 do j=1,5
502 pptl(j,i)=p0(j,i)
503 enddo
504 endif
505 enddo
60ed90b3 506
9ef1c2d9 507 9999 call utprix('utrsph',ish,ishini,4)
508
509 end
510
511cc-----------------------------------------------------------------------
512c double precision function dddlog(xxx)
513cc-----------------------------------------------------------------------
514c double precision xxx
515c dddlog=-1d50
516c if(xxx.gt.0d0)dddlog=dlog(xxx)
517c end
60ed90b3 518c
9ef1c2d9 519ccc-----------------------------------------------------------------------
520c subroutine randfl(jc,iqa0,iflav,ic,isame)
521cc-----------------------------------------------------------------------
522cc returns random flavour ic(2) (iqa0=1:quark,2:antiquark,11:diquark)
523cc-----------------------------------------------------------------------
524c include 'epos.inc'
525c real probab(nflav),probsu(nflav+1)
526c integer jc(nflav,2),jc0(nflav,2),ic(2)
527c if(ish.ge.6)then
528c write(ifch,*)('-',i=1,10)
529c *,' entry sr randfl ',('-',i=1,30)
530c write(ifch,*)'iqa0:',iqa0
531c write(ifch,*)'jc:'
532c write(ifch,*)jc
533c endif
534c iflav=0
535c ic(1)=0
536c ic(2)=0
537c do 10 n=1,nflav
538c do 10 i=1,2
539c10 jc0(n,i)=0
540c iqa1=iqa0*10
541c9999 iqa1=iqa1/10
542c if(iqa1.eq.0)goto9998
543c iqa=mod(iqa1,10)
544c su=0
545c do 20 i=1,nflav
546c probab(i)=jc(i,iqa)-jc0(i,iqa)
547c if(isame.eq.1)probab(i)=probab(i)*(jc(i,3-iqa)-jc0(i,3-iqa))
548c20 su=su+probab(i)
549c if(su.lt..5)then
550c iflav=0
551c ic(1)=0
552c ic(2)=0
553c goto9998
554c endif
555c probsu(1)=0.
556c do 30 i=1,nflav
557c probsu(i+1)=probsu(i)+probab(i)/su
558c if(probsu(i+1)-probsu(i).lt.1e-5)probsu(i+1)=probsu(i)
559c30 continue
560c r=rangen()*probsu(nflav+1)
561c do 50 i=1,nflav
562c if(probsu(i).le.r.and.r.lt.probsu(i+1))iflav=i
563c50 continue
564c jc0(iflav,iqa)=jc0(iflav,iqa)+1
565c if(isame.eq.1)jc0(iflav,3-iqa)=jc0(iflav,3-iqa)+1
566c call idenco(jc0,ic,ireten)
567c if(ireten.eq.1)call utstop('randfl: idenco ret code = 1&')
568c if(ish.ge.6)then
569c write(ifch,*)'probab:'
570c write(ifch,*)probab
571c write(ifch,*)'probsu:'
572c write(ifch,*)probsu
573c write(ifch,*)'ran#:',r,' flav:',iflav
574c endif
575c goto9999
576c9998 continue
577c if(ish.ge.6)write(ifch,*)('-',i=1,30)
578c *,' exit sr randfl ',('-',i=1,10)
579c return
580c end
60ed90b3 581c
582c
9ef1c2d9 583cc-----------------------------------------------------------------------
584c subroutine ranhvy(x,eps)
585cc-----------------------------------------------------------------------
586cc generates x for heavy particle fragmentation according to
587cc the peterson form
588cc d(x)=1/(x*(1-1/x-eps/(1-x))**2)
589cc =d0(x)*d1(x)*d2(x)
590cc d0(x)=(1-x)**2/((1-x)**2+eps)**2
591cc d1(x)=x
592cc d2(x)=(((1-x)**2+eps)/((1-x)**2+eps*x))**2
593cc using x=1-y**pow
594cc generates flat in x if eps>1.
595cc-----------------------------------------------------------------------
596c data aln4/1.3863/
597c if(eps.lt.1.) then
598c pow=alog((3.+eps)/eps)/aln4
599c ymx=(eps*(3.*pow-1.)/(pow+1.))**(.5/pow)
600c zmx=1-ymx**pow
601c d0mx=(1-zmx)**2/((1.-zmx)**2+eps)**2*pow*ymx**(pow-1.)
602c d2mx=2./(2.-sqrt(eps))
603c else
604c pow=1.
605c zmx=0.
606c d0mx=(1.-zmx)**2/((1.-zmx)**2+eps)**2
607c d2mx=1.+eps
608c endif
609cc
610cc generate z according to (1-z)**2/((1-z)**2+eps*z)**2
611c1 continue
612c y=rangen()
613c z=1.-y**pow
614cc
615c d0z=(1.-z)**2/((1.-z)**2+eps)**2*pow*y**(pow-1.)
616c if(d0z.lt.rangen()*d0mx) goto1
617cc
618cc check remaining factors
619c d1=z
620c d2=(((1.-z)**2+eps)/((1.-z)**2+eps*z))**2
621c if(d1*d2.lt.rangen()*d2mx) goto1
622cc
623cc good x
624c x=z
625c return
626c end
60ed90b3 627c
9ef1c2d9 628c-----------------------------------------------------------------------
629 function ransig()
630c-----------------------------------------------------------------------
631c returns randomly +1 or -1
632c-----------------------------------------------------------------------
633 ransig=1
634 if(rangen().gt.0.5)ransig=-1
635 return
636 end
60ed90b3 637
9ef1c2d9 638cc-----------------------------------------------------------------------
639c function ranxq(n,x,q,xmin)
640cc-----------------------------------------------------------------------
641cc returns random number according to x(i) q(i) with x>=xmin
642cc-----------------------------------------------------------------------
643c include 'epos.inc'
644c real x(n),q(n)
645c imin=1
646c if(xmin.eq.0.)goto3
647c i1=1
648c i2=n
649c1 i=i1+(i2-i1)/2
650c if(x(i).lt.xmin)then
651c i1=i
652c elseif(x(i).gt.xmin)then
653c i2=i
654c else
655c imin=i
656c goto3
657c endif
658c if(i2-i1.gt.1)goto1
659c imin=i2
660c3 continue
661c if(q(imin).gt.q(n)*.9999)then
662c ranxq=xmin
663c goto4
664c endif
665c qran=q(imin)+rangen()*(q(n)-q(imin))
666c ranxq=utinvt(n,x,q,qran)
667c4 continue
60ed90b3 668c
9ef1c2d9 669c if(ranxq.lt.xmin)then
670c call utmsg('ranxq ')
671c write(ifch,*)'***** ranxq=',ranxq,' < xmin=',xmin
672c write(ifch,*)'q(imin) q q(n):',q(imin),qran,q(n)
673c write(ifch,*)'x(imin) x x(n):',x(imin),ranxq,x(n)
674c call utmsgf
675c ranxq=xmin
676c endif
60ed90b3 677c
9ef1c2d9 678c return
679c end
60ed90b3 680c
9ef1c2d9 681cc ***** end r-routines
60ed90b3 682cc ***** beg s-routines
9ef1c2d9 683c
684cc-----------------------------------------------------------------------
685c function sbet(z,w)
686cc-----------------------------------------------------------------------
687c sbet=utgam1(z)*utgam1(w)/utgam1(z+w)
688c return
689c end
60ed90b3 690c
9ef1c2d9 691cc-----------------------------------------------------------------------
692c function smass(a,y,z)
693cc-----------------------------------------------------------------------
694cc returns droplet mass (in gev) (per droplet, not (!) per nucleon)
695cc according to berger/jaffe mass formula, prc35(1987)213 eq.2.31,
696cc see also c. dover, BNL-46322, intersections-meeting, tucson, 91.
697cc a: massnr, y: hypercharge, z: charge,
698cc-----------------------------------------------------------------------
699c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
700c ymin=ym*a
701c zmin=cz/(dz/a+zm/a**(1./3.))
702c smass=epsi*a+as*a**(2./3.)+(ac/a**(1./3.)+dz/a/2.)*(z-zmin)**2
703c *+dy/a/2.*(y-ymin)**2
704c return
705c end
60ed90b3 706c
9ef1c2d9 707cc-----------------------------------------------------------------------
708c subroutine smassi(theta)
709cc-----------------------------------------------------------------------
710cc initialization for smass.
711cc calculates parameters for berger/jaffe mass formula
712cc (prc35(1987)213 eq.2.31, see also c. dover, BNL-46322).
713cc theta: parameter that determines all parameters in mass formula.
714cc-----------------------------------------------------------------------
715c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
716c thet=theta
60ed90b3 717c
9ef1c2d9 718c astr=.150
719c pi=3.14159
720c alp=1./137.
60ed90b3 721c
9ef1c2d9 722c co=cos(theta)
723c si=sin(theta)
724c bet=(1+co**3)/2.
725c rzero=si/astr/( 2./3./pi*(1+co**3) )**(1./3.)
726cctp060829 cs=astr/si
727c cz=-astr/si*( ( .5*(1+co**3) )**(1./3.)-1 )
728c sigma=6./8./pi*(astr/si)**3*(co**2/6.-si**2*(1-si)/3.-
729c *1./3./pi*(pi/2.-theta-sin(2*theta)+si**3*alog((1+co)/si)))
60ed90b3 730c
9ef1c2d9 731c epsi=astr*((.5*(1+co**3))**(1./3.)+2)/si
732c as=4*pi*sigma*rzero**2
733c ac=3./5.*alp/rzero
734c dz=astr/si*bet**(1./3.)*co**2*
735c *(co**4*(1+bet**(2./3.))+(1+bet)**2)/
736c *( (2*co**2+bet**(1./3.))*(co**4*(1+bet**(2./3.))+(1+bet)**2)-
737c *(co**4+bet**(1./3.)*(1+bet))*((2*bet**(2./3.)-1)*co**2+1+bet) )
738c dy=astr/6.*(1+co**3)**3/si*
739c *( 1+(1+co)/(4*(1+co**3))**(2./3.) )/
740c *(co**6+co+co*(.5*(1+co**3))**(4./3.))
741c zm=6*alp/(5*rzero)
742c ym=(1-co**3)/(1+co**3)
60ed90b3 743c
9ef1c2d9 744c return
745c end
60ed90b3 746c
9ef1c2d9 747cc-----------------------------------------------------------------------
748c subroutine smassp
749cc-----------------------------------------------------------------------
750cc prints smass.
751cc-----------------------------------------------------------------------
752c include 'epos.inc'
753c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
754c real eng(14),ymi(14),zmi(14)
755c pi=3.14159
756c write(ifch,*)'parameters of mass formula:'
757c write(ifch,*)'theta=',thet,' epsi=',epsi
758c write(ifch,*)'as=',as,' ac=',ac
759c write(ifch,*)'dy=',dy,' dz=',dz
760c write(ifch,*)'ym=',ym
761c write(ifch,*)'cz dz zm=',cz,dz,zm
762c write(ifch,*)'sigma**1/3=',sigma**(1./3.),' rzero=',rzero
763c write(ifch,*)'mass:'
764c write(ifch,5000)(j,j=1,14)
765c5000 format(5x,'a:',14i5)
766c do 4 j=1,14
767c a=j
768c ymi(j)=ym*a
769c4 zmi(j)=cz/(dz/a+zm/a**(1./3.))
770c write(ifch,5002)(ymi(j),j=1,14)
771c5002 format(1x,'ymin: ',14f5.2)
772c write(ifch,5003)(zmi(j),j=1,14)
773c5003 format(1x,'zmin: ',14f5.2)
774c do 2 i=1,15
775c ns=11-i
776c do 3 j=1,14
777c a=j
778c y=a-ns
779c z=0.
780c3 eng(j)=smass(a,y,z)/a
781c write(ifch,5001)ns,(eng(j),j=1,14)
782c5001 format(1x,'s=',i2,2x,14f5.2)
783c2 continue
784c write(ifch,*)'mass-mass(free):'
785c write(ifch,5000)(j,j=1,14)
786c do 5 i=1,15
787c ns=11-i
788c do 6 j=1,14
789c a=j
790c y=a-ns
791c z=0.
792c call smassu(a,y,z,ku,kd,ks,kc)
793c6 eng(j)=(smass(a,y,z)-utamnu(ku,kd,ks,kc,0,0,3))/a
794c write(ifch,5001)ns,(eng(j),j=1,14)
795c5 continue
60ed90b3 796c
9ef1c2d9 797c stop
798c end
60ed90b3 799c
9ef1c2d9 800cc-----------------------------------------------------------------------
801c subroutine smasst(kux,kdx,ksx,kcx,a,y,z)
802cc-----------------------------------------------------------------------
803cc input: kux,kdx,ksx,kcx = net quark numbers (for u,d,s,c quarks).
804cc output: massnr a, hypercharge y and charge z.
805cc-----------------------------------------------------------------------
806c sg=1
807c if(kux+kdx+ksx+kcx.lt.0.)sg=-1
808c ku=sg*kux
809c kd=sg*kdx
810c ks=sg*ksx
811c kc=sg*kcx
812c k=ku+kd+ks+kc
813c if(mod(k,3).ne.0)stop'noninteger baryon number'
814c a=k/3
815c y=a-ks
816c nz=2*ku-kd-ks+2*kc
817c if(mod(nz,3).ne.0)stop'noninteger charge'
818c z=nz/3
819c return
820c end
60ed90b3 821c
9ef1c2d9 822cc-----------------------------------------------------------------------
823c subroutine smassu(ax,yx,zx,ku,kd,ks,kc)
824cc-----------------------------------------------------------------------
825cc input: massnr ax, hypercharge yx and charge zx.
826cc output: ku,kd,ks,kc = net quark numbers (for u,d,s,c quarks).
827cc-----------------------------------------------------------------------
828c sg=1
829c if(ax.lt.0.)sg=-1
830c a=sg*ax
831c y=sg*yx
832c z=sg*zx
833c ku=nint(a+z)
834c kd=nint(a-z+y)
835c ks=nint(a-y)
836c kc=0
837c return
838c end
60ed90b3 839c
9ef1c2d9 840cc-----------------------------------------------------------------------
841c function spoc(a,b,c,d,x)
842cc-----------------------------------------------------------------------
843cc power fctn with cutoff
844cc-----------------------------------------------------------------------
845c spoc=0
846c if(a.eq.0..and.b.eq.0.)return
847c spoc =a+b*x**c
848c spoc0=a+b*d**c
849c spoc=amin1(spoc,spoc0)
850c spoc=amax1(0.,spoc)
851c return
852c end
853c
854c-----------------------------------------------------------------------
855 function utacos(x)
856c-----------------------------------------------------------------------
857c returns acos(x) for -1 <= x <= 1 , acos(+-1) else
858c-----------------------------------------------------------------------
859 include 'epos.inc'
860 argum=x
861 if(x.lt.-1.)then
862 if(ish.ge.1)then
863 call utmsg('utacos')
864 write(ifch,*)'***** argum = ',argum,' set -1'
865 call utmsgf
866 endif
867 argum=-1.
868 elseif(x.gt.1.)then
869 if(ish.ge.1)then
870 call utmsg('utacos')
871 write(ifch,*)'***** argum = ',argum,' set 1'
872 call utmsgf
873 endif
874 argum=1.
875 endif
876 utacos=acos(argum)
877 return
878 end
60ed90b3 879
9ef1c2d9 880c----------------------------------------------------------------------
881 function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus)
882c----------------------------------------------------------------------
883c returns min mass of droplet with given u,d,s,c content
884c keux: net u quark number
885c kedx: net d quark number
886c kesx: net s quark number
887c kecx: net c quark number
888c kebx: net b quark number
889c ketx: net t quark number
890c modus: 4=two lowest multiplets; 5=lowest multiplet
891c----------------------------------------------------------------------
892 common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr
893 common/csjcga/amnull,asuha(7)
894 common/drop4/asuhax(7),asuhay(7)
60ed90b3 895
9ef1c2d9 896 if(modus.lt.4.or.modus.gt.5)stop'UTAMNU: not supported'
897 1 format(' flavours:',6i5 )
898 100 format(' flavours+mass:',6i5,f8.2 )
899c write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx c write(ifmt,*)'wrong mass in gethad ',damss
900
901 amnull=0.
60ed90b3 902
9ef1c2d9 903 do i=1,7
904 if(modus.eq.4)asuha(i)=asuhax(i) !two lowest multiplets
905 if(modus.eq.5)asuha(i)=asuhay(i) !lowest multiplet
906 enddo
907
908 ke=iabs(keux+kedx+kesx+kecx+kebx+ketx)
60ed90b3 909
9ef1c2d9 910 if(keux+kedx+kesx+kecx+kebx+ketx.ge.0)then
911 keu=keux
912 ked=kedx
913 kes=kesx
914 kec=kecx
915 keb=kebx
916 ket=ketx
917 else
918 keu=-keux
919 ked=-kedx
920 kes=-kesx
921 kec=-kecx
922 keb=-kebx
923 ket=-ketx
924 endif
925
60ed90b3 926c write(ifch,*)keu,ked,kes,kec,keb,ket
9ef1c2d9 927
60ed90b3 928c removing top mesons to remove t quarks or antiquarks
9ef1c2d9 929 if(ket.ne.0)then
93012 continue
931 ii=sign(1,ket)
932 ket=ket-ii
933 if(ii*keu.le.ii*ked)then
934 keu=keu+ii
935 else
936 ked=ked+ii
937 endif
938 amnull=amnull+200. ! ???????
939 if(ket.ne.0)goto12
940 endif
941
60ed90b3 942c removing bottom mesons to remove b quarks or antiquarks
9ef1c2d9 943 if(keb.ne.0)then
94411 continue
945 ii=sign(1,keb)
946 keb=keb-ii
947 if(ii*keu.le.ii*ked)then
948 keu=keu+ii
949 else
950 ked=ked+ii
951 endif
952 amnull=amnull+5.28 ! (B-meson)
953 if(keb.ne.0)goto11
954 endif
955
956c removing charm mesons to remove c quarks or antiquarks
957 if(kec.ne.0)then
95810 continue
959 ii=sign(1,kec)
960 kec=kec-ii
961 if(keu*ii.le.ked*ii)then
962 keu=keu+ii
963 else
964 ked=ked+ii
965 endif
966 amnull=amnull+1.87 ! (D-meson)
967 if(kec.ne.0)goto10
968 endif
9ef1c2d9 969
60ed90b3 970c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
971
972c removing mesons to remove s antiquarks
9ef1c2d9 9735 continue
974 if(kes.lt.0)then
975 amnull=amnull+asuha(6)
976 if(keu.ge.ked)then
977 keu=keu-1
978 else
979 ked=ked-1
980 endif
981 kes=kes+1
982 goto5
983 endif
60ed90b3 984
985c removing mesons to remove d antiquarks
9ef1c2d9 9866 continue
987 if(ked.lt.0)then
988 if(keu.ge.kes)then
989 amnull=amnull+asuha(5)
990 keu=keu-1
991 else
992 amnull=amnull+asuha(6)
993 kes=kes-1
994 endif
995 ked=ked+1
996 goto6
997 endif
60ed90b3 998
999c removing mesons to remove u antiquarks
9ef1c2d9 10007 continue
1001 if(keu.lt.0)then
1002 if(ked.ge.kes)then
1003 amnull=amnull+asuha(5)
1004 ked=ked-1
1005 else
1006 amnull=amnull+asuha(6)
1007 kes=kes-1
1008 endif
1009 keu=keu+1
1010 goto7
1011 endif
60ed90b3 1012
1013c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
1014c print*,keu,ked,kes,kec,keb,ket,amnull
1015
9ef1c2d9 1016 if(keu+ked+kes+kec+keb+ket.ne.ke)
1017 *call utstop('utamnu: sum_kei /= ke&')
1018 keq=keu+ked
1019 keqx=keq
1020 amnux=0
60ed90b3 1021
1022c removing strange baryons
9ef1c2d9 1023 i=4
10242 i=i-1
10253 continue
1026 if((4-i)*kes.gt.(i-1)*keq)then
1027 amnux=amnux+asuha(1+i)
1028 kes=kes-i
1029 keq=keq-3+i
1030 if(kes.lt.0)call utstop('utamnu: negative kes&')
1031 if(keq.lt.0)call utstop('utamnu: negative keq&')
1032 goto3
1033 endif
1034 if(i.gt.1)goto2
1035 if(keqx.gt.keq)then
1036 do 8 k=1,keqx-keq
1037 if(keu.ge.ked)then
1038 keu=keu-1
1039 else
1040 ked=ked-1
1041 endif
10428 continue
1043 endif
60ed90b3 1044
9ef1c2d9 1045 if(keu+ked.ne.keq)call utstop('utamnu: keu+ked /= keq&')
60ed90b3 1046c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1047c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1048
1049c removing nonstrange baryons
9ef1c2d9 10509 continue
1051 if(keu.gt.2*ked)then
1052 amnux=amnux+asuha(7)
1053 keu=keu-3
1054 if(keu.lt.0)call utstop('utamnu: negative keu&')
1055 goto9
1056 endif
1057 if(ked.gt.2*keu)then
1058 amnux=amnux+asuha(7)
1059 ked=ked-3
1060 if(ked.lt.0)call utstop('utamnu: negative ked&')
1061 goto9
1062 endif
1063 keq=keu+ked
60ed90b3 1064
1065c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1066c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1067
9ef1c2d9 1068 if(mod(keq,3).ne.0)call utstop('utamnu: mod(keq,3) /= 0&')
1069 amnux=amnux+asuha(1)*keq/3
60ed90b3 1070
1071c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1072c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1073
9ef1c2d9 1074 amnull=amnull+amnux
60ed90b3 1075
9ef1c2d9 1076 if(amnull.eq.0)amnull=asuha(5)
60ed90b3 1077
9ef1c2d9 10781000 utamnu=amnull
1079 return
1080 end
60ed90b3 1081
9ef1c2d9 1082c-----------------------------------------------------------------------
1083 function utamnx(jcp,jcm)
1084c-----------------------------------------------------------------------
1085c returns minimum mass for the decay of jcp---jcm (by calling utamnu).
1086c-----------------------------------------------------------------------
1087 parameter (nflav=6)
1088 integer jcp(nflav,2),jcm(nflav,2)
1089
1090 do i=1,nflav
1091 do j=1,2
1092 if(jcp(i,j).ne.0)goto1
1093 enddo
1094 enddo
1095 keu=jcm(1,1)-jcm(1,2)
1096 ked=jcm(2,1)-jcm(2,2)
1097 kes=jcm(3,1)-jcm(3,2)
1098 kec=jcm(4,1)-jcm(4,2)
1099 keb=jcm(5,1)-jcm(5,2)
1100 ket=jcm(6,1)-jcm(6,2)
1101 utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1102 return
11031 continue
1104
1105 do i=1,nflav
1106 do j=1,2
1107 if(jcm(i,j).ne.0)goto2
1108 enddo
1109 enddo
1110 keu=jcp(1,1)-jcp(1,2)
1111 ked=jcp(2,1)-jcp(2,2)
1112 kes=jcp(3,1)-jcp(3,2)
1113 kec=jcp(4,1)-jcp(4,2)
1114 keb=jcp(5,1)-jcp(5,2)
1115 ket=jcp(6,1)-jcp(6,2)
1116 utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1117 return
11182 continue
1119
1120 keu=jcp(1,1)-jcp(1,2)
1121 ked=jcp(2,1)-jcp(2,2)
1122 kes=jcp(3,1)-jcp(3,2)
1123 kec=jcp(4,1)-jcp(4,2)
1124 keb=jcp(5,1)-jcp(5,2)
1125 ket=jcp(6,1)-jcp(6,2)
1126 ke=keu+ked+kes+kec+keb+ket
1127 if(mod(ke+1,3).eq.0)then
1128 keu=keu+1
1129 amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1130 keu=keu-1
1131 ked=ked+1
1132 amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1133 elseif(mod(ke-1,3).eq.0)then
1134 keu=keu-1
1135 amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1136 keu=keu+1
1137 ked=ked-1
1138 amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1139 else
1140 call utstop('utamnx: no singlet possible (1)&')
1141 endif
1142 keu=jcm(1,1)-jcm(1,2)
1143 ked=jcm(2,1)-jcm(2,2)
1144 kes=jcm(3,1)-jcm(3,2)
1145 kec=jcm(4,1)-jcm(4,2)
1146 keb=jcm(5,1)-jcm(5,2)
1147 ket=jcm(6,1)-jcm(6,2)
1148 ke=keu+ked+kes+kec+keb+ket
1149 if(mod(ke+1,3).eq.0)then
1150 keu=keu+1
1151 amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1152 keu=keu-1
1153 ked=ked+1
1154 amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1155 elseif(mod(ke-1,3).eq.0)then
1156 keu=keu-1
1157 amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1158 keu=keu+1
1159 ked=ked-1
1160 amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1161 else
1162 call utstop('utamnx: no singlet possible (2)&')
1163 endif
1164 utamnx=min(amms1+amms3,amms2+amms4)
1165c print *,amms1,amms3,amms2,amms4,jcp,jcm
1166 return
1167 end
1168
1169
60ed90b3 1170
9ef1c2d9 1171cc-----------------------------------------------------------------------
1172c function utamny(jcp,jcm)
1173cc-----------------------------------------------------------------------
1174cc returns minimum mass of jcp+jcm (by calling utamnu).
1175cc-----------------------------------------------------------------------
1176c parameter (nflav=6)
1177c integer jcp(nflav,2),jcm(nflav,2),jc(nflav,2)
1178c do 7 nf=1,nflav
1179c jc(nf,1)=jcp(nf,1)+jcm(nf,1)
1180c7 jc(nf,2)=jcp(nf,2)+jcm(nf,2)
1181c keu=jc(1,1)-jc(1,2)
1182c ked=jc(2,1)-jc(2,2)
1183c kes=jc(3,1)-jc(3,2)
1184c kec=jc(4,1)-jc(4,2)
1185c keb=jc(5,1)-jc(5,2)
1186c ket=jc(6,1)-jc(6,2)
1187c utamny=utamnu(keu,ked,kes,kec,keb,ket,5)
1188c return
1189c end
1190c
1191c-----------------------------------------------------------------------
1192 function utamnz(jc,modus)
1193c-----------------------------------------------------------------------
1194c returns minimum mass of jc (by calling utamnu).
1195c-----------------------------------------------------------------------
1196 parameter (nflav=6)
1197 integer jc(nflav,2)
1198 keu=jc(1,1)-jc(1,2)
1199 ked=jc(2,1)-jc(2,2)
1200 kes=jc(3,1)-jc(3,2)
1201 kec=jc(4,1)-jc(4,2)
1202 keb=jc(5,1)-jc(5,2)
1203 ket=jc(6,1)-jc(6,2)
1204 utamnz=utamnu(keu,ked,kes,kec,keb,ket,modus)
1205 return
1206 end
60ed90b3 1207
9ef1c2d9 1208c-----------------------------------------------------------------------
1209 subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx)
1210c-----------------------------------------------------------------------
1211c returns the array xx with xx(1)=x0 <= xx(i) <= xx(i3)=x3
1212c-----------------------------------------------------------------------
1213 real xx(i3)
1214 do 1 i=1,i1-1
1215 1 xx(i)=x0+(i-1.)/(i1-1.)*(x1-x0)
1216 do 2 i=i1,i2-1
1217 2 xx(i)=x1+(i-i1*1.)/(i2-i1*1.)*(x2-x1)
1218 do 3 i=i2,i3
1219 3 xx(i)=x2+(i-i2*1.)/(i3-i2*1.)*(x3-x2)
1220 return
1221 end
1222
1223cc---------------------------------------------------------------------
1224c subroutine utaxis(i,j,a1,a2,a3)
1225cc-----------------------------------------------------------------------
1226cc calculates the axis defined by the ptls i,j in the i,j cm system
1227cc---------------------------------------------------------------------
1228c include 'epos.inc'
1229c double precision pi1,pi2,pi3,pi4,pj1,pj2,pj3,pj4,p1,p2,p3,p4,p5
60ed90b3 1230c *,err,a
9ef1c2d9 1231c a1=0
1232c a2=0
1233c a3=1
1234c pi1=dble(pptl(1,i))
60ed90b3 1235c pi2=dble(pptl(2,i))
1236c pi3=dble(pptl(3,i))
1237c pi4=dble(pptl(4,i))
1238c pj1=dble(pptl(1,j))
1239c pj2=dble(pptl(2,j))
1240c pj3=dble(pptl(3,j))
1241c pj4=dble(pptl(4,j))
9ef1c2d9 1242c p1=pi1+pj1
1243c p2=pi2+pj2
1244c p3=pi3+pj3
1245c p4=pi4+pj4
1246c p5=dsqrt(p4**2-p3**2-p2**2-p1**2)
60ed90b3 1247c call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50)
1248c call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51)
1249c err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2
9ef1c2d9 1250c if(err.gt.1d-3)then
60ed90b3 1251c call utmsg('utaxis')
9ef1c2d9 1252c write(ifch,*)'***** err=',err
1253c write(ifch,*)'pi:',pi1,pi2,pi3,pi4
1254c write(ifch,*)'pj:',pj1,pj2,pj3,pj4
60ed90b3 1255c call utmsgf
9ef1c2d9 1256c endif
1257c a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 )
1258c if(a.eq.0.d0)return
1259c a1=sngl((pi1-pj1)/a)
1260c a2=sngl((pi2-pj2)/a)
1261c a3=sngl((pi3-pj3)/a)
1262c return
1263c end
60ed90b3 1264c
9ef1c2d9 1265cc-----------------------------------------------------------------------
1266c subroutine utchm(arp,arm,ii)
1267cc-----------------------------------------------------------------------
1268cc checks whether arp**2=0 and arm**2=0.
1269cc-----------------------------------------------------------------------
1270c include 'epos.inc'
1271c double precision arp(4),arm(4),difp,difm
1272c difp=arp(4)**2-arp(1)**2-arp(2)**2-arp(3)**2
1273c difm=arm(4)**2-arm(1)**2-arm(2)**2-arm(3)**2
1274c if(dabs(difp).gt.1e-3*arp(4)**2
1275c *.or.dabs(difm).gt.1e-3*arm(4)**2)then
1276c call utmsg('utchm ')
1277c write(ifch,*)'***** mass non zero - ',ii
1278c write(ifch,*)'jet-mass**2`s: ',difp,difm
1279c write(ifch,*)'energy**2`s: ',arp(4)**2,arm(4)**2
1280c write(ifch,*)(sngl(arp(i)),i=1,4)
1281c write(ifch,*)(sngl(arm(i)),i=1,4)
1282c call utmsgf
1283c endif
1284c return
1285c end
60ed90b3 1286c
9ef1c2d9 1287c-----------------------------------------------------------------------
1288 subroutine utclea(nptlii,nptl0)
1289c-----------------------------------------------------------------------
1290c starting from nptlii
60ed90b3 1291c overwrites istptl=99 particles in /cptl/, reduces so nptl
9ef1c2d9 1292c and update minfra and maxfra
1293c-----------------------------------------------------------------------
1294 include 'epos.inc'
1295 integer newptl(mxptl)!,oldptl(mxptl),ii(mxptl)
60ed90b3 1296
9ef1c2d9 1297 ish0=ish
1298 if(ishsub/100.eq.18)ish=mod(ishsub,100)
60ed90b3 1299
9ef1c2d9 1300 call utpri('utclea',ish,ishini,2)
1301
1302 nptli=max(maproj+matarg+1,nptlii)
1303 minfra0=minfra
1304 maxfra0=maxfra
1305 minfra1=maxfra
1306 maxfra1=minfra
1307 if(ish.ge.2)write(ifch,*)'entering subr utclea:',nptl
1308 & ,minfra,maxfra
1309 if(ish.ge.7)then
1310 write(ifch,*)('-',l=1,68)
1311 write(ifch,*)'sr utclea. initial.'
1312 write(ifch,*)('-',l=1,68)
1313 do 34 n=nptli,nptl
1314 write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1315 *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1316 *,istptl(n),ityptl(n)
131734 continue
1318116 format(1x,i6,i6,4x,i6,4x,i6,i6,i12,3x,3(e8.2,1x),i3,i3)
1319 endif
1320
1321c ish=ish0
1322c ish0=ish
1323c if(ishsub/100.eq.18)ish=mod(ishsub,100)
1324
1325 i=nptli-1
13261 i=i+1
1327 if(i.gt.nptl)goto 1000
1328 if(istptl(i).eq.99)goto 2
1329 newptl(i)=i
1330c oldptl(i)=i
1331 goto 1
60ed90b3 1332
9ef1c2d9 13332 i=i-1
1334 j=i
13353 i=i+1
13364 j=j+1
1337 if(j.gt.nptl)goto 5
1338 newptl(j)=0
1339 if(istptl(j).eq.99)goto 4
1340 newptl(j)=i
1341c oldptl(i)=j
1342c write(ifch,*)'move',j,' to ',i
1343c write(ifch,*)idptl(i),ityptl(i),idptl(j),ityptl(j),minfra,maxfra
60ed90b3 1344 call utrepl(i,j)
9ef1c2d9 1345 if(j.ge.minfra0.and.j.le.maxfra0)then
1346 minfra1=min(minfra1,i)
1347 maxfra1=max(maxfra1,i)
1348 endif
1349 goto 3
60ed90b3 1350
9ef1c2d9 13515 nptl=i-1
1352 if(nptl.eq.0)then
60ed90b3 1353 nptl0=0
9ef1c2d9 1354 goto 1000
1355 endif
1356
135720 n0=newptl(nptl0)
1358 if(n0.gt.0)then
1359 nptl0=n0
1360 else
1361 nptl0=nptl0-1
1362 if(nptl0.gt.0)goto 20
1363 endif
1364
1365
1366c do 11 k=1,nptl
1367c io=iorptl(k)
1368c if(io.le.0)ii(k)=io
1369c if(io.gt.0)ii(k)=newptl(io)
1370c11 continue
1371c do 12 k=1,nptl
1372c12 iorptl(k)=ii(k)
60ed90b3 1373c
9ef1c2d9 1374c do 13 k=1,nptl
1375c jo=jorptl(k)
1376c if(jo.le.0)ii(k)=jo
1377c if(jo.gt.0)ii(k)=newptl(jo)
1378c13 continue
1379c do 14 k=1,nptl
1380c14 jorptl(k)=ii(k)
60ed90b3 1381c
9ef1c2d9 1382c do 15 k=1,nptl
1383c if1=ifrptl(1,k)
1384c if(if1.le.0)ii(k)=if1
1385c if(if1.gt.0)ii(k)=newptl(if1)
1386c15 continue
1387c do 16 k=1,nptl
1388c16 ifrptl(1,k)=ii(k)
60ed90b3 1389c
9ef1c2d9 1390c do 17 k=1,nptl
1391c if2=ifrptl(2,k)
1392c if(if2.le.0)ii(k)=if2
1393c if(if2.gt.0)ii(k)=newptl(if2)
1394c17 continue
1395c do 18 k=1,nptl
1396c18 ifrptl(2,k)=ii(k)
60ed90b3 1397c
9ef1c2d9 1398c do 19 k=1,nptl
1399c if(ifrptl(1,k).eq.0.and.ifrptl(2,k).gt.0)ifrptl(1,k)=ifrptl(2,k)
1400c if(ifrptl(2,k).eq.0.and.ifrptl(1,k).gt.0)ifrptl(2,k)=ifrptl(1,k)
1401c19 continue
60ed90b3 1402
9ef1c2d9 14031000 continue
60ed90b3 1404
9ef1c2d9 1405 if(minfra1.lt.minfra0)minfra=minfra1
1406 if(maxfra1.ge.minfra1)maxfra=maxfra1
1407
1408 if(ish.ge.2)then
1409 write(ifch,*)'before exiting subr utclea:'
1410 do 35 n=1,nptl
1411 write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1412 *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1413 *,istptl(n),ityptl(n)
141435 continue
1415 endif
60ed90b3 1416
9ef1c2d9 1417 if(ish.ge.2)write(ifch,*)'exiting subr utclea:',nptl
1418 & ,minfra,maxfra
1419
1420 call utprix('utclea',ish,ishini,2)
1421 ish=ish0
1422 return
1423 end
1424
1425c---------------------------------------------------------------------
1426 subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
1427c---------------------------------------------------------------------
60ed90b3 1428c linear fit to data
9ef1c2d9 1429c input:
1430c ndata: nr of data points
1431c x(),y(),sig(): data
1432c mwt: unweighted (0) or weighted (else) data points
1433c output:
1434c a,b: parameters of linear fit a+b*x
60ed90b3 1435c---------------------------------------------------------------------
9ef1c2d9 1436 INTEGER mwt,ndata
1437 REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
60ed90b3 1438CU USES utgmq
9ef1c2d9 1439 INTEGER i
1440 REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,utgmq
1441 sx=0.
1442 sy=0.
1443 st2=0.
1444 b=0.
1445 if(mwt.ne.0) then
1446 ss=0.
1447 do 11 i=1,ndata
1448 wt=1./(sig(i)**2)
1449 ss=ss+wt
1450 sx=sx+x(i)*wt
1451 sy=sy+y(i)*wt
145211 continue
1453 else
1454 do 12 i=1,ndata
1455 sx=sx+x(i)
1456 sy=sy+y(i)
145712 continue
1458 ss=float(ndata)
1459 endif
1460 sxoss=sx/ss
1461 if(mwt.ne.0) then
1462 do 13 i=1,ndata
1463 t=(x(i)-sxoss)/sig(i)
1464 st2=st2+t*t
1465 b=b+t*y(i)/sig(i)
146613 continue
1467 else
1468 do 14 i=1,ndata
1469 t=x(i)-sxoss
1470 st2=st2+t*t
1471 b=b+t*y(i)
147214 continue
1473 endif
1474 b=b/st2
1475 a=(sy-sx*b)/ss
1476 siga=sqrt((1.+sx*sx/(ss*st2))/ss)
1477 sigb=sqrt(1./st2)
1478 chi2=0.
1479 if(mwt.eq.0) then
1480 do 15 i=1,ndata
1481 chi2=chi2+(y(i)-a-b*x(i))**2
148215 continue
1483 q=1.
1484 sigdat=sqrt(chi2/(ndata-2))
1485 siga=siga*sigdat
1486 sigb=sigb*sigdat
1487 else
1488 do 16 i=1,ndata
1489 chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
149016 continue
1491 q=utgmq(0.5*(ndata-2),0.5*chi2)
1492 endif
1493 return
1494 END
60ed90b3 1495
9ef1c2d9 1496c-----------------------------------------------------------------------
1497 function utgam1(x)
1498c-----------------------------------------------------------------------
1499c gamma fctn tabulated
1500c single precision
1501c-----------------------------------------------------------------------
1502 double precision utgamtab,utgam,al,dl
1503 common/gamtab/utgamtab(10000)
1504
1505 if(x.gt.0.01.and.x.lt.99.99)then
1506 al=100.d0*dble(x)
1507 k1=int(al)
1508 k2=k1+1
1509 dl =al-dble(k1)
1510 utgam1=real(utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl))
1511 elseif(x.eq.0.)then
1512 utgam1=0.
1513 else
1514 utgam1=real(utgam(dble(x)))
1515 endif
1516
1517 end
1518
1519c-----------------------------------------------------------------------
1520 double precision function utgam2(x)
1521c-----------------------------------------------------------------------
1522c gamma fctn tabulated
1523c double precision
1524c-----------------------------------------------------------------------
1525 double precision utgamtab,x,al,dl,utgam
1526 common/gamtab/utgamtab(10000)
1527
1528 if(x.gt.0.01d0.and.x.le.99.99d0)then
1529 al=100.d0*x
1530 k1=int(al)
1531 k2=k1+1
1532 dl =al-dble(k1)
1533 utgam2=utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl)
1534 elseif(x.eq.0.d0)then
1535 utgam2=0.d0
1536 else
1537 utgam2=utgam(x)
1538 endif
1539
1540 end
1541
1542c-----------------------------------------------------------------------
1543 double precision function utgam(x)
1544c-----------------------------------------------------------------------
1545c gamma fctn
1546c double precision
1547c-----------------------------------------------------------------------
1548 include 'epos.inc'
1549 double precision c(13),x,z,f
1550 data c
1551 1/ 0.00053 96989 58808, 0.00261 93072 82746, 0.02044 96308 23590,
1552 2 0.07309 48364 14370, 0.27964 36915 78538, 0.55338 76923 85769,
1553 3 0.99999 99999 99998,-0.00083 27247 08684, 0.00469 86580 79622,
1554 4 0.02252 38347 47260,-0.17044 79328 74746,-0.05681 03350 86194,
1555 5 1.13060 33572 86556/
1556 utgam=0d0
1557 z=x
1558 if(x .gt. 170.d0) goto6
1559 if(x .gt. 0.0d0) goto1
1560 if(x .eq. int(x)) goto5
1561 z=1.0d0-z
1562 1 f=1.0d0/z
1563 if(z .le. 1.0d0) goto4
1564 f=1.0d0
1565 2 continue
1566 if(z .lt. 2.0d0) goto3
1567 z=z-1.0d0
1568 f=f*z
1569 goto2
1570 3 z=z-1.0d0
1571 4 utgam=
1572 1 f*((((((c(1)*z+c(2))*z+c(3))*z+c(4))*z+c(5))*z+c(6))*z+c(7))/
60ed90b3 1573 2 ((((((c(8)*z+c(9))*z+c(10))*z+c(11))*z+c(12))*z+c(13))*z+1.0d0)
9ef1c2d9 1574 if(x .gt. 0.0d0) return
1575 utgam=3.141592653589793d0/(sin(3.141592653589793d0*x)*utgam)
1576 return
1577 5 write(ifch,10)sngl(x)
1578 10 format(1x,'argument of gamma fctn = ',e20.5)
1579 call utstop('utgam : negative integer argument&')
1580 6 write(ifch,11)sngl(x)
1581 11 format(1x,'argument of gamma fctn = ',e20.5)
1582 call utstop('utgam : argument too large&')
1583 end
1584
1585c---------------------------------------------------------------------
1586 subroutine utgcf(gammcf,a,x,gln)
1587c---------------------------------------------------------------------
1588 INTEGER ITMAX
1589 REAL a,gammcf,gln,x,EPS,FPMIN
1590 PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
1591CU USES utgmln
1592 INTEGER i
1593 REAL an,b,c,d,del,h,utgmln
1594 gln=utgmln(a)
1595 b=x+1.-a
1596 c=1./FPMIN
1597 d=1./b
1598 h=d
1599 do 11 i=1,ITMAX
1600 an=-i*(i-a)
1601 b=b+2.
1602 d=an*d+b
1603 if(abs(d).lt.FPMIN)d=FPMIN
1604 c=b+an/c
1605 if(abs(c).lt.FPMIN)c=FPMIN
1606 d=1./d
1607 del=d*c
1608 h=h*del
1609 if(abs(del-1.).lt.EPS)goto 1
161011 continue
1611 pause 'a too large, ITMAX too small in utgcf'
16121 gammcf=exp(-x+a*log(x)-gln)*h
1613 return
1614 END
1615
1616c---------------------------------------------------------------------
1617 function utgmln(xx)
1618c---------------------------------------------------------------------
1619 REAL utgmln,xx
1620 INTEGER j
1621 DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
1622 SAVE cof,stp
1623 DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
1624 *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
1625 *-.5395239384953d-5,2.5066282746310005d0/
1626 x=xx
1627 y=x
1628 tmp=x+5.5d0
1629 tmp=(x+0.5d0)*log(tmp)-tmp
1630 ser=1.000000000190015d0
1631 do 11 j=1,6
1632 y=y+1.d0
1633 ser=ser+cof(j)/y
163411 continue
1635 utgmln=tmp+log(stp*ser/x)
1636 return
1637 END
1638
1639c---------------------------------------------------------------------
1640 function utgmq(a,x)
1641c---------------------------------------------------------------------
1642 REAL a,utgmq,x
1643CU USES utgcf,utgser
1644 REAL gammcf,gamser,gln
1645 if(x.lt.0..or.a.le.0.)pause 'bad arguments in utgmq'
1646 if(x.lt.a+1.)then
1647 call utgser(gamser,a,x,gln)
1648 utgmq=1.-gamser
1649 else
1650 call utgcf(gammcf,a,x,gln)
1651 utgmq=gammcf
1652 endif
1653 return
1654 END
1655
1656c---------------------------------------------------------------------
1657 subroutine utgser(gamser,a,x,gln)
1658c---------------------------------------------------------------------
1659 INTEGER ITMAX
1660 REAL a,gamser,gln,x,EPS
1661 PARAMETER (ITMAX=100,EPS=3.e-7)
1662CU USES utgmln
1663 INTEGER n
1664 REAL ap,del,sum,utgmln
1665 gln=utgmln(a)
1666 if(x.le.0.)then
1667 if(x.lt.0.)pause 'x < 0 in utgser'
1668 gamser=0.
1669 return
1670 endif
1671 ap=a
1672 sum=1./a
1673 del=sum
1674 do 11 n=1,ITMAX
1675 ap=ap+1.
1676 del=del*x/ap
1677 sum=sum+del
1678 if(abs(del).lt.abs(sum)*EPS)goto 1
167911 continue
1680 pause 'a too large, ITMAX too small in utgser'
16811 gamser=sum*exp(-x+a*log(x)-gln)
1682 return
1683 END
60ed90b3 1684
9ef1c2d9 1685c-------------------------------------------------------------------------
1686 subroutine uticpl(ic,ifla,iqaq,iret)
60ed90b3 1687c-------------------------------------------------------------------------
9ef1c2d9 1688c adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla
1689c to 2-id ic
60ed90b3 1690c-------------------------------------------------------------------------
9ef1c2d9 1691 include 'epos.inc'
1692 integer jc(nflav,2),ic(2)
1693 iret=0
1694 if(ifla.eq.0)return
1695 call iddeco(ic,jc)
1696 if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
1697 jqaq=3-iqaq
1698 if(jc(ifla,jqaq).gt.0)then
1699 jc(ifla,jqaq)=jc(ifla,jqaq)-1
1700 else
1701 jc(ifla,iqaq)=jc(ifla,iqaq)+1
1702 endif
1703 call idcomj(jc)
1704 call idenco(jc,ic,ireten)
1705 if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
1706 if(ireten.eq.1)iret=1
1707 if(ic(1).eq.0.and.ic(2).eq.0.and.ireten.eq.0)then
1708 ic(1)=100000
1709 ic(2)=100000
1710 endif
1711 return
1712 end
1713
1714cc-----------------------------------------------------------------------
1715c subroutine utindx(n,xar,x,i)
1716cc-----------------------------------------------------------------------
60ed90b3 1717cc input: dimension n
1718cc array xar(n) with xar(i) > xar(i-1)
9ef1c2d9 1719cc some number x between xar(1) and xar(n)
1720cc output: the index i such that x is between xar(i) and xar(i+1)
1721cc-----------------------------------------------------------------------
1722c include 'epos.inc'
1723c real xar(n)
1724c if(x.lt.xar(1))then
1725c if(ish.ge.5)then
1726c call utmsg('utindx')
1727c write(ifch,*)'***** x=',x,' < xar(1)=',xar(1)
1728c call utmsgf
1729c endif
1730c i=1
1731c return
1732c elseif(x.gt.xar(n))then
1733c if(ish.ge.5)then
1734c call utmsg('utindx')
1735c write(ifch,*)'***** x=',x,' > xar(n)=',xar(n)
1736c call utmsgf
1737c endif
1738c i=n
1739c return
1740c endif
1741c lu=1
1742c lo=n
1743c1 lz=(lo+lu)/2
1744c if((xar(lu).le.x).and.(x.le.xar(lz)))then
1745c lo=lz
1746c elseif((xar(lz).lt.x).and.(x.le.xar(lo)))then
1747c lu=lz
1748c else
1749c call utstop('utindx: no interval found&')
1750c endif
1751c if((lo-lu).ge.2) goto1
1752c if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
1753c i=lu
1754c return
1755c end
60ed90b3 1756c
9ef1c2d9 1757c-----------------------------------------------------------------------
1758 function utinvt(n,x,q,y)
1759c-----------------------------------------------------------------------
1760c returns x with y=q(x)
1761c-----------------------------------------------------------------------
1762 include 'epos.inc'
1763 real x(n),q(n)
1764 if(q(n).eq.0.)call utstop('utinvt: q(n)=0&')
1765 if(y.lt.0.)then
1766 if(ish.ge.1)then
1767 call utmsg('utinvt')
1768 write(ifch,*)'***** y=',y,' < 0'
1769 call utmsgf
1770 endif
1771 y=0.
1772 elseif(y.gt.q(n))then
1773 if(ish.ge.1)then
1774 call utmsg('utinvt')
1775 write(ifch,*)'***** y=',y,' > ',q(n)
1776 call utmsgf
1777 endif
1778 y=q(n)
1779 endif
1780 lu=1
1781 lo=n
17821 lz=(lo+lu)/2
1783 if((q(lu).le.y).and.(y.le.q(lz)))then
1784 lo=lz
1785 elseif((q(lz).lt.y).and.(y.le.q(lo)))then
1786 lu=lz
1787 else
1788 write(ifch,*)'q(1),y,q(n):',q(1),y,q(n)
1789 write(ifch,*)'lu,lz,lo:',lu,lz,lo
1790 write(ifch,*)'q(lu),q(lz),q(lo):',q(lu),q(lz),q(lo)
1791 call utstop('utinvt: no interval found&')
1792 endif
1793 if((lo-lu).ge.2) goto1
1794 if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
1795 utinvt=x(lu)+(y-q(lu))*(x(lo)-x(lu))/(q(lo)-q(lu))
1796 return
1797 end
60ed90b3 1798
9ef1c2d9 1799c-----------------------------------------------------------------------
1800 subroutine utlob2(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4,idi)
1801c-----------------------------------------------------------------------
1802c performs a lorentz boost, double prec.
1803c isig=+1 is to boost the four vector x1,x2,x3,x4 such as to obtain it
1804c in the frame specified by the 5-vector p1...p5 (5-vector=4-vector+mass).
1805c isig=-1: the other way round, that means,
60ed90b3 1806c if the 4-vector x1...x4 is given in some frame characterized by
1807c p1...p5 with respect to to some lab-frame, utlob2 returns the 4-vector
9ef1c2d9 1808c x1...x4 in the lab frame.
1809c idi is a call identifyer (integer) to identify the call in case of problem
1810c-----------------------------------------------------------------------
1811 include 'epos.inc'
1812 double precision beta(4),z(4),p1,p2,p3,p4,p5,pp,bp,x1,x2,x3,x4
1813 *,xx0,x10,x20,x30,x40,x4x,x0123
1814 if(ish.ge.2)then
1815 if(ish.ge.9)then
1816 write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1817 write(ifch,301)p1,p2,p3,p4,p5,(p4-p3)*(p4+p3)-p2*p2-p1*p1
1818101 format(' utlob2: x = ',5e13.5)
1819301 format(' p = ',6e13.5)
1820 endif
1821 pp=(p4-p3)*(p4+p3)-p2*p2-p1*p1
1822 if(dabs(pp-p5*p5).gt.1e-3*p4*p4.and.dabs(pp-p5*p5).gt.1e-3)then
1823 call utmsg('utlob2')
1824 write(ifch,*)'***** p**2 .ne. p5**2'
1825 write(ifch,*)'call identifyer:',idi
1826 write(ifch,*)'p**2,p5**2: ',pp,p5*p5
1827 write(ifch,*)'p: ',p1,p2,p3,p4,p5
1828 call utmsgf
1829 endif
1830 x10=x1
1831 x20=x2
1832 x30=x3
1833 x40=x4
1834 endif
1835 xx0=(x4-x3)*(x4+x3)-x2*x2-x1*x1
1836 if(p5.le.0.)then
1837 call utmsg('utlob2')
1838 write(ifch,*)'***** p5 negative.'
1839 write(ifch,*)'call identifyer:',idi
1840 write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
1841 write(ifmt,*)'call identifyer:',idi
1842 write(ifmt,*)'p(5): ',p1,p2,p3,p4,p5
1843 call utmsgf
1844 call utstop('utlob2: p5 negative.&')
1845 endif
1846 z(1)=x1
1847 z(2)=x2
1848 z(3)=x3
1849 z(4)=x4
1850 beta(1)=-p1/p5
1851 beta(2)=-p2/p5
1852 beta(3)=-p3/p5
1853 beta(4)= p4/p5
1854 bp=0.
1855 do 220 k=1,3
1856220 bp=bp+z(k)*isig*beta(k)
1857 do 230 k=1,3
1858230 z(k)=z(k)+isig*beta(k)*z(4)
1859 *+isig*beta(k)*bp/(beta(4)+1.)
1860 z(4)=beta(4)*z(4)+bp
1861 x1=z(1)
1862 x2=z(2)
1863 x3=z(3)
1864 x4=z(4)
1865 if(ish.ge.9)
1866 *write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1867 x4x=x4
1868 x0123=xx0+x1*x1+x2*x2+x3*x3
1869 if(x0123.gt.0.)then
1870 x4=sign( dsqrt(x0123) , x4x )
1871 else
1872 x4=0
1873 endif
1874 if(ish.ge.9)then
1875 write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1876 endif
1877 if(ish.ge.2)then
1878 if(ish.ge.9)write(ifch,*)'check x**2_ini -- x**2_fin'
1879 if(dabs(x4-x4x).gt.1d-2*dabs(x4).and.dabs(x4-x4x).gt.1d-2)then
1880 call utmsg('utlob2')
1881 write(ifch,*)'***** x**2_ini .ne. x**2_fin.'
1882 write(ifch,*)'call identifyer:',idi
1883 write(ifch,*)'x1 x2 x3 x4 x**2 (initial/final/corrected):'
1884102 format(5e13.5)
1885 write(ifch,102)x10,x20,x30,x40,(x40-x30)*(x40+x30)-x20*x20-x10*x10
1886 write(ifch,102)x1,x2,x3,x4x,(x4x-x3)*(x4x+x3)-x2*x2-x1*x1
1887 write(ifch,102)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1888 call utmsgf
1889 endif
1890 endif
1891 if(ish.ge.9)write(ifch,*)'return from utlob2'
1892 return
1893 end
1894
1895c-----------------------------------------------------------------------
1896 subroutine utlob3(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
1897c-----------------------------------------------------------------------
1898c performs a lorentz boost, double prec.
1899c but arguments are single precision
1900c-----------------------------------------------------------------------
1901 double precision xx1,xx2,xx3,xx4
1902 xx1=dble(x1)
1903 xx2=dble(x2)
1904 xx3=dble(x3)
1905 xx4=dble(x4)
1906 call utlob2(isig
1907 *,dble(p1),dble(p2),dble(p3),dble(p4),dble(p5)
1908 *,xx1,xx2,xx3,xx4,52)
1909 x1=sngl(xx1)
1910 x2=sngl(xx2)
1911 x3=sngl(xx3)
1912 x4=sngl(xx4)
1913 return
1914 end
1915
1916c-----------------------------------------------------------------------
1917 subroutine utlob5(yboost,x1,x2,x3,x4,x5)
1918c-----------------------------------------------------------------------
1919 amt=sqrt(x5**2+x1**2+x2**2)
1920 y=sign(1.,x3)*alog((x4+abs(x3))/amt)
1921 y=y-yboost
1922 x4=amt*cosh(y)
1923 x3=amt*sinh(y)
1924 return
1925 end
1926
1927c-----------------------------------------------------------------------
1928 subroutine utlob4(isig,pp1,pp2,pp3,pp4,pp5,x1,x2,x3,x4)
1929c-----------------------------------------------------------------------
1930c performs a lorentz boost, double prec.
1931c but arguments are partly single precision
1932c-----------------------------------------------------------------------
1933 double precision xx1,xx2,xx3,xx4,pp1,pp2,pp3,pp4,pp5
1934 xx1=dble(x1)
1935 xx2=dble(x2)
1936 xx3=dble(x3)
1937 xx4=dble(x4)
1938 call utlob2(isig,pp1,pp2,pp3,pp4,pp5,xx1,xx2,xx3,xx4,53)
1939 x1=sngl(xx1)
1940 x2=sngl(xx2)
1941 x3=sngl(xx3)
1942 x4=sngl(xx4)
1943 return
1944 end
1945
60ed90b3 1946
9ef1c2d9 1947c-----------------------------------------------------------------------
1948 subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
1949c-----------------------------------------------------------------------
1950c performs a lorentz boost
1951c-----------------------------------------------------------------------
1952 include 'epos.inc'
1953 real beta(4),z(4)
1954 if(p5.le.0.)then
1955 call utmsg('utlobo')
1956 write(ifch,*)'***** mass <= 0.'
1957 write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
1958 call utmsgf
1959 call utstop('utlobo: mass <= 0.&')
1960 endif
1961 z(1)=x1
1962 z(2)=x2
1963 z(3)=x3
1964 z(4)=x4
1965 beta(1)=-p1/p5
1966 beta(2)=-p2/p5
1967 beta(3)=-p3/p5
1968 beta(4)= p4/p5
1969 bp=0.
1970 do 220 k=1,3
1971220 bp=bp+z(k)*isig*beta(k)
1972 do 230 k=1,3
1973230 z(k)=z(k)+isig*beta(k)*z(4)
1974 *+isig*beta(k)*bp/(beta(4)+1.)
1975 z(4)=beta(4)*z(4)+bp
1976 x1=z(1)
1977 x2=z(2)
1978 x3=z(3)
1979 x4=z(4)
1980 return
1981 end
60ed90b3 1982
9ef1c2d9 1983c-----------------------------------------------------------------------
1984 subroutine utloc(ar,n,a,l)
1985c-----------------------------------------------------------------------
1986 real ar(n)
1987 do 1 i=1,n
1988 l=i-1
1989 if(a.lt.ar(i))return
19901 continue
1991 l=n
1992 return
1993 end
60ed90b3 1994
9ef1c2d9 1995cc-----------------------------------------------------------------------
1996c subroutine utlow(cone)
1997cc-----------------------------------------------------------------------
1998c character*1 cone
1999c if(cone.eq.'A')cone='a'
2000c if(cone.eq.'B')cone='b'
2001c if(cone.eq.'C')cone='c'
2002c if(cone.eq.'D')cone='d'
2003c if(cone.eq.'E')cone='e'
2004c if(cone.eq.'F')cone='f'
2005c if(cone.eq.'G')cone='g'
2006c if(cone.eq.'H')cone='h'
2007c if(cone.eq.'I')cone='i'
2008c if(cone.eq.'J')cone='j'
2009c if(cone.eq.'K')cone='k'
2010c if(cone.eq.'L')cone='l'
2011c if(cone.eq.'M')cone='m'
2012c if(cone.eq.'N')cone='n'
2013c if(cone.eq.'O')cone='o'
2014c if(cone.eq.'P')cone='p'
2015c if(cone.eq.'Q')cone='q'
2016c if(cone.eq.'R')cone='r'
2017c if(cone.eq.'S')cone='s'
2018c if(cone.eq.'T')cone='t'
2019c if(cone.eq.'U')cone='u'
2020c if(cone.eq.'V')cone='v'
2021c if(cone.eq.'W')cone='w'
2022c if(cone.eq.'X')cone='x'
2023c if(cone.eq.'Y')cone='y'
2024c if(cone.eq.'Z')cone='z'
2025c return
2026c end
60ed90b3 2027c
9ef1c2d9 2028cc-----------------------------------------------------------------------
2029c subroutine utlow3(cthree)
2030cc-----------------------------------------------------------------------
2031c character cthree*3
2032c do 1 i=1,3
2033c1 call utlow(cthree(i:i))
2034c return
2035c end
60ed90b3 2036c
9ef1c2d9 2037cc-----------------------------------------------------------------------
2038c subroutine utlow6(csix)
2039cc-----------------------------------------------------------------------
2040c character csix*6
2041c do 1 i=1,6
2042c1 call utlow(csix(i:i))
2043c return
2044c end
60ed90b3 2045c
9ef1c2d9 2046cc-----------------------------------------------------------------------
2047c function utmom(k,n,x,q)
2048cc-----------------------------------------------------------------------
2049cc calculates kth moment for f(x) with q(i)=int[0,x(i)]f(z)dz
2050cc-----------------------------------------------------------------------
2051c real x(n),q(n)
2052c if(n.lt.2)call utstop('utmom : dimension too small&')
2053c utmom=0
2054c do 1 i=2,n
2055c1 utmom=utmom+((x(i)+x(i-1))/2)**k*(q(i)-q(i-1))
2056c utmom=utmom/q(n)
2057c return
2058c end
60ed90b3 2059c
9ef1c2d9 2060c-----------------------------------------------------------------------
2061 function utpcm(a,b,c)
2062c-----------------------------------------------------------------------
2063c calculates cm momentum for a-->b+c
2064c-----------------------------------------------------------------------
2065 val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2066 if(val.lt.0..and.val.gt.-1e-4)then
2067 utpcm=0
2068 return
2069 endif
2070 utpcm=sqrt(val)/(2.*a)
2071 return
2072 end
2073
2074c-----------------------------------------------------------------------
2075 double precision function utpcmd(a,b,c)
2076c-----------------------------------------------------------------------
2077c calculates cm momentum for a-->b+c
2078c-----------------------------------------------------------------------
2079 double precision a,b,c,val
2080 val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2081 if(val.lt.0d0.and.val.gt.-1d-4)then
2082 utpcmd=0d0
2083 return
2084 elseif(val.lt.0d0)then
2085 call utstop('utpcmd: problem with masses&')
2086 endif
2087 utpcmd=sqrt(val)/(2.d0*a)
2088 return
2089 end
2090
2091c-----------------------------------------------------------------------
2092 subroutine utpri(text,idmmy,ishini,ishx)
2093c-----------------------------------------------------------------------
2094 include 'epos.inc'
2095 character*6 text
2096c double precision seedx !!!
2097 ishini=ish
2098 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2099 if(nrpri.gt.0)then
2100 do nr=1,nrpri
2101 if(subpri(nr)(1:6).eq.text)then
2102 ish=ishpri(nr)
60ed90b3 2103 endif
9ef1c2d9 2104 enddo
2105 endif
2106 if(ish.ge.ishx)then
2107 write(ifch,'(1x,43a)')
2108 * ('-',i=1,10),' entry ',text,' ',('-',i=1,30)
2109c call ranfgt(seedx) !!!
2110c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!!
2111 endif
2112 return
2113 end
60ed90b3 2114
9ef1c2d9 2115c-----------------------------------------------------------------------
2116 subroutine utprix(text,idmmy,ishini,ishx)
2117c-----------------------------------------------------------------------
2118 include 'epos.inc'
2119 character*6 text
2120 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2121 if(ish.ge.ishx)write(ifch,'(1x,44a)')
2122 *('-',i=1,30),' exit ',text,' ',('-',i=1,11)
2123 ish=ishini
2124 return
2125 end
60ed90b3 2126
9ef1c2d9 2127c-----------------------------------------------------------------------
2128 subroutine utprj(text,idmmy,ishini,ishx)
2129c-----------------------------------------------------------------------
2130 include 'epos.inc'
2131 character*20 text
2132c double precision seedx !!!
2133 idx=index(text,' ')-1
2134 ishini=ish
2135 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2136 if(nrpri.gt.0)then
2137 do nr=1,nrpri
2138 if(subpri(nr)(1:idx).eq.text(1:idx))then
2139 ish=ishpri(nr)
60ed90b3 2140 endif
9ef1c2d9 2141 enddo
2142 endif
2143 if(ish.ge.ishx)then
2144 write(ifch,'(1x,43a)')
2145 * ('-',i=1,10),' entry ',text(1:idx),' ',('-',i=1,30)
2146c call ranfgt(seedx) !!!
2147c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!!
2148 endif
2149 return
2150 end
60ed90b3 2151
9ef1c2d9 2152c-----------------------------------------------------------------------
2153 subroutine utprjx(text,idmmy,ishini,ishx)
2154c-----------------------------------------------------------------------
2155 include 'epos.inc'
2156 character*20 text
2157 idx=index(text,' ')-1
2158 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2159 if(ish.ge.ishx)write(ifch,'(1x,44a)')
2160 *('-',i=1,30),' exit ',text(1:idx),' ',('-',i=1,11)
2161 ish=ishini
2162 return
2163 end
2164
2165c-----------------------------------------------------------------------
2166 function utquad(m,x,f,k)
2167c-----------------------------------------------------------------------
2168c performs an integration according to simpson
2169c-----------------------------------------------------------------------
2170 real x(m),f(m)
2171 utquad=0
2172 do 1 i=1,k-1
2173 1 utquad=utquad+(f(i)+f(i+1))/2*(x(i+1)-x(i))
2174 return
2175 end
60ed90b3 2176
9ef1c2d9 2177c-----------------------------------------------------------------------
2178 subroutine utquaf(fu,n,x,q,wo,x0,x1,x2,x3)
2179c-----------------------------------------------------------------------
2180c returns q(i) = integral [x(1)->x(i)] fu(x) dx
2181c-----------------------------------------------------------------------
2182 include 'epos.inc'
2183 real x(n),q(n),wo(n)
2184 parameter (m=10)
2185 real xa(m),fa(m)
2186 external fu
2187 if(x1.lt.x0.or.x2.lt.x1.or.x3.lt.x2)then
2188 if(ish.ge.1)then
2189 call utmsg('utquaf')
2190 write(ifch,*)' xi=',x0,x1,x2,x3
2191 call utmsgf
2192 endif
2193 endif
2194 call utar(n/3,n*2/3,n,x0,x1,x2,x3,x)
2195 q(1)=0
2196 do 2 i=2,n
2197 do 3 k=1,m
2198 z=x(i-1)+(k-1.)/(m-1.)*(x(i)-x(i-1))
2199 xa(k)=z
22003 fa(k)=fu(z)
2201 q(i)=q(i-1)+utquad(m,xa,fa,m)
22022 continue
2203 return
2204 end
2205
2206c-----------------------------------------------------------------------
2207 subroutine utrepl(i,j)
2208c-----------------------------------------------------------------------
2209c i is replaced by j in /cptl/
2210c-----------------------------------------------------------------------
2211 include 'epos.inc'
2212 do 1 k=1,5
22131 pptl(k,i) =pptl(k,j)
2214 iorptl(i) = 0 !iorptl(j)
2215 idptl(i) =idptl(j)
2216 istptl(i) =istptl(j)
2217 do 2 k=1,2
22182 tivptl(k,i)=tivptl(k,j)
2219 do 3 k=1,2
22203 ifrptl(k,i)= 0 !ifrptl(k,j)
2221 jorptl(i) = 0 !jorptl(j)
2222 do 4 k=1,4
22234 xorptl(k,i)=xorptl(k,j)
2224 do 5 k=1,4
22255 ibptl(k,i) =ibptl(k,j)
2226 ityptl(i) =ityptl(j)
2227 iaaptl(i) =iaaptl(j)
2228 radptl(i) =radptl(j)
2229 desptl(i) =desptl(j)
2230 dezptl(i) =dezptl(j)
2231 qsqptl(i) =qsqptl(j)
2232 zpaptl(i) =zpaptl(j)
2233 itsptl(i) =itsptl(j)
2234 rinptl(i) =rinptl(j)
2235 return
2236 end
60ed90b3 2237
9ef1c2d9 2238cc-----------------------------------------------------------------------
2239c subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten)
2240cc-----------------------------------------------------------------------
2241c parameter (nflav=6)
2242c integer icm(2),icp(2),jcm(nflav,2),jcp(nflav,2)
2243c icm(1)=icm1
2244c icm(2)=icm2
2245c icp(1)=icp1
2246c icp(2)=icp2
2247c CALL IDDECO(ICM,JCM)
2248c CALL IDDECO(ICP,JCP)
2249c do 37 nf=1,nflav
2250c do 37 k=1,2
2251c37 jcP(nf,k)=jcp(nf,k)+jcm(nf,k)
2252c CALL IDENCO(JCP,ICP,IRETEN)
2253c IDP=IDTRA(ICP,0,0,3)
2254c call idres(idp,amp,idpr,iadj,0)
2255c return
2256c end
2257c
2258cc-----------------------------------------------------------------------
2259c subroutine utroa1(phi,a1,a2,a3,x1,x2,x3)
2260cc-----------------------------------------------------------------------
2261cc rotates x by angle phi around axis a.
2262cc normalization of a is irrelevant.
2263cc-----------------------------------------------------------------------
2264c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2265c dphi=phi
2266c xx(1)=x1
2267c xx(2)=x2
2268c xx(3)=x3
2269c aa(1)=a1
2270c aa(2)=a2
2271c aa(3)=a3
2272c aaa=0
2273c xxx=0
2274c do i=1,3
2275c aaa=aaa+aa(i)**2
2276c xxx=xxx+xx(i)**2
2277c enddo
2278c if(xxx.eq.0d0)return
2279c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2280c aaa=dsqrt(aaa)
2281c xxx=dsqrt(xxx)
2282cc e3 = a / !a!
2283c do i=1,3
2284c e3(i)=aa(i)/aaa
2285c enddo
2286cc x_parallel
2287c xp=0
2288c do i=1,3
2289c xp=xp+xx(i)*e3(i)
2290c enddo
2291cc x_transverse
2292c if(xxx**2-xp**2.le.0.)return
2293c xt=dsqrt(xxx**2-xp**2)
2294cc e1 = vector x_transverse / absolute value x_transverse
2295c do i=1,3
2296c e1(i)=(xx(i)-e3(i)*xp)/xt
2297c enddo
2298cc e2 orthogonal e3,e1
2299c call utvec2(e3,e1,e2)
2300cc rotate x
2301c do i=1,3
2302c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2303c enddo
2304cc back to single precision
2305c x1=xx(1)
2306c x2=xx(2)
2307c x3=xx(3)
2308c return
2309c end
2310c
2311cc-----------------------------------------------------------------------
2312c subroutine utroa2(phi,a1,a2,a3,x1,x2,x3)
2313cc-----------------------------------------------------------------------
2314cc rotates x by angle phi around axis a.
2315cc normalization of a is irrelevant.
2316cc double precision phi,a1,a2,a3,x1,x2,x3
2317cc-----------------------------------------------------------------------
2318c double precision phi,a1,a2,a3,x1,x2,x3
2319c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2320c dphi=phi
2321c xx(1)=x1
2322c xx(2)=x2
2323c xx(3)=x3
2324c aa(1)=a1
2325c aa(2)=a2
2326c aa(3)=a3
2327c aaa=0
2328c xxx=0
2329c do i=1,3
2330c aaa=aaa+aa(i)**2
2331c xxx=xxx+xx(i)**2
2332c enddo
2333c if(xxx.eq.0d0)return
2334c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2335c aaa=1.0/dsqrt(aaa)
2336cc e3 = a / !a!
2337c do i=1,3
2338c e3(i)=aa(i)*aaa
2339c enddo
2340cc x_parallel
2341c xp=0
2342c do i=1,3
2343c xp=xp+xx(i)*e3(i)
2344c enddo
2345cc x_transverse
2346c if(xxx-xp**2.le.0.)return
2347c xt=dsqrt(xxx-xp**2)
2348cc e1 = vector x_transverse / absolute value x_transverse
2349c do i=1,3
2350c e1(i)=(xx(i)-e3(i)*xp)/xt
2351c enddo
2352cc e2 orthogonal e3,e1
2353c call utvec2(e3,e1,e2)
2354cc rotate x
2355c do i=1,3
2356c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2357c enddo
2358cc back to single precision
2359c x1=xx(1)
2360c x2=xx(2)
2361c x3=xx(3)
2362c return
2363c end
2364c
2365cc--------------------------------------------------------------------
2366c function utroot(funcd,x1,x2,xacc)
2367cc--------------------------------------------------------------------
2368cc combination of newton-raphson and bisection method for root finding
2369cc input:
2370cc funcd: subr returning fctn value and first derivative
2371cc x1,x2: x-interval
2372cc xacc: accuracy
60ed90b3 2373cc output:
9ef1c2d9 2374cc utroot: root
2375cc--------------------------------------------------------------------
2376c include 'epos.inc'
2377c INTEGER MAXIT
2378c REAL utroot,x1,x2,xacc
2379c EXTERNAL funcd
2380c PARAMETER (MAXIT=100)
2381c INTEGER j
2382c REAL df,dx,dxold,f,fh,fl,temp,xh,xl
2383c call funcd(x1,fl,df)
2384c call funcd(x2,fh,df)
2385c if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))
2386c *call utstop('utroot: root must be bracketed&')
2387c if(fl.eq.0.)then
2388c utroot=x1
2389c return
2390c else if(fh.eq.0.)then
2391c utroot=x2
2392c return
2393c else if(fl.lt.0.)then
2394c xl=x1
2395c xh=x2
2396c else
2397c xh=x1
2398c xl=x2
2399c endif
2400c utroot=.5*(x1+x2)
2401c dxold=abs(x2-x1)
2402c dx=dxold
2403c call funcd(utroot,f,df)
2404c do 11 j=1,MAXIT
2405c if(((utroot-xh)*df-f)*((utroot-xl)*df-f).ge.0..or. abs(2.*
2406c *f).gt.abs(dxold*df) ) then
2407c dxold=dx
2408c dx=0.5*(xh-xl)
2409c utroot=xl+dx
2410c if(xl.eq.utroot)return
2411c else
2412c dxold=dx
2413c dx=f/df
2414c temp=utroot
2415c utroot=utroot-dx
2416c if(temp.eq.utroot)return
2417c endif
2418c if(abs(dx).lt.xacc) return
2419c call funcd(utroot,f,df)
2420c if(f.lt.0.) then
2421c xl=utroot
2422c else
2423c xh=utroot
2424c endif
2425c11 continue
2426c call utmsg('utroot')
2427c write(ifch,*)'***** exceeding maximum iterations'
2428c write(ifch,*)'dx:',dx
2429c call utmsgf
2430c return
2431c END
2432c
2433c-----------------------------------------------------------------------
2434 subroutine utrot2(isig,ax,ay,az,x,y,z)
2435c-----------------------------------------------------------------------
2436c performs a rotation, double prec.
2437c-----------------------------------------------------------------------
2438 include 'epos.inc'
2439 double precision ax,ay,az,x,y,z,rx,ry,rz
2440 *,alp,bet,cosa,sina,cosb,sinb,xs,ys,zs
2441 if(ax**2.eq.0.and.ay**2.eq.0.and.az**2.eq.0.)then
2442 write(ifch,*)'ax**2,ay**2,az**2:',ax**2,ay**2,az**2
2443 write(ifch,*)'ax,ay,az:',ax,ay,az
2444 call utstop('utrot2: zero vector.&')
2445 endif
2446 if(az.ge.0.)then
2447 rx=ax
2448 ry=ay
2449 rz=az
2450 else
2451 rx=-ax
2452 ry=-ay
2453 rz=-az
2454 endif
2455 if(rz**2+ry**2.ne.0.)then
2456 alp=dabs(dacos(rz/dsqrt(rz**2+ry**2)))*sign(1.,sngl(ry))
2457 bet=
2458 *dabs(dacos(dsqrt(rz**2+ry**2)/dsqrt(rz**2+ry**2+rx**2)))*
2459 *sign(1.,sngl(rx))
2460 else
2461 alp=3.1415927d0/2d0
2462 bet=3.1415927d0/2d0
2463 endif
2464 cosa=dcos(alp)
2465 sina=dsin(alp)
2466 cosb=dcos(bet)
2467 sinb=dsin(bet)
2468 if(isig.gt.0)then
2469 xs=x*cosb-y*sina*sinb-z*cosa*sinb
2470 ys= y*cosa -z*sina
2471 zs=x*sinb+y*sina*cosb+z*cosa*cosb
2472 elseif(isig.lt.0)then
2473 xs= x*cosb +z*sinb
2474 ys=-x*sinb*sina+y*cosa+z*cosb*sina
2475 zs=-x*sinb*cosa-y*sina+z*cosb*cosa
2476 endif
2477 x=xs
2478 y=ys
2479 z=zs
2480 return
2481 end
60ed90b3 2482
9ef1c2d9 2483
2484c-----------------------------------------------------------------------
2485 subroutine utrot4(isig,ax,ay,az,x,y,z)
2486c-----------------------------------------------------------------------
2487c performs a rotation, double prec.
2488c arguments partly single
2489c-----------------------------------------------------------------------
2490 double precision ax,ay,az,xx,yy,zz
2491 xx=dble(x)
2492 yy=dble(y)
2493 zz=dble(z)
2494 call utrot2(isig,ax,ay,az,xx,yy,zz)
2495 x=sngl(xx)
2496 y=sngl(yy)
2497 z=sngl(zz)
2498 return
2499 end
2500
2501c-----------------------------------------------------------------------
2502 subroutine utrota(isig,ax,ay,az,x,y,z)
2503c-----------------------------------------------------------------------
2504c performs a rotation
60ed90b3 2505c-----------------------------------------------------------------------
9ef1c2d9 2506 if(az.ge.0.)then
2507 rx=ax
2508 ry=ay
2509 rz=az
2510 else
2511 rx=-ax
2512 ry=-ay
2513 rz=-az
2514 endif
2515 if(rz.eq.0..and.ry.eq.0.)then
2516 alp=0.
2517 stop
2518 else
2519 alp=abs(utacos(rz/sqrt(rz**2+ry**2)))*sign(1.,ry)
2520 endif
2521 bet=
2522 *abs(utacos(sqrt(rz**2+ry**2)/sqrt(rz**2+ry**2+rx**2)))*sign(1.,rx)
2523 cosa=cos(alp)
2524 sina=sin(alp)
2525 cosb=cos(bet)
2526 sinb=sin(bet)
2527 if(isig.gt.0)then
2528 xs=x*cosb-y*sina*sinb-z*cosa*sinb
2529 ys= y*cosa -z*sina
2530 zs=x*sinb+y*sina*cosb+z*cosa*cosb
2531 elseif(isig.lt.0)then
2532 xs= x*cosb +z*sinb
2533 ys=-x*sinb*sina+y*cosa+z*cosb*sina
2534 zs=-x*sinb*cosa-y*sina+z*cosb*cosa
2535 endif
2536 x=xs
2537 y=ys
2538 z=zs
2539 return
2540 end
60ed90b3 2541
9ef1c2d9 2542c-----------------------------------------------------------------------
2543 subroutine utstop(text)
2544c-----------------------------------------------------------------------
2545c returns error message and stops execution.
2546c text is an optonal text to appear in the error message.
60ed90b3 2547c text is a character string of length 40;
9ef1c2d9 2548c for shorter text, it has to be terminated by &;
2549c example: call utstop('error in subr xyz&')
2550c-----------------------------------------------------------------------
2551 include 'epos.inc'
2552 parameter(itext=40)
2553 character text*40 ,txt*6
2554 imax=itext+1
2555 do i=itext,1,-1
2556 if(text(i:i).eq.'&')imax=i
2557 enddo
2558 do 1 j=1,2
2559 if(j.eq.1)ifi=ifch
2560 if(j.eq.2)ifi=ifmt
2561 write(ifi,101)('*',k=1,72),text(1:imax-1)
2562 *,nrevt+1,seedc,('*',k=1,72)
2563101 format(
2564 *1x,72a1
2565 */1x,'***** stop in ',a
2566 */1x,'***** current event number: ',i12
2567 */1x,'***** initial seed for current event:',d25.15
2568 */1x,72a1)
25691 continue
2570 stop
2571 entry utmsg(txt)
2572 imsg=imsg+1
2573 write(ifch,'(1x,74a1)')('*',j=1,72)
2574 write(ifch,100)txt,nrevt+1,seedc
2575100 format(1x,'***** msg from ',a6,'. es:',i7,2x,d23.17)
2576 return
2577 entry utmsgf
2578 if(ish.eq.1)return
2579 write(ifch,'(1x,74a1)')('*',j=1,72)
2580 end
60ed90b3 2581
9ef1c2d9 2582c-----------------------------------------------------------------
2583 subroutine uttrap(func,a,b,s)
2584c-----------------------------------------------------------------
2585c trapezoidal method for integration.
2586c input: fctn func and limits a,b
2587c output: value s of the integral
2588c-----------------------------------------------------------------
2589 include 'epos.inc'
2590
2591 INTEGER JMAX
2592 REAL a,b,func,s,epsr
2593 EXTERNAL func
2594 PARAMETER (JMAX=10)
2595CU USES uttras
2596 INTEGER j
2597 REAL olds
2598 olds=-1.e30
2599 do 11 j=1,JMAX
2600 if(ish.ge.9)write(ifch,*)'sr uttrap: j:',j
2601 call uttras(func,a,b,s,j)
2602 ds=abs(s-olds)
2603 if (ds.lt.epsr*abs(olds)) return
2604 olds=s
260511 continue
60ed90b3 2606
9ef1c2d9 2607c-c nepsr=nepsr+1
2608 if(ish.ge.9)then
2609 call utmsg('uttrap')
2610 write(ifch,*)
2611 *'***** requested accuracy could not be achieved'
2612 write(ifch,*)'achieved accuracy: ',ds/abs(olds)
2613 write(ifch,*)'requested accuracy:',epsr
2614 call utmsgf
2615 endif
2616
2617 END
2618
2619c-----------------------------------------------------------------
2620 subroutine uttraq(func,a,b,s)
2621c-----------------------------------------------------------------
2622c trapezoidal method for integration.
2623c input: function func and limits a,b
2624c output: value s of the integral
2625c-----------------------------------------------------------------
2626
2627 REAL a,b,func,s
2628 EXTERNAL func
2629 PARAMETER (eps=1.e-6)
2630CU USES uttras
2631 INTEGER j
2632 REAL olds
2633 olds=-1.e30
2634 j=1
263510 call uttras(func,a,b,s,j)
2636 ds=abs(s-olds)
2637 if (ds.le.eps*abs(olds)) return
2638 olds=s
2639 if(j.ge.15)return
2640 j=j+1
2641 goto10
2642 END
2643
2644c-----------------------------------------------------------------
2645 subroutine uttras(func,a,b,s,n)
2646c-----------------------------------------------------------------
2647c performs one iteration of the trapezoidal method for integration
2648c-----------------------------------------------------------------
2649 INTEGER n
2650 REAL a,b,s,func
2651 EXTERNAL func
2652 INTEGER it,j
2653 REAL del,sum,tnm,x
2654 if (n.eq.1) then
2655 s=0.5*(b-a)*(func(a)+func(b))
2656 else
2657 it=2**(n-2)
2658 tnm=it
2659 del=(b-a)/tnm
2660 x=a+0.5*del
2661 sum=0.
2662 do 11 j=1,it
2663 sum=sum+func(x)
2664 x=x+del
266511 continue
2666 s=0.5*(s+(b-a)*sum/tnm)
2667 endif
2668 return
2669 END
2670
2671cc-----------------------------------------------------------------------
2672c subroutine utvec1(a,b,c)
2673cc-----------------------------------------------------------------------
2674cc returns vector product c = a x b .
2675cc-----------------------------------------------------------------------
2676c real a(3),b(3),c(3)
2677c c(1)=a(2)*b(3)-a(3)*b(2)
2678c c(2)=a(3)*b(1)-a(1)*b(3)
2679c c(3)=a(1)*b(2)-a(2)*b(1)
2680c return
2681c end
2682c
2683cc-----------------------------------------------------------------------
2684c subroutine utvec2(a,b,c)
2685cc-----------------------------------------------------------------------
2686cc returns vector product c = a x b .
2687cc a,b,c double precision.
2688cc-----------------------------------------------------------------------
2689c double precision a(3),b(3),c(3)
2690c c(1)=a(2)*b(3)-a(3)*b(2)
2691c c(2)=a(3)*b(1)-a(1)*b(3)
2692c c(3)=a(1)*b(2)-a(2)*b(1)
2693c return
2694c end
2695c
2696c-------------------------------------------------------------------
2697 subroutine utword(line,i,j,iqu)
2698c-------------------------------------------------------------------
2699c finds the first word of the character string line(j+1:160).
2700c the word is line(i:j) (with new i and j).
2701c if j<0 or if no word found --> new line read.
60ed90b3 2702c a text between quotes "..." is considered one word;
9ef1c2d9 2703c stronger: a text between double quotes ""..."" is consid one word
2704c stronger: a text between "{ and }" is considered one word
2705c-------------------------------------------------------------------
60ed90b3 2706c input:
9ef1c2d9 2707c line: character string (*160)
2708c i: integer between 1 and 160
2709c iqu: for iqu=1 a "?" is written to output before reading a line,
2710c otherwise (iqu/=1) nothing is typed
60ed90b3 2711c output:
9ef1c2d9 2712c i,j: left and right end of word (word=line(i:j))
2713c-------------------------------------------------------------------
2714 include 'epos.inc'
2715 parameter(mempty=2)
2716 character*1 empty(mempty),mk
2717 character line*160
2718 character*2 mrk
2719 data empty/' ',','/
2720 parameter(mxdefine=40)
2721 character w1define*100,w2define*100
2722 common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine)
2723 & ,w1define(mxdefine),w2define(mxdefine)
60ed90b3 2724
9ef1c2d9 2725 if(j.ge.0)then
2726 i=j
2727 goto1
2728 endif
2729
2730 5 continue
2731 if(iqu.eq.1.and.iprmpt.gt.0)write(ifmt,'(a)')'?'
2732 if(nopen.eq.0)ifopx=ifop
2733 if(nopen.gt.0)ifopx=20+nopen
2734 if(nopen.lt.0)ifopx=ifcp
2735 read(ifopx,'(a160)',end=9999)line
2736 if(iecho.eq.1.or.(nopen.ge.0.and.kcpopen.eq.1))then
2737 kmax=2
2738 do k=3,160
2739 if(line(k:k).ne.' ')kmax=k
2740 enddo
2741 endif
2742 if(nopen.ge.0.and.kcpopen.eq.1)
2743 & write(ifcp,'(a)')line(1:kmax)
2744 if(iecho.eq.1)
2745 & write(ifmt,'(a)')line(1:kmax)
2746 i=0
2747
2748 1 i=i+1
6801ca52 2749 if(i.gt.160)goto5
2750 if(line(i:i).eq.'!')goto5
9ef1c2d9 2751 do ne=1,mempty
2752 if(line(i:i).eq.empty(ne))goto1
2753 enddo
2754
2755 mrk=' '
2756 mk=' '
2757 if(line(i:i).eq.'~')mk='~'
2758 if(line(i:i+1).eq.'"{')mrk='}"'
2759 if(line(i:i+1).eq.'""')mrk='""'
2760 if(mrk.ne.' ')goto10
2761 if(line(i:i).eq.'"')mk='"'
2762 if(mk.ne.' ')goto8
2763 j=i-1
2764 6 j=j+1
2765 if(j.gt.160)goto7
2766 if(line(j:j).eq.'!')goto7
2767 do ne=1,mempty
2768 if(line(j:j).eq.empty(ne))goto7
2769 enddo
2770 goto6
2771
2772 8 continue
2773 if(i.ge.160-1)stop'utword: make line shorter!!! '
2774 i=i+1
2775 j=i
2776 if(line(j:j).eq.mk)stop'utword: empty string!!! '
2777 9 j=j+1
2778 if(j.gt.160)goto7
2779 if(line(j:j).eq.mk)then
2780 line(i-1:i-1)=' '
2781 line(j:j)=' '
2782 goto7
2783 endif
2784 goto9
2785
2786 10 continue
2787 if(i.ge.160-3)stop'utword: make line shorter!!!! '
2788 i=i+2
2789 j=i
2790 if(line(j:j+1).eq.mrk)stop'utword: empty string!!!! '
2791 11 j=j+1
2792 if(j+1.gt.160)goto7
2793 if(line(j:j+1).eq.mrk)then
2794 line(i-2:i-1)=' '
2795 line(j:j+1)=' '
2796 goto7
2797 endif
2798 goto11
2799
2800 7 j=j-1
2801 !--------#define---------------
2802 if(ndefine.gt.0)then
2803 do ndf=1,ndefine
2804 l1=l1define(ndf)
2805 l2=l2define(ndf)
2806 do i0=i,j+1-l1
2807 if(line(i0:i0-1+l1).eq.w1define(ndf)(1:l1))then
2808 if(l2.eq.l1)then
2809 line(i0:i0-1+l1)=w2define(ndf)(1:l2)
60ed90b3 2810 elseif(l2.lt.l1)then
9ef1c2d9 2811 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
2812 do k=i0+l2,i0-1+l1
2813 line(k:k)=' '
2814 enddo
60ed90b3 2815 elseif(l2.gt.l1)then
9ef1c2d9 2816 do k=i0+l1,i0+l2-1
2817 if(line(k:k).ne.' ')
2818 & stop'utword: no space for `define` replacement. '
2819 enddo
2820 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
2821 j=i0+l2-1
60ed90b3 2822 endif
9ef1c2d9 2823 endif
2824 enddo
2825 enddo
2826 do k=i,j
2827 if(line(k:k).ne.' ')j0=j
2828 enddo
2829 j=j0
2830 endif
2831 !--------
2832 return
2833
28349999 close(ifopx)
2835 nopen=nopen-1
2836 if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
2837 goto5
2838 end
2839
2840c--------------------------------------------------------------------
2841 subroutine utworn(line,j,ne)
2842c--------------------------------------------------------------------
2843c returns number ne of nonempty characters of line(j+1:160)
2844c--------------------------------------------------------------------
2845 character line*160
2846 ne=0
2847 do l=j+1,160
2848 if(line(l:l).ne.' ')ne=ne+1
2849 enddo
2850 return
2851 end
2852
60ed90b3 2853c-----------------------------------------------------------------------
9ef1c2d9 2854 subroutine getairmol(iz,ia)
60ed90b3 2855c-----------------------------------------------------------------------
9ef1c2d9 2856 include 'epos.inc'
60ed90b3 2857 i=0
2858 r=rangen()
9ef1c2d9 2859 do while(r.gt.0.) ! choose air-molecule
60ed90b3 2860 i=i+1
2861 r=r-airwnxs(i)
2862 enddo
2863 iz = nint(airznxs(i))
2864 ia = nint(airanxs(i))
9ef1c2d9 2865 end
2866
2867c----------------------------------------------------------------------
2868
2869 subroutine factoriel
60ed90b3 2870
9ef1c2d9 2871c----------------------------------------------------------------------
2872c tabulation of fctrl(n)=n!, facto(n)=1/n! and utgamtab(x) for x=0 to 50
2873c----------------------------------------------------------------------
2874 include 'epos.incems'
2875 double precision utgamtab,utgam,x
2876 common/gamtab/utgamtab(10000)
2877
2878 nfctrl=100
2879 fctrl(0)=1.D0
2880 facto(0)=1.D0
2881 do i=1,min(npommx,nfctrl)
2882 fctrl(i)=fctrl(i-1)*dble(i)
2883 facto(i)=1.d0/fctrl(i)
2884 enddo
60ed90b3 2885
9ef1c2d9 2886 do k=1,10000
2887 x=dble(k)/100.d0
2888 utgamtab(k)=utgam(x)
2889 enddo
60ed90b3 2890
9ef1c2d9 2891 return
2892 end
60ed90b3 2893
9ef1c2d9 2894c-----------------------------------------------------------------------
2895 subroutine fremnu(amin,ca,cb,ca0,cb0,ic1,ic2,ic3,ic4)
2896c-----------------------------------------------------------------------
2897 common/hadr2/iomodl,idproj,idtarg,wexcit
2898 real pnll,ptq
2899 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2900 parameter (nflav=6)
2901 integer ic(2),jc(nflav,2)
2902 ic(1)=ca
2903 ic(2)=cb
2904 call iddeco(ic,jc)
2905 keu=jc(1,1)-jc(1,2)
2906 ked=jc(2,1)-jc(2,2)
2907 kes=jc(3,1)-jc(3,2)
2908 kec=jc(4,1)-jc(4,2)
2909 keb=jc(5,1)-jc(5,2)
2910 ket=jc(6,1)-jc(6,2)
2911 amin=utamnu(keu,ked,kes,kec,keb,ket,5) !???4=2mults, 5=1mult
2912 if(ca-ca0.eq.0..and.cb-cb0.eq.0..and.rangen().gt.wexcit)then
2913 ic3=0
2914 ic4=0
2915 ic1=ca
2916 ic2=cb
2917 else
2918 amin=amin+exmass
2919 n=0
2920 do i=1,4
2921 do j=1,2
2922 n=n+jc(i,j)
2923 enddo
2924 enddo
2925 k=1+rangen()*n
2926 do i=1,4
2927 do j=1,2
2928 k=k-jc(i,j)
2929 if(k.le.0)goto 1
2930 enddo
2931 enddo
29321 if(j.eq.1)then
2933 ic3=10**(6-i)
2934 ic4=0
2935 else
2936 ic3=0
2937 ic4=10**(6-i)
2938 endif
2939 ic1=int(ca)-ic3
2940 ic2=int(cb)-ic4
2941 endif
2942 return
2943 end
2944
2945
2946c-----------------------------------------------------------------------
2947 function fremnux(ic)
2948c-----------------------------------------------------------------------
2949 real pnll,ptq
2950 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2951 parameter (nflav=6)
2952 integer ic(2),jc(nflav,2)
2953c ic(1)=210000
2954c ic(2)=0
2955 call iddeco(ic,jc)
2956 keu=jc(1,1)-jc(1,2)
2957 ked=jc(2,1)-jc(2,2)
2958 kes=jc(3,1)-jc(3,2)
2959 kec=jc(4,1)-jc(4,2)
2960 keb=jc(5,1)-jc(5,2)
2961 ket=jc(6,1)-jc(6,2)
2962 fremnux=utamnu(keu,ked,kes,kec,keb,ket,4) !+exmass !???4=2mults, 5=1mult
2963 return
2964 end
2965
2966c-----------------------------------------------------------------------
2967 function fremnux2(ic)
2968c-----------------------------------------------------------------------
2969 real pnll,ptq
2970 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2971 parameter (nflav=6)
2972 integer ic(2),jc(nflav,2)
2973c ic(1)=210000
2974c ic(2)=0
2975 call iddeco(ic,jc)
2976 keu=jc(1,1)-jc(1,2)
2977 ked=jc(2,1)-jc(2,2)
2978 kes=jc(3,1)-jc(3,2)
2979 kec=jc(4,1)-jc(4,2)
2980 keb=jc(5,1)-jc(5,2)
2981 ket=jc(6,1)-jc(6,2)
2982 fremnux2=utamnu(keu,ked,kes,kec,keb,ket,5) !+exmass !???4=2mults, 5=1mult
2983 return
2984 end
2985
2986c-----------------------------------------------------------------------
2987 subroutine fremnx(ammax,amin,sm,ic3,ic4,iret)
2988c-----------------------------------------------------------------------
2989 common/psar9/ alpr
2990 include 'epos.inc'
2991 iret=0
2992 if(ic3.eq.0.and.ic4.eq.0)then
2993 if(ammax.lt.amin**2)then
2994 iret=1
2995 return
2996 endif
2997 sm=amin**2
2998 else
2999c ammax1=min(ammax,(engy/4.)**2)
3000 ammax1=ammax
3001 if(ammax1.lt.amin**2)then
3002 iret=1
3003 return
3004 endif
3005 if(alpr.eq.-1.)then
3006 sm=amin**2*(ammax1/amin**2)**rangen()
3007 else
3008 sm=amin**2*(1.+((ammax1/amin**2)**(1.+alpr)-1.)
3009 * *rangen())**(1./(1.+alpr))
3010 endif
3011 endif
3012 return
3013 end
3014
3015