File indexing completed on 2024-04-06 12:14:05
0001
0002 subroutine timann
0003
0004
0005
0006 include 'epos.inc'
0007 include 'epos.incsem'
0008 common/nfla/nfla
0009 parameter (ntim=1000)
0010 common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0011 &,iorprt(ntim),jorprt(ntim),nprtj
0012 integer jcp(6,2),jcm(6,2)
0013 dimension p1(5),jorprt2(ntim),ng(ntim)
0014
0015 call utpri('timann',ish,ishini,5)
0016 egyevt=engy
0017 qsqevt=engy**2
0018
0019 123 nprtj=1
0020 en=engy
0021 q20=0.
0022 nptl=0
0023 nptl=nptl+1 !electron
0024 pptl(1,nptl)=0.
0025 pptl(2,nptl)=0.
0026 pptl(3,nptl)=sqrt(en**2/4-2.61121e-7)
0027 pptl(4,nptl)=en/2.
0028 pptl(5,nptl)=0.000511
0029 idptl(nptl)=12
0030 istptl(nptl)=1
0031
0032 nptl=nptl+1 !positron
0033 pptl(1,nptl)=0.
0034 pptl(2,nptl)=0.
0035 pptl(3,nptl)=-sqrt(en**2/4-2.61121e-7)
0036 pptl(4,nptl)=en/2.
0037 pptl(5,nptl)=0.000511
0038 idptl(nptl)=-12
0039 istptl(nptl)=1
0040
0041 nptl=nptl+1 !virtual gamma
0042 pptl(1,nptl)=0.
0043 pptl(2,nptl)=0.
0044 pptl(3,nptl)=0.
0045 pptl(4,nptl)=en
0046 pptl(5,nptl)=en
0047 idptl(nptl)=10
0048 iorptl(nptl)=nptl-2
0049 jorptl(nptl)=nptl-1
0050 istptl(nptl)=1
0051
0052 pprt(1,nprtj)=0.
0053 pprt(2,nprtj)=0.
0054 pprt(3,nprtj)=0.
0055 pprt(4,nprtj)=en
0056 pprt(5,nprtj)=en
0057 idaprt(1,nprtj)=2
0058 idaprt(2,nprtj)=3
0059 if(q20.gt.0.)then
0060 pprt(5,nprtj)=sqrt(q20)
0061 else
0062 pprt(5,nprtj)=en
0063 endif
0064 nfla=0
0065 do i=1,naflav
0066 call idmass(i,am)
0067 if (2.*am.lt.en) nfla=i
0068 enddo
0069 if(itflav.eq.0)then
0070 s=engy**2
0071 dlz=2.4
0072 amz=91.1885
0073 al=1./real(137.035989d0)
0074 gf=1.16639e-5
0075 ak=sqrt(2.)*gf*amz**2/(16*pi*al)
0076 chi1=ak*s*(s-amz**2)/((s-amz**2)**2+dlz**2*amz**2)
0077 chi2=ak**2*s**2/((s-amz**2)**2+dlz**2*amz**2)
0078 qe=-1.
0079 ve=0.25-2.*qe*0.232
0080 ae=-0.5
0081 qf=2./3.
0082 vf=sign(.5,qf)-2.*qf*0.232
0083 af=sign(.5,qf)
0084 dsmax1=
0085 $ 2.*(qf**2-2.*qf*ve*vf*chi1
0086 $ +(ae**2+ve**2)*(af**2+vf**2)*chi2)
0087 $ + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
0088
0089 qf=-1./3.
0090 vf=sign(.5,qf)-2.*qf*0.232
0091 af=sign(.5,qf)
0092 dsmax2=
0093 $ 2.*(qf**2-2.*qf*ve*vf*chi1
0094 $ +(ae**2+ve**2)*(af**2+vf**2)*chi2)
0095 $ + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
0096
0097 100 iq1=1+INT(nfla*rangen())
0098 call idchrg(iq1,qf)
0099 ct=-1.+2.*rangen()
0100 vf=sign(.5,qf)-2.*qf*0.232
0101 af=sign(.5,qf)
0102 dsigma=
0103 $ (1.+ct**2)*(qf**2-2.*qf*ve*vf*chi1
0104 $ +(ae**2+ve**2)*(af**2+vf**2)*chi2)
0105 $ + ct*(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
0106 if(rangen().gt.dsigma/max(dsmax1,dsmax2)) goto 100
0107 else
0108 iq1=itflav
0109 endif
0110 if(rangen().lt.0.5)iq1=-iq1
0111
0112
0113 nprtj=nprtj+1
0114 idprt(nprtj)=iq1
0115 pprt(4,nprtj)=en/2.
0116 pprt(5,nprtj)=en
0117 iorprt(nprtj)=1
0118 jorprt(nprtj)=iq1
0119 q2prt(nprtj)=en**2
0120
0121 nprtj=nprtj+1
0122 idprt(nprtj)=-iq1
0123 pprt(4,nprtj)=en/2.
0124 pprt(5,nprtj)=en
0125 iorprt(nprtj)=1
0126 jorprt(nprtj)=-iq1
0127 q2prt(nprtj)=en**2
0128
0129 call timsho(2,3)
0130
0131 jorprt2(1)=0 !!color-connection, no origin!!
0132 jt=1
0133 if(idprt(idaprt(1,1)).lt.0)jt=2
0134 do i=1,nprtj
0135 ng(i)=0
0136 if(idaprt(1,i).ne.0) then
0137 js=jt
0138 if(idprt(i).lt.0.and.
0139 & ((idprt(idaprt(2,i)).eq.9.and.jt.eq.1).or.
0140 & (idprt(idaprt(1,i)).eq.9.and.jt.eq.2)))then
0141 js=3-jt
0142 elseif(idprt(i).gt.0.and.idprt(i).ne.9.and.
0143 & ((idprt(idaprt(2,i)).eq.9.and.jt.eq.2).or.
0144 & (idprt(idaprt(1,i)).eq.9.and.jt.eq.1)))then
0145 js=3-jt
0146 elseif(idprt(i).eq.9.and.idprt(idaprt(1,i)).ne.9.and.
0147 & ((idprt(idaprt(1,i)).lt.0.and.jt.eq.2).or.
0148 & (idprt(idaprt(1,i)).gt.0.and.jt.eq.1)))then
0149 js=3-jt
0150 endif
0151 jorprt2(idaprt(3-js,i))=jorprt2(i)
0152 jorprt2(i)=idaprt(js,i)
0153 jorprt2(idaprt(js,i))=idaprt(3-js,i)
0154 else
0155 j=i
0156 do while(iorprt(j).ne.0)
0157 ng(i)=ng(i)+1
0158 j=iorprt(j)
0159 enddo
0160 endif
0161 enddo
0162
0163 if(ish.ge.5)then
0164 i=1
0165 do while(i.ne.0)
0166 if(idaprt(1,i) .eq. 0 ) then
0167 write(ifch,*) idprt(i)
0168 write(ifch,*) '|'
0169 endif
0170 i=jorprt2(i)
0171 enddo
0172 endif
0173
0174 iptl=nptl
0175 i=1
0176 do while(i.gt.0)
0177 if(idaprt(1,i) .eq. 0) then
0178 nptl=nptl+1
0179 do j=1,5
0180 pptl(j,nptl)=pprt(j,i)
0181 enddo
0182 idptl(nptl)=idprt(i)
0183 istptl(nptl)=20
0184 ityptl(nptl)=30
0185 iorptl(nptl)=iptl
0186 jorptl(nptl)=jorprt(i) !type of mother parton
0187 qsqptl(nptl)=q2prt(i)
0188 endif
0189 i=jorprt2(i)
0190 enddo
0191 ifrptl(1,iptl)=iptl+1
0192 ifrptl(2,iptl)=nptl
0193
0194 nk1=iptl+1
0195 441 nk2=nk1+1
0196 do while (idptl(nk2).eq.9)
0197 nk2=nk2+1
0198 enddo
0199 do i=1,4
0200 p1(i)=0.
0201 do j=nk1,nk2
0202 p1(i)=p1(i)+pptl(i,j)
0203 enddo
0204 enddo
0205 p1(5)=sqrt(max(0.,p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2))
0206 do i=1,2
0207 do j=1,6
0208 jcm(j,i)=0
0209 jcp(j,i)=0
0210 enddo
0211 enddo
0212 ii=1
0213 if(idptl(nk1).lt.0)ii=2
0214 jcp(abs(idptl(nk1)),ii)=1
0215 jcm(abs(idptl(nk2)),3-ii)=1
0216 amm= utamnx(jcp,jcm)
0217 if(amm.gt.p1(5))goto 123
0218 nk1=nk2+1
0219 if(nk1.lt.nptl)goto 441
0220
0221
0222 if(ish.ge.5)then
0223 write(ifch,*)
0224 do i=1,nprtj
0225 write(ifch,98) i,(pprt(j,i),j=1,5),idprt(i)
0226 & ,iorprt(i),idaprt(1,i),idaprt(2,i),jorprt2(i),ng(i)
0227 enddo
0228 write(ifch,*)
0229 do i=1,nprtj
0230 if(pprt(5,i).eq.0.)
0231 & write(ifch,99) i,(pprt(j,i),j=1,5),idprt(i),ng(i)
0232 enddo
0233 write(ifch,*)
0234 write(ifch,*)
0235 endif
0236
0237 99 format(i4,5g10.3,2i4)
0238 98 format(i4,5g10.3,6i4)
0239
0240 call utprix('timann',ish,ishini,5)
0241 return
0242 end
0243
0244
0245 subroutine timsh1(q20,en,idfla,jo)
0246
0247 include 'epos.inc'
0248 include 'epos.incsem'
0249 common/nfla/nfla
0250 parameter (ntim=1000)
0251 common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0252 &,iorprt(ntim),jorprt(ntim),nprtj
0253 common/cprtx/nprtjx,pprtx(5,2)
0254 nfla=naflav
0255 nprtj=1
0256 pprt(1,nprtj)=0.
0257 pprt(2,nprtj)=0.
0258 pprt(3,nprtj)=0.
0259 pprt(4,nprtj)=en
0260 pprt(5,nprtj)=sqrt(q20)
0261 q2prt(nprtj)=q20
0262 idprt(nprtj)=idfla
0263 iorprt(nprtj)=nprtj
0264 jorprt(nprtj)=jo
0265 call timsho(1,0)
0266 nprtjx=0
0267 return
0268 end
0269
0270
0271
0272 subroutine timsh2(q20,q21,en,idfla1,idfla2,jo1,jo2)
0273
0274 include 'epos.inc'
0275 include 'epos.incsem'
0276 common/nfla/nfla
0277 parameter (ntim=1000)
0278 common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0279 &,iorprt(ntim),jorprt(ntim),nprtj
0280 common/cprtx/nprtjx,pprtx(5,2)
0281 nfla=naflav
0282
0283 nprtj=1
0284 pprt(1,nprtj)=0.
0285 pprt(2,nprtj)=0.
0286 pprt(3,nprtj)=0.
0287 pprt(4,nprtj)=en/2.
0288 pprt(5,nprtj)=sqrt(q20)
0289 q2prt(nprtj)=q20
0290 idprt(nprtj)=idfla1
0291 iorprt(nprtj)=nprtj
0292 jorprt(nprtj)=jo1
0293
0294 nprtj=2
0295 pprt(1,nprtj)=0.
0296 pprt(2,nprtj)=0.
0297 pprt(3,nprtj)=0.
0298 pprt(4,nprtj)=en/2.
0299 pprt(5,nprtj)=sqrt(q21)
0300 q2prt(nprtj)=q21
0301 idprt(nprtj)=idfla2
0302 iorprt(nprtj)=nprtj
0303 jorprt(nprtj)=jo2
0304
0305 call timsho(1,2)
0306
0307 nprtjx=1
0308 pprtx(1,nprtjx)=0.
0309 pprtx(2,nprtjx)=0.
0310 pprtx(3,nprtjx)=en/2.
0311 pprtx(4,nprtjx)=en/2.
0312 pprtx(5,nprtjx)=0
0313 nprtjx=2
0314 pprtx(1,nprtjx)=0.
0315 pprtx(2,nprtjx)=0.
0316 pprtx(3,nprtjx)=-en/2.
0317 pprtx(4,nprtjx)=en/2.
0318 pprtx(5,nprtjx)=0
0319
0320 return
0321 end
0322
0323
0324 subroutine timsho(j1,j2)
0325
0326
0327
0328 include 'epos.inc'
0329 include 'epos.incsem'
0330 parameter (ntim=1000)
0331 common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0332 &,iorprt(ntim),jorprt(ntim),nprtj
0333 dimension pz(ntim),id(2,ntim)
0334 dimension ij(2),jo(2,2),ee(2),amm2(-6:10),q2s(2)
0335 logical first
0336 common/nfla/nfla
0337 call utpri('timsho',ish,ishini,5)
0338 if(iappl.eq.6)then
0339 do i=1,6
0340 call idmass(i,am)
0341 amm2(i)=am**2
0342 amm2(-i)=am**2
0343 enddo
0344 else
0345 do i=1,6
0346 amm2(i)=0.
0347 amm2(-i)=0.
0348 enddo
0349 endif
0350 amm2(9)=0.
0351 amm2(0)=0.
0352 amm2(10)=0.
0353 nn=nprtj
0354 ij(1)=j1
0355 ij(2)=j2
0356 ii2=2
0357 if(j2.eq.0)ii2=1
0358
0359 ii=1
0360 first=.true.
0361 n10=0
0362 10 n10=n10+1
0363 if(float(n10).gt.1e7)then
0364 goto9999
0365 endif
0366 io=iorprt(ij(ii))
0367 idfl=idprt(ij(ii))
0368 q2start=pprt(5,ij(ii))**2
0369 E=pprt(4,ij(ii))
0370 if(ij(2).eq.j2.and.ii2.eq.2)E=pprt(4,ij(1))+pprt(4,ij(2))
0371 zetamx=0.
0372 if(ij(1).ne.j1)then
0373 zetamx=pprt(5,io)/pprt(4,io)/sqrt(pz(io)*(1.-pz(io)))
0374 endif
0375
0376
0377
0378
0379 q2=q2start
0380 z=0.5
0381 PT2MIN=max(qcdlam*1.1,q2fin)
0382 ALFM=LOG(PT2MIN/qcdlam)
0383 if(ish.ge.9)then
0384 write (ifch,*) '---------------------',ii
0385 $ ,pprt(5,ij(1)) !,pprt(5,ij(2))
0386 write(ifch,*) ' idfl,q2start,zetamx:',idfl,q2start,zetamx
0387 endif
0388 if (q2.lt.4.*q2fin+2.*amm2(idfl) )then
0389 if(ish.ge.9)then
0390 write(ifch,'(a,f3.0)') 'null:',0.
0391 endif
0392 q2=amm2(idfl)
0393 idfla=0
0394 idflb=0
0395 goto 999
0396 endif
0397
0398
0399 390 zc=.5*(1.-sqrt(max(0.000001,1.-4.*q2fin/q2)))
0400 if(ish.ge.9)then
0401 write(ifch,*) 'zc=',zc
0402 endif
0403 IF(idfl.EQ.9) THEN
0404 FBRqqb=nfla*(0.5-ZC)
0405 FBR=6.*LOG((1.-ZC)/ZC)+FBRqqb
0406 ELSE
0407 FBRqqb=0.
0408 FBR=(8./3.)*LOG((1.-ZC)/ZC)
0409 endif
0410
0411 B0=(33.-2.*nfla)/6.
0412
0413
0414
0415 r=rangen()
0416 q2=q2*exp(log(r)*B0*ALFM/FBR)
0417 if(ish.ge.9)then
0418 write(ifch,*) 'q^2=',q2
0419 endif
0420 if (q2.lt.4.*q2fin+2.*amm2(idfl))then
0421 q2=amm2(idfl)
0422 if(ish.ge.9)then
0423 write(ifch,'(a,f3.0)') 'null:',0.
0424 endif
0425 idfla=0
0426 idflb=0
0427 goto 999
0428 endif
0429
0430
0431
0432 IF(idfl.EQ.9) THEN
0433 if(rangen()*FBR.lt.FBRqqb)then
0434 ! .................g -> qqbar
0435 Z=ZC+(1.-2.*ZC)*rangen()
0436 IF(Z**2+(1.-Z)**2.LT.rangen()) GOTO 390
0437 idfla=int(1.+rangen()*real(nfla))
0438 idflb=-idfla
0439 else !..................g -> gg
0440 Z=(1.-ZC)*(ZC/(1.-ZC))**rangen()
0441 IF(rangen().GT.0.5) Z=1.-Z
0442 IF((1.-Z*(1.-Z))**2.lt.rangen()) GOTO 390
0443 idfla=9
0444 idflb=9
0445 endif
0446 ELSE
0447 Z=1.-(1.-ZC)*(ZC/(1.-ZC))**rangen() !!........q -> qg
0448 IF(1.+Z**2.LT.2.*rangen()) GOTO 390
0449 idfla=idfl
0450 idflb=9
0451 endif
0452
0453 if(q2*z*(1.-z).le.qcdlam) goto 390
0454 if(alfm.lt.rangen()*log(q2*z*(1.-z)/qcdlam)) goto 390
0455
0456
0457 if(ij(1).ne.j1.or.ij(2).eq.0)then
0458 if(E.le.sqrt(q2))goto 390
0459 pzz=sqrt((E-sqrt(q2))*(E+sqrt(q2)))
0460 pt2=(E**2*(z*(1.-z)*q2-z*amm2(idflb)-(1.-z)*amm2(idfla))
0461 $ -.25*(amm2(idflb)-amm2(idfla)-q2)**2+q2*amm2(idfla))/pzz**2
0462 if (pt2.le.0.) then
0463 if(ish.ge.9)then
0464 write(ifch,*) 'z not good for pt2:',z,pt2
0465 endif
0466 goto 390
0467 endif
0468 endif
0469
0470 zeta = sqrt(q2)/E/sqrt(z*(1.-z))
0471 if(zetamx.gt.0)then
0472 if (zeta.gt.zetamx)then
0473 if(ish.ge.9)then
0474 write(ifch,*) zeta,' > ',zetamx,'zeta-Ablehnung'
0475 endif
0476 goto 390
0477 endif
0478 endif
0479 999 continue
0480
0481
0482 if(zeta.gt.zetacut)then
0483 q2s(ii)=q2
0484 if(idfla.eq.9)then !anything -> gluon jet
0485 jo(1,ii)=-9
0486 elseif(jorprt(ij(ii)).eq.-9)then !gluon -> quark jet
0487 jo(1,ii)=idfla
0488 else !no change
0489 jo(1,ii)=jorprt(ij(ii))
0490 endif
0491 if(idflb.eq.9)then !anything -> gluon jet
0492 jo(2,ii)=-9
0493 elseif(jorprt(ij(ii)).eq.-9)then !gluon -> quark jet
0494 jo(2,ii)=idflb
0495 else
0496 jo(2,ii)=jorprt(ij(ii))
0497 endif
0498 else !no change
0499 jo(1,ii)=jorprt(ij(ii))
0500 jo(2,ii)=jorprt(ij(ii))
0501 q2s(ii)=q2prt(ij(ii))
0502 endif
0503
0504 pprt(5,ij(ii))=sqrt(q2)
0505 id(1,ii)=idfla
0506 id(2,ii)=idflb
0507 pz(ij(ii))=z
0508 if(first)then
0509 ii=2
0510 first=.false.
0511 if(ii2.eq.2)goto 10
0512 endif
0513
0514 if(ij(2).eq.0)then
0515 z1=1.
0516 z2=0.
0517 pt=0.
0518 alpha=0.
0519 pprt(3,ij(1))=sqrt(E**2-pprt(5,ij(1))**2)
0520 elseif(ij(1).eq.j1.and.ij(2).eq.j2)then
0521 E=pprt(4,ij(1))+pprt(4,ij(2))
0522 ee(1)=E*.5+(pprt(5,ij(1))**2-pprt(5,ij(2))**2)/2./E
0523 ee(2)=E-ee(1)
0524 ii=1
0525 do while (ii.le.2)
0526 if(ee(ii)-pprt(5,ij(ii)).lt.0.) then
0527 if(ish.ge.5)then
0528 write(ifch,*) 'goto 11'
0529 endif
0530 ii=1
0531 if ( pprt(5,ij(1))**2-amm2(idprt(ij(1))).lt.
0532 $ pprt(5,ij(2))**2-amm2(idprt(ij(2))) ) ii=2
0533 goto 10
0534 endif
0535
0536
0537
0538
0539
0540
0541
0542 z=pz(ij(ii))
0543 q2=pprt(5,ij(ii))**2
0544 pzz=sqrt((ee(ii)-pprt(5,ij(ii)))*(ee(ii)+pprt(5,ij(ii))))
0545 if(pzz.gt.0.)then
0546 pt2=(ee(ii)**2*(z*(1.-z)*q2-z*amm2(id(2,ii))
0547 $ -(1.-z)*amm2(id(1,ii)))
0548 $ -.25*(amm2(id(2,ii))-amm2(id(1,ii))-q2)**2
0549 $ +q2*amm2(id(1,ii)))/pzz**2
0550 else
0551 pt2=0.
0552 endif
0553 if(id(1,ii).ne.0.and.pt2.le.0.)then
0554 if(ish.ge.7)then
0555 write(ifch,*) 'first branching rejected for pt2',ii
0556 $ ,z1,q2,ee(ii),pprt(5,ij(ii)),id(1,ii),id(2,ii)
0557 endif
0558 goto 10
0559 endif
0560 ii=ii+1
0561 enddo
0562 z1=ee(1)/E
0563 z2=ee(2)/E
0564 if(ish.ge.7)then
0565 write(ifch,*) 'z of first branching',z1
0566 endif
0567 pprt(3,ij(1))= sqrt(max(0.,ee(1)**2-pprt(5,ij(1))**2))
0568 pprt(3,ij(2))=-sqrt(max(0.,ee(2)**2-pprt(5,ij(2))**2))
0569 pt=0.
0570 alpha=0.
0571 else
0572 E=pprt(4,io)
0573 z1=pz(io)
0574 z2=1.-z1
0575 am0=pprt(5,io)
0576 am1=pprt(5,ij(1))
0577 am2=pprt(5,ij(2))
0578 aM=am2**2-am0**2-am1**2
0579 pzz=sqrt((E-am0)*(E+am0))
0580 pprt(3,ij(1))=.5*(aM+2.*z1*E**2)/pzz
0581 pprt(3,ij(2))=pzz-pprt(3,ij(1))
0582 pt2=(E**2*(z1*z2*am0**2-z1*am2**2-z2*am1**2)
0583 $ -.25*aM**2+am0**2*am1**2)/pzz**2
0584 if(ish.ge.8)then
0585 write(ifch,*) 'pt2,pzz=',pt2,pzz,z1**2*E*pprt(5,io),z1,E
0586 $ ,pprt(5,io)
0587 endif
0588 if(pt2.lt.0.) then
0589 ii=1
0590 if ( pprt(5,ij(1))**2-amm2(idprt(ij(1))).lt.
0591 $ pprt(5,ij(2))**2-amm2(idprt(ij(2))) ) ii=2
0592 goto 10
0593 endif
0594 pt=sqrt(pt2)
0595 alpha=2.*pi*rangen()
0596 endif
0597 pprt(4,ij(1))=z1*E
0598 pprt(1,ij(1))=cos(alpha)*pt
0599 pprt(2,ij(1))=sin(alpha)*pt
0600 if(ii2.eq.2)then
0601 pprt(4,ij(2))=z2*E
0602 pprt(1,ij(2))=-cos(alpha)*pt
0603 pprt(2,ij(2))=-sin(alpha)*pt
0604 endif
0605 if(ij(1).ne.j1)then
0606 if(pprt(1,io).ne.0..or.pprt(2,io).ne.0.)then
0607 do ii=1,2
0608 call utrota(-1,pprt(1,io),pprt(2,io),pprt(3,io)
0609 & ,pprt(1,ij(ii)),pprt(2,ij(ii)),pprt(3,ij(ii)))
0610 enddo
0611 endif
0612 endif
0613 if(ij(1).ne.j1)then
0614 if(pprt(3,io).lt.0.)then
0615 do k=1,3
0616 do ii=1,2
0617 pprt(k,ij(ii)) = -pprt(k,ij(ii))
0618 enddo
0619 enddo
0620 endif
0621 endif
0622 do ii=1,ii2
0623 if(id(1,ii).ne.0)then
0624 idprt(nprtj+1)=id(1,ii)
0625 idprt(nprtj+2)=id(2,ii)
0626 pprt(4,nprtj+1)=pz(ij(ii))*pprt(4,ij(ii))
0627 pprt(5,nprtj+1)=pz(ij(ii))*pprt(5,ij(ii))
0628 pprt(4,nprtj+2)=(1.-pz(ij(ii)))*pprt(4,ij(ii))
0629 pprt(5,nprtj+2)=(1.-pz(ij(ii)))*pprt(5,ij(ii))
0630 q2prt(nprtj+1)=q2s(ii)
0631 q2prt(nprtj+2)=q2s(ii)
0632 iorprt(nprtj+1)=ij(ii)
0633 iorprt(nprtj+2)=ij(ii)
0634 jorprt(nprtj+1)=jo(1,ii)
0635 jorprt(nprtj+2)=jo(2,ii)
0636 idaprt(1,ij(ii))=nprtj+1
0637 idaprt(2,ij(ii))=nprtj+2
0638 idaprt(1,nprtj+1)=0
0639 idaprt(2,nprtj+1)=0
0640 idaprt(1,nprtj+2)=0
0641 idaprt(2,nprtj+2)=0
0642 nprtj=nprtj+2
0643 else
0644 idaprt(1,ij(ii))=0
0645 idaprt(2,ij(ii))=0
0646 endif
0647 enddo
0648 if(ish.ge.5)then
0649 if(ij(1).ne.j1)then
0650 write(ifch,98) io,(pprt(j,io),j=1,5),pz(io),jorprt(io)
0651 & ,idprt(io)
0652 write(ifch,*) '->'
0653 endif
0654 write(ifch,99) ij(1),(pprt(j,ij(1)),j=1,5),pz(ij(1))
0655 & ,jo(1,1),jo(2,1),idprt(ij(1)),'->',id(1,1),id(2,1)
0656 if(ij(2).ne.0)write(ifch,99) ij(2),
0657 & (pprt(j,ij(2)),j=1,5),pz(ij(2))
0658 & ,jo(1,2),jo(2,2),idprt(ij(2)),'->',id(1,2),id(2,2)
0659 write(ifch,*)
0660 endif
0661
0662 98 format(i4,6g10.3,3i3)
0663 99 format(i4,6g10.3,3i3,a,2i4)
0664
0665
0666 if(ij(1).le.nn)ij(1)=nn-1
0667 ij(1)=ij(1)+2
0668 ij(2)=ij(1)+1
0669 ii=1
0670 ii2=2
0671 first=.true.
0672 if(ij(1).le.nprtj)goto 10
0673
0674 if(ish.ge.5)then
0675 write(ifch,*)
0676 do i=1,nprtj
0677 write(ifch,'(i4,1x,5g10.4,2i4,a,2i4)')
0678 & i,(pprt(j,i),j=1,5),idprt(i)
0679 & ,iorprt(i),' -->',idaprt(1,i),idaprt(2,i)
0680 enddo
0681 write(ifch,*)
0682 endif
0683 9999 call utprix('timsho',ish,ishini,5)
0684 return
0685 end
0686
