Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:29

0001 c---------------------------------------------------------------------
0002       subroutine timann
0003 c---------------------------------------------------------------------
0004 c      electron-positron
0005 c---------------------------------------------------------------------
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 c---------------------------------------------------------------------
0245       subroutine timsh1(q20,en,idfla,jo)
0246 c---------------------------------------------------------------------
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 c---------------------------------------------------------------------
0272       subroutine timsh2(q20,q21,en,idfla1,idfla2,jo1,jo2)
0273 c---------------------------------------------------------------------
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 c---------------------------------------------------------------------
0324       subroutine timsho(j1,j2)
0325 c---------------------------------------------------------------------
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 c      call timdev(idfl,q2start,E,zetamx,idfla,idflb,qa2,z)
0378 c......................................................................
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 c.....select new q2
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 c.....select flavor and z-value .....................................
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 c......................................................................
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 c          zc=.5*(1.-sqrt(1.-pprt(5,ij(ii))**2/ee(ii)**2))
0536 c          if(pz(ij(ii)).lt.zc.or.pz(ij(ii)).gt.1.-zc)then
0537 c            if(ish.ge.7)then
0538 c              write(ifch,*) 'first branching rejected'
0539 c            endif
0540 c            goto 10
0541 c          endif
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