| 9ef1c2d9 |
1 | c--------------------------------------------------------------------- |
| 60ed90b3 |
2 | subroutine timann |
| 9ef1c2d9 |
3 | c--------------------------------------------------------------------- |
| 4 | c electron-positron |
| 5 | c--------------------------------------------------------------------- |
| 6 | include 'epos.inc' |
| 7 | include 'epos.incsem' |
| 60ed90b3 |
8 | common/nfla/nfla |
| 9ef1c2d9 |
9 | parameter (ntim=1000) |
| 10 | common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim) |
| 11 | &,idaprt(2,ntim) |
| 12 | integer jcp(6,2),jcm(6,2) |
| 13 | dimension jorprt(ntim),p1(5) |
| 14 | |
| 15 | call utpri('timann',ish,ishini,5) |
| 16 | egyevt=engy |
| 17 | qsqevt=engy**2 |
| 18 | |
| 19 | 123 nprtj=1 |
| 20 | en=engy |
| 21 | q20=0. |
| 22 | nptl=0 |
| 23 | nptl=nptl+1 !electron |
| 24 | pptl(1,nptl)=0. |
| 25 | pptl(2,nptl)=0. |
| 26 | pptl(3,nptl)=sqrt(en**2/4-2.61121e-7) |
| 27 | pptl(4,nptl)=en/2. |
| 28 | pptl(5,nptl)=0.000511 |
| 29 | idptl(nptl)=12 |
| 30 | istptl(nptl)=1 |
| 31 | |
| 32 | nptl=nptl+1 !positron |
| 33 | pptl(1,nptl)=0. |
| 34 | pptl(2,nptl)=0. |
| 35 | pptl(3,nptl)=-sqrt(en**2/4-2.61121e-7) |
| 36 | pptl(4,nptl)=en/2. |
| 37 | pptl(5,nptl)=0.000511 |
| 38 | idptl(nptl)=-12 |
| 39 | istptl(nptl)=1 |
| 40 | |
| 41 | nptl=nptl+1 !virtual gamma |
| 42 | pptl(1,nptl)=0. |
| 43 | pptl(2,nptl)=0. |
| 44 | pptl(3,nptl)=0. |
| 45 | pptl(4,nptl)=en |
| 46 | pptl(5,nptl)=en |
| 47 | idptl(nptl)=10 |
| 48 | iorptl(nptl)=nptl-2 |
| 49 | jorptl(nptl)=nptl-1 |
| 50 | istptl(nptl)=1 |
| 51 | |
| 52 | pprt(1,nprtj)=0. |
| 53 | pprt(2,nprtj)=0. |
| 54 | pprt(3,nprtj)=0. |
| 60ed90b3 |
55 | pprt(4,nprtj)=en |
| 56 | pprt(5,nprtj)=en |
| 9ef1c2d9 |
57 | idaprt(1,nprtj)=2 |
| 58 | idaprt(2,nprtj)=3 |
| 59 | if(q20.gt.0.)then |
| 60 | pprt(5,nprtj)=sqrt(q20) |
| 61 | else |
| 62 | pprt(5,nprtj)=en |
| 63 | endif |
| 64 | nfla=0 |
| 65 | do i=1,naflav |
| 66 | call idmass(i,am) |
| 67 | if (2.*am.lt.en) nfla=i |
| 68 | enddo |
| 69 | if(itflav.eq.0)then |
| 60ed90b3 |
70 | s=engy**2 |
| 9ef1c2d9 |
71 | dlz=2.4 |
| 72 | amz=91.1885 |
| 73 | al=1./real(137.035989d0) |
| 74 | gf=1.16639e-5 |
| 75 | ak=sqrt(2.)*gf*amz**2/(16*pi*al) |
| 76 | chi1=ak*s*(s-amz**2)/((s-amz**2)**2+dlz**2*amz**2) |
| 77 | chi2=ak**2*s**2/((s-amz**2)**2+dlz**2*amz**2) |
| 78 | qe=-1. |
| 79 | ve=0.25-2.*qe*0.232 |
| 80 | ae=-0.5 |
| 81 | qf=2./3. |
| 60ed90b3 |
82 | vf=sign(.5,qf)-2.*qf*0.232 |
| 9ef1c2d9 |
83 | af=sign(.5,qf) |
| 84 | dsmax1= |
| 60ed90b3 |
85 | $ 2.*(qf**2-2.*qf*ve*vf*chi1 |
| 9ef1c2d9 |
86 | $ +(ae**2+ve**2)*(af**2+vf**2)*chi2) |
| 87 | $ + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2) |
| 60ed90b3 |
88 | |
| 9ef1c2d9 |
89 | qf=-1./3. |
| 60ed90b3 |
90 | vf=sign(.5,qf)-2.*qf*0.232 |
| 9ef1c2d9 |
91 | af=sign(.5,qf) |
| 92 | dsmax2= |
| 60ed90b3 |
93 | $ 2.*(qf**2-2.*qf*ve*vf*chi1 |
| 9ef1c2d9 |
94 | $ +(ae**2+ve**2)*(af**2+vf**2)*chi2) |
| 95 | $ + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2) |
| 96 | |
| 60ed90b3 |
97 | 100 iq1=1+INT(nfla*rangen()) |
| 9ef1c2d9 |
98 | call idchrg(iq1,qf) |
| 99 | ct=-1.+2.*rangen() |
| 60ed90b3 |
100 | vf=sign(.5,qf)-2.*qf*0.232 |
| 9ef1c2d9 |
101 | af=sign(.5,qf) |
| 102 | |
| 103 | dsigma= |
| 60ed90b3 |
104 | $ (1.+ct**2)*(qf**2-2.*qf*ve*vf*chi1 |
| 9ef1c2d9 |
105 | $ +(ae**2+ve**2)*(af**2+vf**2)*chi2) |
| 106 | $ + ct*(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2) |
| 107 | if(rangen().gt.dsigma/max(dsmax1,dsmax2)) goto 100 |
| 108 | else |
| 109 | iq1=itflav |
| 110 | endif |
| 111 | if(rangen().lt.0.5)iq1=-iq1 |
| 60ed90b3 |
112 | |
| 9ef1c2d9 |
113 | nprtj=nprtj+1 |
| 114 | idprt(nprtj)=iq1 |
| 115 | pprt(4,nprtj)=en/2. |
| 116 | pprt(5,nprtj)=en |
| 117 | iorprt(nprtj)=1 |
| 60ed90b3 |
118 | |
| 9ef1c2d9 |
119 | nprtj=nprtj+1 |
| 120 | idprt(nprtj)=-iq1 |
| 121 | pprt(4,nprtj)=en/2. |
| 122 | pprt(5,nprtj)=en |
| 123 | iorprt(nprtj)=1 |
| 124 | |
| 125 | call timsho(2,3) |
| 126 | |
| 127 | jorprt(1)=0 !!color-connection, no origin!! |
| 128 | jt=1 |
| 129 | if(idprt(idaprt(1,1)).lt.0)jt=2 |
| 60ed90b3 |
130 | do i=1,nprtj |
| 9ef1c2d9 |
131 | if(idaprt(1,i).ne.0) then |
| 132 | js=jt |
| 133 | if(idprt(i).lt.0.and. |
| 134 | & ((idprt(idaprt(2,i)).eq.9.and.jt.eq.1).or. |
| 60ed90b3 |
135 | & (idprt(idaprt(1,i)).eq.9.and.jt.eq.2)))then |
| 9ef1c2d9 |
136 | js=3-jt |
| 137 | elseif(idprt(i).gt.0.and.idprt(i).ne.9.and. |
| 138 | & ((idprt(idaprt(2,i)).eq.9.and.jt.eq.2).or. |
| 60ed90b3 |
139 | & (idprt(idaprt(1,i)).eq.9.and.jt.eq.1)))then |
| 9ef1c2d9 |
140 | js=3-jt |
| 141 | elseif(idprt(i).eq.9.and.idprt(idaprt(1,i)).ne.9.and. |
| 142 | & ((idprt(idaprt(1,i)).lt.0.and.jt.eq.2).or. |
| 143 | & (idprt(idaprt(1,i)).gt.0.and.jt.eq.1)))then |
| 144 | js=3-jt |
| 145 | endif |
| 146 | jorprt(idaprt(3-js,i))=jorprt(i) |
| 147 | jorprt(i)=idaprt(js,i) |
| 148 | jorprt(idaprt(js,i))=idaprt(3-js,i) |
| 149 | endif |
| 150 | enddo |
| 60ed90b3 |
151 | |
| 9ef1c2d9 |
152 | if(ish.ge.5)then |
| 60ed90b3 |
153 | i=1 |
| 9ef1c2d9 |
154 | do while(i.ne.0) |
| 60ed90b3 |
155 | if(idaprt(1,i) .eq. 0 ) then |
| 9ef1c2d9 |
156 | write(ifch,*) idprt(i) |
| 157 | write(ifch,*) '|' |
| 158 | endif |
| 159 | i=jorprt(i) |
| 160 | enddo |
| 161 | endif |
| 162 | |
| 163 | iptl=nptl |
| 60ed90b3 |
164 | i=1 |
| 9ef1c2d9 |
165 | do while(i.gt.0) |
| 60ed90b3 |
166 | if(idaprt(1,i) .eq. 0) then |
| 9ef1c2d9 |
167 | nptl=nptl+1 |
| 168 | do j=1,5 |
| 169 | pptl(j,nptl)=pprt(j,i) |
| 170 | enddo |
| 171 | idptl(nptl)=idprt(i) |
| 172 | istptl(nptl)=20 |
| 173 | iorptl(nptl)=iptl |
| 174 | qsqptl(nptl)=en**2 |
| 175 | endif |
| 176 | i=jorprt(i) |
| 177 | enddo |
| 178 | ifrptl(1,iptl)=iptl+1 |
| 179 | ifrptl(2,iptl)=nptl |
| 180 | |
| 181 | nk1=iptl+1 |
| 182 | 441 nk2=nk1+1 |
| 183 | do while (idptl(nk2).eq.9) |
| 184 | nk2=nk2+1 |
| 185 | enddo |
| 186 | do i=1,4 |
| 187 | p1(i)=0. |
| 188 | do j=nk1,nk2 |
| 189 | p1(i)=p1(i)+pptl(i,j) |
| 190 | enddo |
| 191 | enddo |
| 192 | p1(5)=sqrt(max(0.,p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2)) |
| 193 | do i=1,2 |
| 194 | do j=1,6 |
| 195 | jcm(j,i)=0 |
| 196 | jcp(j,i)=0 |
| 197 | enddo |
| 198 | enddo |
| 199 | ii=1 |
| 200 | if(idptl(nk1).lt.0)ii=2 |
| 201 | jcp(abs(idptl(nk1)),ii)=1 |
| 202 | jcm(abs(idptl(nk2)),3-ii)=1 |
| 203 | amm= utamnx(jcp,jcm) |
| 204 | if(amm.gt.p1(5))goto 123 |
| 205 | nk1=nk2+1 |
| 206 | if(nk1.lt.nptl)goto 441 |
| 207 | |
| 208 | |
| 209 | if(ish.ge.5)then |
| 210 | write(ifch,*) |
| 211 | do i=1,nprtj |
| 212 | write(ifch,98) i,(pprt(j,i),j=1,5),idprt(i) |
| 213 | & ,iorprt(i),idaprt(1,i),idaprt(2,i),jorprt(i) |
| 214 | enddo |
| 215 | write(ifch,*) |
| 216 | do i=1,nprtj |
| 60ed90b3 |
217 | if(pprt(5,i).eq.0.) |
| 9ef1c2d9 |
218 | & write(ifch,99) i,(pprt(j,i),j=1,5),idprt(i) |
| 219 | enddo |
| 220 | write(ifch,*) |
| 221 | write(ifch,*) |
| 222 | endif |
| 60ed90b3 |
223 | |
| 9ef1c2d9 |
224 | 99 format(i4,5g10.3,1i4) |
| 225 | 98 format(i4,5g10.3,5i4) |
| 60ed90b3 |
226 | |
| 9ef1c2d9 |
227 | call utprix('timann',ish,ishini,5) |
| 228 | return |
| 229 | end |
| 60ed90b3 |
230 | |
| 9ef1c2d9 |
231 | c--------------------------------------------------------------------- |
| 232 | subroutine timsh1(q20,en,idfla) |
| 233 | c--------------------------------------------------------------------- |
| 234 | include 'epos.inc' |
| 235 | include 'epos.incsem' |
| 236 | common/nfla/nfla |
| 237 | parameter (ntim=1000) |
| 238 | common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim) |
| 239 | &,idaprt(2,ntim) |
| 240 | common/cprtx/nprtjx,pprtx(5,2) |
| 241 | nfla=naflav |
| 242 | nprtj=1 |
| 243 | pprt(1,nprtj)=0. |
| 244 | pprt(2,nprtj)=0. |
| 245 | pprt(3,nprtj)=0. |
| 60ed90b3 |
246 | pprt(4,nprtj)=en |
| 9ef1c2d9 |
247 | pprt(5,nprtj)=sqrt(q20) |
| 248 | idprt(nprtj)=idfla |
| 249 | call timsho(1,0) |
| 250 | nprtjx=0 |
| 251 | return |
| 252 | end |
| 253 | |
| 254 | |
| 255 | c--------------------------------------------------------------------- |
| 60ed90b3 |
256 | subroutine timsh2(q20,q21,en,idfla1,idfla2) |
| 9ef1c2d9 |
257 | c--------------------------------------------------------------------- |
| 258 | include 'epos.inc' |
| 259 | include 'epos.incsem' |
| 260 | common/nfla/nfla |
| 261 | parameter (ntim=1000) |
| 262 | common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim) |
| 263 | &,idaprt(2,ntim) |
| 264 | common/cprtx/nprtjx,pprtx(5,2) |
| 265 | nfla=naflav |
| 266 | |
| 267 | nprtj=1 |
| 268 | pprt(1,nprtj)=0. |
| 269 | pprt(2,nprtj)=0. |
| 270 | pprt(3,nprtj)=0. |
| 60ed90b3 |
271 | pprt(4,nprtj)=en/2. |
| 9ef1c2d9 |
272 | pprt(5,nprtj)=sqrt(q20) |
| 273 | idprt(nprtj)=idfla1 |
| 274 | |
| 275 | nprtj=2 |
| 276 | pprt(1,nprtj)=0. |
| 277 | pprt(2,nprtj)=0. |
| 278 | pprt(3,nprtj)=0. |
| 279 | pprt(4,nprtj)=en/2. |
| 280 | pprt(5,nprtj)=sqrt(q21) |
| 281 | idprt(nprtj)=idfla2 |
| 282 | |
| 283 | call timsho(1,2) |
| 284 | |
| 285 | nprtjx=1 |
| 286 | pprtx(1,nprtjx)=0. |
| 287 | pprtx(2,nprtjx)=0. |
| 60ed90b3 |
288 | pprtx(3,nprtjx)=en/2. |
| 289 | pprtx(4,nprtjx)=en/2. |
| 9ef1c2d9 |
290 | pprtx(5,nprtjx)=0 |
| 291 | nprtjx=2 |
| 292 | pprtx(1,nprtjx)=0. |
| 293 | pprtx(2,nprtjx)=0. |
| 60ed90b3 |
294 | pprtx(3,nprtjx)=-en/2. |
| 9ef1c2d9 |
295 | pprtx(4,nprtjx)=en/2. |
| 296 | pprtx(5,nprtjx)=0 |
| 297 | |
| 298 | return |
| 299 | end |
| 300 | |
| 301 | c--------------------------------------------------------------------- |
| 302 | subroutine timsho(j1,j2) |
| 303 | c--------------------------------------------------------------------- |
| 304 | |
| 305 | |
| 306 | include 'epos.inc' |
| 307 | include 'epos.incsem' |
| 308 | parameter (ntim=1000) |
| 309 | common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim) |
| 310 | &,idaprt(2,ntim) |
| 311 | dimension pz(ntim),id(2,ntim) |
| 312 | dimension ij(2),ee(2),amm2(-6:10) |
| 313 | logical first |
| 314 | common/nfla/nfla |
| 315 | call utpri('timsho',ish,ishini,5) |
| 316 | if(iappl.eq.6)then |
| 317 | do i=1,6 |
| 318 | call idmass(i,am) |
| 319 | amm2(i)=am**2 |
| 320 | amm2(-i)=am**2 |
| 321 | enddo |
| 322 | else |
| 323 | do i=1,6 |
| 324 | amm2(i)=0. |
| 325 | amm2(-i)=0. |
| 326 | enddo |
| 327 | endif |
| 328 | amm2(9)=0. |
| 329 | amm2(0)=0. |
| 330 | amm2(10)=0. |
| 331 | nn=nprtj |
| 332 | ij(1)=j1 |
| 333 | ij(2)=j2 |
| 334 | ii2=2 |
| 335 | if(j2.eq.0)ii2=1 |
| 336 | |
| 337 | ii=1 |
| 338 | first=.true. |
| 339 | n10=0 |
| 340 | 10 n10=n10+1 |
| 341 | if(float(n10).gt.1e7)then |
| 342 | goto9999 |
| 343 | endif |
| 344 | io=iorprt(ij(ii)) |
| 345 | idfl=idprt(ij(ii)) |
| 346 | q2start=pprt(5,ij(ii))**2 |
| 347 | E=pprt(4,ij(ii)) |
| 348 | if(ij(2).eq.j2.and.ii2.eq.2)E=pprt(4,ij(1))+pprt(4,ij(2)) |
| 349 | zetamx=0. |
| 350 | if(ij(1).ne.j1)then |
| 60ed90b3 |
351 | zetamx=pprt(5,io)/pprt(4,io)/sqrt(pz(io)*(1.-pz(io))) |
| 9ef1c2d9 |
352 | endif |
| 353 | |
| 354 | c call timdev(idfl,q2start,E,zetamx,idfla,idflb,qa2,z) |
| 355 | c...................................................................... |
| 356 | 380 q2=q2start |
| 357 | z=0.5 |
| 358 | PT2MIN=max(qcdlam*1.1,q2fin) |
| 60ed90b3 |
359 | ALFM=LOG(PT2MIN/qcdlam) |
| 9ef1c2d9 |
360 | if(ish.ge.9)then |
| 361 | write (ifch,*) '---------------------',ii |
| 362 | $ ,pprt(5,ij(1)) !,pprt(5,ij(2)) |
| 363 | write(ifch,*) ' idfl,q2start,zetamx:',idfl,q2start,zetamx |
| 364 | endif |
| 365 | if (q2.lt.4.*q2fin+2.*amm2(idfl) )then |
| 366 | if(ish.ge.9)then |
| 6a15fbe7 |
367 | write(ifch,'(a,i4,i4,2f15.8)') 'null:',0 |
| 9ef1c2d9 |
368 | endif |
| 369 | q2=amm2(idfl) |
| 370 | idfla=0 |
| 371 | idflb=0 |
| 372 | goto 999 |
| 373 | endif |
| 60ed90b3 |
374 | |
| 9ef1c2d9 |
375 | 390 zc=.5*(1.-sqrt(max(0.000001,1.-4.*q2fin/q2))) |
| 376 | if(ish.ge.9)then |
| 377 | write(ifch,*) 'zc=',zc |
| 378 | endif |
| 60ed90b3 |
379 | IF(idfl.EQ.9) THEN |
| 380 | FBR=6.*LOG((1.-ZC)/ZC)+nfla*(0.5-ZC) |
| 9ef1c2d9 |
381 | ELSE |
| 60ed90b3 |
382 | FBR=(8./3.)*LOG((1.-ZC)/ZC) |
| 9ef1c2d9 |
383 | endif |
| 60ed90b3 |
384 | |
| 9ef1c2d9 |
385 | B0=(33.-2.*nfla)/6. |
| 60ed90b3 |
386 | |
| 387 | |
| 9ef1c2d9 |
388 | c.....select new q2 |
| 389 | r=rangen() |
| 390 | q2=q2*exp(log(r)*B0*ALFM/FBR) |
| 391 | if(ish.ge.9)then |
| 392 | write(ifch,*) 'q^2=',q2 |
| 393 | endif |
| 394 | if (q2.lt.4.*q2fin+2.*amm2(idfl))then |
| 395 | q2=amm2(idfl) |
| 396 | if(ish.ge.9)then |
| 1733af6a |
397 | write(ifch,'(a,i4,i4,2f15.8)') 'null:',0 |
| 9ef1c2d9 |
398 | endif |
| 399 | idfla=0 |
| 400 | idflb=0 |
| 401 | goto 999 |
| 402 | endif |
| 60ed90b3 |
403 | |
| 9ef1c2d9 |
404 | c.....select flavor and z-value ..................................... |
| 60ed90b3 |
405 | IF(idfl.EQ.9) THEN |
| 9ef1c2d9 |
406 | if(rangen()*FBR.lt.nfla*(0.5-ZC))then |
| 407 | ! .................g -> qqbar |
| 408 | Z=ZC+(1.-2.*ZC)*rangen() |
| 60ed90b3 |
409 | IF(Z**2+(1.-Z)**2.LT.rangen()) GOTO 390 |
| 9ef1c2d9 |
410 | idfla=int(1.+rangen()*real(nfla)) |
| 411 | idflb=-idfla |
| 412 | else !..................g -> gg |
| 413 | Z=(1.-ZC)*(ZC/(1.-ZC))**rangen() |
| 60ed90b3 |
414 | IF(rangen().GT.0.5) Z=1.-Z |
| 415 | IF((1.-Z*(1.-Z))**2.lt.rangen()) GOTO 390 |
| 9ef1c2d9 |
416 | idfla=9 |
| 417 | idflb=9 |
| 418 | endif |
| 419 | ELSE |
| 420 | Z=1.-(1.-ZC)*(ZC/(1.-ZC))**rangen() !!........q -> qg |
| 60ed90b3 |
421 | IF(1.+Z**2.LT.2.*rangen()) GOTO 390 |
| 9ef1c2d9 |
422 | idfla=idfl |
| 423 | idflb=9 |
| 424 | endif |
| 425 | |
| 8c9e69d9 |
426 | if(alfm.lt.rangen()*log(q2*z*(1.-z)/qcdlam)) goto 390 |
| 60ed90b3 |
427 | |
| 9ef1c2d9 |
428 | if(ij(1).ne.j1.or.ij(2).eq.0)then |
| f1d31e44 |
429 | if(E.le.sqrt(q2))goto 390 |
| 9ef1c2d9 |
430 | pzz=sqrt((E-sqrt(q2))*(E+sqrt(q2))) |
| 431 | pt2=(E**2*(z*(1.-z)*q2-z*amm2(idflb)-(1.-z)*amm2(idfla)) |
| 432 | $ -.25*(amm2(idflb)-amm2(idfla)-q2)**2+q2*amm2(idfla))/pzz**2 |
| 433 | if (pt2.le.0.) then |
| 434 | if(ish.ge.9)then |
| 435 | write(ifch,*) 'z not good for pt2:',z,pt2 |
| 436 | endif |
| 437 | goto 390 |
| 438 | endif |
| 439 | endif |
| 440 | |
| 441 | zeta = sqrt(q2)/E/sqrt(z*(1.-z)) |
| 442 | if(zetamx.gt.0)then |
| 443 | if (zeta.gt.zetamx)then |
| 444 | if(ish.ge.9)then |
| 445 | write(ifch,*) zeta,' > ',zetamx,'zeta-Ablehnung' |
| 446 | endif |
| 447 | goto 390 |
| 448 | endif |
| 449 | endif |
| 450 | 999 continue |
| 451 | c...................................................................... |
| 60ed90b3 |
452 | |
| 9ef1c2d9 |
453 | pprt(5,ij(ii))=sqrt(q2) |
| 454 | id(1,ii)=idfla |
| 455 | id(2,ii)=idflb |
| 456 | pz(ij(ii))=z |
| 457 | if(first)then |
| 458 | ii=2 |
| 459 | first=.false. |
| 460 | if(ii2.eq.2)goto 10 |
| 461 | endif |
| 462 | |
| 463 | if(ij(2).eq.0)then |
| 464 | z1=1. |
| 465 | pt=0. |
| 466 | alpha=0. |
| 467 | pprt(3,ij(1))=sqrt(E**2-pprt(5,ij(1))**2) |
| 468 | elseif(ij(1).eq.j1.and.ij(2).eq.j2)then |
| 469 | E=pprt(4,ij(1))+pprt(4,ij(2)) |
| 470 | ee(1)=E*.5+(pprt(5,ij(1))**2-pprt(5,ij(2))**2)/2./E |
| 471 | ee(2)=E-ee(1) |
| 472 | ii=1 |
| 473 | do while (ii.le.2) |
| 474 | if(ee(ii)-pprt(5,ij(ii)).lt.0.) then |
| 475 | if(ish.ge.5)then |
| 476 | write(ifch,*) 'goto 11' |
| 477 | endif |
| 478 | ii=1 |
| 479 | if ( pprt(5,ij(1))**2-amm2(idprt(ij(1))).lt. |
| 480 | $ pprt(5,ij(2))**2-amm2(idprt(ij(2))) ) ii=2 |
| 481 | goto 10 |
| 482 | endif |
| 483 | c zc=.5*(1.-sqrt(1.-pprt(5,ij(ii))**2/ee(ii)**2)) |
| 484 | c if(pz(ij(ii)).lt.zc.or.pz(ij(ii)).gt.1.-zc)then |
| 485 | c if(ish.ge.7)then |
| 486 | c write(ifch,*) 'first branching rejected' |
| 487 | c endif |
| 488 | c goto 10 |
| 489 | c endif |
| 490 | z=pz(ij(ii)) |
| 491 | q2=pprt(5,ij(ii))**2 |
| 492 | pzz=sqrt((ee(ii)-pprt(5,ij(ii)))*(ee(ii)+pprt(5,ij(ii)))) |
| 493 | if(pzz.gt.0.)then |
| 494 | pt2=(ee(ii)**2*(z*(1.-z)*q2-z*amm2(id(2,ii)) |
| 495 | $ -(1.-z)*amm2(id(1,ii))) |
| 496 | $ -.25*(amm2(id(2,ii))-amm2(id(1,ii))-q2)**2 |
| 497 | $ +q2*amm2(id(1,ii)))/pzz**2 |
| 498 | else |
| 499 | pt2=0. |
| 500 | endif |
| 501 | if(id(1,ii).ne.0.and.pt2.le.0.)then |
| 502 | if(ish.ge.7)then |
| 503 | write(ifch,*) 'first branching rejected for pt2',ii |
| 504 | $ ,z1,q2,ee(ii),pprt(5,ij(ii)),id(1,ii),id(2,ii) |
| 505 | endif |
| 506 | goto 10 |
| 507 | endif |
| 508 | ii=ii+1 |
| 509 | enddo |
| 510 | z1=ee(1)/E |
| 511 | z2=ee(2)/E |
| 512 | if(ish.ge.7)then |
| 513 | write(ifch,*) 'z of first branching',z1 |
| 514 | endif |
| 515 | pprt(3,ij(1))= sqrt(max(0.,ee(1)**2-pprt(5,ij(1))**2)) |
| 516 | pprt(3,ij(2))=-sqrt(max(0.,ee(2)**2-pprt(5,ij(2))**2)) |
| 517 | pt=0. |
| 518 | alpha=0. |
| 519 | else |
| 520 | E=pprt(4,io) |
| 521 | z1=pz(io) |
| 522 | z2=1.-z1 |
| 523 | am0=pprt(5,io) |
| 524 | am1=pprt(5,ij(1)) |
| 525 | am2=pprt(5,ij(2)) |
| 526 | aM=am2**2-am0**2-am1**2 |
| 527 | pzz=sqrt((E-am0)*(E+am0)) |
| 528 | pprt(3,ij(1))=.5*(aM+2.*z1*E**2)/pzz |
| 529 | pprt(3,ij(2))=pzz-pprt(3,ij(1)) |
| 530 | pt2=(E**2*(z1*z2*am0**2-z1*am2**2-z2*am1**2) |
| 531 | $ -.25*aM**2+am0**2*am1**2)/pzz**2 |
| 532 | if(ish.ge.8)then |
| 533 | write(ifch,*) 'pt2,pzz=',pt2,pzz,z1**2*E*pprt(5,io),z1,E |
| 534 | $ ,pprt(5,io) |
| 535 | endif |
| 60ed90b3 |
536 | 111 if(pt2.lt.0.) then |
| 9ef1c2d9 |
537 | ii=1 |
| 538 | if ( pprt(5,ij(1))**2-amm2(idprt(ij(1))).lt. |
| 539 | $ pprt(5,ij(2))**2-amm2(idprt(ij(2))) ) ii=2 |
| 60ed90b3 |
540 | goto 10 |
| 9ef1c2d9 |
541 | endif |
| 542 | pt=sqrt(pt2) |
| 543 | alpha=2.*pi*rangen() |
| 544 | endif |
| 545 | pprt(4,ij(1))=z1*E |
| 546 | pprt(1,ij(1))=cos(alpha)*pt |
| 547 | pprt(2,ij(1))=sin(alpha)*pt |
| 548 | if(ii2.eq.2)then |
| 549 | pprt(4,ij(2))=z2*E |
| 550 | pprt(1,ij(2))=-cos(alpha)*pt |
| 551 | pprt(2,ij(2))=-sin(alpha)*pt |
| 552 | endif |
| 6801ca52 |
553 | if(ij(1).ne.j1)then |
| 554 | if(pprt(1,io).ne.0..or.pprt(2,io).ne.0.)then |
| 9ef1c2d9 |
555 | do ii=1,2 |
| 556 | call utrota(-1,pprt(1,io),pprt(2,io),pprt(3,io) |
| 557 | & ,pprt(1,ij(ii)),pprt(2,ij(ii)),pprt(3,ij(ii))) |
| 558 | enddo |
| 559 | endif |
| 6801ca52 |
560 | endif |
| 561 | if(ij(1).ne.j1)then |
| 562 | if(pprt(3,io).lt.0.)then |
| 9ef1c2d9 |
563 | do k=1,3 |
| 564 | do ii=1,2 |
| 565 | pprt(k,ij(ii)) = -pprt(k,ij(ii)) |
| 566 | enddo |
| 567 | enddo |
| 568 | endif |
| 6801ca52 |
569 | endif |
| 9ef1c2d9 |
570 | do ii=1,ii2 |
| 571 | if(id(1,ii).ne.0)then |
| 572 | idprt(nprtj+1)=id(1,ii) |
| 573 | idprt(nprtj+2)=id(2,ii) |
| 574 | pprt(4,nprtj+1)=pz(ij(ii))*pprt(4,ij(ii)) |
| 575 | pprt(5,nprtj+1)=pz(ij(ii))*pprt(5,ij(ii)) |
| 576 | pprt(4,nprtj+2)=(1.-pz(ij(ii)))*pprt(4,ij(ii)) |
| 577 | pprt(5,nprtj+2)=(1.-pz(ij(ii)))*pprt(5,ij(ii)) |
| 578 | iorprt(nprtj+1)=ij(ii) |
| 579 | iorprt(nprtj+2)=ij(ii) |
| 580 | idaprt(1,ij(ii))=nprtj+1 |
| 581 | idaprt(2,ij(ii))=nprtj+2 |
| 582 | idaprt(1,nprtj+1)=0 |
| 583 | idaprt(2,nprtj+1)=0 |
| 584 | idaprt(1,nprtj+2)=0 |
| 585 | idaprt(2,nprtj+2)=0 |
| 586 | nprtj=nprtj+2 |
| 587 | else |
| 588 | idaprt(1,ij(ii))=0 |
| 589 | idaprt(2,ij(ii))=0 |
| 590 | endif |
| 591 | enddo |
| 592 | if(ish.ge.5)then |
| 593 | if(ij(1).ne.j1)then |
| 594 | write(ifch,98) io,(pprt(j,io),j=1,5),pz(io) |
| 595 | & ,idprt(io) |
| 596 | write(ifch,*) '->' |
| 597 | endif |
| 598 | write(ifch,99) ij(1),(pprt(j,ij(1)),j=1,5),pz(ij(1)) |
| 599 | & ,idprt(ij(1)),'->',id(1,1),id(2,1) |
| 600 | if(ij(2).ne.0)write(ifch,99) ij(2), |
| 601 | & (pprt(j,ij(2)),j=1,5),pz(ij(2)) |
| 602 | & ,idprt(ij(2)),'->',id(1,2),id(2,2) |
| 603 | write(ifch,*) |
| 604 | endif |
| 605 | |
| 606 | 98 format(i4,6g10.3,1i4) |
| 607 | 99 format(i4,6g10.3,1i4,a,2i4) |
| 608 | |
| 609 | |
| 610 | if(ij(1).le.nn)ij(1)=nn-1 |
| 611 | ij(1)=ij(1)+2 |
| 612 | ij(2)=ij(1)+1 |
| 613 | ii=1 |
| 614 | ii2=2 |
| 615 | first=.true. |
| 616 | if(ij(1).le.nprtj)goto10 |
| 617 | |
| 618 | if(ish.ge.5)then |
| 619 | write(ifch,*) |
| 620 | do i=1,nprtj |
| 60ed90b3 |
621 | write(ifch,'(i4,1x,5g10.4,2i4,a,2i4)') |
| 9ef1c2d9 |
622 | & i,(pprt(j,i),j=1,5),idprt(i) |
| 623 | & ,iorprt(i),' -->',idaprt(1,i),idaprt(2,i) |
| 624 | enddo |
| 625 | write(ifch,*) |
| 626 | endif |
| 627 | 9999 call utprix('timsho',ish,ishini,5) |
| 628 | return |
| 629 | end |
| 630 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame