]>
Commit | Line | Data |
---|---|---|
1 | c----------------------------------------------------------------------- | |
2 | subroutine utresc(iret) | |
3 | c----------------------------------------------------------------------- | |
4 | c if irescl=1 rescaling is done, otherwise the purpose of going through | |
5 | c this routine is to not change the seed in case of no rescaling | |
6 | c----------------------------------------------------------------------- | |
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) | |
17 | ||
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 | |
23 | ||
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 | |
40 | c 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 | |
44 | c 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. | |
57 | c write(ifch,*)'nolead',i | |
58 | else | |
59 | nolead(i)=.false. | |
60 | c 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' | |
72 | c call utstop('utresc&') | |
73 | goto 9999 | |
74 | endif | |
75 | esoll=p1(5) | |
76 | if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector: ',p1 | |
77 | ||
78 | c trafo | |
79 | c ----- | |
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&') | |
91 | endif | |
92 | ||
93 | if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------' | |
94 | ||
95 | c rescale momenta in rest frame | |
96 | c ----------------------------- | |
97 | scal=1. | |
98 | diff0=0. | |
99 | c 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 | |
113 | c modified particles | |
114 | if(nolead(j))then | |
115 | c if(j.gt.nptlpt)then | |
116 | c 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. | |
122 | if( force .or.( | |
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)), | |
127 | & rangen()*abs(difft)),difft) | |
128 | pptl3new=scal*(pptl(3,j)-diff) | |
129 | c write(ifch,*)'par',j,pptl3new,pptl(3,j),diff,difft | |
130 | c & ,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 | |
134 | c 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 | |
144 | c 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 | |
155 | if(abs(scal-1.).le.errlim.and.abs(diff/pnullx).lt.errlim) | |
156 | $ goto 300 | |
157 | if(ndif.gt.0.and.(force.or.ipass.lt.150))then | |
158 | c ndif0=ndif | |
159 | diff0=diff | |
160 | elseif(abs(scal-1.).le.1e-2.and.abs(diff/pnullx).lt.5e-2)then | |
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 | |
183 | ||
184 | c trafo | |
185 | c ----- | |
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 | |
199 | ||
200 | if(ish.ge.4)call alist('list after rescaling&',1,nptl) | |
201 | ||
202 | 9999 continue | |
203 | if(ish.ge.2)then | |
204 | scalmean=scalmean+scal | |
205 | write(ifch,*)' average rescaling factor: ',scalmean | |
206 | & /float(nrevt+1) | |
207 | endif | |
208 | call ranfst(seedp) ! ... after this subroutine | |
209 | call utprix('utresc',ish,ishini,4) | |
210 | ||
211 | end | |
212 | ||
213 | c----------------------------------------------------------------------- | |
214 | subroutine utghost(iret) | |
215 | c----------------------------------------------------------------------- | |
216 | c if irescl=1 make particle on-shell if not | |
217 | c----------------------------------------------------------------------- | |
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 | |
227 | ||
228 | if(ish.ge.5)write(ifch,'(a)')'---------mark ghosts---------' | |
229 | ||
230 | c mark ghosts | |
231 | c ----------- | |
232 | do j=nptlpt+1,nptl | |
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) | |
256 | if(ish.ge.5)then | |
257 | write(ifch,'(a,$)')'ghost:' | |
258 | call alistc("&",j,j) | |
259 | endif | |
260 | ityptl(j)=100+ityptl(j)/10 | |
261 | endif | |
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---------' | |
273 | ||
274 | c treat ghosts | |
275 | c ------------ | |
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 | |
285 | do j=nptlpt+1,nptl | |
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 | |
289 | if(ifirst.eq.1)then | |
290 | pfif=pfif+pptl(3,j) | |
291 | if(pptl(4,j).gt.0.)efif=efif+pptl(4,j) | |
292 | endif | |
293 | if(irescl.ge.1) then | |
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 | |
317 | c 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 | |
330 | endif | |
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 | ||
371 | c Check Ghost list | |
372 | ||
373 | if(ish.ge.5)then | |
374 | do j=nptlpt+1,nptl | |
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 | ||
384 | ||
385 | 9999 continue | |
386 | call ranfst(seedp) ! ... after this subroutine | |
387 | call utprix('ughost',ish,ishini,4) | |
388 | ||
389 | end | |
390 | ||
391 | c----------------------------------------------------------------------- | |
392 | subroutine utrsph(iret) | |
393 | c----------------------------------------------------------------------- | |
394 | c if irescl=1 and ispherio=1 rescaling is done for particle used by | |
395 | c spherio as initial condition. | |
396 | c----------------------------------------------------------------------- | |
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 | |
404 | ||
405 | iret=0 | |
406 | nptlpt=iabs(maproj)+iabs(matarg) | |
407 | if(nptl.le.nptlpt) goto 9999 | |
408 | ||
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)) | |
417 | $ .or.istptl(i).eq.20.or.istptl(i).eq.21)then | |
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 | ||
435 | c trafo | |
436 | c ----- | |
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 | |
441 | $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then | |
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 | ||
450 | c rescale momenta in rest frame | |
451 | c ----------------------------- | |
452 | ||
453 | scal=1. | |
454 | diff=0. | |
455 | do ipass=1,1000 | |
456 | sum=0. | |
457 | sum3=0. | |
458 | ndif=0 | |
459 | do j=1,nptl | |
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 | |
463 | $ .or.(istptl(j).eq.0.and.j.le.nptlpt))then | |
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 | ||
487 | ||
488 | c trafo | |
489 | c ----- | |
490 | 300 continue | |
491 | c 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 | |
496 | $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then | |
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 | |
506 | ||
507 | 9999 call utprix('utrsph',ish,ishini,4) | |
508 | ||
509 | end | |
510 | ||
511 | cc----------------------------------------------------------------------- | |
512 | c double precision function dddlog(xxx) | |
513 | cc----------------------------------------------------------------------- | |
514 | c double precision xxx | |
515 | c dddlog=-1d50 | |
516 | c if(xxx.gt.0d0)dddlog=dlog(xxx) | |
517 | c end | |
518 | c | |
519 | ccc----------------------------------------------------------------------- | |
520 | c subroutine randfl(jc,iqa0,iflav,ic,isame) | |
521 | cc----------------------------------------------------------------------- | |
522 | cc returns random flavour ic(2) (iqa0=1:quark,2:antiquark,11:diquark) | |
523 | cc----------------------------------------------------------------------- | |
524 | c include 'epos.inc' | |
525 | c real probab(nflav),probsu(nflav+1) | |
526 | c integer jc(nflav,2),jc0(nflav,2),ic(2) | |
527 | c if(ish.ge.6)then | |
528 | c write(ifch,*)('-',i=1,10) | |
529 | c *,' entry sr randfl ',('-',i=1,30) | |
530 | c write(ifch,*)'iqa0:',iqa0 | |
531 | c write(ifch,*)'jc:' | |
532 | c write(ifch,*)jc | |
533 | c endif | |
534 | c iflav=0 | |
535 | c ic(1)=0 | |
536 | c ic(2)=0 | |
537 | c do 10 n=1,nflav | |
538 | c do 10 i=1,2 | |
539 | c10 jc0(n,i)=0 | |
540 | c iqa1=iqa0*10 | |
541 | c9999 iqa1=iqa1/10 | |
542 | c if(iqa1.eq.0)goto9998 | |
543 | c iqa=mod(iqa1,10) | |
544 | c su=0 | |
545 | c do 20 i=1,nflav | |
546 | c probab(i)=jc(i,iqa)-jc0(i,iqa) | |
547 | c if(isame.eq.1)probab(i)=probab(i)*(jc(i,3-iqa)-jc0(i,3-iqa)) | |
548 | c20 su=su+probab(i) | |
549 | c if(su.lt..5)then | |
550 | c iflav=0 | |
551 | c ic(1)=0 | |
552 | c ic(2)=0 | |
553 | c goto9998 | |
554 | c endif | |
555 | c probsu(1)=0. | |
556 | c do 30 i=1,nflav | |
557 | c probsu(i+1)=probsu(i)+probab(i)/su | |
558 | c if(probsu(i+1)-probsu(i).lt.1e-5)probsu(i+1)=probsu(i) | |
559 | c30 continue | |
560 | c r=rangen()*probsu(nflav+1) | |
561 | c do 50 i=1,nflav | |
562 | c if(probsu(i).le.r.and.r.lt.probsu(i+1))iflav=i | |
563 | c50 continue | |
564 | c jc0(iflav,iqa)=jc0(iflav,iqa)+1 | |
565 | c if(isame.eq.1)jc0(iflav,3-iqa)=jc0(iflav,3-iqa)+1 | |
566 | c call idenco(jc0,ic,ireten) | |
567 | c if(ireten.eq.1)call utstop('randfl: idenco ret code = 1&') | |
568 | c if(ish.ge.6)then | |
569 | c write(ifch,*)'probab:' | |
570 | c write(ifch,*)probab | |
571 | c write(ifch,*)'probsu:' | |
572 | c write(ifch,*)probsu | |
573 | c write(ifch,*)'ran#:',r,' flav:',iflav | |
574 | c endif | |
575 | c goto9999 | |
576 | c9998 continue | |
577 | c if(ish.ge.6)write(ifch,*)('-',i=1,30) | |
578 | c *,' exit sr randfl ',('-',i=1,10) | |
579 | c return | |
580 | c end | |
581 | c | |
582 | c | |
583 | cc----------------------------------------------------------------------- | |
584 | c subroutine ranhvy(x,eps) | |
585 | cc----------------------------------------------------------------------- | |
586 | cc generates x for heavy particle fragmentation according to | |
587 | cc the peterson form | |
588 | cc d(x)=1/(x*(1-1/x-eps/(1-x))**2) | |
589 | cc =d0(x)*d1(x)*d2(x) | |
590 | cc d0(x)=(1-x)**2/((1-x)**2+eps)**2 | |
591 | cc d1(x)=x | |
592 | cc d2(x)=(((1-x)**2+eps)/((1-x)**2+eps*x))**2 | |
593 | cc using x=1-y**pow | |
594 | cc generates flat in x if eps>1. | |
595 | cc----------------------------------------------------------------------- | |
596 | c data aln4/1.3863/ | |
597 | c if(eps.lt.1.) then | |
598 | c pow=alog((3.+eps)/eps)/aln4 | |
599 | c ymx=(eps*(3.*pow-1.)/(pow+1.))**(.5/pow) | |
600 | c zmx=1-ymx**pow | |
601 | c d0mx=(1-zmx)**2/((1.-zmx)**2+eps)**2*pow*ymx**(pow-1.) | |
602 | c d2mx=2./(2.-sqrt(eps)) | |
603 | c else | |
604 | c pow=1. | |
605 | c zmx=0. | |
606 | c d0mx=(1.-zmx)**2/((1.-zmx)**2+eps)**2 | |
607 | c d2mx=1.+eps | |
608 | c endif | |
609 | cc | |
610 | cc generate z according to (1-z)**2/((1-z)**2+eps*z)**2 | |
611 | c1 continue | |
612 | c y=rangen() | |
613 | c z=1.-y**pow | |
614 | cc | |
615 | c d0z=(1.-z)**2/((1.-z)**2+eps)**2*pow*y**(pow-1.) | |
616 | c if(d0z.lt.rangen()*d0mx) goto1 | |
617 | cc | |
618 | cc check remaining factors | |
619 | c d1=z | |
620 | c d2=(((1.-z)**2+eps)/((1.-z)**2+eps*z))**2 | |
621 | c if(d1*d2.lt.rangen()*d2mx) goto1 | |
622 | cc | |
623 | cc good x | |
624 | c x=z | |
625 | c return | |
626 | c end | |
627 | c | |
628 | c----------------------------------------------------------------------- | |
629 | function ransig() | |
630 | c----------------------------------------------------------------------- | |
631 | c returns randomly +1 or -1 | |
632 | c----------------------------------------------------------------------- | |
633 | ransig=1 | |
634 | if(rangen().gt.0.5)ransig=-1 | |
635 | return | |
636 | end | |
637 | ||
638 | cc----------------------------------------------------------------------- | |
639 | c function ranxq(n,x,q,xmin) | |
640 | cc----------------------------------------------------------------------- | |
641 | cc returns random number according to x(i) q(i) with x>=xmin | |
642 | cc----------------------------------------------------------------------- | |
643 | c include 'epos.inc' | |
644 | c real x(n),q(n) | |
645 | c imin=1 | |
646 | c if(xmin.eq.0.)goto3 | |
647 | c i1=1 | |
648 | c i2=n | |
649 | c1 i=i1+(i2-i1)/2 | |
650 | c if(x(i).lt.xmin)then | |
651 | c i1=i | |
652 | c elseif(x(i).gt.xmin)then | |
653 | c i2=i | |
654 | c else | |
655 | c imin=i | |
656 | c goto3 | |
657 | c endif | |
658 | c if(i2-i1.gt.1)goto1 | |
659 | c imin=i2 | |
660 | c3 continue | |
661 | c if(q(imin).gt.q(n)*.9999)then | |
662 | c ranxq=xmin | |
663 | c goto4 | |
664 | c endif | |
665 | c qran=q(imin)+rangen()*(q(n)-q(imin)) | |
666 | c ranxq=utinvt(n,x,q,qran) | |
667 | c4 continue | |
668 | c | |
669 | c if(ranxq.lt.xmin)then | |
670 | c call utmsg('ranxq ') | |
671 | c write(ifch,*)'***** ranxq=',ranxq,' < xmin=',xmin | |
672 | c write(ifch,*)'q(imin) q q(n):',q(imin),qran,q(n) | |
673 | c write(ifch,*)'x(imin) x x(n):',x(imin),ranxq,x(n) | |
674 | c call utmsgf | |
675 | c ranxq=xmin | |
676 | c endif | |
677 | c | |
678 | c return | |
679 | c end | |
680 | c | |
681 | cc ***** end r-routines | |
682 | cc ***** beg s-routines | |
683 | c | |
684 | cc----------------------------------------------------------------------- | |
685 | c function sbet(z,w) | |
686 | cc----------------------------------------------------------------------- | |
687 | c sbet=utgam1(z)*utgam1(w)/utgam1(z+w) | |
688 | c return | |
689 | c end | |
690 | c | |
691 | cc----------------------------------------------------------------------- | |
692 | c function smass(a,y,z) | |
693 | cc----------------------------------------------------------------------- | |
694 | cc returns droplet mass (in gev) (per droplet, not (!) per nucleon) | |
695 | cc according to berger/jaffe mass formula, prc35(1987)213 eq.2.31, | |
696 | cc see also c. dover, BNL-46322, intersections-meeting, tucson, 91. | |
697 | cc a: massnr, y: hypercharge, z: charge, | |
698 | cc----------------------------------------------------------------------- | |
699 | c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero | |
700 | c ymin=ym*a | |
701 | c zmin=cz/(dz/a+zm/a**(1./3.)) | |
702 | c smass=epsi*a+as*a**(2./3.)+(ac/a**(1./3.)+dz/a/2.)*(z-zmin)**2 | |
703 | c *+dy/a/2.*(y-ymin)**2 | |
704 | c return | |
705 | c end | |
706 | c | |
707 | cc----------------------------------------------------------------------- | |
708 | c subroutine smassi(theta) | |
709 | cc----------------------------------------------------------------------- | |
710 | cc initialization for smass. | |
711 | cc calculates parameters for berger/jaffe mass formula | |
712 | cc (prc35(1987)213 eq.2.31, see also c. dover, BNL-46322). | |
713 | cc theta: parameter that determines all parameters in mass formula. | |
714 | cc----------------------------------------------------------------------- | |
715 | c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero | |
716 | c thet=theta | |
717 | c | |
718 | c astr=.150 | |
719 | c pi=3.14159 | |
720 | c alp=1./137. | |
721 | c | |
722 | c co=cos(theta) | |
723 | c si=sin(theta) | |
724 | c bet=(1+co**3)/2. | |
725 | c rzero=si/astr/( 2./3./pi*(1+co**3) )**(1./3.) | |
726 | cctp060829 cs=astr/si | |
727 | c cz=-astr/si*( ( .5*(1+co**3) )**(1./3.)-1 ) | |
728 | c sigma=6./8./pi*(astr/si)**3*(co**2/6.-si**2*(1-si)/3.- | |
729 | c *1./3./pi*(pi/2.-theta-sin(2*theta)+si**3*alog((1+co)/si))) | |
730 | c | |
731 | c epsi=astr*((.5*(1+co**3))**(1./3.)+2)/si | |
732 | c as=4*pi*sigma*rzero**2 | |
733 | c ac=3./5.*alp/rzero | |
734 | c dz=astr/si*bet**(1./3.)*co**2* | |
735 | c *(co**4*(1+bet**(2./3.))+(1+bet)**2)/ | |
736 | c *( (2*co**2+bet**(1./3.))*(co**4*(1+bet**(2./3.))+(1+bet)**2)- | |
737 | c *(co**4+bet**(1./3.)*(1+bet))*((2*bet**(2./3.)-1)*co**2+1+bet) ) | |
738 | c dy=astr/6.*(1+co**3)**3/si* | |
739 | c *( 1+(1+co)/(4*(1+co**3))**(2./3.) )/ | |
740 | c *(co**6+co+co*(.5*(1+co**3))**(4./3.)) | |
741 | c zm=6*alp/(5*rzero) | |
742 | c ym=(1-co**3)/(1+co**3) | |
743 | c | |
744 | c return | |
745 | c end | |
746 | c | |
747 | cc----------------------------------------------------------------------- | |
748 | c subroutine smassp | |
749 | cc----------------------------------------------------------------------- | |
750 | cc prints smass. | |
751 | cc----------------------------------------------------------------------- | |
752 | c include 'epos.inc' | |
753 | c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero | |
754 | c real eng(14),ymi(14),zmi(14) | |
755 | c pi=3.14159 | |
756 | c write(ifch,*)'parameters of mass formula:' | |
757 | c write(ifch,*)'theta=',thet,' epsi=',epsi | |
758 | c write(ifch,*)'as=',as,' ac=',ac | |
759 | c write(ifch,*)'dy=',dy,' dz=',dz | |
760 | c write(ifch,*)'ym=',ym | |
761 | c write(ifch,*)'cz dz zm=',cz,dz,zm | |
762 | c write(ifch,*)'sigma**1/3=',sigma**(1./3.),' rzero=',rzero | |
763 | c write(ifch,*)'mass:' | |
764 | c write(ifch,5000)(j,j=1,14) | |
765 | c5000 format(5x,'a:',14i5) | |
766 | c do 4 j=1,14 | |
767 | c a=j | |
768 | c ymi(j)=ym*a | |
769 | c4 zmi(j)=cz/(dz/a+zm/a**(1./3.)) | |
770 | c write(ifch,5002)(ymi(j),j=1,14) | |
771 | c5002 format(1x,'ymin: ',14f5.2) | |
772 | c write(ifch,5003)(zmi(j),j=1,14) | |
773 | c5003 format(1x,'zmin: ',14f5.2) | |
774 | c do 2 i=1,15 | |
775 | c ns=11-i | |
776 | c do 3 j=1,14 | |
777 | c a=j | |
778 | c y=a-ns | |
779 | c z=0. | |
780 | c3 eng(j)=smass(a,y,z)/a | |
781 | c write(ifch,5001)ns,(eng(j),j=1,14) | |
782 | c5001 format(1x,'s=',i2,2x,14f5.2) | |
783 | c2 continue | |
784 | c write(ifch,*)'mass-mass(free):' | |
785 | c write(ifch,5000)(j,j=1,14) | |
786 | c do 5 i=1,15 | |
787 | c ns=11-i | |
788 | c do 6 j=1,14 | |
789 | c a=j | |
790 | c y=a-ns | |
791 | c z=0. | |
792 | c call smassu(a,y,z,ku,kd,ks,kc) | |
793 | c6 eng(j)=(smass(a,y,z)-utamnu(ku,kd,ks,kc,0,0,3))/a | |
794 | c write(ifch,5001)ns,(eng(j),j=1,14) | |
795 | c5 continue | |
796 | c | |
797 | c stop | |
798 | c end | |
799 | c | |
800 | cc----------------------------------------------------------------------- | |
801 | c subroutine smasst(kux,kdx,ksx,kcx,a,y,z) | |
802 | cc----------------------------------------------------------------------- | |
803 | cc input: kux,kdx,ksx,kcx = net quark numbers (for u,d,s,c quarks). | |
804 | cc output: massnr a, hypercharge y and charge z. | |
805 | cc----------------------------------------------------------------------- | |
806 | c sg=1 | |
807 | c if(kux+kdx+ksx+kcx.lt.0.)sg=-1 | |
808 | c ku=sg*kux | |
809 | c kd=sg*kdx | |
810 | c ks=sg*ksx | |
811 | c kc=sg*kcx | |
812 | c k=ku+kd+ks+kc | |
813 | c if(mod(k,3).ne.0)stop'noninteger baryon number' | |
814 | c a=k/3 | |
815 | c y=a-ks | |
816 | c nz=2*ku-kd-ks+2*kc | |
817 | c if(mod(nz,3).ne.0)stop'noninteger charge' | |
818 | c z=nz/3 | |
819 | c return | |
820 | c end | |
821 | c | |
822 | cc----------------------------------------------------------------------- | |
823 | c subroutine smassu(ax,yx,zx,ku,kd,ks,kc) | |
824 | cc----------------------------------------------------------------------- | |
825 | cc input: massnr ax, hypercharge yx and charge zx. | |
826 | cc output: ku,kd,ks,kc = net quark numbers (for u,d,s,c quarks). | |
827 | cc----------------------------------------------------------------------- | |
828 | c sg=1 | |
829 | c if(ax.lt.0.)sg=-1 | |
830 | c a=sg*ax | |
831 | c y=sg*yx | |
832 | c z=sg*zx | |
833 | c ku=nint(a+z) | |
834 | c kd=nint(a-z+y) | |
835 | c ks=nint(a-y) | |
836 | c kc=0 | |
837 | c return | |
838 | c end | |
839 | c | |
840 | cc----------------------------------------------------------------------- | |
841 | c function spoc(a,b,c,d,x) | |
842 | cc----------------------------------------------------------------------- | |
843 | cc power fctn with cutoff | |
844 | cc----------------------------------------------------------------------- | |
845 | c spoc=0 | |
846 | c if(a.eq.0..and.b.eq.0.)return | |
847 | c spoc =a+b*x**c | |
848 | c spoc0=a+b*d**c | |
849 | c spoc=amin1(spoc,spoc0) | |
850 | c spoc=amax1(0.,spoc) | |
851 | c return | |
852 | c end | |
853 | c | |
854 | c----------------------------------------------------------------------- | |
855 | function utacos(x) | |
856 | c----------------------------------------------------------------------- | |
857 | c returns acos(x) for -1 <= x <= 1 , acos(+-1) else | |
858 | c----------------------------------------------------------------------- | |
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 | |
879 | ||
880 | c---------------------------------------------------------------------- | |
881 | function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus) | |
882 | c---------------------------------------------------------------------- | |
883 | c returns min mass of droplet with given u,d,s,c content | |
884 | c keux: net u quark number | |
885 | c kedx: net d quark number | |
886 | c kesx: net s quark number | |
887 | c kecx: net c quark number | |
888 | c kebx: net b quark number | |
889 | c ketx: net t quark number | |
890 | c modus: 4=two lowest multiplets; 5=lowest multiplet | |
891 | c---------------------------------------------------------------------- | |
892 | common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr | |
893 | common/csjcga/amnull,asuha(7) | |
894 | common/drop4/asuhax(7),asuhay(7) | |
895 | ||
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 ) | |
899 | c write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx c write(ifmt,*)'wrong mass in gethad ',damss | |
900 | ||
901 | amnull=0. | |
902 | ||
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) | |
909 | ||
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 | ||
926 | c write(ifch,*)keu,ked,kes,kec,keb,ket | |
927 | ||
928 | c removing top mesons to remove t quarks or antiquarks | |
929 | if(ket.ne.0)then | |
930 | 12 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 | ||
942 | c removing bottom mesons to remove b quarks or antiquarks | |
943 | if(keb.ne.0)then | |
944 | 11 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 | ||
956 | c removing charm mesons to remove c quarks or antiquarks | |
957 | if(kec.ne.0)then | |
958 | 10 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 | |
969 | ||
970 | c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull | |
971 | ||
972 | c removing mesons to remove s antiquarks | |
973 | 5 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 | |
984 | ||
985 | c removing mesons to remove d antiquarks | |
986 | 6 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 | |
998 | ||
999 | c removing mesons to remove u antiquarks | |
1000 | 7 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 | |
1012 | ||
1013 | c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull | |
1014 | c print*,keu,ked,kes,kec,keb,ket,amnull | |
1015 | ||
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 | |
1021 | ||
1022 | c removing strange baryons | |
1023 | i=4 | |
1024 | 2 i=i-1 | |
1025 | 3 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 | |
1042 | 8 continue | |
1043 | endif | |
1044 | ||
1045 | if(keu+ked.ne.keq)call utstop('utamnu: keu+ked /= keq&') | |
1046 | c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux | |
1047 | c print*,keu,ked,kes,kec,keb,ket,amnull+amnux | |
1048 | ||
1049 | c removing nonstrange baryons | |
1050 | 9 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 | |
1064 | ||
1065 | c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux | |
1066 | c print*,keu,ked,kes,kec,keb,ket,amnull+amnux | |
1067 | ||
1068 | if(mod(keq,3).ne.0)call utstop('utamnu: mod(keq,3) /= 0&') | |
1069 | amnux=amnux+asuha(1)*keq/3 | |
1070 | ||
1071 | c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux | |
1072 | c print*,keu,ked,kes,kec,keb,ket,amnull+amnux | |
1073 | ||
1074 | amnull=amnull+amnux | |
1075 | ||
1076 | if(amnull.eq.0)amnull=asuha(5) | |
1077 | ||
1078 | 1000 utamnu=amnull | |
1079 | return | |
1080 | end | |
1081 | ||
1082 | c----------------------------------------------------------------------- | |
1083 | function utamnx(jcp,jcm) | |
1084 | c----------------------------------------------------------------------- | |
1085 | c returns minimum mass for the decay of jcp---jcm (by calling utamnu). | |
1086 | c----------------------------------------------------------------------- | |
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 | |
1103 | 1 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 | |
1118 | 2 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) | |
1165 | c print *,amms1,amms3,amms2,amms4,jcp,jcm | |
1166 | return | |
1167 | end | |
1168 | ||
1169 | ||
1170 | ||
1171 | cc----------------------------------------------------------------------- | |
1172 | c function utamny(jcp,jcm) | |
1173 | cc----------------------------------------------------------------------- | |
1174 | cc returns minimum mass of jcp+jcm (by calling utamnu). | |
1175 | cc----------------------------------------------------------------------- | |
1176 | c parameter (nflav=6) | |
1177 | c integer jcp(nflav,2),jcm(nflav,2),jc(nflav,2) | |
1178 | c do 7 nf=1,nflav | |
1179 | c jc(nf,1)=jcp(nf,1)+jcm(nf,1) | |
1180 | c7 jc(nf,2)=jcp(nf,2)+jcm(nf,2) | |
1181 | c keu=jc(1,1)-jc(1,2) | |
1182 | c ked=jc(2,1)-jc(2,2) | |
1183 | c kes=jc(3,1)-jc(3,2) | |
1184 | c kec=jc(4,1)-jc(4,2) | |
1185 | c keb=jc(5,1)-jc(5,2) | |
1186 | c ket=jc(6,1)-jc(6,2) | |
1187 | c utamny=utamnu(keu,ked,kes,kec,keb,ket,5) | |
1188 | c return | |
1189 | c end | |
1190 | c | |
1191 | c----------------------------------------------------------------------- | |
1192 | function utamnz(jc,modus) | |
1193 | c----------------------------------------------------------------------- | |
1194 | c returns minimum mass of jc (by calling utamnu). | |
1195 | c----------------------------------------------------------------------- | |
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 | |
1207 | ||
1208 | c----------------------------------------------------------------------- | |
1209 | subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx) | |
1210 | c----------------------------------------------------------------------- | |
1211 | c returns the array xx with xx(1)=x0 <= xx(i) <= xx(i3)=x3 | |
1212 | c----------------------------------------------------------------------- | |
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 | ||
1223 | cc--------------------------------------------------------------------- | |
1224 | c subroutine utaxis(i,j,a1,a2,a3) | |
1225 | cc----------------------------------------------------------------------- | |
1226 | cc calculates the axis defined by the ptls i,j in the i,j cm system | |
1227 | cc--------------------------------------------------------------------- | |
1228 | c include 'epos.inc' | |
1229 | c double precision pi1,pi2,pi3,pi4,pj1,pj2,pj3,pj4,p1,p2,p3,p4,p5 | |
1230 | c *,err,a | |
1231 | c a1=0 | |
1232 | c a2=0 | |
1233 | c a3=1 | |
1234 | c pi1=dble(pptl(1,i)) | |
1235 | c pi2=dble(pptl(2,i)) | |
1236 | c pi3=dble(pptl(3,i)) | |
1237 | c pi4=dble(pptl(4,i)) | |
1238 | c pj1=dble(pptl(1,j)) | |
1239 | c pj2=dble(pptl(2,j)) | |
1240 | c pj3=dble(pptl(3,j)) | |
1241 | c pj4=dble(pptl(4,j)) | |
1242 | c p1=pi1+pj1 | |
1243 | c p2=pi2+pj2 | |
1244 | c p3=pi3+pj3 | |
1245 | c p4=pi4+pj4 | |
1246 | c p5=dsqrt(p4**2-p3**2-p2**2-p1**2) | |
1247 | c call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50) | |
1248 | c call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51) | |
1249 | c err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2 | |
1250 | c if(err.gt.1d-3)then | |
1251 | c call utmsg('utaxis') | |
1252 | c write(ifch,*)'***** err=',err | |
1253 | c write(ifch,*)'pi:',pi1,pi2,pi3,pi4 | |
1254 | c write(ifch,*)'pj:',pj1,pj2,pj3,pj4 | |
1255 | c call utmsgf | |
1256 | c endif | |
1257 | c a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 ) | |
1258 | c if(a.eq.0.d0)return | |
1259 | c a1=sngl((pi1-pj1)/a) | |
1260 | c a2=sngl((pi2-pj2)/a) | |
1261 | c a3=sngl((pi3-pj3)/a) | |
1262 | c return | |
1263 | c end | |
1264 | c | |
1265 | cc----------------------------------------------------------------------- | |
1266 | c subroutine utchm(arp,arm,ii) | |
1267 | cc----------------------------------------------------------------------- | |
1268 | cc checks whether arp**2=0 and arm**2=0. | |
1269 | cc----------------------------------------------------------------------- | |
1270 | c include 'epos.inc' | |
1271 | c double precision arp(4),arm(4),difp,difm | |
1272 | c difp=arp(4)**2-arp(1)**2-arp(2)**2-arp(3)**2 | |
1273 | c difm=arm(4)**2-arm(1)**2-arm(2)**2-arm(3)**2 | |
1274 | c if(dabs(difp).gt.1e-3*arp(4)**2 | |
1275 | c *.or.dabs(difm).gt.1e-3*arm(4)**2)then | |
1276 | c call utmsg('utchm ') | |
1277 | c write(ifch,*)'***** mass non zero - ',ii | |
1278 | c write(ifch,*)'jet-mass**2`s: ',difp,difm | |
1279 | c write(ifch,*)'energy**2`s: ',arp(4)**2,arm(4)**2 | |
1280 | c write(ifch,*)(sngl(arp(i)),i=1,4) | |
1281 | c write(ifch,*)(sngl(arm(i)),i=1,4) | |
1282 | c call utmsgf | |
1283 | c endif | |
1284 | c return | |
1285 | c end | |
1286 | c | |
1287 | c----------------------------------------------------------------------- | |
1288 | subroutine utclea(nptlii,nptl0) | |
1289 | c----------------------------------------------------------------------- | |
1290 | c starting from nptlii | |
1291 | c overwrites istptl=99 particles in /cptl/, reduces so nptl | |
1292 | c and update minfra and maxfra | |
1293 | c----------------------------------------------------------------------- | |
1294 | include 'epos.inc' | |
1295 | integer newptl(mxptl)!,oldptl(mxptl),ii(mxptl) | |
1296 | ||
1297 | ish0=ish | |
1298 | if(ishsub/100.eq.18)ish=mod(ishsub,100) | |
1299 | ||
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) | |
1317 | 34 continue | |
1318 | 116 format(1x,i6,i6,4x,i6,4x,i6,i6,i12,3x,3(e8.2,1x),i3,i3) | |
1319 | endif | |
1320 | ||
1321 | c ish=ish0 | |
1322 | c ish0=ish | |
1323 | c if(ishsub/100.eq.18)ish=mod(ishsub,100) | |
1324 | ||
1325 | i=nptli-1 | |
1326 | 1 i=i+1 | |
1327 | if(i.gt.nptl)goto 1000 | |
1328 | if(istptl(i).eq.99)goto 2 | |
1329 | newptl(i)=i | |
1330 | c oldptl(i)=i | |
1331 | goto 1 | |
1332 | ||
1333 | 2 i=i-1 | |
1334 | j=i | |
1335 | 3 i=i+1 | |
1336 | 4 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 | |
1341 | c oldptl(i)=j | |
1342 | c write(ifch,*)'move',j,' to ',i | |
1343 | c write(ifch,*)idptl(i),ityptl(i),idptl(j),ityptl(j),minfra,maxfra | |
1344 | call utrepl(i,j) | |
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 | |
1350 | ||
1351 | 5 nptl=i-1 | |
1352 | if(nptl.eq.0)then | |
1353 | nptl0=0 | |
1354 | goto 1000 | |
1355 | endif | |
1356 | ||
1357 | 20 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 | ||
1366 | c do 11 k=1,nptl | |
1367 | c io=iorptl(k) | |
1368 | c if(io.le.0)ii(k)=io | |
1369 | c if(io.gt.0)ii(k)=newptl(io) | |
1370 | c11 continue | |
1371 | c do 12 k=1,nptl | |
1372 | c12 iorptl(k)=ii(k) | |
1373 | c | |
1374 | c do 13 k=1,nptl | |
1375 | c jo=jorptl(k) | |
1376 | c if(jo.le.0)ii(k)=jo | |
1377 | c if(jo.gt.0)ii(k)=newptl(jo) | |
1378 | c13 continue | |
1379 | c do 14 k=1,nptl | |
1380 | c14 jorptl(k)=ii(k) | |
1381 | c | |
1382 | c do 15 k=1,nptl | |
1383 | c if1=ifrptl(1,k) | |
1384 | c if(if1.le.0)ii(k)=if1 | |
1385 | c if(if1.gt.0)ii(k)=newptl(if1) | |
1386 | c15 continue | |
1387 | c do 16 k=1,nptl | |
1388 | c16 ifrptl(1,k)=ii(k) | |
1389 | c | |
1390 | c do 17 k=1,nptl | |
1391 | c if2=ifrptl(2,k) | |
1392 | c if(if2.le.0)ii(k)=if2 | |
1393 | c if(if2.gt.0)ii(k)=newptl(if2) | |
1394 | c17 continue | |
1395 | c do 18 k=1,nptl | |
1396 | c18 ifrptl(2,k)=ii(k) | |
1397 | c | |
1398 | c do 19 k=1,nptl | |
1399 | c if(ifrptl(1,k).eq.0.and.ifrptl(2,k).gt.0)ifrptl(1,k)=ifrptl(2,k) | |
1400 | c if(ifrptl(2,k).eq.0.and.ifrptl(1,k).gt.0)ifrptl(2,k)=ifrptl(1,k) | |
1401 | c19 continue | |
1402 | ||
1403 | 1000 continue | |
1404 | ||
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) | |
1414 | 35 continue | |
1415 | endif | |
1416 | ||
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 | ||
1425 | c--------------------------------------------------------------------- | |
1426 | subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q) | |
1427 | c--------------------------------------------------------------------- | |
1428 | c linear fit to data | |
1429 | c input: | |
1430 | c ndata: nr of data points | |
1431 | c x(),y(),sig(): data | |
1432 | c mwt: unweighted (0) or weighted (else) data points | |
1433 | c output: | |
1434 | c a,b: parameters of linear fit a+b*x | |
1435 | c--------------------------------------------------------------------- | |
1436 | INTEGER mwt,ndata | |
1437 | REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata) | |
1438 | CU USES utgmq | |
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 | |
1452 | 11 continue | |
1453 | else | |
1454 | do 12 i=1,ndata | |
1455 | sx=sx+x(i) | |
1456 | sy=sy+y(i) | |
1457 | 12 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) | |
1466 | 13 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) | |
1472 | 14 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 | |
1482 | 15 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 | |
1490 | 16 continue | |
1491 | q=utgmq(0.5*(ndata-2),0.5*chi2) | |
1492 | endif | |
1493 | return | |
1494 | END | |
1495 | ||
1496 | c----------------------------------------------------------------------- | |
1497 | function utgam1(x) | |
1498 | c----------------------------------------------------------------------- | |
1499 | c gamma fctn tabulated | |
1500 | c single precision | |
1501 | c----------------------------------------------------------------------- | |
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 | ||
1519 | c----------------------------------------------------------------------- | |
1520 | double precision function utgam2(x) | |
1521 | c----------------------------------------------------------------------- | |
1522 | c gamma fctn tabulated | |
1523 | c double precision | |
1524 | c----------------------------------------------------------------------- | |
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 | ||
1542 | c----------------------------------------------------------------------- | |
1543 | double precision function utgam(x) | |
1544 | c----------------------------------------------------------------------- | |
1545 | c gamma fctn | |
1546 | c double precision | |
1547 | c----------------------------------------------------------------------- | |
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))/ | |
1573 | 2 ((((((c(8)*z+c(9))*z+c(10))*z+c(11))*z+c(12))*z+c(13))*z+1.0d0) | |
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 | ||
1585 | c--------------------------------------------------------------------- | |
1586 | subroutine utgcf(gammcf,a,x,gln) | |
1587 | c--------------------------------------------------------------------- | |
1588 | INTEGER ITMAX | |
1589 | REAL a,gammcf,gln,x,EPS,FPMIN | |
1590 | PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30) | |
1591 | CU 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 | |
1610 | 11 continue | |
1611 | pause 'a too large, ITMAX too small in utgcf' | |
1612 | 1 gammcf=exp(-x+a*log(x)-gln)*h | |
1613 | return | |
1614 | END | |
1615 | ||
1616 | c--------------------------------------------------------------------- | |
1617 | function utgmln(xx) | |
1618 | c--------------------------------------------------------------------- | |
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 | |
1634 | 11 continue | |
1635 | utgmln=tmp+log(stp*ser/x) | |
1636 | return | |
1637 | END | |
1638 | ||
1639 | c--------------------------------------------------------------------- | |
1640 | function utgmq(a,x) | |
1641 | c--------------------------------------------------------------------- | |
1642 | REAL a,utgmq,x | |
1643 | CU 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 | ||
1656 | c--------------------------------------------------------------------- | |
1657 | subroutine utgser(gamser,a,x,gln) | |
1658 | c--------------------------------------------------------------------- | |
1659 | INTEGER ITMAX | |
1660 | REAL a,gamser,gln,x,EPS | |
1661 | PARAMETER (ITMAX=100,EPS=3.e-7) | |
1662 | CU 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 | |
1679 | 11 continue | |
1680 | pause 'a too large, ITMAX too small in utgser' | |
1681 | 1 gamser=sum*exp(-x+a*log(x)-gln) | |
1682 | return | |
1683 | END | |
1684 | ||
1685 | c------------------------------------------------------------------------- | |
1686 | subroutine uticpl(ic,ifla,iqaq,iret) | |
1687 | c------------------------------------------------------------------------- | |
1688 | c adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla | |
1689 | c to 2-id ic | |
1690 | c------------------------------------------------------------------------- | |
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 | ||
1714 | cc----------------------------------------------------------------------- | |
1715 | c subroutine utindx(n,xar,x,i) | |
1716 | cc----------------------------------------------------------------------- | |
1717 | cc input: dimension n | |
1718 | cc array xar(n) with xar(i) > xar(i-1) | |
1719 | cc some number x between xar(1) and xar(n) | |
1720 | cc output: the index i such that x is between xar(i) and xar(i+1) | |
1721 | cc----------------------------------------------------------------------- | |
1722 | c include 'epos.inc' | |
1723 | c real xar(n) | |
1724 | c if(x.lt.xar(1))then | |
1725 | c if(ish.ge.5)then | |
1726 | c call utmsg('utindx') | |
1727 | c write(ifch,*)'***** x=',x,' < xar(1)=',xar(1) | |
1728 | c call utmsgf | |
1729 | c endif | |
1730 | c i=1 | |
1731 | c return | |
1732 | c elseif(x.gt.xar(n))then | |
1733 | c if(ish.ge.5)then | |
1734 | c call utmsg('utindx') | |
1735 | c write(ifch,*)'***** x=',x,' > xar(n)=',xar(n) | |
1736 | c call utmsgf | |
1737 | c endif | |
1738 | c i=n | |
1739 | c return | |
1740 | c endif | |
1741 | c lu=1 | |
1742 | c lo=n | |
1743 | c1 lz=(lo+lu)/2 | |
1744 | c if((xar(lu).le.x).and.(x.le.xar(lz)))then | |
1745 | c lo=lz | |
1746 | c elseif((xar(lz).lt.x).and.(x.le.xar(lo)))then | |
1747 | c lu=lz | |
1748 | c else | |
1749 | c call utstop('utindx: no interval found&') | |
1750 | c endif | |
1751 | c if((lo-lu).ge.2) goto1 | |
1752 | c if(lo.le.lu)call utstop('utinvt: lo.le.lu&') | |
1753 | c i=lu | |
1754 | c return | |
1755 | c end | |
1756 | c | |
1757 | c----------------------------------------------------------------------- | |
1758 | function utinvt(n,x,q,y) | |
1759 | c----------------------------------------------------------------------- | |
1760 | c returns x with y=q(x) | |
1761 | c----------------------------------------------------------------------- | |
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 | |
1782 | 1 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 | |
1798 | ||
1799 | c----------------------------------------------------------------------- | |
1800 | subroutine utlob2(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4,idi) | |
1801 | c----------------------------------------------------------------------- | |
1802 | c performs a lorentz boost, double prec. | |
1803 | c isig=+1 is to boost the four vector x1,x2,x3,x4 such as to obtain it | |
1804 | c in the frame specified by the 5-vector p1...p5 (5-vector=4-vector+mass). | |
1805 | c isig=-1: the other way round, that means, | |
1806 | c if the 4-vector x1...x4 is given in some frame characterized by | |
1807 | c p1...p5 with respect to to some lab-frame, utlob2 returns the 4-vector | |
1808 | c x1...x4 in the lab frame. | |
1809 | c idi is a call identifyer (integer) to identify the call in case of problem | |
1810 | c----------------------------------------------------------------------- | |
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 | |
1818 | 101 format(' utlob2: x = ',5e13.5) | |
1819 | 301 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 | |
1856 | 220 bp=bp+z(k)*isig*beta(k) | |
1857 | do 230 k=1,3 | |
1858 | 230 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):' | |
1884 | 102 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 | ||
1895 | c----------------------------------------------------------------------- | |
1896 | subroutine utlob3(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4) | |
1897 | c----------------------------------------------------------------------- | |
1898 | c performs a lorentz boost, double prec. | |
1899 | c but arguments are single precision | |
1900 | c----------------------------------------------------------------------- | |
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 | ||
1916 | c----------------------------------------------------------------------- | |
1917 | subroutine utlob5(yboost,x1,x2,x3,x4,x5) | |
1918 | c----------------------------------------------------------------------- | |
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 | ||
1927 | c----------------------------------------------------------------------- | |
1928 | subroutine utlob4(isig,pp1,pp2,pp3,pp4,pp5,x1,x2,x3,x4) | |
1929 | c----------------------------------------------------------------------- | |
1930 | c performs a lorentz boost, double prec. | |
1931 | c but arguments are partly single precision | |
1932 | c----------------------------------------------------------------------- | |
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 | ||
1946 | ||
1947 | c----------------------------------------------------------------------- | |
1948 | subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4) | |
1949 | c----------------------------------------------------------------------- | |
1950 | c performs a lorentz boost | |
1951 | c----------------------------------------------------------------------- | |
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 | |
1971 | 220 bp=bp+z(k)*isig*beta(k) | |
1972 | do 230 k=1,3 | |
1973 | 230 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 | |
1982 | ||
1983 | c----------------------------------------------------------------------- | |
1984 | subroutine utloc(ar,n,a,l) | |
1985 | c----------------------------------------------------------------------- | |
1986 | real ar(n) | |
1987 | do 1 i=1,n | |
1988 | l=i-1 | |
1989 | if(a.lt.ar(i))return | |
1990 | 1 continue | |
1991 | l=n | |
1992 | return | |
1993 | end | |
1994 | ||
1995 | cc----------------------------------------------------------------------- | |
1996 | c subroutine utlow(cone) | |
1997 | cc----------------------------------------------------------------------- | |
1998 | c character*1 cone | |
1999 | c if(cone.eq.'A')cone='a' | |
2000 | c if(cone.eq.'B')cone='b' | |
2001 | c if(cone.eq.'C')cone='c' | |
2002 | c if(cone.eq.'D')cone='d' | |
2003 | c if(cone.eq.'E')cone='e' | |
2004 | c if(cone.eq.'F')cone='f' | |
2005 | c if(cone.eq.'G')cone='g' | |
2006 | c if(cone.eq.'H')cone='h' | |
2007 | c if(cone.eq.'I')cone='i' | |
2008 | c if(cone.eq.'J')cone='j' | |
2009 | c if(cone.eq.'K')cone='k' | |
2010 | c if(cone.eq.'L')cone='l' | |
2011 | c if(cone.eq.'M')cone='m' | |
2012 | c if(cone.eq.'N')cone='n' | |
2013 | c if(cone.eq.'O')cone='o' | |
2014 | c if(cone.eq.'P')cone='p' | |
2015 | c if(cone.eq.'Q')cone='q' | |
2016 | c if(cone.eq.'R')cone='r' | |
2017 | c if(cone.eq.'S')cone='s' | |
2018 | c if(cone.eq.'T')cone='t' | |
2019 | c if(cone.eq.'U')cone='u' | |
2020 | c if(cone.eq.'V')cone='v' | |
2021 | c if(cone.eq.'W')cone='w' | |
2022 | c if(cone.eq.'X')cone='x' | |
2023 | c if(cone.eq.'Y')cone='y' | |
2024 | c if(cone.eq.'Z')cone='z' | |
2025 | c return | |
2026 | c end | |
2027 | c | |
2028 | cc----------------------------------------------------------------------- | |
2029 | c subroutine utlow3(cthree) | |
2030 | cc----------------------------------------------------------------------- | |
2031 | c character cthree*3 | |
2032 | c do 1 i=1,3 | |
2033 | c1 call utlow(cthree(i:i)) | |
2034 | c return | |
2035 | c end | |
2036 | c | |
2037 | cc----------------------------------------------------------------------- | |
2038 | c subroutine utlow6(csix) | |
2039 | cc----------------------------------------------------------------------- | |
2040 | c character csix*6 | |
2041 | c do 1 i=1,6 | |
2042 | c1 call utlow(csix(i:i)) | |
2043 | c return | |
2044 | c end | |
2045 | c | |
2046 | cc----------------------------------------------------------------------- | |
2047 | c function utmom(k,n,x,q) | |
2048 | cc----------------------------------------------------------------------- | |
2049 | cc calculates kth moment for f(x) with q(i)=int[0,x(i)]f(z)dz | |
2050 | cc----------------------------------------------------------------------- | |
2051 | c real x(n),q(n) | |
2052 | c if(n.lt.2)call utstop('utmom : dimension too small&') | |
2053 | c utmom=0 | |
2054 | c do 1 i=2,n | |
2055 | c1 utmom=utmom+((x(i)+x(i-1))/2)**k*(q(i)-q(i-1)) | |
2056 | c utmom=utmom/q(n) | |
2057 | c return | |
2058 | c end | |
2059 | c | |
2060 | c----------------------------------------------------------------------- | |
2061 | function utpcm(a,b,c) | |
2062 | c----------------------------------------------------------------------- | |
2063 | c calculates cm momentum for a-->b+c | |
2064 | c----------------------------------------------------------------------- | |
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 | ||
2074 | c----------------------------------------------------------------------- | |
2075 | double precision function utpcmd(a,b,c) | |
2076 | c----------------------------------------------------------------------- | |
2077 | c calculates cm momentum for a-->b+c | |
2078 | c----------------------------------------------------------------------- | |
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 | ||
2091 | c----------------------------------------------------------------------- | |
2092 | subroutine utpri(text,idmmy,ishini,ishx) | |
2093 | c----------------------------------------------------------------------- | |
2094 | include 'epos.inc' | |
2095 | character*6 text | |
2096 | c 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) | |
2103 | endif | |
2104 | enddo | |
2105 | endif | |
2106 | if(ish.ge.ishx)then | |
2107 | write(ifch,'(1x,43a)') | |
2108 | * ('-',i=1,10),' entry ',text,' ',('-',i=1,30) | |
2109 | c call ranfgt(seedx) !!! | |
2110 | c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!! | |
2111 | endif | |
2112 | return | |
2113 | end | |
2114 | ||
2115 | c----------------------------------------------------------------------- | |
2116 | subroutine utprix(text,idmmy,ishini,ishx) | |
2117 | c----------------------------------------------------------------------- | |
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 | |
2126 | ||
2127 | c----------------------------------------------------------------------- | |
2128 | subroutine utprj(text,idmmy,ishini,ishx) | |
2129 | c----------------------------------------------------------------------- | |
2130 | include 'epos.inc' | |
2131 | character*20 text | |
2132 | c 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) | |
2140 | endif | |
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) | |
2146 | c call ranfgt(seedx) !!! | |
2147 | c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!! | |
2148 | endif | |
2149 | return | |
2150 | end | |
2151 | ||
2152 | c----------------------------------------------------------------------- | |
2153 | subroutine utprjx(text,idmmy,ishini,ishx) | |
2154 | c----------------------------------------------------------------------- | |
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 | ||
2165 | c----------------------------------------------------------------------- | |
2166 | function utquad(m,x,f,k) | |
2167 | c----------------------------------------------------------------------- | |
2168 | c performs an integration according to simpson | |
2169 | c----------------------------------------------------------------------- | |
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 | |
2176 | ||
2177 | c----------------------------------------------------------------------- | |
2178 | subroutine utquaf(fu,n,x,q,wo,x0,x1,x2,x3) | |
2179 | c----------------------------------------------------------------------- | |
2180 | c returns q(i) = integral [x(1)->x(i)] fu(x) dx | |
2181 | c----------------------------------------------------------------------- | |
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 | |
2200 | 3 fa(k)=fu(z) | |
2201 | q(i)=q(i-1)+utquad(m,xa,fa,m) | |
2202 | 2 continue | |
2203 | return | |
2204 | end | |
2205 | ||
2206 | c----------------------------------------------------------------------- | |
2207 | subroutine utrepl(i,j) | |
2208 | c----------------------------------------------------------------------- | |
2209 | c i is replaced by j in /cptl/ | |
2210 | c----------------------------------------------------------------------- | |
2211 | include 'epos.inc' | |
2212 | do 1 k=1,5 | |
2213 | 1 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 | |
2218 | 2 tivptl(k,i)=tivptl(k,j) | |
2219 | do 3 k=1,2 | |
2220 | 3 ifrptl(k,i)= 0 !ifrptl(k,j) | |
2221 | jorptl(i) = 0 !jorptl(j) | |
2222 | do 4 k=1,4 | |
2223 | 4 xorptl(k,i)=xorptl(k,j) | |
2224 | do 5 k=1,4 | |
2225 | 5 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 | |
2237 | ||
2238 | cc----------------------------------------------------------------------- | |
2239 | c subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten) | |
2240 | cc----------------------------------------------------------------------- | |
2241 | c parameter (nflav=6) | |
2242 | c integer icm(2),icp(2),jcm(nflav,2),jcp(nflav,2) | |
2243 | c icm(1)=icm1 | |
2244 | c icm(2)=icm2 | |
2245 | c icp(1)=icp1 | |
2246 | c icp(2)=icp2 | |
2247 | c CALL IDDECO(ICM,JCM) | |
2248 | c CALL IDDECO(ICP,JCP) | |
2249 | c do 37 nf=1,nflav | |
2250 | c do 37 k=1,2 | |
2251 | c37 jcP(nf,k)=jcp(nf,k)+jcm(nf,k) | |
2252 | c CALL IDENCO(JCP,ICP,IRETEN) | |
2253 | c IDP=IDTRA(ICP,0,0,3) | |
2254 | c call idres(idp,amp,idpr,iadj,0) | |
2255 | c return | |
2256 | c end | |
2257 | c | |
2258 | cc----------------------------------------------------------------------- | |
2259 | c subroutine utroa1(phi,a1,a2,a3,x1,x2,x3) | |
2260 | cc----------------------------------------------------------------------- | |
2261 | cc rotates x by angle phi around axis a. | |
2262 | cc normalization of a is irrelevant. | |
2263 | cc----------------------------------------------------------------------- | |
2264 | c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi | |
2265 | c dphi=phi | |
2266 | c xx(1)=x1 | |
2267 | c xx(2)=x2 | |
2268 | c xx(3)=x3 | |
2269 | c aa(1)=a1 | |
2270 | c aa(2)=a2 | |
2271 | c aa(3)=a3 | |
2272 | c aaa=0 | |
2273 | c xxx=0 | |
2274 | c do i=1,3 | |
2275 | c aaa=aaa+aa(i)**2 | |
2276 | c xxx=xxx+xx(i)**2 | |
2277 | c enddo | |
2278 | c if(xxx.eq.0d0)return | |
2279 | c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&') | |
2280 | c aaa=dsqrt(aaa) | |
2281 | c xxx=dsqrt(xxx) | |
2282 | cc e3 = a / !a! | |
2283 | c do i=1,3 | |
2284 | c e3(i)=aa(i)/aaa | |
2285 | c enddo | |
2286 | cc x_parallel | |
2287 | c xp=0 | |
2288 | c do i=1,3 | |
2289 | c xp=xp+xx(i)*e3(i) | |
2290 | c enddo | |
2291 | cc x_transverse | |
2292 | c if(xxx**2-xp**2.le.0.)return | |
2293 | c xt=dsqrt(xxx**2-xp**2) | |
2294 | cc e1 = vector x_transverse / absolute value x_transverse | |
2295 | c do i=1,3 | |
2296 | c e1(i)=(xx(i)-e3(i)*xp)/xt | |
2297 | c enddo | |
2298 | cc e2 orthogonal e3,e1 | |
2299 | c call utvec2(e3,e1,e2) | |
2300 | cc rotate x | |
2301 | c do i=1,3 | |
2302 | c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i) | |
2303 | c enddo | |
2304 | cc back to single precision | |
2305 | c x1=xx(1) | |
2306 | c x2=xx(2) | |
2307 | c x3=xx(3) | |
2308 | c return | |
2309 | c end | |
2310 | c | |
2311 | cc----------------------------------------------------------------------- | |
2312 | c subroutine utroa2(phi,a1,a2,a3,x1,x2,x3) | |
2313 | cc----------------------------------------------------------------------- | |
2314 | cc rotates x by angle phi around axis a. | |
2315 | cc normalization of a is irrelevant. | |
2316 | cc double precision phi,a1,a2,a3,x1,x2,x3 | |
2317 | cc----------------------------------------------------------------------- | |
2318 | c double precision phi,a1,a2,a3,x1,x2,x3 | |
2319 | c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi | |
2320 | c dphi=phi | |
2321 | c xx(1)=x1 | |
2322 | c xx(2)=x2 | |
2323 | c xx(3)=x3 | |
2324 | c aa(1)=a1 | |
2325 | c aa(2)=a2 | |
2326 | c aa(3)=a3 | |
2327 | c aaa=0 | |
2328 | c xxx=0 | |
2329 | c do i=1,3 | |
2330 | c aaa=aaa+aa(i)**2 | |
2331 | c xxx=xxx+xx(i)**2 | |
2332 | c enddo | |
2333 | c if(xxx.eq.0d0)return | |
2334 | c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&') | |
2335 | c aaa=1.0/dsqrt(aaa) | |
2336 | cc e3 = a / !a! | |
2337 | c do i=1,3 | |
2338 | c e3(i)=aa(i)*aaa | |
2339 | c enddo | |
2340 | cc x_parallel | |
2341 | c xp=0 | |
2342 | c do i=1,3 | |
2343 | c xp=xp+xx(i)*e3(i) | |
2344 | c enddo | |
2345 | cc x_transverse | |
2346 | c if(xxx-xp**2.le.0.)return | |
2347 | c xt=dsqrt(xxx-xp**2) | |
2348 | cc e1 = vector x_transverse / absolute value x_transverse | |
2349 | c do i=1,3 | |
2350 | c e1(i)=(xx(i)-e3(i)*xp)/xt | |
2351 | c enddo | |
2352 | cc e2 orthogonal e3,e1 | |
2353 | c call utvec2(e3,e1,e2) | |
2354 | cc rotate x | |
2355 | c do i=1,3 | |
2356 | c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i) | |
2357 | c enddo | |
2358 | cc back to single precision | |
2359 | c x1=xx(1) | |
2360 | c x2=xx(2) | |
2361 | c x3=xx(3) | |
2362 | c return | |
2363 | c end | |
2364 | c | |
2365 | cc-------------------------------------------------------------------- | |
2366 | c function utroot(funcd,x1,x2,xacc) | |
2367 | cc-------------------------------------------------------------------- | |
2368 | cc combination of newton-raphson and bisection method for root finding | |
2369 | cc input: | |
2370 | cc funcd: subr returning fctn value and first derivative | |
2371 | cc x1,x2: x-interval | |
2372 | cc xacc: accuracy | |
2373 | cc output: | |
2374 | cc utroot: root | |
2375 | cc-------------------------------------------------------------------- | |
2376 | c include 'epos.inc' | |
2377 | c INTEGER MAXIT | |
2378 | c REAL utroot,x1,x2,xacc | |
2379 | c EXTERNAL funcd | |
2380 | c PARAMETER (MAXIT=100) | |
2381 | c INTEGER j | |
2382 | c REAL df,dx,dxold,f,fh,fl,temp,xh,xl | |
2383 | c call funcd(x1,fl,df) | |
2384 | c call funcd(x2,fh,df) | |
2385 | c if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.)) | |
2386 | c *call utstop('utroot: root must be bracketed&') | |
2387 | c if(fl.eq.0.)then | |
2388 | c utroot=x1 | |
2389 | c return | |
2390 | c else if(fh.eq.0.)then | |
2391 | c utroot=x2 | |
2392 | c return | |
2393 | c else if(fl.lt.0.)then | |
2394 | c xl=x1 | |
2395 | c xh=x2 | |
2396 | c else | |
2397 | c xh=x1 | |
2398 | c xl=x2 | |
2399 | c endif | |
2400 | c utroot=.5*(x1+x2) | |
2401 | c dxold=abs(x2-x1) | |
2402 | c dx=dxold | |
2403 | c call funcd(utroot,f,df) | |
2404 | c do 11 j=1,MAXIT | |
2405 | c if(((utroot-xh)*df-f)*((utroot-xl)*df-f).ge.0..or. abs(2.* | |
2406 | c *f).gt.abs(dxold*df) ) then | |
2407 | c dxold=dx | |
2408 | c dx=0.5*(xh-xl) | |
2409 | c utroot=xl+dx | |
2410 | c if(xl.eq.utroot)return | |
2411 | c else | |
2412 | c dxold=dx | |
2413 | c dx=f/df | |
2414 | c temp=utroot | |
2415 | c utroot=utroot-dx | |
2416 | c if(temp.eq.utroot)return | |
2417 | c endif | |
2418 | c if(abs(dx).lt.xacc) return | |
2419 | c call funcd(utroot,f,df) | |
2420 | c if(f.lt.0.) then | |
2421 | c xl=utroot | |
2422 | c else | |
2423 | c xh=utroot | |
2424 | c endif | |
2425 | c11 continue | |
2426 | c call utmsg('utroot') | |
2427 | c write(ifch,*)'***** exceeding maximum iterations' | |
2428 | c write(ifch,*)'dx:',dx | |
2429 | c call utmsgf | |
2430 | c return | |
2431 | c END | |
2432 | c | |
2433 | c----------------------------------------------------------------------- | |
2434 | subroutine utrot2(isig,ax,ay,az,x,y,z) | |
2435 | c----------------------------------------------------------------------- | |
2436 | c performs a rotation, double prec. | |
2437 | c----------------------------------------------------------------------- | |
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 | |
2482 | ||
2483 | ||
2484 | c----------------------------------------------------------------------- | |
2485 | subroutine utrot4(isig,ax,ay,az,x,y,z) | |
2486 | c----------------------------------------------------------------------- | |
2487 | c performs a rotation, double prec. | |
2488 | c arguments partly single | |
2489 | c----------------------------------------------------------------------- | |
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 | ||
2501 | c----------------------------------------------------------------------- | |
2502 | subroutine utrota(isig,ax,ay,az,x,y,z) | |
2503 | c----------------------------------------------------------------------- | |
2504 | c performs a rotation | |
2505 | c----------------------------------------------------------------------- | |
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 | |
2541 | ||
2542 | c----------------------------------------------------------------------- | |
2543 | subroutine utstop(text) | |
2544 | c----------------------------------------------------------------------- | |
2545 | c returns error message and stops execution. | |
2546 | c text is an optonal text to appear in the error message. | |
2547 | c text is a character string of length 40; | |
2548 | c for shorter text, it has to be terminated by &; | |
2549 | c example: call utstop('error in subr xyz&') | |
2550 | c----------------------------------------------------------------------- | |
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) | |
2563 | 101 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) | |
2569 | 1 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 | |
2575 | 100 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 | |
2581 | ||
2582 | c----------------------------------------------------------------- | |
2583 | subroutine uttrap(func,a,b,s) | |
2584 | c----------------------------------------------------------------- | |
2585 | c trapezoidal method for integration. | |
2586 | c input: fctn func and limits a,b | |
2587 | c output: value s of the integral | |
2588 | c----------------------------------------------------------------- | |
2589 | include 'epos.inc' | |
2590 | ||
2591 | INTEGER JMAX | |
2592 | REAL a,b,func,s,epsr | |
2593 | EXTERNAL func | |
2594 | PARAMETER (JMAX=10) | |
2595 | CU 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 | |
2605 | 11 continue | |
2606 | ||
2607 | c-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 | ||
2619 | c----------------------------------------------------------------- | |
2620 | subroutine uttraq(func,a,b,s) | |
2621 | c----------------------------------------------------------------- | |
2622 | c trapezoidal method for integration. | |
2623 | c input: function func and limits a,b | |
2624 | c output: value s of the integral | |
2625 | c----------------------------------------------------------------- | |
2626 | ||
2627 | REAL a,b,func,s | |
2628 | EXTERNAL func | |
2629 | PARAMETER (eps=1.e-6) | |
2630 | CU USES uttras | |
2631 | INTEGER j | |
2632 | REAL olds | |
2633 | olds=-1.e30 | |
2634 | j=1 | |
2635 | 10 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 | ||
2644 | c----------------------------------------------------------------- | |
2645 | subroutine uttras(func,a,b,s,n) | |
2646 | c----------------------------------------------------------------- | |
2647 | c performs one iteration of the trapezoidal method for integration | |
2648 | c----------------------------------------------------------------- | |
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 | |
2665 | 11 continue | |
2666 | s=0.5*(s+(b-a)*sum/tnm) | |
2667 | endif | |
2668 | return | |
2669 | END | |
2670 | ||
2671 | cc----------------------------------------------------------------------- | |
2672 | c subroutine utvec1(a,b,c) | |
2673 | cc----------------------------------------------------------------------- | |
2674 | cc returns vector product c = a x b . | |
2675 | cc----------------------------------------------------------------------- | |
2676 | c real a(3),b(3),c(3) | |
2677 | c c(1)=a(2)*b(3)-a(3)*b(2) | |
2678 | c c(2)=a(3)*b(1)-a(1)*b(3) | |
2679 | c c(3)=a(1)*b(2)-a(2)*b(1) | |
2680 | c return | |
2681 | c end | |
2682 | c | |
2683 | cc----------------------------------------------------------------------- | |
2684 | c subroutine utvec2(a,b,c) | |
2685 | cc----------------------------------------------------------------------- | |
2686 | cc returns vector product c = a x b . | |
2687 | cc a,b,c double precision. | |
2688 | cc----------------------------------------------------------------------- | |
2689 | c double precision a(3),b(3),c(3) | |
2690 | c c(1)=a(2)*b(3)-a(3)*b(2) | |
2691 | c c(2)=a(3)*b(1)-a(1)*b(3) | |
2692 | c c(3)=a(1)*b(2)-a(2)*b(1) | |
2693 | c return | |
2694 | c end | |
2695 | c | |
2696 | c------------------------------------------------------------------- | |
2697 | subroutine utword(line,i,j,iqu) | |
2698 | c------------------------------------------------------------------- | |
2699 | c finds the first word of the character string line(j+1:160). | |
2700 | c the word is line(i:j) (with new i and j). | |
2701 | c if j<0 or if no word found --> new line read. | |
2702 | c a text between quotes "..." is considered one word; | |
2703 | c stronger: a text between double quotes ""..."" is consid one word | |
2704 | c stronger: a text between "{ and }" is considered one word | |
2705 | c------------------------------------------------------------------- | |
2706 | c input: | |
2707 | c line: character string (*160) | |
2708 | c i: integer between 1 and 160 | |
2709 | c iqu: for iqu=1 a "?" is written to output before reading a line, | |
2710 | c otherwise (iqu/=1) nothing is typed | |
2711 | c output: | |
2712 | c i,j: left and right end of word (word=line(i:j)) | |
2713 | c------------------------------------------------------------------- | |
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) | |
2724 | ||
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 | |
2749 | if(i.gt.160)goto5 | |
2750 | if(line(i:i).eq.'!')goto5 | |
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) | |
2810 | elseif(l2.lt.l1)then | |
2811 | line(i0:i0+l2-1)=w2define(ndf)(1:l2) | |
2812 | do k=i0+l2,i0-1+l1 | |
2813 | line(k:k)=' ' | |
2814 | enddo | |
2815 | elseif(l2.gt.l1)then | |
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 | |
2822 | endif | |
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 | ||
2834 | 9999 close(ifopx) | |
2835 | nopen=nopen-1 | |
2836 | if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1 | |
2837 | goto5 | |
2838 | end | |
2839 | ||
2840 | c-------------------------------------------------------------------- | |
2841 | subroutine utworn(line,j,ne) | |
2842 | c-------------------------------------------------------------------- | |
2843 | c returns number ne of nonempty characters of line(j+1:160) | |
2844 | c-------------------------------------------------------------------- | |
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 | ||
2853 | c----------------------------------------------------------------------- | |
2854 | subroutine getairmol(iz,ia) | |
2855 | c----------------------------------------------------------------------- | |
2856 | include 'epos.inc' | |
2857 | i=0 | |
2858 | r=rangen() | |
2859 | do while(r.gt.0.) ! choose air-molecule | |
2860 | i=i+1 | |
2861 | r=r-airwnxs(i) | |
2862 | enddo | |
2863 | iz = nint(airznxs(i)) | |
2864 | ia = nint(airanxs(i)) | |
2865 | end | |
2866 | ||
2867 | c---------------------------------------------------------------------- | |
2868 | ||
2869 | subroutine factoriel | |
2870 | ||
2871 | c---------------------------------------------------------------------- | |
2872 | c tabulation of fctrl(n)=n!, facto(n)=1/n! and utgamtab(x) for x=0 to 50 | |
2873 | c---------------------------------------------------------------------- | |
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 | |
2885 | ||
2886 | do k=1,10000 | |
2887 | x=dble(k)/100.d0 | |
2888 | utgamtab(k)=utgam(x) | |
2889 | enddo | |
2890 | ||
2891 | return | |
2892 | end | |
2893 | ||
2894 | c----------------------------------------------------------------------- | |
2895 | subroutine fremnu(amin,ca,cb,ca0,cb0,ic1,ic2,ic3,ic4) | |
2896 | c----------------------------------------------------------------------- | |
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 | |
2932 | 1 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 | ||
2946 | c----------------------------------------------------------------------- | |
2947 | function fremnux(ic) | |
2948 | c----------------------------------------------------------------------- | |
2949 | real pnll,ptq | |
2950 | common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg | |
2951 | parameter (nflav=6) | |
2952 | integer ic(2),jc(nflav,2) | |
2953 | c ic(1)=210000 | |
2954 | c 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 | ||
2966 | c----------------------------------------------------------------------- | |
2967 | function fremnux2(ic) | |
2968 | c----------------------------------------------------------------------- | |
2969 | real pnll,ptq | |
2970 | common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg | |
2971 | parameter (nflav=6) | |
2972 | integer ic(2),jc(nflav,2) | |
2973 | c ic(1)=210000 | |
2974 | c 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 | ||
2986 | c----------------------------------------------------------------------- | |
2987 | subroutine fremnx(ammax,amin,sm,ic3,ic4,iret) | |
2988 | c----------------------------------------------------------------------- | |
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 | |
2999 | c 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 |