EPOS/epos167/epos-tim-155.f

   1 c---------------------------------------------------------------------
   2       subroutine timann
   3 c---------------------------------------------------------------------
   4 c      electron-positron
   5 c---------------------------------------------------------------------
   6       include 'epos.inc'
   7       include 'epos.incsem'
   8       common/nfla/nfla
   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.
  55       pprt(4,nprtj)=en
  56       pprt(5,nprtj)=en
  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
  70         s=engy**2
  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.
  82         vf=sign(.5,qf)-2.*qf*0.232
  83         af=sign(.5,qf)
  84         dsmax1=
  85      $       2.*(qf**2-2.*qf*ve*vf*chi1
  86      $       +(ae**2+ve**2)*(af**2+vf**2)*chi2)
  87      $       + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
  88
  89         qf=-1./3.
  90         vf=sign(.5,qf)-2.*qf*0.232
  91         af=sign(.5,qf)
  92         dsmax2=
  93      $       2.*(qf**2-2.*qf*ve*vf*chi1
  94      $       +(ae**2+ve**2)*(af**2+vf**2)*chi2)
  95      $       + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
  96
  97  100    iq1=1+INT(nfla*rangen())
  98         call idchrg(iq1,qf)
  99         ct=-1.+2.*rangen()
 100         vf=sign(.5,qf)-2.*qf*0.232
 101         af=sign(.5,qf)
 102
 103         dsigma=
 104      $       (1.+ct**2)*(qf**2-2.*qf*ve*vf*chi1
 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
 112
 113       nprtj=nprtj+1
 114       idprt(nprtj)=iq1
 115       pprt(4,nprtj)=en/2.
 116       pprt(5,nprtj)=en
 117       iorprt(nprtj)=1
 118
 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
 130       do i=1,nprtj
 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.
 135      &      (idprt(idaprt(1,i)).eq.9.and.jt.eq.2)))then
 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.
 139      &        (idprt(idaprt(1,i)).eq.9.and.jt.eq.1)))then
 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
 151
 152       if(ish.ge.5)then
 153         i=1
 154         do while(i.ne.0)
 155           if(idaprt(1,i) .eq. 0 ) then
 156             write(ifch,*) idprt(i)
 157             write(ifch,*) '|'
 158           endif
 159           i=jorprt(i)
 160         enddo
 161       endif
 162
 163       iptl=nptl
 164       i=1
 165       do while(i.gt.0)
 166         if(idaprt(1,i) .eq. 0) then
 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
 217           if(pprt(5,i).eq.0.)
 218      &    write(ifch,99) i,(pprt(j,i),j=1,5),idprt(i)
 219         enddo
 220         write(ifch,*)
 221         write(ifch,*)
 222       endif
 223
 224  99   format(i4,5g10.3,1i4)
 225  98   format(i4,5g10.3,5i4)
 226
 227       call utprix('timann',ish,ishini,5)
 228       return
 229       end
 230
 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.
 246       pprt(4,nprtj)=en
 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---------------------------------------------------------------------
 256       subroutine timsh2(q20,q21,en,idfla1,idfla2)
 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.
 271       pprt(4,nprtj)=en/2.
 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.
 288       pprtx(3,nprtjx)=en/2.
 289       pprtx(4,nprtjx)=en/2.
 290       pprtx(5,nprtjx)=0
 291       nprtjx=2
 292       pprtx(1,nprtjx)=0.
 293       pprtx(2,nprtjx)=0.
 294       pprtx(3,nprtjx)=-en/2.
 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
 351         zetamx=pprt(5,io)/pprt(4,io)/sqrt(pz(io)*(1.-pz(io)))
 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)
 359       ALFM=LOG(PT2MIN/qcdlam)
 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
 367           write(ifch,'(a,i4,i4,2f15.8)') 'null:',0.
 368         endif
 369         q2=amm2(idfl)
 370         idfla=0
 371         idflb=0
 372         goto 999
 373       endif
 374
 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
 379       IF(idfl.EQ.9) THEN
 380         FBR=6.*LOG((1.-ZC)/ZC)+nfla*(0.5-ZC)
 381       ELSE
 382         FBR=(8./3.)*LOG((1.-ZC)/ZC)
 383       endif
 384
 385       B0=(33.-2.*nfla)/6.
 386
 387
 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
 397           write(ifch,'(a,i4,i4,2f15.8)') 'null:',0.
 398         endif
 399         idfla=0
 400         idflb=0
 401         goto 999
 402       endif
 403
 404 c.....select flavor and z-value .....................................
 405       IF(idfl.EQ.9) THEN
 406         if(rangen()*FBR.lt.nfla*(0.5-ZC))then
 407                                 ! .................g -> qqbar
 408           Z=ZC+(1.-2.*ZC)*rangen()
 409           IF(Z**2+(1.-Z)**2.LT.rangen()) GOTO 390
 410           idfla=int(1.+rangen()*real(nfla))
 411           idflb=-idfla
 412         else                    !..................g -> gg
 413           Z=(1.-ZC)*(ZC/(1.-ZC))**rangen()
 414           IF(rangen().GT.0.5) Z=1.-Z
 415           IF((1.-Z*(1.-Z))**2.lt.rangen()) GOTO 390
 416           idfla=9
 417           idflb=9
 418         endif
 419       ELSE
 420         Z=1.-(1.-ZC)*(ZC/(1.-ZC))**rangen() !!........q -> qg
 421         IF(1.+Z**2.LT.2.*rangen()) GOTO 390
 422           idfla=idfl
 423           idflb=9
 424       endif
 425
 426       if(alfm/log(q2*z*(1.-z)/qcdlam).lt.rangen()) goto 390
 427
 428       if(ij(1).ne.j1.or.ij(2).eq.0)then
 429         if(E.lt.sqrt(q2))goto 390
 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......................................................................
 452
 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
 536  111    if(pt2.lt.0.) then
 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
 540           goto 10
 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
 553       if(ij(1).ne.j1.and.(pprt(1,io).ne.0..or.pprt(2,io).ne.0.))then
 554         do ii=1,2
 555           call utrota(-1,pprt(1,io),pprt(2,io),pprt(3,io)
 556      &    ,pprt(1,ij(ii)),pprt(2,ij(ii)),pprt(3,ij(ii)))
 557         enddo
 558       endif
 559       if(ij(1).ne.j1.and.pprt(3,io).lt.0.)then
 560         do k=1,3
 561           do ii=1,2
 562             pprt(k,ij(ii)) = -pprt(k,ij(ii))
 563           enddo
 564         enddo
 565       endif
 566       do ii=1,ii2
 567         if(id(1,ii).ne.0)then
 568           idprt(nprtj+1)=id(1,ii)
 569           idprt(nprtj+2)=id(2,ii)
 570           pprt(4,nprtj+1)=pz(ij(ii))*pprt(4,ij(ii))
 571           pprt(5,nprtj+1)=pz(ij(ii))*pprt(5,ij(ii))
 572           pprt(4,nprtj+2)=(1.-pz(ij(ii)))*pprt(4,ij(ii))
 573           pprt(5,nprtj+2)=(1.-pz(ij(ii)))*pprt(5,ij(ii))
 574           iorprt(nprtj+1)=ij(ii)
 575           iorprt(nprtj+2)=ij(ii)
 576           idaprt(1,ij(ii))=nprtj+1
 577           idaprt(2,ij(ii))=nprtj+2
 578           idaprt(1,nprtj+1)=0
 579           idaprt(2,nprtj+1)=0
 580           idaprt(1,nprtj+2)=0
 581           idaprt(2,nprtj+2)=0
 582           nprtj=nprtj+2
 583         else
 584           idaprt(1,ij(ii))=0
 585           idaprt(2,ij(ii))=0
 586         endif
 587       enddo
 588       if(ish.ge.5)then
 589         if(ij(1).ne.j1)then
 590           write(ifch,98) io,(pprt(j,io),j=1,5),pz(io)
 591      &    ,idprt(io)
 592           write(ifch,*) '->'
 593         endif
 594         write(ifch,99) ij(1),(pprt(j,ij(1)),j=1,5),pz(ij(1))
 595      &  ,idprt(ij(1)),'->',id(1,1),id(2,1)
 596         if(ij(2).ne.0)write(ifch,99) ij(2),
 597      &  (pprt(j,ij(2)),j=1,5),pz(ij(2))
 598      &  ,idprt(ij(2)),'->',id(1,2),id(2,2)
 599         write(ifch,*)
 600       endif
 601
 602  98   format(i4,6g10.3,1i4)
 603  99   format(i4,6g10.3,1i4,a,2i4)
 604
 605
 606       if(ij(1).le.nn)ij(1)=nn-1
 607       ij(1)=ij(1)+2
 608       ij(2)=ij(1)+1
 609       ii=1
 610       ii2=2
 611       first=.true.
 612       if(ij(1).le.nprtj)goto10
 613
 614       if(ish.ge.5)then
 615         write(ifch,*)
 616         do i=1,nprtj
 617           write(ifch,'(i4,1x,5g10.4,2i4,a,2i4)')
 618      &    i,(pprt(j,i),j=1,5),idprt(i)
 619      &    ,iorprt(i),'  -->',idaprt(1,i),idaprt(2,i)
 620         enddo
 621         write(ifch,*)
 622       endif
 623  9999 call utprix('timsho',ish,ishini,5)
 624       return
 625       end
 626