C This file consists of subroutines that modify Pythia so as to make C it work with partly developed showers from beowulf. The code C is mostly by S. Mrenna. C C The file pythiamore.f conttains two subroutines that were removed C from Pythia to make pythialess.f: C C SUBROUTINE PYSHOW(IP1,IP2,QMAX) C SUBROUTINE PYADSH(NFIN) C C There are also a number of utility subroutines and functions. C C------------------------------------------------------------------------ C Modifications to Pythia to make it work with beowulf. S. Mrenna C Some lines commented out by DES, 7 April 2005. C More lines relating to nup0 = 6 or 7 commented out by DES, 16 May 2005 C IPASS common modified, 21 May 2005 DES C KTLIMIT used, 24 May 2005 DES C unused variables removed 3 September 2005, DES SUBROUTINE USRINIT(tmpstr) IMPLICIT NONE INTEGER HBSIZE, HPAW, INTUPLE PARAMETER (HBSIZE=500000) COMMON /PAWC/HPAW(HBSIZE) COMMON/MYANLC/CHNAME, TOPDIR,FRZNAM COMMON/MYANLI/LUNRZ INTEGER LUNRZ CHARACTER*6 CHNAME CHARACTER*8 TOPDIR CDES INTEGER IOERR CHARACTER*60 FRZNAM,tmpstr INTEGER NPASS,NTOT COMMON/NPASS/NPASS,NTOT integer idvect(200),nout common/ibook2/idvect,nout NTOT=0 NPASS=0 frznam=tmpstr if(frznam.eq.' ') FRZNAM='output.hb' c4now c DES CALL HLIMIT(HBSIZE) LUNRZ = 66 CHNAME = 'ONETOP' TOPDIR = CHNAME C OPEN NEW RZ FILE c DES CALL HROPEN(LUNRZ,CHNAME,FRZNAM,'N',1024,IOERR) c DES OPEN(LUNRZ,file=FRZNAM,STATUS='unknown') c IF(IOERR.NE.0) WRITE(*,*) 'HROPEN STATUS = ',IOERR C INTUPLE=0 C IF(INTUPLE.NE.0) THEN ELSE c DES call hbook1(1,'Njets (raw) ',10,-0.5,9.5,0.0) c DES call hbook1(2,' log10(y3) b4 pythia ',100,-5.0,0.0,0.0) c DES call hbook1(3,' log10(y4) b4 pythia ',100,-5.0,0.0,0.0) c DES call hbook1(14,' NN ',100,0.0,100.0,0.0) c DES call hbook1(15,' NUP ',20,0.0,20.0,0.0) c DES call hbook1(20,' max mass ',35,0.0,35.0,0.0) c DES call hbook1(21,' mid mass ',35,0.0,35.0,0.0) c DES call hbook1(22,' min mass ',35,0.0,35.0,0.0) c DES call hbook1(30,' max E ',50,0.0,50.0,0.0) c DES call hbook1(31,' mid E ',50,0.0,50.0,0.0) c DES call hbook1(32,' min E ',50,0.0,50.0,0.0) c DES call hbook1(40,' system mass ',100,0.0,100.0,0.0) c DES call hbook1(41,' Thrust ',50,0.6,1.1,0.0) c DES call hbook1(42,' logThrust ',50,-.3,.1,0.0) c DES call hbook1(43,' dmin ',60,0.0,30.0,0.0) C C......some "inclusive" quantities c DES call hbook1(60,' log10(y1) ',100,10.0,20.0,0.0) c DES call hbook1(61,' log10(y2) ',100,-5.0,0.0,0.0) c DES call hbook1(62,' log10(y3) ',100,-5.0,0.0,0.0) c DES call hbook1(63,' log10(y4) ',100,-5.0,0.0,0.0) c DES call hbook1(64,' log10(y5) ',100,-5.0,0.0,0.0) c DES call hbook1(65,' log10(y6) ',100,-5.0,0.0,0.0) c DES call hbook1(66,' log10(y7) ',100,-5.0,0.0,0.0) c DES call hbook1(67,' log10(y8) ',100,-5.0,0.0,0.0) c DES call hbook1(68,' log10(y9) ',100,-5.0,0.0,0.0) c DES call hbook1(69,' log10(y10) ',100,-5.0,0.0,0.0) c DES call hbook1(160,' log10(y1) ',100,10.0,20.0,0.0) c DES call hbook1(161,' log10(y2) ',100,-5.0,0.0,0.0) c DES call hbook1(162,' log10(y3) ',100,-5.0,0.0,0.0) c DES call hbook1(163,' log10(y4) ',100,-5.0,0.0,0.0) c DES call hbook1(164,' log10(y5) ',100,-5.0,0.0,0.0) c DES call hbook1(165,' log10(y6) ',100,-5.0,0.0,0.0) c DES call hbook1(166,' log10(y7) ',100,-5.0,0.0,0.0) c DES call hbook1(167,' log10(y8) ',100,-5.0,0.0,0.0) c DES call hbook1(168,' log10(y9) ',100,-5.0,0.0,0.0) c DES call hbook1(169,' log10(y10) ',100,-5.0,0.0,0.0) ENDIF RETURN END C---------------------------------------------------------------------- SUBROUTINE USREVNT(WT) IMPLICIT NONE C New standard event common REAL*8 PHEP,VHEP INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(4000),IDHEP(4000), & JMOHEP(2,4000),JDAHEP(2,4000),PHEP(5,4000),VHEP(4,4000) SAVE /HEPEVT/ C--GUP Event common block INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP, & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP), & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP), & SPINUP(MAXNUP) SAVE /HEPEUP/ integer mstu,mstj real*8 paru,parj COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) save /pydat1/ INTEGER NJETS CDES REAL*8 PT1,PT2,WT,XMZ,temp REAL*8 WT,temp CDES INTEGER IG1,IG2,ii INTEGER ii REAL WT4 CDES REAL DUMB,DUMX REAL DUMB CDES INTEGER I,J,IPROC,ILSP,I1,I2,IGRV,IO INTEGER I,J CDES INTEGER IDE,IDP,IME,IMP,ICUT,IDE0,IDP0 INTEGER ICUT CDES DOUBLE PRECISION PYP, PYR CDES INTEGER IS,IJ,IDJET(5),IJE,IM1,IM2,IP,IH,IN2 CDES INTEGER NPASS,NTOT,ITIME,IHIGGS,ITIME2,IKNT,IQ1,IQ2 INTEGER NPASS,NTOT,ITIME,ITIME2 CDES INTEGER IT1,IT2,IPS DATA ICUT/1/ DATA ITIME/0/,ITIME2/0/ COMMON/NPASS/NPASS,NTOT SAVE ITIME CDES real*8 ppair(4) CDES integer idabs,ncjet integer idabs CDES real*8 pcjet integer njmax PARAMETER (NJMAX=512) CDES dimension PCJET(4,NJMAX) CDES integer iphix,iyx CDES real*8 px integer n,npad,k real*8 p,v COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) SAVE /PYJETS/ INTEGER kMATCH CDES Note that there is a separate NMATCH defined here. CDES REAL*8 Q0RES,QXRES,QHARD,PASSKT(20),scaleq(20),thord(20) CDES COMMON/IPASS/Q0RES,QXRES,QHARD,PASSKT,scaleq,thord,kMATCH REAL*8 Q0RES,QXRES,QHARD,PASSKT(20),SCALEQ(20),THORD(20) REAL*8 KTLIMIT COMMON/IPASS/Q0RES,QXRES,QHARD,PASSKT,SCALEQ,THORD,KTLIMIT,kMATCH SAVE /IPASS/ C C.....My variables CDES real*8 ecut,y(njmax),pp(4,njmax),drmin,etmax real*8 ecut,y(njmax),pp(4,njmax) CDES integer nn,id,ilep integer nn CDES integer imatch,nmatch,iok(10),ijet,idl,idn integer idl,idn CDES real*8 cosj,cosi,etai,etaj,phii,phij,deta,dphi,dr,dely CDES real*8 p12(3),p34(3),parsave CDES real*8 a(3),thrustnew,thrustold,obl real*8 a(3) CDES real*8 scalkt,ptw real*8 scalkt integer NUPMAX common /scales/ scalkt,NUPMAX save /scales/ real*8 eccut real*8 total_cross_section common/statistics/total_cross_section CDES integer isuccess integer ifirst save ifirst INTEGER TYPE,ANGL,MONO,RECO COMMON/CLSTR/TYPE,ANGL,MONO,RECO SAVE /CLSTR/ INTEGER IKTTOPT COMMON/KTSWITCH/IKTTOPT SAVE /KTSWITCH/ data ifirst/0/ ifirst=ifirst+1 C/* Convert HEPEVT listing of 1st event to PYTHIA for analysis c print*,' ** getting ready for hadronization ** ' csm call myhepc(2) c call pylist(1) c print*,' ** calling pyexec ** ' c call pyexec c call lunhep(1) C call pylist(7) C call pylist(2) nevhep=ifirst total_cross_section = total_cross_section + wt C/* I still want to PLOT kt-clustered variables. ikttopt=0 wt4=sngl(wt) C.....M(1,2) at the parton level C eccut=0.1D0 298 continue idl=0 idn=0 NN=0 DO 299 I=3,N if(k(i,1).ne.1.and.k(i,1).ne.2) goto 299 idabs=abs(k(i,2)) if(idabs.eq.22.and.p(i,4).lt.0.0001D0) goto 299 if(sqrt(p(i,1)**2+p(i,2)**2).lt.eccut) GOTO 299 NN=NN+1 IF(NN.GT.512) THEN eccut=eccut*1.5D0 goto 298 ENDIF DO J=1,4 PP(J,NN)=P(I,J) ENDDO 299 CONTINUE c call hf1(14,sngl(nn),wt4) ! That's wrong. It should be "float." DES. IF(TYPE.EQ.1) THEN DO I=1,10 Y(I)=0D0 ENDDO IF(NN.GT.1) THEN ECUT=0D0 DO I=1,MAX(NN,12) Y(I)=0D0 ENDDO IKTTOPT=0 CALL KTCLUS(1111,PP,NN,ECUT,Y,*999) C DO II=1,min(NN,10) dumb=sngl(y(ii)) if(ii.eq.4 .and. y(3).LT..05D0) RETURN if(dumb.ne.0d0) then dumb=log10(dumb) c DES call hf1(59+ii,dumb,wt4) endif enddo IKTTOPT=1 CALL KTCLUS(1111,PP,NN,ECUT,Y,*999) DO II=1,min(NN,10) dumb=sngl(y(ii)) if(dumb.ne.0d0) then dumb=log10(dumb) c DES call hf1(159+ii,dumb,wt4) endif enddo IKTTOPT=0 ENDIF ENDIF njets=0 call pyclus(njets) dumb=njets c DES call hf1(1,dumb,wt4) c DES call hf1(40,sngl(paru(61)),wt4) c DES call hf1(41,sngl(paru(62)),wt4) c DES call hf1(42,sngl(log10(paru(62))),wt4) c DES call hf1(43,sngl(paru(63)),wt4) a(1)=p(n+1,5) a(2)=p(n+2,5) a(3)=p(n+3,5) i=2 if(a(3).gt.a(2)) i=3 j=2 if(a(1).lt.a(2)) j=1 if(i.gt.j) then temp=a(i) a(i)=a(j) a(j)=temp i=2 if(a(3).gt.a(2)) i=3 j=2 if(a(1).lt.a(2)) j=1 if(i.gt.j) then temp=a(i) a(i)=a(j) a(j)=temp endif endif if(a(1).lt.a(2).or.a(3).gt.a(2).or.a(3).gt.a(1)) then print*,' error in sorting masses ' endif c DES call hf1(20,sngl(a(1)),wt4) c DES call hf1(21,sngl(a(2)),wt4) c DES call hf1(22,sngl(a(3)),wt4) C c DES call hf1(30,sngl(p(n+1,4)),wt4) c DES call hf1(31,sngl(p(n+2,4)),wt4) c DES call hf1(32,sngl(p(n+3,4)),wt4) C..... 999 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE USRDONE(rnorm) IMPLICIT NONE COMMON/MYANLC/CHNAME, TOPDIR,FRZNAM COMMON/MYANLI/LUNRZ INTEGER LUNRZ CHARACTER*6 CHNAME CHARACTER*8 TOPDIR CHARACTER*60 FRZNAM CDES INTEGER ICYCLE INTEGER NPASS,NTOT COMMON/NPASS/NPASS,NTOT REAL*8 RNORM INTEGER I,IH integer idvect(200),nout common/ibook2/idvect,nout integer nevlast common /evpass/ nevlast C4now c DES CALL HID1(IDVECT,NOUT) RNORM=1D0/RNORM DO 100 I=1,NOUT IH=IDVECT(I) c DES CALL HOPERA(IH,'+',IH,IH,sngl(RNORM),0.0) 100 continue c4now CALL PYDUMP(1,LUNRZ,0,IH) C c CALL HPRNT(1) c CALL HCDIR('//PAWC', ' ') c CALL HCDIR(TOPDIR,' ') c DES CALL HROUT(0, ICYCLE, ' ') c print*,' ICYCLE = ',ICYCLE c CALL HLDIR(TOPDIR, ' ') c CALL HPRINT(100) c CALL HPRINT(200) c DES CALL HREND(TOPDIR) C CLOSE(LUNRZ) RETURN END C/* C23456789012345678901234567890123456789012345678901234567890123456789012x4567890 C Extra subroutines to help Pythia work with partly developed showers. C Contains subroutines PREP_EVENT, FULLCOPY(IIN,IOUT), C PYSHOW(IP1,IP2,QMAX), PYADSH(NFIN), PYPTFS(NPART,IPART,PTMAX,PTMIN,PTGEN) C C S. Mrenna, 29 June 2004 (with some print statements removed by DES). C Modified 14 August 2004, DES. C KTLIMIT added to common /IPASS/. C This is used to require kt < KTLIMIT in PREP_EVENT. C/* SUBROUTINE PREP_EVENT IMPLICIT NONE CDES INTEGER J,IJ,I,NLO,NHI INTEGER J,IJ,I INTEGER idabs CDES integer ic1,ic2,jc1,jc2,jdabs,iconn,nup0,idm integer ic1,ic2,jc1,jc2,jdabs,iconn,nup0 CDES REAL*8 thang,thtmp,Zi,rmmin,qold,FACT REAL*8 thang,thtmp CDES real*8 xmass C Common block for extra info. INTEGER NMATCH REAL*8 Q0RES,QXRES,QHARD,PASSKT(20),SCALEQ(20),THORD(20) REAL*8 KTLIMIT COMMON/IPASS/Q0RES,QXRES,QHARD,PASSKT,SCALEQ,THORD,KTLIMIT,NMATCH SAVE /IPASS/ C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ INTEGER IKTTOPT COMMON/KTSWITCH/IKTTOPT SAVE /KTSWITCH/ REAL*8 RESKT(0:10) COMMON/RES/RESKT SAVE /RES/ NUP0=NUP C/* qhard=pup(4,1)+pup(4,2) c$$$ rmmin=qhard c$$$ if(nup.ge.8) then c$$$ rmmin=min(rmmin,xmass(4,5)) c$$$ rmmin=min(rmmin,xmass(4,6)) c$$$ rmmin=min(rmmin,xmass(5,6)) c$$$ endif ij=101 DO I=1,NUP call fullcopy(i,ij) ij=ij+1 ENDDO C/* NUP0=NUP NUP=3 do ij=4,nup0 IF(istup(ij).eq.1.or.istup(ij).eq.7) THEN nup=nup+1 scaleq(nup)=scaleq(ij) scaleq(ij)=0D0 passkt(nup)=passkt(ij) passkt(ij)=0D0 thord(nup)=thord(ij) thord(ij)=0D0 call fullcopy(ij,nup) if(istup(nup).eq.1) then C------- End DES addition else scaleq(nup)=1D-6 ! small scale prevents showering for istup = 7 istup(nup)=1 ! now reset istup to 1. endif mothup(1,nup)=3 ! set mother to photon. mothup(2,nup)=3 c passkt(i)=qold*sqrt(Zi*(1D0-Zi)) ! limits kt subsequent emissions. C New DES addition IF (passkt(i).GT.KTLIMIT) then passkt(i) = KTLIMIT END IF C End new DES addition spinup(i)=0D0 ENDIF enddo C/* C/******* C/* Set the angles for angular ordering. DO I=4,NUP thord(i)=1D6 ENDDO DO 20 I=4,NUP IF(ABS(ISTUP(I)).NE.1) GOTO 20 IDABS=ABS(IDUP(I)) IF(IDABS.GE.11.AND.IDABS.LE.16) GOTO 20 IF(IDABS.GT.21) GOTO 20 IC1=ICOLUP(1,I) IC2=ICOLUP(2,I) C /* gluons have two color connections DO 30 J=4,NUP IF(J.EQ.I) GOTO 30 IF(ABS(ISTUP(J)).NE.1) GOTO 30 JDABS=ABS(IDUP(J)) IF(JDABS.GE.11.AND.JDABS.LE.16) GOTO 30 IF(JDABS.GT.21) GOTO 30 JC1=ICOLUP(1,J) JC2=ICOLUP(2,J) ICONN=0 C /* FSR IF(ISTUP(I)*ISTUP(J).GT.0) THEN IF(IC1.EQ.0.AND.IC2.NE.0.AND.JC1.EQ.IC2) THEN ICONN=1 ELSEIF(IC1.NE.0.AND.IC2.EQ.0.AND.JC2.EQ.IC1) THEN ICONN=1 ELSEIF(IC1*IC2.NE.0.AND.(JC1.EQ.IC2.OR.JC2.EQ.IC1)) THEN ICONN=1 ENDIF C /* ISR ELSE IF(IC2.EQ.0.AND.IC1.NE.0.AND.JC1.EQ.IC1) THEN ICONN=1 ELSEIF(IC2.NE.0.AND.IC1.EQ.0.AND.JC2.EQ.IC2) THEN ICONN=1 ELSEIF(IC1*IC2.NE.0.AND.(JC1.EQ.IC1.OR.JC2.EQ.IC2)) THEN ICONN=1 ENDIF ENDIF IF(ICONN.EQ.1) THEN thtmp=thang(pup(1,I),pup(1,J)) !calculate angle if(thtmp.lt.thord(i)) thord(i)=thtmp if(thtmp.lt.thord(j)) thord(j)=thtmp ENDIF 30 CONTINUE 20 CONTINUE C DES IF(NUP0.EQ.6 .OR. NUP0.EQ.7) THEN C DES NMATCH=NUP-5 C DES NLO=101 C DES IKTTOPT=1 C DES DO I=0,10 C DES RESKT(I)=0D0 C DES ENDDO C DES CALL CLUSTER(NLO,NHI,NMATCH,0) C DES if(.false.) then C DES print*,' nmatch = ',nmatch C DES do i=nlo,nhi C DES write(*,588) i,istup(i),idup(i),mothup(1,i),mothup(2,i), C DES $ icolup(1,i),icolup(2,i), C DES $ pup(5,i),pup(1,i) C DES enddo C DES 588 format(7I6,2(3X,E12.6)) C DES print*,'reskt',nup C DES print*,(reskt(i),i=0,min(nup,10)) C DES endif C DES CALL REWEIGHT(NLO,NHI,FACT,0,2) C DES IKTTOPT=0 C DES ENDIF RETURN END SUBROUTINE FULLCOPY(IIN,IOUT) IMPLICIT NONE INTEGER IIN,IOUT INTEGER J C/* COPY PARTON FROM ONE COMMON TO ANOTHER C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ IDUP(IOUT)=IDUP(IIN) ISTUP(IOUT)=ISTUP(IIN) MOTHUP(1,IOUT)=MOTHUP(1,IIN) MOTHUP(2,IOUT)=MOTHUP(2,IIN) ICOLUP(1,IOUT)=ICOLUP(1,IIN) ICOLUP(2,IOUT)=ICOLUP(2,IIN) DO J=1,5 PUP(J,IOUT)=PUP(J,IIN) ENDDO RETURN END C********************************************************************* C...PYSHOW C...Generates timelike parton showers from given partons. SUBROUTINE PYSHOW(IP1,IP2,QMAX) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) CDES INTEGER PYK,PYCHGE,PYCOMP INTEGER PYCOMP C...Parameter statement to help give large particle numbers. PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000, &KEXCIT=4000000,KDIMEN=5000000) C...Commonblocks. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) COMMON/PYINT1/MINT(400),VINT(400) SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYINT1/ C...Local arrays. DIMENSION PMTH(5,50),PS(5),PMA(19),PMSD(10),IEP(10),IPA(10), &KFLA(10),KFLD(10),KFL(10),ITRY(10),ISI(10),ISL(10),DP(10), &DPT(5,4),KSH(0:40),KCII(2),NIIS(2),IIIS(2,2),THEIIS(2,2), &PHIIIS(2,2),ISII(2),ISSET(10),ISCOL(0:40),ISCHG(0:40), &IREF(1000) C Common block for extra info. INTEGER NMATCH REAL*8 Q0RES,QXRES,QHARD,PASSKT(20),SCALEQ(20),THORD(20) REAL*8 KTLIMIT COMMON/IPASS/Q0RES,QXRES,QHARD,PASSKT,SCALEQ,THORD,KTLIMIT,NMATCH SAVE /IPASS/ PARAMETER (MAXNUP=500) DIMENSION IPART(MAXNUP) real*8 tempq(20) integer itime data itime/0/ save itime do i=1,20 tempq(i)=scaleq(i) enddo itime=itime+1 C...Check that QMAX not too low. IF(MSTJ(41).LE.0) THEN RETURN ELSEIF(MSTJ(41).EQ.1) THEN IF(QMAX.LE.PARJ(82).AND.IP2.GT.-8) RETURN ELSE IF(QMAX.LE.MIN(PARJ(82),PARJ(83),PARJ(90)).AND.IP2.GT.-8) & RETURN ENDIF C...Initialization of cutoff masses etc. DO 100 IFL=0,40 ISCOL(IFL)=0 ISCHG(IFL)=0 KSH(IFL)=0 100 CONTINUE ISCOL(21)=1 KSH(21)=1 PMTH(1,21)=PYMASS(21) PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25D0*PARJ(82)**2) PMTH(3,21)=2D0*PMTH(2,21) PMTH(4,21)=PMTH(3,21) PMTH(5,21)=PMTH(3,21) PMTH(1,22)=PYMASS(22) PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25D0*PARJ(83)**2) PMTH(3,22)=2D0*PMTH(2,22) PMTH(4,22)=PMTH(3,22) PMTH(5,22)=PMTH(3,22) PMQTH1=PARJ(82) IF(MSTJ(41).GE.2) PMQTH1=MIN(PARJ(82),PARJ(83)) PMQT1E=MIN(PMQTH1,PARJ(90)) PMQTH2=PMTH(2,21) IF(MSTJ(41).GE.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22)) PMQT2E=MIN(PMQTH2,0.5D0*PARJ(90)) DO 110 IFL=1,5 ISCOL(IFL)=1 IF(MSTJ(41).GE.2) ISCHG(IFL)=1 KSH(IFL)=1 PMTH(1,IFL)=PYMASS(IFL) PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PMQTH1**2) PMTH(3,IFL)=PMTH(2,IFL)+PMQTH2 PMTH(4,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(82)**2)+PMTH(2,21) PMTH(5,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(83)**2)+PMTH(2,22) 110 CONTINUE DO 120 IFL=11,15,2 IF(MSTJ(41).GE.2) ISCHG(IFL)=1 IF(MSTJ(41).GE.2) KSH(IFL)=1 PMTH(1,IFL)=PYMASS(IFL) PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(90)**2) PMTH(3,IFL)=PMTH(2,IFL)+0.5D0*PARJ(90) PMTH(4,IFL)=PMTH(3,IFL) PMTH(5,IFL)=PMTH(3,IFL) 120 CONTINUE PT2MIN=MAX(0.5D0*PARJ(82),1.1D0*PARJ(81))**2 ALAMS=PARJ(81)**2 ALFM=LOG(PT2MIN/ALAMS) C...Store positions of shower initiating partons. MPSPD=0 IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN NPA=1 IPA(1)=IP1 ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)- & MSTU(32))) THEN NPA=2 IPA(1)=IP1 IPA(2)=IP2 ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0 & .AND.IP2.GE.-7) THEN NPA=IABS(IP2) DO 130 I=1,NPA IPA(I)=IP1+I-1 130 CONTINUE ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND. &IP2.EQ.-8) THEN MPSPD=1 NPA=2 IPA(1)=IP1+6 IPA(2)=IP1+7 ELSE CALL PYERRM(12, & '(PYSHOW:) failed to reconstruct showering system') IF(MSTU(21).GE.1) RETURN ENDIF C...Send off to PYPTFS for pT-ordered evolution if requested, C...if at least 2 partons, and without predefined shower branchings. IF((MSTJ(41).EQ.11.OR.MSTJ(41).EQ.12).AND.NPA.GE.2.AND. &MPSPD.EQ.0) THEN NPART=NPA DO 135 II=1,NPART IPART(II)=IPA(II) 135 CONTINUE CALL PYPTFS(NPART,IPART,0.5D0*QMAX,0D0,PTGEN) RETURN ENDIF C...Check on phase space available for emission. IREJ=0 DO 140 J=1,5 PS(J)=0D0 140 CONTINUE PM=0D0 DO 160 I=1,NPA KFLA(I)=IABS(K(IPA(I),2)) PMA(I)=P(IPA(I),5) C...Special cutoff masses for initial partons (may be a heavy quark, C...squark, ..., and need not be on the mass shell). IR=30+I IF(NPA.LE.1) IREF(I)=IR IF(NPA.GE.2) IREF(I+1)=IR IF(KFLA(I).LE.8) THEN ISCOL(IR)=1 IF(MSTJ(41).GE.2) ISCHG(IR)=1 ELSEIF(KFLA(I).EQ.11.OR.KFLA(I).EQ.13.OR.KFLA(I).EQ.15.OR. & KFLA(I).EQ.17) THEN IF(MSTJ(41).GE.2) ISCHG(IR)=1 ELSEIF(KFLA(I).EQ.21) THEN ISCOL(IR)=1 ELSEIF((KFLA(I).GE.KSUSY1+1.AND.KFLA(I).LE.KSUSY1+8).OR. & (KFLA(I).GE.KSUSY2+1.AND.KFLA(I).LE.KSUSY2+8)) THEN ISCOL(IR)=1 ELSEIF(KFLA(I).EQ.KSUSY1+21) THEN ISCOL(IR)=1 ENDIF IF(ISCOL(IR).EQ.1.OR.ISCHG(IR).EQ.1) KSH(IR)=1 PMTH(1,IR)=PMA(I) IF(ISCOL(IR).EQ.1.AND.ISCHG(IR).EQ.1) THEN PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PMQTH1**2) PMTH(3,IR)=PMTH(2,IR)+PMQTH2 PMTH(4,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2)+PMTH(2,21) PMTH(5,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(83)**2)+PMTH(2,22) ELSEIF(ISCOL(IR).EQ.1) THEN PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2) PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(82) PMTH(4,IR)=PMTH(3,IR) PMTH(5,IR)=PMTH(3,IR) ELSEIF(ISCHG(IR).EQ.1) THEN PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(90)**2) PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(90) PMTH(4,IR)=PMTH(3,IR) PMTH(5,IR)=PMTH(3,IR) ENDIF IF(KSH(IR).EQ.1) PMA(I)=PMTH(3,IR) PM=PM+PMA(I) IF(KSH(IR).EQ.0.OR.PMA(I).GT.10D0*QMAX) IREJ=IREJ+1 DO 150 J=1,4 PS(J)=PS(J)+P(IPA(I),J) 150 CONTINUE 160 CONTINUE IF(IREJ.EQ.NPA.AND.IP2.GE.-7) RETURN PS(5)=SQRT(MAX(0D0,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2)) IF(NPA.EQ.1) PS(5)=PS(4) IF(PS(5).LE.PM+PMQT1E) RETURN C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0). KFSRCE=0 IF(IP2.LE.0) THEN ELSEIF(K(IP1,3).EQ.K(IP2,3).AND.K(IP1,3).GT.0) THEN KFSRCE=IABS(K(K(IP1,3),2)) ELSE IPAR1=MAX(1,K(IP1,3)) IPAR2=MAX(1,K(IP2,3)) IF(K(IPAR1,3).EQ.K(IPAR2,3).AND.K(IPAR1,3).GT.0) & KFSRCE=IABS(K(K(IPAR1,3),2)) ENDIF ITYPES=0 IF(KFSRCE.GE.1.AND.KFSRCE.LE.8) ITYPES=1 IF(KFSRCE.GE.KSUSY1+1.AND.KFSRCE.LE.KSUSY1+8) ITYPES=2 IF(KFSRCE.GE.KSUSY2+1.AND.KFSRCE.LE.KSUSY2+8) ITYPES=2 IF(KFSRCE.GE.21.AND.KFSRCE.LE.24) ITYPES=3 IF(KFSRCE.GE.32.AND.KFSRCE.LE.34) ITYPES=3 IF(KFSRCE.EQ.25.OR.(KFSRCE.GE.35.AND.KFSRCE.LE.37)) ITYPES=4 IF(KFSRCE.GE.KSUSY1+22.AND.KFSRCE.LE.KSUSY1+37) ITYPES=5 IF(KFSRCE.EQ.KSUSY1+21) ITYPES=6 C...Identify two primary showerers. ITYPE1=0 IF(KFLA(1).GE.1.AND.KFLA(1).LE.8) ITYPE1=1 IF(KFLA(1).GE.KSUSY1+1.AND.KFLA(1).LE.KSUSY1+8) ITYPE1=2 IF(KFLA(1).GE.KSUSY2+1.AND.KFLA(1).LE.KSUSY2+8) ITYPE1=2 IF(KFLA(1).GE.21.AND.KFLA(1).LE.24) ITYPE1=3 IF(KFLA(1).GE.32.AND.KFLA(1).LE.34) ITYPE1=3 IF(KFLA(1).EQ.25.OR.(KFLA(1).GE.35.AND.KFLA(1).LE.37)) ITYPE1=4 IF(KFLA(1).GE.KSUSY1+22.AND.KFLA(1).LE.KSUSY1+37) ITYPE1=5 IF(KFLA(1).EQ.KSUSY1+21) ITYPE1=6 ITYPE2=0 IF(KFLA(2).GE.1.AND.KFLA(2).LE.8) ITYPE2=1 IF(KFLA(2).GE.KSUSY1+1.AND.KFLA(2).LE.KSUSY1+8) ITYPE2=2 IF(KFLA(2).GE.KSUSY2+1.AND.KFLA(2).LE.KSUSY2+8) ITYPE2=2 IF(KFLA(2).GE.21.AND.KFLA(2).LE.24) ITYPE2=3 IF(KFLA(2).GE.32.AND.KFLA(2).LE.34) ITYPE2=3 IF(KFLA(2).EQ.25.OR.(KFLA(2).GE.35.AND.KFLA(2).LE.37)) ITYPE2=4 IF(KFLA(2).GE.KSUSY1+22.AND.KFLA(2).LE.KSUSY1+37) ITYPE2=5 IF(KFLA(2).EQ.KSUSY1+21) ITYPE2=6 C...Order of showerers. Presence of gluino. ITYPMN=MIN(ITYPE1,ITYPE2) ITYPMX=MAX(ITYPE1,ITYPE2) IORD=1 IF(ITYPE1.GT.ITYPE2) IORD=2 IGLUI=0 IF(ITYPE1.EQ.6.OR.ITYPE2.EQ.6) IGLUI=1 C...Check if 3-jet matrix elements to be used. M3JC=0 ALPHA=0.5D0 IF(NPA.EQ.2.AND.MSTJ(47).GE.1.AND.MPSPD.EQ.0) THEN IF(MSTJ(38).NE.0) THEN M3JC=MSTJ(38) ALPHA=PARJ(80) MSTJ(38)=0 ELSEIF(MSTJ(47).GE.6) THEN M3JC=MSTJ(47) ELSE ICLASS=1 ICOMBI=4 C...Vector/axial vector -> q + qbar; q -> q + V. IF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.3)) THEN ICLASS=2 IF(KFSRCE.EQ.21.OR.KFSRCE.EQ.22) THEN ICOMBI=1 ELSEIF(KFSRCE.EQ.23.OR.(KFSRCE.EQ.0.AND. & K(IP1,2)+K(IP2,2).EQ.0)) THEN C...gamma*/Z0: assume e+e- initial state if unknown. EI=-1D0 IF(KFSRCE.EQ.23) THEN IANNFL=K(K(IP1,3),3) IF(IANNFL.NE.0) THEN KANNFL=IABS(K(IANNFL,2)) IF(KANNFL.GE.1.AND.KANNFL.LE.18) EI=KCHG(KANNFL,1)/3D0 ENDIF ENDIF AI=SIGN(1D0,EI+0.1D0) VI=AI-4D0*EI*PARU(102) EF=KCHG(KFLA(1),1)/3D0 AF=SIGN(1D0,EF+0.1D0) VF=AF-4D0*EF*PARU(102) XWC=1D0/(16D0*PARU(102)*(1D0-PARU(102))) SH=PS(5)**2 SQMZ=PMAS(23,1)**2 SQWZ=PS(5)*PMAS(23,2) SBWZ=1D0/((SH-SQMZ)**2+SQWZ**2) VECT=EI**2*EF**2+2D0*EI*VI*EF*VF*XWC*SH*(SH-SQMZ)*SBWZ+ & (VI**2+AI**2)*VF**2*XWC**2*SH**2*SBWZ AXIV=(VI**2+AI**2)*AF**2*XWC**2*SH**2*SBWZ ICOMBI=3 ALPHA=VECT/(VECT+AXIV) ELSEIF(KFSRCE.EQ.24.OR.KFSRCE.EQ.0) THEN ICOMBI=4 ENDIF C...For chi -> chi q qbar, use V/A -> q qbar as first approximation. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.5) THEN ICLASS=2 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=3 C...Scalar/pseudoscalar -> q + qbar; q -> q + S. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.4) THEN ICLASS=4 IF(KFSRCE.EQ.25.OR.KFSRCE.EQ.35.OR.KFSRCE.EQ.37) THEN ICOMBI=1 ELSEIF(KFSRCE.EQ.36) THEN ICOMBI=2 ENDIF ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=5 C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S. ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.3)) THEN ICLASS=6 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=7 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.4) THEN ICLASS=8 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=9 C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.5)) THEN ICLASS=10 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=11 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=12 C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.6) THEN ICLASS=13 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=14 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=15 C...g -> ~g + ~g (eikonal approximation). ELSEIF(ITYPMN.EQ.6.AND.ITYPMX.EQ.6.AND.ITYPES.EQ.0) THEN ICLASS=16 ENDIF M3JC=5*ICLASS+ICOMBI ENDIF ENDIF C...Find if interference with initial state partons. MIIS=0 IF(MSTJ(50).GE.1.AND.MSTJ(50).LE.3.AND.NPA.EQ.2.AND.KFSRCE.EQ.0 &.AND.MPSPD.EQ.0) MIIS=MSTJ(50) IF(MSTJ(50).GE.4.AND.MSTJ(50).LE.6.AND.NPA.EQ.2.AND.MPSPD.EQ.0) &MIIS=MSTJ(50)-3 IF(MIIS.NE.0) THEN DO 180 I=1,2 KCII(I)=0 KCA=PYCOMP(KFLA(I)) IF(KCA.NE.0) KCII(I)=KCHG(KCA,2)*ISIGN(1,K(IPA(I),2)) NIIS(I)=0 IF(KCII(I).NE.0) THEN DO 170 J=1,2 ICSI=MOD(K(IPA(I),3+J)/MSTU(5),MSTU(5)) IF(ICSI.GT.0.AND.ICSI.NE.IPA(1).AND.ICSI.NE.IPA(2).AND. & (KCII(I).EQ.(-1)**(J+1).OR.KCII(I).EQ.2)) THEN NIIS(I)=NIIS(I)+1 IIIS(I,NIIS(I))=ICSI ENDIF 170 CONTINUE ENDIF 180 CONTINUE IF(NIIS(1)+NIIS(2).EQ.0) MIIS=0 ENDIF C...Boost interfering initial partons to rest frame C...and reconstruct their polar and azimuthal angles. IF(MIIS.NE.0) THEN DO 200 I=1,2 DO 190 J=1,5 K(N+I,J)=K(IPA(I),J) P(N+I,J)=P(IPA(I),J) V(N+I,J)=0D0 190 CONTINUE 200 CONTINUE DO 220 I=3,2+NIIS(1) DO 210 J=1,5 K(N+I,J)=K(IIIS(1,I-2),J) P(N+I,J)=P(IIIS(1,I-2),J) V(N+I,J)=0D0 210 CONTINUE 220 CONTINUE DO 240 I=3+NIIS(1),2+NIIS(1)+NIIS(2) DO 230 J=1,5 K(N+I,J)=K(IIIS(2,I-2-NIIS(1)),J) P(N+I,J)=P(IIIS(2,I-2-NIIS(1)),J) V(N+I,J)=0D0 230 CONTINUE 240 CONTINUE CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,0D0,-PS(1)/PS(4), & -PS(2)/PS(4),-PS(3)/PS(4)) PHI=PYANGL(P(N+1,1),P(N+1,2)) CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,-PHI,0D0,0D0,0D0) THE=PYANGL(P(N+1,3),P(N+1,1)) CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),-THE,0D0,0D0,0D0,0D0) DO 250 I=3,2+NIIS(1) THEIIS(1,I-2)=PYANGL(P(N+I,3),SQRT(P(N+I,1)**2+P(N+I,2)**2)) PHIIIS(1,I-2)=PYANGL(P(N+I,1),P(N+I,2)) 250 CONTINUE DO 260 I=3+NIIS(1),2+NIIS(1)+NIIS(2) THEIIS(2,I-2-NIIS(1))=PARU(1)-PYANGL(P(N+I,3), & SQRT(P(N+I,1)**2+P(N+I,2)**2)) PHIIIS(2,I-2-NIIS(1))=PYANGL(P(N+I,1),P(N+I,2)) 260 CONTINUE ENDIF C...Boost 3 or more partons to their rest frame. IF(NPA.GE.3) CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4), &-PS(2)/PS(4),-PS(3)/PS(4)) C...Define imagined single initiator of shower for parton system. NS=N IF(N.GT.MSTU(4)-MSTU(32)-10) THEN CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS') IF(MSTU(21).GE.1) RETURN ENDIF 270 N=NS IF(NPA.GE.2) THEN K(N+1,1)=11 K(N+1,2)=21 K(N+1,3)=0 K(N+1,4)=0 K(N+1,5)=0 P(N+1,1)=0D0 P(N+1,2)=0D0 P(N+1,3)=0D0 P(N+1,4)=PS(5) P(N+1,5)=PS(5) V(N+1,5)=PS(5)**2 N=N+1 IREF(1)=21 ENDIF C...Loop over partons that may branch. NEP=NPA IM=NS IF(NPA.EQ.1) IM=NS-1 280 IM=IM+1 IF(N.GT.NS) THEN IF(IM.GT.N) GOTO 590 KFLM=IABS(K(IM,2)) IR=IREF(IM-NS) IF(KSH(IR).EQ.0) GOTO 280 IF(P(IM,5).LT.PMTH(2,IR)) GOTO 280 IGM=K(IM,3) ELSE IGM=-1 ENDIF IF(N+NEP.GT.MSTU(4)-MSTU(32)-10) THEN CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS') IF(MSTU(21).GE.1) RETURN ENDIF C...Position of aunt (sister to branching parton). C...Origin and flavour of daughters. IAU=0 IF(IGM.GT.0) THEN IF(K(IM-1,3).EQ.IGM) IAU=IM-1 IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1 ENDIF IF(IGM.GE.0) THEN K(IM,4)=N+1 DO 290 I=1,NEP K(N+I,3)=IM 290 CONTINUE ELSE K(N+1,3)=IPA(1) ENDIF IF(IGM.LE.0) THEN DO 300 I=1,NEP K(N+I,2)=K(IPA(I),2) 300 CONTINUE ELSEIF(KFLM.NE.21) THEN K(N+1,2)=K(IM,2) K(N+2,2)=K(IM,5) IREF(N+1-NS)=IREF(IM-NS) IREF(N+2-NS)=IABS(K(N+2,2)) ELSEIF(K(IM,5).EQ.21) THEN K(N+1,2)=21 K(N+2,2)=21 IREF(N+1-NS)=21 IREF(N+2-NS)=21 ELSE K(N+1,2)=K(IM,5) K(N+2,2)=-K(IM,5) IREF(N+1-NS)=IABS(K(N+1,2)) IREF(N+2-NS)=IABS(K(N+2,2)) ENDIF C...Reset flags on daughters and tries made. DO 310 IP=1,NEP K(N+IP,1)=3 K(N+IP,4)=0 K(N+IP,5)=0 KFLD(IP)=IABS(K(N+IP,2)) IF(KCHG(PYCOMP(KFLD(IP)),2).EQ.0) K(N+IP,1)=1 ITRY(IP)=0 ISL(IP)=0 ISI(IP)=0 IF(KSH(IREF(N+IP-NS)).EQ.1) ISI(IP)=1 310 CONTINUE ISLM=0 C...Maximum virtuality of daughters. IF(IGM.LE.0) THEN DO 320 I=1,NPA IF(NPA.GE.3) P(N+I,4)=P(IPA(I),4) P(N+I,5)=MIN(QMAX,PS(5)) IR=IREF(N+I-NS) IF(IP2.LE.-8) P(N+I,5)=MAX(P(N+I,5),2D0*PMTH(3,IR)) IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5) C.....Select upper scale for parton showers IF(MINT(400).EQ.1.and.i.le.17) THEN C/* if(itime.le.50) print*,tempq c DES if(itime.le.50) print*,' i ',i,k(n+i,2),k(n+i,1) C/* P(N+I,5)=TEMPQ(I+2) IF(ABS(K(N+I,2)).LE.6.or.K(N+I,2).eq.21) THEN P(N+I,5)=TEMPQ(I+3) c if(itime.le.50) print*,' max virt ',n+i,k(N+i,2),p(n+i,5) ENDIF c DES if(itime.le.50) print*,' max virt ',n+i,k(N+i,2),p(n+i,5) ENDIF 320 CONTINUE ELSE IF(MSTJ(43).LE.2) PEM=V(IM,2) IF(MSTJ(43).GE.3) PEM=P(IM,4) P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM) P(N+2,5)=MIN(P(IM,5),(1D0-V(IM,1))*PEM) IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22) ENDIF DO 330 I=1,NEP PMSD(I)=P(N+I,5) IF(ISI(I).EQ.1) THEN IR=IREF(N+I-NS) IF(P(N+I,5).LE.PMTH(3,IR)) P(N+I,5)=PMTH(1,IR) ENDIF V(N+I,5)=P(N+I,5)**2 330 CONTINUE C...Choose one of the daughters for evolution. 340 INUM=0 IF(NEP.EQ.1) INUM=1 DO 350 I=1,NEP IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I 350 CONTINUE DO 360 I=1,NEP IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN IR=IREF(N+I-NS) IF(P(N+I,5).GE.PMTH(2,IR)) INUM=I ENDIF 360 CONTINUE IF(INUM.EQ.0) THEN RMAX=0D0 DO 370 I=1,NEP IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQT2E) THEN RPM=P(N+I,5)/PMSD(I) IR=IREF(N+I-NS) IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,IR)) THEN RMAX=RPM INUM=I ENDIF ENDIF 370 CONTINUE ENDIF C...Cancel choice of predetermined daughter already treated. INUM=MAX(1,INUM) INUMT=INUM IF(MPSPD.EQ.1.AND.IGM.EQ.0.AND.ITRY(INUMT).GE.1) THEN IF(K(IP1-1+INUM,4).GT.0) INUM=3-INUM ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2.AND.ITRY(INUMT).GE.1) THEN IF(KFLD(INUMT).NE.21.AND.K(IP1+2,4).GT.0) INUM=3-INUM IF(KFLD(INUMT).EQ.21.AND.K(IP1+3,4).GT.0) INUM=3-INUM ENDIF C...Store information on choice of evolving daughter. IEP(1)=N+INUM DO 380 I=2,NEP IEP(I)=IEP(I-1)+1 IF(IEP(I).GT.N+NEP) IEP(I)=N+1 380 CONTINUE DO 390 I=1,NEP KFL(I)=IABS(K(IEP(I),2)) 390 CONTINUE ITRY(INUM)=ITRY(INUM)+1 IF(ITRY(INUM).GT.200) THEN CALL PYERRM(14,'(PYSHOW:) caught in infinite loop') IF(MSTU(21).GE.1) RETURN ENDIF Z=0.5D0 IR=IREF(IEP(1)-NS) IF(KSH(IR).EQ.0) GOTO 440 IF(P(IEP(1),5).LT.PMTH(2,IR)) GOTO 440 C...Check if evolution already predetermined for daughter. IPSPD=0 IF(MPSPD.EQ.1.AND.IGM.EQ.0) THEN IF(K(IP1-1+INUM,4).GT.0) IPSPD=IP1-1+INUM ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2) THEN IF(KFL(1).NE.21.AND.K(IP1+2,4).GT.0) IPSPD=IP1+2 IF(KFL(1).EQ.21.AND.K(IP1+3,4).GT.0) IPSPD=IP1+3 ENDIF ISSET(INUM)=0 IF(IPSPD.NE.0) ISSET(INUM)=1 C...Select side for interference with initial state partons. IF(MIIS.GE.1.AND.IEP(1).LE.NS+3) THEN III=IEP(1)-NS-1 ISII(III)=0 IF(IABS(KCII(III)).EQ.1.AND.NIIS(III).EQ.1) THEN ISII(III)=1 ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.1) THEN IF(PYR(0).GT.0.5D0) ISII(III)=1 ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.2) THEN ISII(III)=1 IF(PYR(0).GT.0.5D0) ISII(III)=2 ENDIF ENDIF C...Calculate allowed z range. IF(NEP.EQ.1) THEN PMED=PS(4) ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN PMED=P(IM,5) ELSE IF(INUM.EQ.1) PMED=V(IM,1)*PEM IF(INUM.EQ.2) PMED=(1D0-V(IM,1))*PEM ENDIF IF(MOD(MSTJ(43),2).EQ.1) THEN ZC=PMTH(2,21)/PMED ZCE=PMTH(2,22)/PMED IF(ISCOL(IR).EQ.0) ZCE=0.5D0*PARJ(90)/PMED ELSE ZC=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTH(2,21)/PMED)**2))) IF(ZC.LT.1D-6) ZC=(PMTH(2,21)/PMED)**2 PMTMPE=PMTH(2,22) IF(ISCOL(IR).EQ.0) PMTMPE=0.5D0*PARJ(90) ZCE=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTMPE/PMED)**2))) IF(ZCE.LT.1D-6) ZCE=(PMTMPE/PMED)**2 ENDIF ZC=MIN(ZC,0.491D0) ZCE=MIN(ZCE,0.49991D0) IF(((MSTJ(41).EQ.1.AND.ZC.GT.0.49D0).OR.(MSTJ(41).GE.2.AND. &MIN(ZC,ZCE).GT.0.4999D0)).AND.IPSPD.EQ.0) THEN P(IEP(1),5)=PMTH(1,IR) V(IEP(1),5)=P(IEP(1),5)**2 GOTO 440 ENDIF C...Integral of Altarelli-Parisi z kernel for QCD. C...(Includes squark and gluino; with factor N_C/C_F extra for latter). IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN FBR=6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0 ELSEIF(MSTJ(49).EQ.0) THEN FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC) IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0) C...Integral of Altarelli-Parisi z kernel for scalar gluon. ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1D0-2D0*ZC) ELSEIF(MSTJ(49).EQ.1) THEN FBR=(1D0-2D0*ZC)/3D0 IF(IGM.EQ.0.AND.M3JC.GE.1) FBR=4D0*FBR C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon. ELSEIF(KFL(1).EQ.21) THEN FBR=6D0*MSTJ(45)*(0.5D0-ZC) ELSE FBR=2D0*LOG((1D0-ZC)/ZC) ENDIF C...Reset QCD probability for colourless. IF(ISCOL(IR).EQ.0) FBR=0D0 C...Integral of Altarelli-Parisi kernel for photon emission. FBRE=0D0 IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1) THEN IF(KFL(1).LE.18) THEN FBRE=(KCHG(KFL(1),1)/3D0)**2*2D0*LOG((1D0-ZCE)/ZCE) ENDIF IF(MSTJ(41).EQ.10) FBRE=PARJ(84)*FBRE ENDIF C/* C.....Select which partons to veto IF(MINT(400).EQ.1.and.IEP(1).LE.NS+NEP+1) THEN IARG=IEP(1)-N IF(IARG.LE.17) THEN VETOKT=PASSKT(IARG+3) VETOTH=THORD(IARG+3) ENDIF C.....This should correspond to no veto ELSE C/* IF(VETOKT.LE.0D0) VETOKT=QMAX VETOKT=QMAX VETOTH=10.0 ENDIF c DES if(itime.lt.50.and.mint(400).eq.1) then c DES write(*,111) ' vetokt = ',vetokt,vetoth,iarg,iep(1),ns+nep+1 c DES 111 format(A10,1X,2G12.6,3I6) c DES endif C*/ C...Inner veto algorithm starts. Find maximum mass for evolution. 400 PMS=V(IEP(1),5) IF(IGM.GE.0) THEN PM2=0D0 DO 410 I=2,NEP PM=P(IEP(I),5) IRI=IREF(IEP(I)-NS) IF(KSH(IRI).EQ.1) PM=PMTH(2,IRI) PM2=PM2+PM 410 CONTINUE PMS=MIN(PMS,(P(IM,5)-PM2)**2) ENDIF C...Select mass for daughter in QCD evolution. B0=27D0/6D0 B1=(153D0-19D0*3D0)/6D0 DO 420 IFF=4,MSTJ(45) IF(PMS.GT.4D0*PMTH(2,IFF)**2) B0=(33D0-2D0*IFF)/6D0 IF(PMS.GT.4D0*PMTH(2,IFF)**2) B1=(153D0-19D0*IFF)/6D0 420 CONTINUE B1=0D0 C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2. PMSC=MAX(0.5D0*PARJ(82),PMS-PMTH(1,IR)**2) C...Already predetermined choice. IF(IPSPD.NE.0) THEN PMSQCD=P(IPSPD,5)**2 ELSEIF(FBR.LT.1D-3) THEN PMSQCD=0D0 ELSEIF(MSTJ(44).LE.0) THEN PMSQCD=PMSC*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/(PARU(111)*FBR))) ELSEIF(MSTJ(44).EQ.1) THEN PMSQCD=4D0*ALAMS*(0.25D0*PMSC/ALAMS)**(PYR(0)**(B0/FBR)) ELSE ccccc PMSQCD=PMSC*EXP(MAX(-50D0,ALFM*B0*LOG(PYR(0))/FBR)) ATEMP=ALFM/(1D0-B1*LOG(ALFM)/(B0**2*ALFM)) PMSQCD=PMSC*EXP(MAX(-50D0,ATEMP*B0*LOG(PYR(0))/FBR)) ENDIF C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2. IF(IPSPD.EQ.0) PMSQCD=PMSQCD+PMTH(1,IR)**2 IF(ZC.GT.0.49D0.OR.PMSQCD.LE.PMTH(4,IR)**2) PMSQCD=PMTH(2,IR)**2 V(IEP(1),5)=PMSQCD MCE=1 C...Select mass for daughter in QED evolution. IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1.AND.IPSPD.EQ.0) THEN C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2. PMSE=MAX(0.5D0*PARJ(83),PMS-PMTH(1,IR)**2) IF(FBRE.LT.1D-3) THEN PMSQED=0D0 ELSE PMSQED=PMSE*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/ & (PARU(101)*FBRE))) ENDIF C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2. PMSQED=PMSQED+PMTH(1,IR)**2 IF(ZCE.GT.0.4999D0.OR.PMSQED.LE.PMTH(5,IR)**2) PMSQED= & PMTH(2,IR)**2 IF(PMSQED.GT.PMSQCD) THEN V(IEP(1),5)=PMSQED MCE=2 ENDIF ENDIF C...Check whether daughter mass below cutoff. P(IEP(1),5)=SQRT(V(IEP(1),5)) IF(P(IEP(1),5).LE.PMTH(3,IR)) THEN P(IEP(1),5)=PMTH(1,IR) V(IEP(1),5)=P(IEP(1),5)**2 GOTO 440 ENDIF C...Already predetermined choice of z, and flavour in g -> qqbar. IF(IPSPD.NE.0) THEN IPSGD1=K(IPSPD,4) IPSGD2=K(IPSPD,5) PMSGD1=P(IPSGD1,5)**2 PMSGD2=P(IPSGD2,5)**2 ALAMPS=SQRT(MAX(1D-10,(PMSQCD-PMSGD1-PMSGD2)**2- & 4D0*PMSGD1*PMSGD2)) Z=0.5D0*(PMSQCD*(2D0*P(IPSGD1,4)/P(IPSPD,4)-1D0)+ALAMPS- & PMSGD1+PMSGD2)/ALAMPS Z=MAX(0.00001D0,MIN(0.99999D0,Z)) IF(KFL(1).NE.21) THEN K(IEP(1),5)=21 ELSE K(IEP(1),5)=IABS(K(IPSGD1,2)) ENDIF C...Select z value of branching: q -> qgamma. ELSEIF(MCE.EQ.2) THEN Z=1D0-(1D0-ZCE)*(ZCE/(1D0-ZCE))**PYR(0) IF(1D0+Z**2.LT.2D0*PYR(0)) GOTO 400 K(IEP(1),5)=22 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar. ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN Z=1D0-(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0) C...Only do z weighting when no ME correction afterwards. IF(M3JC.EQ.0.AND.1D0+Z**2.LT.2D0*PYR(0)) GOTO 400 K(IEP(1),5)=21 ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*0.5D0.LT.PYR(0)*FBR) THEN Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0) IF(PYR(0).GT.0.5D0) Z=1D0-Z IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 400 K(IEP(1),5)=21 ELSEIF(MSTJ(49).NE.1) THEN Z=PYR(0) IF(Z**2+(1D0-Z)**2.LT.PYR(0)) GOTO 400 KFLB=1+INT(MSTJ(45)*PYR(0)) PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5) IF(PMQ.GE.1D0) GOTO 400 IF(MSTJ(44).LE.2.OR.MSTJ(44).EQ.4) THEN IF(Z.LT.ZC.OR.Z.GT.1D0-ZC) GOTO 400 PMQ0=4D0*PMTH(2,21)**2/V(IEP(1),5) IF(MOD(MSTJ(43),2).EQ.0.AND.(1D0+0.5D0*PMQ)*SQRT(1D0-PMQ) & .LT.PYR(0)*(1D0+0.5D0*PMQ0)*SQRT(1D0-PMQ0)) GOTO 400 ELSE IF((1D0+0.5D0*PMQ)*SQRT(1D0-PMQ).LT.PYR(0)) GOTO 400 ENDIF K(IEP(1),5)=KFLB C...Ditto for scalar gluon model. ELSEIF(KFL(1).NE.21) THEN Z=1D0-SQRT(ZC**2+PYR(0)*(1D0-2D0*ZC)) K(IEP(1),5)=21 ELSEIF(PYR(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN Z=ZC+(1D0-2D0*ZC)*PYR(0) K(IEP(1),5)=21 ELSE Z=ZC+(1D0-2D0*ZC)*PYR(0) KFLB=1+INT(MSTJ(45)*PYR(0)) PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5) IF(PMQ.GE.1D0) GOTO 400 K(IEP(1),5)=KFLB ENDIF C...Check if z consistent with chosen m. IF(KFL(1).EQ.21) THEN IRGD1=IABS(K(IEP(1),5)) IRGD2=IRGD1 ELSE IRGD1=IR IRGD2=IABS(K(IEP(1),5)) ENDIF IF(NEP.EQ.1) THEN PED=PS(4) ELSEIF(NEP.GE.3) THEN PED=P(IEP(1),4) ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN PED=0.5D0*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5) ELSE IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM IF(IEP(1).EQ.N+2) PED=(1D0-V(IM,1))*PEM ENDIF IF(MOD(MSTJ(43),2).EQ.1) THEN PMQTH3=0.5D0*PARJ(82) IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83) IF(IRGD2.EQ.22.AND.ISCOL(IR).EQ.0) PMQTH3=0.5D0*PARJ(90) PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(IEP(1),5) PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(IEP(1),5) ZD=SQRT(MAX(0D0,(1D0-V(IEP(1),5)/PED**2)*((1D0-PMQ1-PMQ2)**2- & 4D0*PMQ1*PMQ2))) ZH=1D0+PMQ1-PMQ2 ELSE ZD=SQRT(MAX(0D0,1D0-V(IEP(1),5)/PED**2)) ZH=1D0 ENDIF c DES if(itime.lt.50) then c DES write(*,116) ' im, ped ',im,iep(1),ns,nep,ped,v(iep(1),5),v(im,5) c DES 116 format(A9,4I6,3G12.6) c DES endif IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND. &(MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN ELSEIF(IPSPD.NE.0) THEN ELSE ZL=0.5D0*(ZH-ZD) ZU=0.5D0*(ZH+ZD) IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 400 ENDIF C...Correct to alpha_s(pT^2) (optionally m^2/4 for g -> q qbar). IF(MCE.EQ.1.AND.MSTJ(44).GE.2.AND.IPSPD.EQ.0) THEN IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND. & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN IF(ALFM/LOG(V(IEP(1),5)*0.25D0/ALAMS).LT.PYR(0)) GOTO 400 ELSE PT2APP=Z*(1D0-Z)*V(IEP(1),5) C/* IF( MIN(Z,ZI)*V(IEP(1),5).GT. C/* $ VETOKT**2+GUESSM2*(Z+ZI) ) C/* $ GOTO 400 IF(MSTJ(44).GE.4) PT2APP=PT2APP* & (1D0-PMTH(1,IR)**2/V(IEP(1),5))**2 IF(PT2APP.GT.VETOKT**2) THEN c DES IF(itime.le.50) print*,' fsr veto ',sqrt(pt2app),' > ',vetokt GOTO 400 ENDIF C/* cc print*,paru(2)/(B0*LOG(qxres**2/ALAMS)),qxres,B0,sqrt(alams) C*/ IF(PT2APP.LT.PT2MIN) GOTO 400 BLFM=LOG(MAX(PT2MIN,PT2APP)/ALAMS) BTEMP=BLFM/(1D0-B1*LOG(BLFM)/(B0**2*BLFM)) csmxxx IF(ALFM/LOG(PT2APP/ALAMS).LT.PYR(0)) GOTO 400 IF(ATEMP/BTEMP.LT.PYR(0)) GOTO 400 ENDIF ENDIF IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1D0-ZL)/MAX(1D-20,ZL* &(1D0-ZU))) IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU)) C...Width suppression for q -> q + g. IF(MSTJ(40).NE.0.AND.KFL(1).NE.21.AND.IPSPD.EQ.0) THEN IF(IGM.EQ.0) THEN EGLU=0.5D0*PS(5)*(1D0-Z)*(1D0+V(IEP(1),5)/V(NS+1,5)) ELSE EGLU=PMED*(1D0-Z) ENDIF CHI=PARJ(89)**2/(PARJ(89)**2+EGLU**2) IF(MSTJ(40).EQ.1) THEN IF(CHI.LT.PYR(0)) GOTO 400 ELSEIF(MSTJ(40).EQ.2) THEN IF(1D0-CHI.LT.PYR(0)) GOTO 400 ENDIF ENDIF C...Three-jet matrix element correction. IF(M3JC.GE.1) THEN WME=1D0 WSHOW=1D0 C...QED matrix elements: only for massless case so far. IF(MCE.EQ.2.AND.IGM.EQ.0) THEN X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5)) X2=1D0-V(IEP(1),5)/V(NS+1,5) X3=(1D0-X1)+(1D0-X2) KI1=K(IPA(INUM),2) KI2=K(IPA(3-INUM),2) QF1=KCHG(PYCOMP(KI1),1)*ISIGN(1,KI1)/3D0 QF2=KCHG(PYCOMP(KI2),1)*ISIGN(1,KI2)/3D0 WSHOW=QF1**2*(1D0-X1)/X3*(1D0+(X1/(2D0-X2))**2)+ & QF2**2*(1D0-X2)/X3*(1D0+(X2/(2D0-X1))**2) WME=(QF1*(1D0-X1)/X3-QF2*(1D0-X2)/X3)**2*(X1**2+X2**2) ELSEIF(MCE.EQ.2) THEN C...QCD matrix elements, including mass effects. ELSEIF(MSTJ(49).NE.1.AND.K(IEP(1),2).NE.21) THEN PS1ME=V(IEP(1),5) PM1ME=PMTH(1,IR) M3JCC=M3JC IF(IR.GE.31.AND.IGM.EQ.0) THEN C...QCD ME: original parton, first branching. PM2ME=PMTH(1,63-IR) ECMME=PS(5) ELSEIF(IR.GE.31) THEN C...QCD ME: original parton, subsequent branchings. PM2ME=PMTH(1,63-IR) PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5)) ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2)) ELSEIF(K(IM,2).EQ.21) THEN C...QCD ME: secondary partons, first branching. PM2ME=PM1ME ZMME=V(IM,1) IF(IEP(1).GT.IEP(2)) ZMME=1D0-ZMME PMLME=SQRT(MAX(0D0,(V(IM,5)-PS1ME-PM2ME**2)**2- & 4D0*PS1ME*PM2ME**2)) PEDME=PEM*(0.5D0*(V(IM,5)-PMLME+PS1ME-PM2ME**2)+PMLME*ZMME)/ & V(IM,5) ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2)) M3JCC=66 ELSE C...QCD ME: secondary partons, subsequent branchings. PM2ME=PM1ME PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5)) ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2)) M3JCC=66 ENDIF C...Construct ME variables. R1ME=PM1ME/ECMME R2ME=PM2ME/ECMME X1=(1D0+PS1ME/ECMME**2-R2ME**2)*(Z+(1D0-Z)*PM1ME**2/PS1ME) X2=1D0+R2ME**2-PS1ME/ECMME**2 C...Call ME, with right order important for two inequivalent showerers. IF(IR.EQ.IORD+30) THEN WME=PYMAEL(M3JCC,X1,X2,R1ME,R2ME,ALPHA) ELSE WME=PYMAEL(M3JCC,X2,X1,R2ME,R1ME,ALPHA) ENDIF C...Split up total ME when two radiating partons. ISPRAD=1 IF((M3JCC.GE.16.AND.M3JCC.LE.19).OR. & (M3JCC.GE.26.AND.M3JCC.LE.29).OR. & (M3JCC.GE.36.AND.M3JCC.LE.39).OR. & (M3JCC.GE.46.AND.M3JCC.LE.49).OR. & (M3JCC.GE.56.AND.M3JCC.LE.64)) ISPRAD=0 IF(ISPRAD.EQ.1) WME=WME*MAX(1D-10,1D0+R1ME**2-R2ME**2-X1)/ & MAX(1D-10,2D0-X1-X2) C...Evaluate shower rate to be compared with. WSHOW=2D0/(MAX(1D-10,2D0-X1-X2)* & MAX(1D-10,1D0+R2ME**2-R1ME**2-X2)) IF(IGLUI.EQ.1.AND.IR.GE.31) WSHOW=(9D0/4D0)*WSHOW ELSEIF(MSTJ(49).NE.1) THEN C...Toy model scalar theory matrix elements; no mass effects. ELSE X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5)) X2=1D0-V(IEP(1),5)/V(NS+1,5) X3=(1D0-X1)+(1D0-X2) WSHOW=4D0*X3*((1D0-X1)/(2D0-X2)**2+(1D0-X2)/(2D0-X1)**2) WME=X3**2 IF(MSTJ(102).GE.2) WME=X3**2-2D0*(1D0+X3)*(1D0-X1)*(1D0-X2)* & PARJ(171) ENDIF IF(WME.LT.PYR(0)*WSHOW) GOTO 400 ENDIF C...Impose angular ordering by rejection of nonordered emission. IF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2.AND.IPSPD.EQ.0) THEN PEMAO=V(IM,1)*P(IM,4) IF(IEP(1).EQ.N+2) PEMAO=(1D0-V(IM,1))*P(IM,4) IF(IR.GE.31.AND.MSTJ(42).GE.5) THEN MAOD=0 ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.4 & .OR.MSTJ(42).EQ.7)) THEN MAOD=0 ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.3 & .OR.MSTJ(42).EQ.6)) THEN MAOD=1 PMDAO=PMTH(2,K(IEP(1),5)) THE2ID=Z*(1D0-Z)*PEMAO**2/(V(IEP(1),5)-4D0*PMDAO**2) ELSE MAOD=1 THE2ID=Z*(1D0-Z)*PEMAO**2/V(IEP(1),5) IF(MSTJ(42).GE.3.AND.MSTJ(42).NE.5) THE2ID=THE2ID* & (1D0+PMTH(1,IR)**2*(1D0-Z)/(V(IEP(1),5)*Z))**2 ENDIF ccc CLKT2D=PEMAO**2/THE2ID*MIN(Z,1D0-Z)**2 MAOM=1 IAOM=IM 430 IF(K(IAOM,5).EQ.22) THEN IAOM=K(IAOM,3) IF(K(IAOM,3).LE.NS) MAOM=0 IF(MAOM.EQ.1) GOTO 430 ENDIF IF(MAOM.EQ.1.AND.MAOD.EQ.1) THEN THE2IM=V(IAOM,1)*(1D0-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5) ccc CLKT2M=P(IAOM,4)**2/THE2IM*MIN(V(IAOM,1),1D0-V(IAOM,1))**2 IF(THE2ID.LT.THE2IM) GOTO 400 ccc IF(CLKT2D.GT.CKLT2M) GOTO 400 ENDIF ENDIF C...Impose user-defined maximum angle at first branching. IF(IPSPD.EQ.0.AND.MCE.EQ.1.AND.IEP(1).LT.NS+NEP+1) THEN THE2ID=1D0-.5D0*V(IEP(1),5)/(Z*(1D0-Z)*PED**2) IF(THE2ID.LT.0D0) THE2ID=0D0 IF(THE2ID.GT.2D0) THE2ID=2D0 THE2ID=1D0/ACOS(THE2ID)**2 C/* THE2ID=Z*(1D0-Z)*PED**2/V(IEP(1),5) c DES if(itime.lt.50) then c DES print*,' theta ',vetoth*1D0/sqrt(the2id)-1D0 c DES endif IF(VETOTH**2*THE2ID.LT.1D0) GOTO 400 ENDIF C...Impose user-defined maximum angle at first branching. IF(MSTJ(48).EQ.1.AND.IPSPD.EQ.0) THEN IF(NEP.EQ.1.AND.IM.EQ.NS) THEN THE2ID=Z*(1D0-Z)*PS(4)**2/V(IEP(1),5) IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 400 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5) IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 400 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5) IF(PARJ(86)**2*THE2ID.LT.1D0) GOTO 400 ENDIF ENDIF C...Impose angular constraint in first branching from interference C...with initial state partons. IF(MIIS.GE.2.AND.IEP(1).LE.NS+3) THEN THE2D=MAX((1D0-Z)/Z,Z/(1D0-Z))*V(IEP(1),5)/(0.5D0*P(IM,4))**2 IF(IEP(1).EQ.NS+2.AND.ISII(1).GE.1) THEN IF(THE2D.GT.THEIIS(1,ISII(1))**2) GOTO 400 ELSEIF(IEP(1).EQ.NS+3.AND.ISII(2).GE.1) THEN IF(THE2D.GT.THEIIS(2,ISII(2))**2) GOTO 400 ENDIF ENDIF C...End of inner veto algorithm. Check if only one leg evolved so far. 440 V(IEP(1),1)=Z c$$$ print*,' branching at z = ',z,sqrt(v(iep(1),5)) c$$$ c$$$ print*,' got here ',iep(1),nep c$$$ print*,(v(iep(1),i),i=1,5) c$$$ print*,(p(iep(1),i),i=1,5) c$$$ print*,(v(iep(2),i),i=1,5) c$$$ print*,(p(iep(2),i),i=1,5) c$$$ print*,im,p(im,5) ISL(1)=0 ISL(2)=0 IF(NEP.EQ.1) GOTO 480 IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 340 DO 450 I=1,NEP IR=IREF(N+I-NS) IF(ITRY(I).EQ.0.AND.KSH(IR).EQ.1) THEN IF(P(N+I,5).GE.PMTH(2,IR)) GOTO 340 ENDIF 450 CONTINUE C...Check if chosen multiplet m1,m2,z1,z2 is physical. IF(NEP.GE.3) THEN PMSUM=0D0 DO 460 I=1,NEP PMSUM=PMSUM+P(N+I,5) 460 CONTINUE IF(PMSUM.GE.PS(5)) GOTO 340 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN DO 470 I1=N+1,N+2 IRDA=IREF(I1-NS) IF(KSH(IRDA).EQ.0) GOTO 470 IF(P(I1,5).LT.PMTH(2,IRDA)) GOTO 470 IF(IRDA.EQ.21) THEN IRGD1=IABS(K(I1,5)) IRGD2=IRGD1 ELSE IRGD1=IRDA IRGD2=IABS(K(I1,5)) ENDIF I2=2*N+3-I1 IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN PED=0.5D0*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5) ELSE IF(I1.EQ.N+1) ZM=V(IM,1) IF(I1.EQ.N+2) ZM=1D0-V(IM,1) PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2- & 4D0*V(N+1,5)*V(N+2,5)) PED=PEM*(0.5D0*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/ & V(IM,5) ENDIF IF(MOD(MSTJ(43),2).EQ.1) THEN PMQTH3=0.5D0*PARJ(82) IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83) IF(IRGD2.EQ.22.AND.ISCOL(IRDA).EQ.0) PMQTH3=0.5D0*PARJ(90) PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(I1,5) PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(I1,5) ZD=SQRT(MAX(0D0,(1D0-V(I1,5)/PED**2)*((1D0-PMQ1-PMQ2)**2- & 4D0*PMQ1*PMQ2))) ZH=1D0+PMQ1-PMQ2 ELSE ZD=SQRT(MAX(0D0,1D0-V(I1,5)/PED**2)) ZH=1D0 ENDIF IF(IRDA.EQ.21.AND.IRGD1.LT.10.AND. & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN ELSE ZL=0.5D0*(ZH-ZD) ZU=0.5D0*(ZH+ZD) IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND. & ISSET(1).EQ.0) THEN ISL(1)=1 ELSEIF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND. & ISSET(2).EQ.0) THEN ISL(2)=1 ENDIF ENDIF IF(IRDA.EQ.21) V(I1,4)=LOG(ZU*(1D0-ZL)/MAX(1D-20, & ZL*(1D0-ZU))) IF(IRDA.NE.21) V(I1,4)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU)) 470 CONTINUE IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN ISL(3-ISLM)=0 ISLM=3-ISLM ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN ZDR1=MAX(0D0,V(N+1,3)/MAX(1D-6,V(N+1,4))-1D0) ZDR2=MAX(0D0,V(N+2,3)/MAX(1D-6,V(N+2,4))-1D0) IF(ZDR2.GT.PYR(0)*(ZDR1+ZDR2)) ISL(1)=0 IF(ISL(1).EQ.1) ISL(2)=0 IF(ISL(1).EQ.0) ISLM=1 IF(ISL(2).EQ.0) ISLM=2 ENDIF IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 340 ENDIF IRD1=IREF(N+1-NS) IRD2=IREF(N+2-NS) IF(IGM.GT.0) THEN IF(MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE. & PMTH(2,IRD1).OR.P(N+2,5).GE.PMTH(2,IRD2))) THEN PMQ1=V(N+1,5)/V(IM,5) PMQ2=V(N+2,5)/V(IM,5) ZD=SQRT(MAX(0D0,(1D0-V(IM,5)/PEM**2)*((1D0-PMQ1-PMQ2)**2- & 4D0*PMQ1*PMQ2))) ZH=1D0+PMQ1-PMQ2 ZL=0.5D0*(ZH-ZD) ZU=0.5D0*(ZH+ZD) IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 340 ENDIF ENDIF C...Accepted branch. Construct four-momentum for initial partons. 480 MAZIP=0 MAZIC=0 IF(NEP.EQ.1) THEN P(N+1,1)=0D0 P(N+1,2)=0D0 P(N+1,3)=SQRT(MAX(0D0,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)- & P(N+1,5)))) P(N+1,4)=P(IPA(1),4) V(N+1,2)=P(N+1,4) ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN PED1=0.5D0*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5) P(N+1,1)=0D0 P(N+1,2)=0D0 P(N+1,3)=SQRT(MAX(0D0,(PED1+P(N+1,5))*(PED1-P(N+1,5)))) P(N+1,4)=PED1 P(N+2,1)=0D0 P(N+2,2)=0D0 P(N+2,3)=-P(N+1,3) P(N+2,4)=P(IM,5)-PED1 V(N+1,2)=P(N+1,4) V(N+2,2)=P(N+2,4) ELSEIF(NEP.GE.3) THEN C...Rescale all momenta for energy conservation. LOOP=0 PES=0D0 PQS=0D0 DO 500 I=1,NEP DO 490 J=1,4 P(N+I,J)=P(IPA(I),J) 490 CONTINUE PES=PES+P(N+I,4) PQS=PQS+P(N+I,5)**2/P(N+I,4) 500 CONTINUE 510 LOOP=LOOP+1 FAC=(PS(5)-PQS)/(PES-PQS) PES=0D0 PQS=0D0 DO 530 I=1,NEP DO 520 J=1,3 P(N+I,J)=FAC*P(N+I,J) 520 CONTINUE P(N+I,4)=SQRT(P(N+I,5)**2+P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2) V(N+I,2)=P(N+I,4) PES=PES+P(N+I,4) PQS=PQS+P(N+I,5)**2/P(N+I,4) 530 CONTINUE IF(LOOP.LT.10.AND.ABS(PES-PS(5)).GT.1D-12*PS(5)) GOTO 510 C...Construct transverse momentum for ordinary branching in shower. ELSE ZM=V(IM,1) LOOPPT=0 540 LOOPPT=LOOPPT+1 PZM=SQRT(MAX(0D0,(PEM+P(IM,5))*(PEM-P(IM,5)))) PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4D0*V(N+1,5)*V(N+2,5) IF(PZM.LE.0D0) THEN PTS=0D0 ELSEIF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND. & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN PTS=PMLS*ZM*(1D0-ZM)/V(IM,5) ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN PTS=(PEM**2*(ZM*(1D0-ZM)*V(IM,5)-(1D0-ZM)*V(N+1,5)- & ZM*V(N+2,5))-0.25D0*PMLS)/PZM**2 ELSE PTS=PMLS*(ZM*(1D0-ZM)*PEM**2/V(IM,5)-0.25D0)/PZM**2 ENDIF IF(PTS.LT.0D0.AND.LOOPPT.LT.10) THEN ZM=0.05D0+0.9D0*ZM GOTO 540 ELSEIF(PTS.LT.0D0) THEN GOTO 270 ENDIF PT=SQRT(MAX(0D0,PTS)) C...Find coefficient of azimuthal asymmetry due to gluon polarization. HAZIP=0D0 IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21 & .AND.IAU.NE.0) THEN IF(K(IGM,3).NE.0) MAZIP=1 ZAU=V(IGM,1) IF(IAU.EQ.IM+1) ZAU=1D0-V(IGM,1) IF(MAZIP.EQ.0) ZAU=0D0 IF(K(IGM,2).NE.21) THEN HAZIP=2D0*ZAU/(1D0+ZAU**2) ELSE HAZIP=(ZAU/(1D0-ZAU*(1D0-ZAU)))**2 ENDIF IF(K(N+1,2).NE.21) THEN HAZIP=HAZIP*(-2D0*ZM*(1D0-ZM))/(1D0-2D0*ZM*(1D0-ZM)) ELSE HAZIP=HAZIP*(ZM*(1D0-ZM)/(1D0-ZM*(1D0-ZM)))**2 ENDIF ENDIF C...Find coefficient of azimuthal asymmetry due to soft gluon C...interference. HAZIC=0D0 IF(MSTJ(49).NE.2.AND.MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR. & K(N+2,2).EQ.21).AND.IAU.NE.0) THEN IF(K(IGM,3).NE.0) MAZIC=N+1 IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2 IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND. & ZM.GT.0.5D0) MAZIC=N+2 IF(K(IAU,2).EQ.22) MAZIC=0 ZS=ZM IF(MAZIC.EQ.N+2) ZS=1D0-ZM ZGM=V(IGM,1) IF(IAU.EQ.IM-1) ZGM=1D0-V(IGM,1) IF(MAZIC.EQ.0) ZGM=1D0 IF(MAZIC.NE.0) HAZIC=(P(IM,5)/P(IGM,5))* & SQRT((1D0-ZS)*(1D0-ZGM)/(ZS*ZGM)) HAZIC=MIN(0.95D0,HAZIC) ENDIF ENDIF C...Construct energies for ordinary branching in shower. 550 IF(NEP.EQ.2.AND.IGM.GT.0) THEN IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND. & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN P(N+1,4)=0.5D0*(PEM*(V(IM,5)+V(N+1,5)-V(N+2,5))+ & PZM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5) ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN P(N+1,4)=PEM*V(IM,1) ELSE CMrenna...Modify here for new z P(N+1,4)=PEM*(0.5D0*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+ & SQRT(PMLS)*ZM)/V(IM,5) CX P(N+1,4)=V(IM,1)*PZM+EBSTAR/P(IM,5)*(PEM-PZM) ENDIF C...Already predetermined choice of phi angle or not PHI=PARU(2)*PYR(0) IF(MPSPD.EQ.1.AND.IGM.EQ.NS+1) THEN IPSPD=IP1+IM-NS-2 IF(K(IPSPD,4).GT.0) THEN IPSGD1=K(IPSPD,4) IF(IM.EQ.NS+2) THEN PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2)) ELSE PHI=PYANGL(-P(IPSGD1,1),P(IPSGD1,2)) ENDIF ENDIF ELSEIF(MPSPD.EQ.1.AND.IGM.EQ.NS+2) THEN IPSPD=IP1+IM-NS-2 IF(K(IPSPD,4).GT.0) THEN IPSGD1=K(IPSPD,4) PHIPSM=PYANGL(P(IPSPD,1),P(IPSPD,2)) THEPSM=PYANGL(P(IPSPD,3),SQRT(P(IPSPD,1)**2+P(IPSPD,2)**2)) CALL PYROBO(IPSGD1,IPSGD1,0D0,-PHIPSM,0D0,0D0,0D0) CALL PYROBO(IPSGD1,IPSGD1,-THEPSM,0D0,0D0,0D0,0D0) PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2)) CALL PYROBO(IPSGD1,IPSGD1,THEPSM,PHIPSM,0D0,0D0,0D0) ENDIF ENDIF C...Construct momenta for ordinary branching in shower. P(N+1,1)=PT*COS(PHI) P(N+1,2)=PT*SIN(PHI) IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND. & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN P(N+1,3)=0.5D0*(PZM*(V(IM,5)+V(N+1,5)-V(N+2,5))+ & PEM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5) ELSEIF(PZM.GT.0D0) THEN CMrenna...Modify here for new z P(N+1,3)=0.5D0*(V(N+2,5)-V(N+1,5)-V(IM,5)+ & 2D0*PEM*P(N+1,4))/PZM CX P(N+1,3)=V(IM,1)*PEM-EBSTAR/P(IM,5)*(PEM-PZM) ELSE P(N+1,3)=0D0 ENDIF P(N+2,1)=-P(N+1,1) P(N+2,2)=-P(N+1,2) P(N+2,3)=PZM-P(N+1,3) P(N+2,4)=PEM-P(N+1,4) IF(MSTJ(43).LE.2) THEN V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5) V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5) ENDIF ENDIF C...Rotate and boost daughters. IF(IGM.GT.0) THEN IF(MSTJ(43).LE.2) THEN BEX=P(IGM,1)/P(IGM,4) BEY=P(IGM,2)/P(IGM,4) BEZ=P(IGM,3)/P(IGM,4) GA=P(IGM,4)/P(IGM,5) GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1D0+GA)- & P(IM,4)) ELSE BEX=0D0 BEY=0D0 BEZ=0D0 GA=1D0 GABEP=0D0 ENDIF PTIMB=SQRT((P(IM,1)+GABEP*BEX)**2+(P(IM,2)+GABEP*BEY)**2) THE=PYANGL(P(IM,3)+GABEP*BEZ,PTIMB) IF(PTIMB.GT.1D-4) THEN PHI=PYANGL(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY) ELSE PHI=0D0 ENDIF DO 560 I=N+1,N+2 DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+ & SIN(THE)*COS(PHI)*P(I,3) DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+ & SIN(THE)*SIN(PHI)*P(I,3) DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3) DP(4)=P(I,4) DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3) DGABP=GA*(GA*DBP/(1D0+GA)+DP(4)) P(I,1)=DP(1)+DGABP*BEX P(I,2)=DP(2)+DGABP*BEY P(I,3)=DP(3)+DGABP*BEZ P(I,4)=GA*(DP(4)+DBP) 560 CONTINUE ENDIF C...Weight with azimuthal distribution, if required. IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN DO 570 J=1,3 DPT(1,J)=P(IM,J) DPT(2,J)=P(IAU,J) DPT(3,J)=P(N+1,J) 570 CONTINUE DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3) DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3) DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2 DO 580 J=1,3 DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/MAX(1D-10,DPMM) DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/MAX(1D-10,DPMM) 580 CONTINUE DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2) DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2) IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1D0*PARJ(82)) THEN CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+ & DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4)) IF(MAZIP.NE.0) THEN IF(1D0+HAZIP*(2D0*CAD**2-1D0).LT.PYR(0)*(1D0+ABS(HAZIP))) & GOTO 550 ENDIF IF(MAZIC.NE.0) THEN IF(MAZIC.EQ.N+2) CAD=-CAD IF((1D0-HAZIC)*(1D0-HAZIC*CAD)/(1D0+HAZIC**2-2D0*HAZIC*CAD) & .LT.PYR(0)) GOTO 550 ENDIF ENDIF ENDIF C...Azimuthal anisotropy due to interference with initial state partons. IF(MOD(MIIS,2).EQ.1.AND.IGM.EQ.NS+1.AND.(K(N+1,2).EQ.21.OR. &K(N+2,2).EQ.21)) THEN III=IM-NS-1 IF(ISII(III).GE.1) THEN IAZIID=N+1 IF(K(N+1,2).NE.21) IAZIID=N+2 IF(K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND. & P(N+1,4).GT.P(N+2,4)) IAZIID=N+2 THEIID=PYANGL(P(IAZIID,3),SQRT(P(IAZIID,1)**2+P(IAZIID,2)**2)) IF(III.EQ.2) THEIID=PARU(1)-THEIID PHIIID=PYANGL(P(IAZIID,1),P(IAZIID,2)) HAZII=MIN(0.95D0,THEIID/THEIIS(III,ISII(III))) CAD=COS(PHIIID-PHIIIS(III,ISII(III))) PHIREL=ABS(PHIIID-PHIIIS(III,ISII(III))) IF(PHIREL.GT.PARU(1)) PHIREL=PARU(2)-PHIREL IF((1D0-HAZII)*(1D0-HAZII*CAD)/(1D0+HAZII**2-2D0*HAZII*CAD) & .LT.PYR(0)) GOTO 550 ENDIF ENDIF C...Continue loop over partons that may branch, until none left. IF(IGM.GE.0) K(IM,1)=14 N=N+NEP NEP=2 IF(N.GT.MSTU(4)-MSTU(32)-10) THEN CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS') IF(MSTU(21).GE.1) N=NS IF(MSTU(21).GE.1) RETURN ENDIF GOTO 280 C...Set information on imagined shower initiator. 590 IF(NPA.GE.2) THEN K(NS+1,1)=11 K(NS+1,2)=94 K(NS+1,3)=IP1 IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2 K(NS+1,4)=NS+2 K(NS+1,5)=NS+1+NPA IIM=1 ELSE IIM=0 ENDIF C...Reconstruct string drawing information. DO 600 I=NS+1+IIM,N KQ=KCHG(PYCOMP(K(I,2)),2) IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN K(I,1)=1 ELSEIF(K(I,1).LE.10.AND.IABS(K(I,2)).GE.11.AND. & IABS(K(I,2)).LE.18) THEN K(I,1)=1 ELSEIF(K(I,1).LE.10) THEN K(I,4)=MSTU(5)*(K(I,4)/MSTU(5)) K(I,5)=MSTU(5)*(K(I,5)/MSTU(5)) ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN ID1=MOD(K(I,4),MSTU(5)) IF(KQ.EQ.1.AND.K(I,2).GT.0) ID1=MOD(K(I,4),MSTU(5))+1 IF(KQ.EQ.2.AND.(K(ID1,2).EQ.21.OR.K(ID1+1,2).EQ.21).AND. & PYR(0).GT.0.5D0) ID1=MOD(K(I,4),MSTU(5))+1 ID2=2*MOD(K(I,4),MSTU(5))+1-ID1 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2 K(ID1,4)=K(ID1,4)+MSTU(5)*I K(ID1,5)=K(ID1,5)+MSTU(5)*ID2 K(ID2,4)=K(ID2,4)+MSTU(5)*ID1 K(ID2,5)=K(ID2,5)+MSTU(5)*I ELSE ID1=MOD(K(I,4),MSTU(5)) ID2=ID1+1 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1 IF(KQ.EQ.1.OR.K(ID1,1).GE.11) THEN K(ID1,4)=K(ID1,4)+MSTU(5)*I K(ID1,5)=K(ID1,5)+MSTU(5)*I ELSE K(ID1,4)=0 K(ID1,5)=0 ENDIF K(ID2,4)=0 K(ID2,5)=0 ENDIF 600 CONTINUE C...Transformation from CM frame. IF(NPA.EQ.1) THEN THE=PYANGL(P(IPA(1),3),SQRT(P(IPA(1),1)**2+P(IPA(1),2)**2)) PHI=PYANGL(P(IPA(1),1),P(IPA(1),2)) MSTU(33)=1 CALL PYROBO(NS+1,N,THE,PHI,0D0,0D0,0D0) ELSEIF(NPA.EQ.2) THEN BEX=PS(1)/PS(4) BEY=PS(2)/PS(4) BEZ=PS(3)/PS(4) GA=PS(4)/PS(5) GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3)) & /(1D0+GA)-P(IPA(1),4)) THE=PYANGL(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1) & +GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2)) PHI=PYANGL(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY) MSTU(33)=1 CALL PYROBO(NS+1,N,THE,PHI,BEX,BEY,BEZ) ELSE CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,PS(1)/PS(4),PS(2)/PS(4), & PS(3)/PS(4)) MSTU(33)=1 CALL PYROBO(NS+1,N,0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),PS(3)/PS(4)) ENDIF C...Decay vertex of shower. DO 620 I=NS+1,N DO 610 J=1,5 V(I,J)=V(IP1,J) 610 CONTINUE 620 CONTINUE C...Delete trivial shower, else connect initiators. IF(N.LE.NS+NPA+IIM) THEN N=NS ELSE DO 630 IP=1,NPA K(IPA(IP),1)=14 K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP K(NS+IIM+IP,3)=IPA(IP) IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1 IF(K(NS+IIM+IP,1).NE.1) THEN K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4) K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5) ENDIF 630 CONTINUE ENDIF RETURN END C********************************************************************* C...PYADSH C...Administers the generation of successive final-state showers C...in external processes. SUBROUTINE PYADSH(NFIN) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C INTEGER PYK,PYCHGE,PYCOMP C...Commonblocks. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) COMMON/PYINT1/MINT(400),VINT(400) SAVE /PYJETS/,/PYDAT1/,/PYPARS/,/PYINT1/ C...Local array. C DIMENSION IBEG(100),KSAV(10,5),IORD(10),PSUM(4),BETA(3) DIMENSION IBEG(100),KSAV(10,5),PSUM(4),BETA(3) INTEGER IFIRST DATA IFIRST /0/ SAVE IFIRST ifirst=ifirst+1 C...Set primary vertex. DO 100 J=1,5 V(MINT(83)+5,J)=0D0 V(MINT(83)+6,J)=0D0 V(MINT(84)+1,J)=0D0 V(MINT(84)+2,J)=0D0 100 CONTINUE C...Isolate systems of particles with the same mother. NSYS=0 IMS=-1 DO 140 I=MINT(84)+3,NFIN IM=K(I,3) IF(IM.GT.0.AND.IM.LE.MINT(84)) IM=K(IM,3) IF(IM.NE.IMS) THEN NSYS=NSYS+1 IBEG(NSYS)=I IMS=IM ENDIF C...Set production vertices. IF(IM.LE.MINT(83)+6.OR.(IM.GT.MINT(84).AND.IM.LE.MINT(84)+2)) & THEN DO 110 J=1,4 V(I,J)=0D0 110 CONTINUE ELSE DO 120 J=1,4 V(I,J)=V(IM,J)+V(IM,5)*P(IM,J)/P(IM,5) 120 CONTINUE ENDIF IF(MSTP(125).GE.1) THEN IDOC=I-MSTP(126)+4 DO 130 J=1,5 V(IDOC,J)=V(I,J) 130 CONTINUE ENDIF 140 CONTINUE C...End loop over systems. Return if no showers to be performed. IBEG(NSYS+1)=NFIN+1 IF(MSTP(71).LE.0) RETURN C...Loop through systems of particles; check that sensible size. DO 260 ISYS=1,NSYS NSIZ=IBEG(ISYS+1)-IBEG(ISYS) IF(NSIZ.EQ.1.AND.ISYS.EQ.1) THEN ELSEIF(NSIZ.LE.1) THEN CALL PYERRM(2,'(PYADSH:) only one particle in system') ELSEIF(NSIZ.GT.7) THEN CALL PYERRM(2,'(PYADSH:) more than seven particles in system') ELSE C...Save status codes and daughters of showering pair; reset them. DO 150 J=1,4 PSUM(J)=0D0 150 CONTINUE DO 170 II=1,NSIZ I=IBEG(ISYS)-1+II KSAV(II,1)=K(I,1) IF(K(I,1).GT.10) THEN K(I,1)=1 IF(KSAV(II,1).EQ.14) K(I,1)=3 ENDIF IF(KSAV(II,1).LE.10) THEN ELSEIF(K(I,1).EQ.1) THEN KSAV(II,4)=K(I,4) KSAV(II,5)=K(I,5) K(I,4)=0 K(I,5)=0 ELSE KSAV(II,4)=MOD(K(I,4),MSTU(5)) KSAV(II,5)=MOD(K(I,5),MSTU(5)) K(I,4)=K(I,4)-KSAV(II,4) K(I,5)=K(I,5)-KSAV(II,5) ENDIF DO 160 J=1,4 PSUM(J)=PSUM(J)+P(I,J) 160 CONTINUE 170 CONTINUE C...Perform shower. QMAX=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2- & PSUM(3)**2)) IF(ISYS.EQ.1) QMAX=MIN(QMAX,SQRT(PARP(71))*VINT(55)) NSAV=N MINT(400)=0 IF(ISYS.EQ.2) MINT(400)=1 c DES if(ifirst.lt.10) then c DES print*,' QMAX, ISYS = ',qmax,isys,nsiz,nsys c DES print*,' beg ',ibeg(isys),ibeg(isys)+1 c DES endif IF(NSIZ.EQ.2.and.MINT(400).EQ.0) THEN CALL PYSHOW(IBEG(ISYS),IBEG(ISYS)+1,QMAX) ELSE CALL PYSHOW(IBEG(ISYS),-NSIZ,QMAX) ENDIF MINT(400)=0 C...Look up showered copies of original showering particles. DO 250 II=1,NSIZ I=IBEG(ISYS)-1+II IMV=I IF(N.EQ.NSAV.OR.K(I,1).LE.10) THEN ELSEIF(K(I,1).EQ.11) THEN 180 IMV=MOD(K(IMV,4),MSTU(5)) IF(K(IMV,1).EQ.11) GOTO 180 ELSE KDA1=MOD(K(I,4),MSTU(5)) KDA2=MOD(K(I,5),MSTU(5)) DO 190 I3=I+1,N IF(K(I3,2).EQ.K(I,2).AND.(I3.EQ.KDA1.OR.I3.EQ.KDA2)) & THEN IMV=I3 KDA1=MOD(K(I3,4),MSTU(5)) KDA2=MOD(K(I3,5),MSTU(5)) ENDIF 190 CONTINUE ENDIF C...Restore daughter info of original partons to showered copies. IF(KSAV(II,1).GT.10) K(IMV,1)=KSAV(II,1) IF(KSAV(II,1).LE.10) THEN ELSEIF(K(I,1).EQ.1) THEN K(IMV,4)=KSAV(II,4) K(IMV,5)=KSAV(II,5) ELSE K(IMV,4)=K(IMV,4)+KSAV(II,4) K(IMV,5)=K(IMV,5)+KSAV(II,5) ENDIF C...Reset mother info of existing daughters to showered copies. DO 200 I3=IBEG(ISYS+1),NFIN IF(K(I3,3).EQ.I) K(I3,3)=IMV IF(K(I3,1).EQ.3.OR.K(I3,1).EQ.14) THEN IF(K(I3,4)/MSTU(5).EQ.I) K(I3,4)=K(I3,4)+MSTU(5)*(IMV-I) IF(K(I3,5)/MSTU(5).EQ.I) K(I3,5)=K(I3,5)+MSTU(5)*(IMV-I) ENDIF 200 CONTINUE C...Boost all original daughters to new frame of showered copy. IF(IMV.NE.I) THEN DO 210 J=1,3 BETA(J)=(P(IMV,J)-P(I,J))/(P(IMV,4)+P(I,4)) 210 CONTINUE FAC=2D0/(1D0+BETA(1)**2+BETA(2)**2+BETA(3)**2) DO 220 J=1,3 BETA(J)=FAC*BETA(J) 220 CONTINUE DO 240 I3=IBEG(ISYS+1),NFIN IMO=I3 230 IMO=K(IMO,3) IF(IMO.GT.0.AND.IMO.NE.I.AND.IMO.NE.K(I,3)) GOTO 230 IF(IMO.EQ.I.OR.(K(I,3).LE.MINT(84).AND.IMO.EQ.K(I,3))) & CALL PYROBO(I3,I3,0D0,0D0,BETA(1),BETA(2),BETA(3)) 240 CONTINUE ENDIF 250 CONTINUE C...End of loop over showering systems ENDIF 260 CONTINUE RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE CLUSTER(NL,NH,NMATCH,IVERB) IMPLICIT NONE C...Double precision and integer declarations. cc INTEGER PYK,PYCHGE,PYCOMP C...Commonblocks. C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ C...Local arrays and saved variables. REAL*8 R2,R2MIN,R2MAX,D2(500) INTEGER IMIN,IKNT,ITRY1,ITRY2,K,IVERB INTEGER INDEX(200),JNDEX(200),IDH2 C.....Hard-wired --> may need to change for more partons INTEGER IOKBEG,IOKEND PARAMETER(IOKBEG=95,IOKEND=135) INTEGER IOK(IOKBEG:IOKEND) INTEGER IPICK(IOKBEG:IOKEND) INTEGER IPAIR,NPAIR integer icid INTEGER IC1,IC2,JC1,JC2,ITMP,IGL1,IGL2 INTEGER ID1,ID2,IP1,IP2,Imatch INTEGER NL1,NH1 integer icola,jcola REAL*8 R2M INTEGER I,NL,NH,NMATCH,NSOFAR,IMIN1,IMIN2,J INTEGER ISIG integer iflag,inext save iflag,inext INTEGER ICOLL REAL*8 PUPA(12,MAXNUP) COMMON/HEPA/PUPA,ICOLL SAVE /HEPA/ REAL*8 RESKT(0:10) COMMON/RES/RESKT SAVE /RES/ real*8 pt,ptt,PPLUS,PMINS integer IT,IDA,IDB,NL0,NH0,IST,IOG,NUPP,ICLUS integer NCHG,is1 REAL*8 KTMDPI real*8 etamax,eps INTEGER TYPE,ANGL,MONO,RECO COMMON/CLSTR/TYPE,ANGL,MONO,RECO SAVE /CLSTR/ integer iskip common/skip/iskip save /skip/ real*8 qboson c integer ifam integer lcola,lcolb,n_fermions integer pychge logical check_color data iflag/0/ DATA ETAMAX,EPS/10,1D-6/ DATA R2MAX/1D20/ iskip=0 check_color=.true. c check_color=.false. NL1=NL NH1=NH C/* Flag for when check_color set false INEXT=0 100 CONTINUE NL=NL1 NH=NH1 ICID=500 ICLUS=0 C/* Initialize status and Id do I=IOKBEG,IOKEND ISTUP(I)=0 IDUP(I)=0 enddo C/* Copy over information to 100+ NUPP=0 NCHG=0 QBOSON=0 DO 5 I=1,NUP if(istup(i).eq.2) goto 5 NCHG=NCHG+PYCHGE(IDUP(I))*ISTUP(I) IDA=ABS(IDUP(I)) NUPP=NUPP+1 IT=100+NUPP IF(IDA.GE.11.AND.IDA.LE.16) THEN QBOSON=QBOSON+PYCHGE(IDUP(I)) IF(ISTUP(I).EQ.-1) THEN ISTUP(IT)=0 ELSE ISTUP(IT)=0 ENDIF ELSE ISTUP(IT)=ISTUP(I) ENDIF if(istup(i).eq.7) then istup(it)=1 elseif(istup(i).gt.1) then istup(it)=100 elseif(istup(i).lt.-1) then istup(it)=-100 endif IDUP(IT)=IDUP(I) cc if(i.le.2.and.ida.eq.11) then if(i.le.2) then mothup(1,it)=0 mothup(2,it)=0 else mothup(1,it)=101 mothup(2,it)=102 endif ICOLUP(1,IT)=ICOLUP(1,I) ICOLUP(2,IT)=ICOLUP(2,I) C---COPY FROM A TO B. 5TH=1/(3-MTM), 6TH=PT, 7TH=ETA, 8TH=PHI, 9TH=PT**2 DO J=1,4 PUPA(J,IT)=PUP(J,I) ENDDO PUPA(10,IT)=0D0 PUPA(12,IT)=0D0 if(istup(i).eq.-1.or.istup(i).eq.1) then PUPA(9,IT)=PUPA(1,IT)**2+PUPA(2,IT)**2 PUPA(9,IT)=MAX(1D-6,PUPA(9,IT)) IF(PUPA(9,IT).GT.1D-6) THEN PUPA(8,IT)=DATAN2(pupa(2,it),pupa(1,it)) ELSE PUPA(8,IT)=1D-6 ENDIF PUPA(5,IT)=1D0/DSQRT(PUPA(9,IT)+PUPA(3,IT)**2) PPLUS=(PUPA(4,IT)+PUPA(3,IT)) PMINS=(PUPA(4,IT)-PUPA(3,IT)) IF(ABS(PPLUS).LT.1D-6) THEN PUPA(7,IT)=-1D6 ELSEIF(ABS(PMINS).LT.1D-6) THEN PUPA(7,IT)=1D6 ELSE PUPA(7,IT)=.5D0*DLOG(PPLUS/PMINS) ENDIF PUPA(6,IT)=SQRT(PUPA(9,IT)) ENDIF 5 CONTINUE C/* 2x2 CKM matrix c IFAM=0 c IF(ABS(ABS(QBOSON/3.D0)-1D0).LT.1D-6) IFAM=1 C/* Avoid obvious non-matches DO I=IOKBEG,IOKEND IDA=ABS(ISTUP(I)) IOK(I)=0 IF(IDA.NE.1) IOK(I)=1 IPICK(I)=0 IDA=ABS(IDUP(I)) IF(IDA.GE.11.AND.IDA.LE.16) IOK(I)=1 ENDDO if(iverb.eq.1) print*,' nchg = ',nchg CCCCC NL=101 NH=100+NUPP NL0=NL NH0=NH C/* Check on Number of f-fbar pairs (for e+e-) IPAIR=0 NPAIR=0 DO 10 I=NL,NH-1 IDA=IDUP(I)*ISTUP(I) DO 20 J=I+1,NH IDB=IDUP(J)*ISTUP(J) IF(IPICK(J).NE.0) GOTO 20 IF(IDA.EQ.-IDB.AND.ABS(IDA).LE.6.AND.IDA.NE.0) THEN NPAIR=NPAIR+1 if(iverb.eq.1) then print*,' pairs ',i,j,ida,idb,npair endif IPICK(I)=NPAIR IPICK(J)=NPAIR GOTO 10 ENDIF 20 CONTINUE 10 CONTINUE IDH2=NL+1 C.....Must match to NMATCH partons c.......For Z-> jets, match down to 2 partons NSOFAR=0 ccc R2SFAR=0D0 C...Find two closest jets. 300 R2MIN=R2MAX IF(NSOFAR.EQ.NMATCH) THEN if(iverb.eq.1) then print*,' nosfar = nmatch ',nsofar,nmatch endif GOTO 402 ENDIF n_fermions=0 DO 30 I=NL,NH IF(IOK(I).NE.0) GOTO 30 is1=abs(istup(i)) if(is1.le.6D0) n_fermions=n_fermions+1 30 CONTINUE C/* set up for a core process of q qbar <-> boson IMIN=0 IKNT=0 DO 400 ITRY1=NL,NH-1 IS1=ABS(ISTUP(ITRY1)) ID1=IDUP(ITRY1) IF(IOK(ITRY1).NE.0) GOTO 400 cc IF(NSOFAR.NE.0.and.IS1.GE.NSOFAR+1) GOTO 400 IC1=ICOLUP(1,ITRY1) JC1=ICOLUP(2,ITRY1) IGL1=IC1*JC1 IF(ISTUP(ITRY1).LT.0) THEN C.......make all colors "outgoing" ITMP=IC1 IC1=JC1 JC1=ITMP ENDIF DO 390 ITRY2=MAX(IDH2+1,ITRY1+1),NH cc IS2=ABS(ISTUP(ITRY2)) ID2=IDUP(ITRY2) IF(IOK(ITRY2).NE.0) GOTO 390 cc IF(NSOFAR.NE.0.and.IS2.GE.NSOFAR+1) GOTO 390 IC2=ICOLUP(1,ITRY2) JC2=ICOLUP(2,ITRY2) IGL2=IC2*JC2 IF(ISTUP(ITRY2).LT.0) THEN C.......make all colors "outgoing" ITMP=IC2 IC2=JC2 JC2=ITMP ENDIF isig=1 if(istup(itry1).lt.0 .or. istup(itry2).lt.0) isig=-1 R2=R2M(ITRY1,ITRY2,ISIG) if(iverb.eq.1) then print*,' r2(',itry1,itry2,')=',r2 endif IF(ISTUP(ITRY1).EQ.-1.and.ISTUP(ITRY2).EQ.-1) THEN R2=R2MAX C.......two fermions ELSEIF(IGL1.EQ.0.AND.IGL2.EQ.0) THEN if(iverb.eq.1) then print*,itry1,itry2,check_color,ipair,npair endif C/* do not allow a cluster to only gluons IF(ISIG.EQ.-1.and.(id1.eq.id2)) then if(n_fermions.eq.2) R2=R2MAX ELSEIF(ISIG.EQ.1.and.(id1.eq.-id2)) then if(n_fermions.eq.2) R2=R2MAX ELSEIF(Isig.eq.-1.and.(id1.ne.id2)) then R2=R2MAX ELSEIF(Isig.eq.1.and.(id1.ne.-id2)) then R2=R2MAX ENDIF IF(R2.NE.R2MAX) THEN if(check_color) then IF(IC1.NE.0.AND.JC2.NE.0) THEN JCOLa=IC1 ICOLa=JC2 IF(IC1.EQ.JC2) R2=R2MAX ELSEIF(JC1.NE.0.AND.IC2.NE.0) THEN JCOLa=IC2 ICOLa=JC1 IF(JC1.EQ.IC2) R2=R2MAX ELSE R2=R2MAX ENDIF endif ENDIF IF(R2.NE.R2MAX) THEN C........Check that the emitted gluon makes sense C.........no color singlets IF(check_color) THEN K=NL-1 if(iverb.eq.1) print*,' checking color ',icola,jcola do while(k.lt.nh) k=k+1 if(iverb.eq.1) print*,k,idup(k),iok(k) if(idup(k).eq.21.and.iok(k).eq.0) then lcola=icolup(1,k) lcolb=icolup(2,k) if(istup(k).eq.-1) then lcola=icolup(2,k) lcolb=icolup(1,k) endif if(iverb.eq.1) print*,' * ',lcola,' to ', $ icola,' : ',lcolb,' to ',jcola if(lcola.eq.icola.and.lcolb.eq.jcola) R2=R2MAX endif enddo ENDIF C.........Type=1 is a e+e- collider IF(TYPE.EQ.1.and.IPAIR.EQ.NPAIR-1) R2=R2MAX IF(type.eq.1.and.R2.NE.R2MAX.AND.(NPAIR.GT.1)) THEN C.........Is there another pair with the same flavor??? IST=0 DO 88 I=NL0,NH0 IF(I.EQ.ITRY1) GOTO 88 IF(IOK(I).NE.0) GOTO 88 IF(IDUP(I).EQ.IDUP(ITRY1)) THEN IST=IST+1 IOG=I ENDIF 88 CONTINUE C..........Only concerned if there is just one pair like it IF(IST.EQ.1) THEN DO 77 I=NL0,NH0 IF(I.EQ.ITRY1.OR.I.EQ.ITRY2) GOTO 77 IF(IOK(I).NE.0) GOTO 77 IF(IDUP(I).EQ.-IDUP(IOG)) THEN IF(IDUP(I).GT.0.and.check_color) THEN IF(ICOLUP(1,I).EQ.ICOLUP(2,IOG)) R2=R2MAX ELSE IF(ICOLUP(2,I).EQ.ICOLUP(1,IOG)) R2=R2MAX ENDIF GOTO 78 ENDIF 77 CONTINUE 78 CONTINUE ENDIF ENDIF ENDIF C.......gluon+fermion (case i) ELSEIF(IGL1.NE.0.AND.IGL2.EQ.0.and.check_color) THEN IF(IC2.NE.JC1 .and. JC2.NE.IC1) THEN R2=R2MAX ENDIF C.......gluon+fermion (case ii) ELSEIF(IGL2.NE.0.AND.IGL1.EQ.0.and.check_color) THEN IF(IC2.NE.JC1 .and. JC2.NE.IC1) THEN R2=R2MAX ENDIF C.......gluon+gluon ELSEIF(IGL1.NE.0.AND.IGL2.NE.0.and.check_color) THEN IF(IC1.EQ.JC2.AND.JC1.EQ.IC2) THEN R2=R2MAX ELSEIF(IC1.NE.JC2.AND.JC1.NE.IC2) THEN R2=R2MAX ENDIF IF(R2.NE.R2MAX.and.check_color) THEN IF(ISTUP(ITRY1).GT.0.and.ISTUP(ITRY2).GT.0) THEN C........Don't form a color singlet system with another gluon K=NL-1 do while(k.lt.nh) k=k+1 if(idup(k).eq.21) then if(icolup(1,k).eq.icola.and.icolup(2,k).eq.jcola) R2=R2MAX endif enddo ENDIF ENDIF ENDIF C........No top quarks in PDF IF(ISTUP(ITRY1)*ISTUP(ITRY2).LT.0) THEN IF(ABS(IDUP(ITRY1)).EQ.6.OR.ABS(IDUP(ITRY2)).EQ.6) R2=R2MAX ENDIF IKNT=IKNT+1 D2(IKNT)=R2 R2=ABS(R2) C.......Allow for somewhat non-ordered showers. Cxxxxx IF(R2.LT.R2SFAR) R2=R2MAX IF(R2.LT.R2MIN) THEN R2MIN=R2 IMIN=IKNT ENDIF INDEX(IKNT)=ITRY1 JNDEX(IKNT)=ITRY2 390 CONTINUE 400 CONTINUE NSOFAR=NSOFAR+1 if(imin.eq.0.and.inext.eq.0) then inext=1 if(iflag.lt.5) then print*,' imin = ',imin,iflag,inext call pylist(7) print*,' d2 = ',(d2(j),j=1,iknt) print*,' in = ',(index(j),j=1,iknt) print*,' jn = ',(jndex(j),j=1,iknt) print*,' imin = ',imin,r2min print*,' no color check ' do i=nl,nh write(*,588) i,istup(i),idup(i),mothup(1,i),mothup(2,i), $ icolup(1,i),icolup(2,i), $ pup(5,i),pup(1,i),pup(2,i),pup(3,i),pup(4,i) enddo 588 format(7I4,5(2X,E10.4)) endif check_color=.false. iflag=iflag+1 goto 100 endif if(inext.eq.1.and.imin.eq.0) then inext=0 iskip=1 return endif if(iverb.eq.1) then print*,' d2 = ',(d2(j),j=1,iknt) print*,' in = ',(index(j),j=1,iknt) print*,' jn = ',(jndex(j),j=1,iknt) print*,' imin = ',imin,r2min endif do j=1,iknt if(d2(j).lt.0) print*,' D2 < 0 ??? ' enddo I=IMIN cc R2SFAR=abs(D2(IMIN)) if(i.eq.0) GOTO 402 IMIN2=JNDEX(I) IMIN1=INDEX(I) IOK(IMIN1)=1 IOK(IMIN2)=1 IF(IMIN2.LT.IMIN1) THEN ITMP=IMIN1 IMIN1=IMIN2 IMIN2=ITMP ENDIF IC1=ICOLUP(1,IMIN1) JC1=ICOLUP(2,IMIN1) IC2=ICOLUP(1,IMIN2) JC2=ICOLUP(2,IMIN2) IF(ISTUP(IMIN1).GT.0.and.ISTUP(IMIN2).GT.0) THEN ISTUP(IMIN1)=1+NSOFAR ISTUP(IMIN2)=1+NSOFAR NH=NH+1 IOK(NH)=0 ISTUP(NH)=1 MOTHUP(1,IMIN1)=NH MOTHUP(2,IMIN1)=IMIN2 MOTHUP(1,IMIN2)=NH MOTHUP(2,IMIN2)=IMIN1 MOTHUP(1,NH)=101 MOTHUP(2,NH)=102 DO J=1,4 PUPA(J,NH)=PUPA(J,IMIN1)+PUPA(J,IMIN2) ENDDO PUPA(11,NH)=PUPA(4,IMIN1)*PUPA(4,IMIN2)- $ (PUPA(1,IMIN1)*PUPA(1,IMIN2)+PUPA(2,IMIN1)*PUPA(2,IMIN2)+ $ PUPA(3,IMIN1)*PUPA(3,IMIN2)) C......Add them here IF (RECO.EQ.1) THEN C---VECTOR ADDITION PUPA(5,NH)=SQRT(PUPA(1,NH)**2+PUPA(2,NH)**2+PUPA(3,NH)**2) IF (PUPA(5,NH).EQ.0D0) THEN PUPA(5,NH)=1D0 ELSE PUPA(5,NH)=1D0/PUPA(5,NH) ENDIF PUPA(9,NH)=PUPA(1,NH)**2+PUPA(2,NH)**2 PUPA(7,NH)=PUPA(4,NH)**2-PUPA(3,NH)**2 IF (PUPA(7,NH).LE.EPS*PUPA(4,NH)**2) PUPA(7,NH)=PUPA(9,NH) IF (PUPA(7,NH).GT.0) THEN PUPA(7,NH)=LOG((PUPA(4,NH)+ABS(PUPA(3,NH)))**2/PUPA(7,NH))/2 IF (PUPA(7,NH).GT.ETAMAX) PUPA(7,NH)=ETAMAX+2 ELSE PUPA(7,NH)=ETAMAX+2 ENDIF PUPA(7,NH)=SIGN(PUPA(7,NH),PUPA(3,NH)) IF (PUPA(1,NH).NE.0.AND.PUPA(2,NH).NE.0) THEN PUPA(8,NH)=ATAN2(PUPA(2,NH),PUPA(1,NH)) ELSE PUPA(8,NH)=0 ENDIF ELSEIF (RECO.EQ.2) THEN C---PT WEIGHTED ETA-PHI ADDITION PT=PUPA(6,IMIN1)+PUPA(6,IMIN2) IF (PT.EQ.0) THEN PTT=1D0 ELSE PTT=1D0/PT ENDIF PUPA(7,NH)=(PUPA(6,IMIN1)*PUPA(7,IMIN1)+ $ PUPA(6,IMIN2)*PUPA(7,IMIN2))*PTT PUPA(8,NH)=KTMDPI(PUPA(8,IMIN1)+ $ PUPA(6,IMIN2)*PTT*KTMDPI(PUPA(8,IMIN2)-PUPA(8,IMIN1))) PUPA(6,NH)=PT PUPA(9,NH)=PT**2 ELSEIF (RECO.EQ.3) THEN C---PT**2 WEIGHTED ETA-PHI ADDITION PT=PUPA(9,IMIN1)+PUPA(9,IMIN2) IF (PT.EQ.0) THEN PTT=1D0 ELSE PTT=1D0/PT ENDIF PUPA(7,NH)=(PUPA(9,IMIN1)*PUPA(7,IMIN1)+ $ PUPA(9,IMIN2)*PUPA(7,IMIN2))*PTT PUPA(8,NH)=KTMDPI(PUPA(8,IMIN1)+ $ PUPA(9,IMIN2)*PTT*KTMDPI(PUPA(8,IMIN2)-PUPA(8,IMIN1))) PUPA(6,NH)=PUPA(6,IMIN1)+PUPA(6,IMIN2) PUPA(9,NH)=PUPA(6,NH)**2 ENDIF IF (ANGL.NE.1.AND.RECO.EQ.1) THEN ELSEIF (ANGL.EQ.1.AND.RECO.NE.1) THEN PUPA(1,NH)=PUPA(6,NH)*COS(PUPA(8,NH)) PUPA(2,NH)=PUPA(6,NH)*SIN(PUPA(8,NH)) PUPA(3,NH)=PUPA(6,NH)*SINH(PUPA(7,NH)) PUPA(4,NH)=PUPA(6,NH)*COSH(PUPA(7,NH)) IF (PUPA(4,NH).NE.0) THEN PUPA(5,NH)=1D0/PUPA(4,NH) ELSE PUPA(5,NH)=1D0 ENDIF ENDIF PUPA(10,NH)=sqrt(R2M(IMIN1,IMIN2,1)) ICLUS=ICLUS+1 RESKT(ICLUS)=PUPA(10,NH) C.......gluon -> quark+antiquark IF(IC1*JC1.EQ.0 .and. IC2*JC2.EQ.0) THEN IPAIR=IPAIR+1 IDUP(NH)=21 IF(IC1.EQ.0) THEN ICOLUP(1,NH)=IC2 ICOLUP(2,NH)=JC1 ELSE ICOLUP(1,NH)=IC1 ICOLUP(2,NH)=JC2 ENDIF C.......gluon -> gluon+gluon ELSEIF(IC1*JC1.NE.0 .and. IC2*JC2.NE.0) THEN IDUP(NH)=21 IF(IC1.EQ.JC2) THEN ICOLUP(1,NH)=IC2 ICOLUP(2,NH)=JC1 ELSE ICOLUP(1,NH)=IC1 ICOLUP(2,NH)=JC2 ENDIF C.......quark -> quark+gluon ELSEIF(IC1*JC1.EQ.0) THEN IDUP(NH)=IDUP(IMIN1) IF(JC1.EQ.0) THEN ICOLUP(1,NH)=IC2 ICOLUP(2,NH)=0 ELSE ICOLUP(1,NH)=0 ICOLUP(2,NH)=JC2 ENDIF C.......quark -> quark+gluon ELSEIF(IC2*JC2.EQ.0) THEN IDUP(NH)=IDUP(IMIN2) IF(JC2.EQ.0) THEN ICOLUP(1,NH)=IC1 ICOLUP(2,NH)=0 ELSE ICOLUP(1,NH)=0 ICOLUP(2,NH)=JC1 ENDIF C......bizarre configurations ELSE IDUP(NH)=IDUP(IMIN2) IF(JC2.EQ.0) THEN ICOLUP(1,NH)=IC1 ICOLUP(2,NH)=0 ELSE ICOLUP(1,NH)=0 ICOLUP(2,NH)=JC1 ENDIF ENDIF if(.not.check_color) then if(idup(nh).eq.21) then icolup(1,nh)=1 icolup(2,nh)=1 elseif(abs(idup(nh)).ge.1.and.abs(idup(nh)).le.6) then if(idup(nh).gt.0) then icolup(1,nh)=1 icolup(2,nh)=0 else icolup(2,nh)=1 icolup(1,nh)=0 endif else print*,' problem ' endif endif cc PUP(5,NH)=PUP(5,NH)-PMS(IDUP(NH))**2 ELSE NL=NL-1 IOK(NL)=0 ISTUP(IMIN1)=SIGN((1+NSOFAR),ISTUP(IMIN1)) ISTUP(IMIN2)=SIGN((1+NSOFAR),ISTUP(IMIN2)) ISTUP(NL)=-1 MOTHUP(1,IMIN1)=NL MOTHUP(2,IMIN1)=IMIN2 MOTHUP(1,IMIN2)=NL MOTHUP(2,IMIN2)=IMIN1 MOTHUP(1,NL)=0 MOTHUP(2,NL)=0 IF(ISTUP(IMIN1).LT.0) THEN DO J=1,4 PUPA(J,NL)=PUPA(J,IMIN1)-PUPA(J,IMIN2) ENDDO ELSE DO J=1,4 PUPA(J,NL)=-PUPA(J,IMIN1)+PUPA(J,IMIN2) ENDDO ENDIF PUPA(11,NL)=PUPA(4,IMIN1)*PUPA(4,IMIN2)- $ (PUPA(1,IMIN1)*PUPA(1,IMIN2)+PUPA(2,IMIN1)*PUPA(2,IMIN2)+ $ PUPA(3,IMIN1)*PUPA(3,IMIN2)) C......Add them here IF (RECO.EQ.1) THEN C---VECTOR ADDITION PUPA(5,NL)=SQRT(PUPA(1,NL)**2+PUPA(2,NL)**2+PUPA(3,NL)**2) IF (PUPA(5,NL).EQ.0D0) THEN PUPA(5,NL)=1D0 ELSE PUPA(5,NL)=1D0/PUPA(5,NL) ENDIF ELSEIF (RECO.EQ.2) THEN C---PT WEIGHTED ETA-PHI ADDITION PT=PUPA(6,IMIN1)+PUPA(6,IMIN2) IF (PT.EQ.0) THEN PTT=1D0 ELSE PTT=1D0/PT ENDIF PUPA(7,NL)=(PUPA(6,IMIN1)*PUPA(7,IMIN1)+ $ PUPA(6,IMIN2)*PUPA(7,IMIN2))*PTT PUPA(8,NL)=KTMDPI(PUPA(8,IMIN1)+ $ PUPA(6,IMIN2)*PTT*KTMDPI(PUPA(8,IMIN2)-PUPA(8,IMIN1))) PUPA(6,NL)=PT PUPA(9,NL)=PT**2 ELSEIF (RECO.EQ.3) THEN C---PT**2 WEIGHTED ETA-PHI ADDITION PT=PUPA(9,IMIN1)+PUPA(9,IMIN2) IF (PT.EQ.0) THEN PTT=1D0 ELSE PTT=1D0/PT ENDIF PUPA(7,NL)=(PUPA(9,IMIN1)*PUPA(7,IMIN1)+ $ PUPA(9,IMIN2)*PUPA(7,IMIN2))*PTT PUPA(8,NL)=KTMDPI(PUPA(8,IMIN1)+ $ PUPA(9,IMIN2)*PTT*KTMDPI(PUPA(8,IMIN2)-PUPA(8,IMIN1))) PUPA(6,NL)=PUPA(6,IMIN1)+PUPA(6,IMIN2) PUPA(9,NL)=PUPA(6,NL)**2 ENDIF cc IF (ANGL.NE.1.AND.RECO.EQ.1) THEN IF (RECO.EQ.1) THEN PUPA(9,NL)=PUPA(1,NL)**2+PUPA(2,NL)**2 PUPA(7,NL)=PUPA(4,NL)**2-PUPA(3,NL)**2 IF (PUPA(7,NL).LE.EPS*PUPA(4,NL)**2) PUPA(7,NL)=PUPA(9,NL) IF (PUPA(7,NL).GT.0) THEN PUPA(7,NL)=LOG((PUPA(4,NL)+ABS(PUPA(3,NL)))**2/PUPA(7,NL))/2 IF (PUPA(7,NL).GT.ETAMAX) PUPA(7,NL)=ETAMAX+2D0 ELSE PUPA(7,NL)=ETAMAX+2D0 ENDIF PUPA(7,NL)=SIGN(PUPA(7,NL),PUPA(3,NL)) IF (PUPA(1,NL).NE.0D0.AND.PUPA(2,NL).NE.0D0) THEN PUPA(8,NL)=ATAN2(PUPA(2,NL),PUPA(1,NL)) ELSE PUPA(8,NL)=0D0 ENDIF ELSEIF (ANGL.EQ.1.AND.RECO.NE.1) THEN PUPA(1,NL)=PUPA(6,NL)*COS(PUPA(8,NL)) PUPA(2,NL)=PUPA(6,NL)*SIN(PUPA(8,NL)) PUPA(3,NL)=PUPA(6,NL)*SINH(PUPA(7,NL)) PUPA(4,NL)=PUPA(6,NL)*COSH(PUPA(7,NL)) IF (PUPA(4,NL).NE.0D0) THEN PUPA(5,NL)=1D0/PUPA(4,NL) ELSE PUPA(5,NL)=1D0 ENDIF ENDIF PUPA(10,NL)=sqrt(R2M(IMIN1,IMIN2,-1)) ICLUS=ICLUS+1 RESKT(ICLUS)=PUPA(10,NL) C......Clustering of initial state particles PUPA(9,NL)=0D0 C.......quark -> quark + (gluon) IF(IC1*JC1.EQ.0 .and. IC2*JC2.EQ.0) THEN IDUP(NL)=21 IF(IC1.EQ.0) THEN ICOLUP(1,NL)=JC2 ICOLUP(2,NL)=JC1 ELSE ICOLUP(1,NL)=IC1 ICOLUP(2,NL)=IC2 ENDIF C.......gluon -> gluon+ (gluon) ELSEIF(IC1*JC1.NE.0 .and. IC2*JC2.NE.0) THEN IDUP(NL)=21 IF(IC1.EQ.IC2) THEN ICOLUP(1,NL)=JC2 ICOLUP(2,NL)=JC1 ELSE ICOLUP(1,NL)=IC1 ICOLUP(2,NL)=IC2 ENDIF C.......quark -> gluon + (quark) ELSEIF(IC1*JC1.EQ.0) THEN IDUP(NL)=IDUP(IMIN1) IF(JC1.EQ.0) THEN ICOLUP(1,NL)=JC2 ICOLUP(2,NL)=0 ELSE ICOLUP(2,NL)=IC2 ICOLUP(1,NL)=0 ENDIF C.......gluon -> quark + (antiquark) ELSEIF(IC2*JC2.EQ.0) THEN IDUP(NL)=-IDUP(IMIN2) IF(JC2.EQ.0) THEN ICOLUP(2,NL)=JC1 ICOLUP(1,NL)=0 ELSE ICOLUP(1,NL)=IC1 ICOLUP(2,NL)=0 ENDIF C.......bizarre configurations ELSE IDUP(NL)=-IDUP(IMIN2) IF(JC2.EQ.0) THEN ICOLUP(2,NL)=JC1 ICOLUP(1,NL)=0 ELSE ICOLUP(1,NL)=IC1 ICOLUP(2,NL)=0 ENDIF ENDIF ccc PUP(5,NL)=PUP(5,NL)-PMS(IDUP(NL))**2 if(.not.check_color) then if(idup(nl).eq.21) then icolup(1,nl)=1 icolup(2,nl)=1 elseif(abs(idup(nl)).ge.1.and.abs(idup(nl)).le.6) then if(idup(nl).gt.0) then icolup(1,nl)=1 icolup(2,nl)=0 else icolup(2,nl)=1 icolup(1,nl)=0 endif else print*,' problem ' endif endif ENDIF 401 CONTINUE GOTO 300 402 CONTINUE if(iverb.eq.1) then print*,' exiting here ' endif C/* Pick the color flow for final state partons IF NECESSARY IF(.not.check_color) THEN C/* Initialize all color to 0 do i=nl,nh icolup(1,i)=0 icolup(2,i)=0 enddo C/* Match to the "hardest" scatterers imatch=1 imin1=0 imin2=0 j=nl do while(imin1*imin2.eq.0.and.j.le.nh) if(abs(istup(j)).eq.imatch) then if(imin1.eq.0) then imin1=j else imin2=j endif endif j=j+1 enddo C/* Assumption that the initial state is a color singlet IC2=1 JC2=2 C/* Either a gluon IF(IDUP(IMIN1).EQ.21) THEN ICID=ICID+1 ICOLUP(1,IMIN1)=ICID ICID=ICID+1 ICOLUP(2,IMIN1)=ICID IF(IDUP(IMIN2).EQ.21) THEN ICOLUP(IC2,IMIN2)=ICOLUP(2,IMIN1) ICID=ICID+1 ICOLUP(JC2,IMIN2)=ICID ELSEIF(IDUP(IMIN2).GE.1.and.IDUP(IMIN2).LE.6) THEN ICOLUP(1,IMIN2)=ICOLUP(JC2,IMIN1) ELSEIF(IDUP(IMIN2).GE.-6.and.IDUP(IMIN2).LE.-1) THEN ICOLUP(2,IMIN2)=ICOLUP(IC2,IMIN1) ENDIF C/* Or a quark ELSEIF(IDUP(IMIN1).GE.1.and.IDUP(IMIN1).LE.6) THEN ICID=ICID+1 ICOLUP(1,IMIN1)=ICID IF(IDUP(IMIN2).EQ.21) THEN ICOLUP(JC2,IMIN2)=ICOLUP(1,IMIN1) ICID=ICID+1 ICOLUP(IC2,IMIN2)=ICID ELSEIF(IDUP(IMIN2).GE.-6.and.IDUP(IMIN2).LE.-1) THEN ICOLUP(JC2,IMIN2)=ICID ENDIF C/* Or an anti-quark ELSEIF(IDUP(IMIN1).GE.-6.and.IDUP(IMIN1).LE.-1) THEN ICID=ICID+1 ICOLUP(2,IMIN1)=ICID IF(IDUP(IMIN2).EQ.21) THEN ICOLUP(IC2,IMIN2)=ICOLUP(2,IMIN1) ICID=ICID+1 ICOLUP(JC2,IMIN2)=ICID ELSEIF(IDUP(IMIN2).GE.1.and.IDUP(IMIN2).LE.6) THEN ICOLUP(IC2,IMIN2)=ICID ENDIF ENDIF ID1=IMIN1 ID2=IMIN2 IP1=ID1 IP2=ID2 C/* start a tree iknt=0 imatch=0 do while(imatch.eq.0.and.iknt.lt.50) imin1=0 j=nl do while(imin1.eq.0.and.j.le.nh) C/* change: hard-wired if(mothup(1,j).eq.id1.and.(j.ne.103.and.j.ne.104)) then if(icolup(1,j).eq.0.AND.icolup(2,j).eq.0) THEN imin1=j call mapcolor(id1,imin1,icid) elseif(icolup(1,mothup(2,j)).eq.0.AND. $ icolup(2,mothup(2,j)).EQ.0) then imin1=mothup(2,j) call mapcolor(id1,imin1,icid) endif endif j=j+1 enddo iknt=iknt+1 C/* didn't find a branching if(imin1.eq.0) then if(id1.eq.ip1) imatch=1 id1=mothup(1,id1) if(id1.eq.0) imatch=1 else id1=imin1 endif enddo C/* Finish the 2nd side ip1=ip2 id1=id2 imatch=0 iknt=0 do while(imatch.eq.0.and.iknt.lt.50) imin1=0 j=nl do while(imin1.eq.0.and.j.le.nh) C/* change: hard-wired if(mothup(1,j).eq.id1.and.(j.ne.103.and.j.ne.104)) then if(icolup(1,j).eq.0.AND.icolup(2,j).eq.0) THEN imin1=j call mapcolor(id1,imin1,icid) elseif(icolup(1,mothup(2,j)).eq.0.AND. $ icolup(2,mothup(2,j)).EQ.0) then imin1=mothup(2,j) call mapcolor(id1,imin1,icid) endif endif j=j+1 enddo C/* didn't find a branching if(imin1.eq.0) then if(id1.eq.ip1) imatch=1 id1=mothup(1,id1) if(id1.eq.0) imatch=1 else id1=imin1 endif iknt=iknt+1 enddo ENDIF DO I=NL,NH DO J=1,4 PUP(J,I)=PUPA(J,I) ENDDO PUP(5,I)=PUPA(10,I) PUPA(12,I)=PUPA(4,I)**2-PUPA(1,I)**2-PUPA(2,I)**2- $ PUPA(3,I)**2 PUPA(12,I)=SQRT(ABS(PUPA(12,I))) PUPA(11,I)=SQRT(2D0*ABS(PUPA(11,I))) IF(PUP(5,I).GT.1.01D0*PUPA(12,I)) THEN PUPA(12,I)=PUP(5,I) ENDIF if(.not.check_color.and.(i.ge.101.and.i.le.100+NUP)) then icolup(1,i-100)=icolup(1,i) icolup(2,i-100)=icolup(2,i) endif ENDDO RETURN END C/* subroutine mapcolor(imo,ida,icid) implicit none integer imo,ida,icid,ipa,imin1,imin2 integer ic1,ic2,nl,jc1,jc2 C.....Global C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ C/* Assumption that the initial state is a color singlet ipa=mothup(2,ida) if(icolup(1,ipa).eq.0.and.icolup(2,ipa).eq.0) then ic2=1 jc2=2 imin1=ida imin2=imo IF(IDUP(IMIN1).EQ.21) THEN IF(IDUP(IMIN2).EQ.21) THEN ICOLUP(IC2,IMIN1)=ICOLUP(1,IMIN2) ICID=ICID+1 ICOLUP(JC2,IMIN1)=ICID ELSEIF(IDUP(IMIN2).GE.1.and.IDUP(IMIN2).LE.6) THEN ICOLUP(IC2,IMIN1)=ICOLUP(1,IMIN2) ICID=ICID+1 ICOLUP(JC2,IMIN1)=ICID ELSEIF(IDUP(IMIN2).GE.-6.and.IDUP(IMIN2).LE.-1) THEN ICOLUP(JC2,IMIN1)=ICOLUP(2,IMIN2) ICID=ICID+1 ICOLUP(IC2,IMIN1)=ICID ENDIF C/* Or a quark ELSEIF(IDUP(IMIN1).GE.1.and.IDUP(IMIN1).LE.6) THEN IF(IDUP(IMIN2).EQ.21) THEN ICOLUP(1,IMIN1)=ICOLUP(IC2,IMIN2) ELSEIF(IDUP(IMIN2).GE.-6.and.IDUP(IMIN2).LE.-1) THEN ICID=ICID+1 ICOLUP(1,IMIN1)=ICID ELSEIF(IDUP(IMIN2).GE.1.and.IDUP(imin2).le.6) then ICID=ICID+1 ICOLUP(1,IMIN1)=ICID ENDIF C/* Or an anti-quark ELSEIF(IDUP(IMIN1).GE.-6.and.IDUP(IMIN1).LE.-1) THEN IF(IDUP(IMIN2).EQ.21) THEN ICOLUP(2,IMIN1)=ICOLUP(JC2,IMIN2) ELSEIF(IDUP(IMIN2).GE.-6.and.IDUP(IMIN2).LE.-1) THEN ICID=ICID+1 ICOLUP(2,IMIN1)=ICID ELSEIF(IDUP(IMIN2).GE.1.and.IDUP(imin2).le.6) then ICID=ICID+1 ICOLUP(2,IMIN1)=ICID ENDIF ENDIF C/* color flow is predeterimined else imin1=imo imin2=ipa IC2=ICOLUP(1,IMIN2) JC2=ICOLUP(2,IMIN2) IC1=ICOLUP(1,IMIN1) JC1=ICOLUP(2,IMIN1) IF(.true.) then NL=IDA C.......quark -> quark + (gluon) IF(IC1*JC1.EQ.0 .and. IC2*JC2.EQ.0) THEN C/* e,g IF(IC1.EQ.0) THEN ICOLUP(2,NL)=JC2 ICOLUP(1,NL)=JC1 ELSE ICOLUP(2,NL)=IC1 ICOLUP(1,NL)=IC2 ENDIF C.......gluon -> gluon+ (gluon) ELSEIF(IC1*JC1.NE.0 .and. IC2*JC2.NE.0) THEN C/* a IF(IC1.EQ.IC2) THEN ICOLUP(1,NL)=JC1 ICOLUP(2,NL)=JC2 ELSE ICOLUP(1,NL)=IC2 ICOLUP(2,NL)=IC1 ENDIF C.......quark -> gluon + (quark) ELSEIF(IC1*JC1.EQ.0) THEN C/* d, f IF(JC1.EQ.0) THEN ICOLUP(2,NL)=JC2 ICOLUP(1,NL)=0 ELSE ICOLUP(1,NL)=IC2 ICOLUP(2,NL)=0 ENDIF C.......gluon -> quark + (antiquark) ELSEIF(IC2*JC2.EQ.0) THEN C/* IF(JC2.EQ.0) THEN ICOLUP(1,NL)=JC1 ICOLUP(2,NL)=0 ELSE ICOLUP(2,NL)=IC1 ICOLUP(1,NL)=0 ENDIF ELSE ENDIF IF(ISTUP(IMO).GT.0) THEN IC1=ICOLUP(1,NL) JC1=ICOLUP(2,NL) ICOLUP(1,NL)=JC1 ICOLUP(2,NL)=IC1 ENDIF ENDIF ENDIF RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE REWEIGHT(NLO,NHI,FACT,IVERB,ilast) IMPLICIT NONE REAL*8 FACT INTEGER NLO,NHI,IVERB,ilast C.....Global C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ C Common block for extra info. INTEGER NMATCH REAL*8 Q0RES,QXRES,QHARD,PASSKT(20),SCALEQ(20),THORD(20) REAL*8 KTLIMIT COMMON/IPASS/Q0RES,QXRES,QHARD,PASSKT,SCALEQ,THORD,KTLIMIT,NMATCH SAVE /IPASS/ INTEGER ICOLL REAL*8 PUPA(12,MAXNUP) COMMON/HEPA/PUPA,ICOLL SAVE /HEPA/ INTEGER ISHWVAR COMMON/SHWPSS/ISHWVAR SAVE /SHWPSS/ INTEGER TYPE,ANGL,MONO,RECO COMMON/CLSTR/TYPE,ANGL,MONO,RECO SAVE /CLSTR/ C.....Local integer ilep,imcheck integer im1,im2,i,j,id,ida,ima,ipa real*8 rktmin,qnow integer iskip common/skip/iskip save /skip/ C fact=1D0 Cxxxxx im1=0 im2=0 imcheck=-1 if(type.eq.1) imcheck=1 do i=nlo,nhi if(istup(i).eq.imcheck) then if(im1.eq.0) then im1=i elseif(im2.eq.0) then im2=i else iskip=1 return endif endif enddo CMRENNA if(im1*im2.eq.0) then print*,' found NO mothers ',im1,im2 iskip=1 return endif if((icolup(1,im1).ne.icolup(2,im2)).AND. $ (icolup(2,im1).ne.icolup(1,im2))) then print*,' not a color singlet ' print*,im1,im2,icolup(1,im1),icolup(2,im1),icolup(1,im2), $ icolup(2,im2) endif Cxxxxx C.....Passkt values are used for vetoing ISR and FSR C........Set default values for the kt-clus values do i=1,nup passkt(i)=0D0 scaleq(i)=0D0 enddo scaleq(1)=qhard scaleq(2)=qhard ilep=0 do 55 i=101,100+NUP id=idup(i) ida=abs(id) if(ida.ne.21.and.(ida.lt.1.or.ida.gt.6)) then if(ida.ge.11.and.ida.le.16) then ilep=ilep+1 if(ilep.eq.1) then qnow=0D0 do j=1,4 qnow=qnow+sign((pup(j,i)+pup(j,i+1))**2,dble(j)-3.5D0) enddo else ilep=0 endif scaleq(i-100)=0D0 else scaleq(i-100)=qhard endif passkt(i-100)=1D4 goto 55 endif ima=mothup(1,i) ipa=mothup(2,i) if(ima.eq.0) then scaleq(i-100)=qhard passkt(i-100)=qhard goto 55 endif scaleq(i-100)=pupa(ISHWVAR,ima) if(ima.lt.101) goto 55 if(ida.eq.21) then if(idup(ima).eq.21) then if(pup(4,i).gt.pup(4,ipa)) then scaleq(i-100)=pupa(ISHWVAR,mothup(1,ima)) endif endif else do while(ima.ne.0.and.idup(ima).ne.21 .and. ima.ne.101) ima=mothup(1,ima) enddo if(ima.eq.0) then scaleq(i-100)=qhard elseif(mothup(1,ima).eq.0) then scaleq(i-100)=qhard elseif(idup(ima).eq.21) then scaleq(i-100)=pupa(ISHWVAR,mothup(1,ima)) else scaleq(i-100)=qhard endif endif ipa=mothup(1,i) if(ipa.eq.0) then passkt(i-100)=qhard elseif(mothup(1,ipa).eq.0) then passkt(i-100)=qhard else passkt(i-100)=Pupa(10,ipa) endif if(iverb.eq.1) write(*,56) ' i ',i,ida,ima, $ scaleq(i-100),passkt(i-100) 56 format(A3,3I6,2G12.6) 55 continue if(iverb.eq.1) print*,' scaleq ',scaleq if(iverb.eq.1) print*,' passkt ',passkt RKTMIN=1D6 do i=1,nup passkt(i)=MAX(Passkt(i),QXRES) if(ilast.eq.1) passkt(i)=MAX(Passkt(i),RKTMIN) if(scaleq(i).le.qxres.and. $ (abs(idup(i)).lt.11.or.abs(idup(i)).gt.16)) scaleq(i)=rktmin enddo if(iverb.eq.1) print*,' passkt ',passkt if(iverb.eq.1) print*,' fact before ISR ',fact if(iverb.eq.1) print*,' scaleq ',scaleq return END FUNCTION R2M(I1,I2,IS) IMPLICIT NONE REAL*8 R2M INTEGER I1,I2,I3,IS REAL*8 PA,PB C New standard event common INTEGER MAXNUP PARAMETER (MAXNUP=500) REAL*8 PUPA(12,MAXNUP) INTEGER ICOLL COMMON/HEPA/PUPA,ICOLL SAVE /HEPA/ INTEGER TYPE,ANGL,MONO,RECO COMMON/CLSTR/TYPE,ANGL,MONO,RECO SAVE /CLSTR/ INTEGER IKTTOPT COMMON/KTSWITCH/IKTTOPT SAVE /KTSWITCH/ C---COPY FROM A TO B. 5TH=1/(3-MTM), 6TH=PT, 7TH=ETA, 8TH=PHI, 9TH=PT**2 REAL*8 COSTH REAL*8 R,ESQ,ETA,PHI C cccccc ANGL=ICOLL IF(IS.GT.0) THEN IF (ANGL.EQ.1) THEN PA=1D0/DSQRT(PUPA(1,I1)**2+PUPA(2,I1)**2+PUPA(3,I1)**2) PB=1D0/DSQRT(PUPA(1,I2)**2+PUPA(2,I2)**2+PUPA(3,I2)**2) R=2D0*(1D0-(PUPA(1,I1)*PUPA(1,I2)+PUPA(2,I1)*PUPA(2,I2)+ $ PUPA(3,I1)*PUPA(3,I2))*(PA*PB)) C/* IF(IKTTOPT.EQ.1) THEN ESQ=(PUPA(4,I1)*PUPA(4,I2)/(PUPA(4,I1)+PUPA(4,I2)))**2 ELSE ESQ=MIN(PUPA(4,I1),PUPA(4,I2))**2 ENDIF C/* ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN ETA=PUPA(7,I1)-PUPA(7,I2) PHI=PUPA(8,I1)-PUPA(8,I2) IF (ANGL.EQ.2) THEN R=ETA**2+PHI**2 ELSE R=2D0*(COSH(ETA)-COS(PHI)) ENDIF ESQ=MIN(PUPA(9,I1),PUPA(9,I2)) ENDIF R2M=ESQ*R ELSE IF(MIN(PUPA(9,I2),PUPA(9,I1)).GT.1D-6) THEN print*,' kt measure is wrong for isr ' ENDIF IF(PUPA(9,I2).GT.PUPA(9,I1)) THEN I3=I2 ELSE I3=I1 ENDIF IF (ANGL.EQ.1) THEN COSTH=PUPA(3,I3)*PUPA(5,I3) IF (TYPE.EQ.2) THEN COSTH=-COSTH ELSEIF (TYPE.EQ.4) THEN COSTH=ABS(COSTH) ENDIF R=2D0*(1D0-COSTH) C---IF CLOSE TO BEAM, USE APPROX 2*(1-COS(THETA))=SIN**2(THETA) IF (R.LT.1D-4) R=(PUPA(1,I3)**2+PUPA(2,I3)**2)*PUPA(5,I3)**2 R2M=PUPA(4,I3)**2*R ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN R2M=PUPA(9,I3) ENDIF C......Max, because one of these should be zero! R2M=R2M-SIGN(1D-3,PUPa(3,I1)*PUPa(3,I2)) ENDIF IF(R2M.LT.0D0) THEN print*,' R2M < 0 ',R2M,I1,I2,IS endif C/* R2M=MAX(R2M,Q0RES**2) RETURN END FUNCTION THANG(P1,P2) IMPLICIT NONE REAL*8 THANG,P1(4),P2(4),E1,E2,XI INTEGER I THANG=4D0*ATAN(1D0) C/* E1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2) E2=SQRT(P2(1)**2+P2(2)**2+P2(3)**2) XI=1D0 DO I=1,3 XI=XI-P1(I)*P2(I)/(E1*E2) ENDDO IF(XI.GT.2D0) XI=2D0 IF(XI.LT.0) XI=0D0 THANG=ACOS(1D0-XI) IF(THANG.LT.0D0) THANG=THANG+8D0*ATAN(1D0) C/* RETURN END function xmass(i,j) C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ real*8 xmass,v1(4) integer i,j,k xmass=1D6 if(idup(i).eq.-idup(j)) return if(max(istup(i),istup(j)).eq.7) return xmass=0D0 do k=1,4 v1(k)=pup(k,i)+pup(k,j) xmass=xmass+sign(v1(k)**2,2.d0*dble(k)-7.d0) enddo xmass=sqrt(abs(xmass)) return end C----------------------------------------------------------------------- FUNCTION KTMDPI(PHI) IMPLICIT NONE C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) DOUBLE PRECISION KTMDPI,PHI,PI,TWOPI,THRPI,EPS PARAMETER (PI=3.14159265358979324D0,TWOPI=6.28318530717958648D0, & THRPI=9.42477796076937972D0) PARAMETER (EPS=1D-15) KTMDPI=PHI IF (KTMDPI.LE.PI) THEN IF (KTMDPI.GT.-PI) THEN GOTO 100 ELSEIF (KTMDPI.GT.-THRPI) THEN KTMDPI=KTMDPI+TWOPI ELSE KTMDPI=-MOD(PI-KTMDPI,TWOPI)+PI ENDIF ELSEIF (KTMDPI.LE.THRPI) THEN KTMDPI=KTMDPI-TWOPI ELSE KTMDPI=MOD(PI+KTMDPI,TWOPI)-PI ENDIF 100 IF (ABS(KTMDPI).LT.EPS) KTMDPI=0 END C...PYPTFS C...Generates pT-ordered timelike final-state parton showers. C...Calling parameters: C...NPART: the number of partons in the system to be showered, C... must be at least 2. Is updated at return to be the new C... number of partons. C...IPART: array, of dimension 500 (cf. Les Houches accord), C... in which the positions of the relevant partons are stored. C... To allow the identification of matrix elements, (the showered C... copies of) the original resonance decay products should be C... stored in the first two slots, IPART(1) and IPART(2). C... Is updated at return to be the new list of partons, by C... adding new particles at the end. Thus, for q -> q g the C... quark position is updated and the gluon added, while for C... g -> g g and g -> q qbar the choice of original and new is C... arbitrary. C...PTMAX : upper scale of shower evolution. An absolute limit is C... set by kinematical constraints inside each "dipole". C...PTMIN : lower scale of shower evolution. Some absolute lower C... limit is set by PARP(82)/2 and 1.1*Lambda(3flavours). C...PTGEN : returns the hardest pT generated; if none then PTGEN=0. C...The evolution is factorized, so that a set of successive calls, C...where the PTMIN of one call becomes the PTMAX of the next, gives C...the same result (on the average) as one single call for the full C...pT range. Note that the IPART(1) and IPART(2) entries continue C...to point to (the showered copies of) the original decay products C...of a resonance if they did so to begin with. C...Essentially all the relevant functionality of PYSHOW remains, C...but many non-standard options have not (yet) been implemented. C...Therefore the available options and parameters are: C...MSTJ(41) : (D=2) type of branching allowed in shower. C... = 0 : no showering at all. C... = 1, 11 : QCD branchings only. C... = 2, 12 : QCD branchings and photon emission. C...MSTJ(45) : (D=5) maximum flavour allowed in g -> q qbar branchings. C...MSTJ(46) : (D=3) nonhomogeneous azimuthal distribution in gluon C... branching induced by gluon polarization. C... = 0, 2 : no. C... = 1, 3 : yes C...MSTJ(47) (and MSTJ(38), PARJ(80)) : possibility to pick different C... matrix element corrections; default tries to make best match C... to given list of partons. C...PARJ(81) : (D=0.29 GeV) Lambda_QCD; should be reduced to around C... 0.14 in new shower. C...PARJ(82) : (D=1 GeV) twice pTmin cutoff for QCD shower. C...PARJ(83) : (D=1 GeV) twice pTmin cutoff for photon emission C... off a quark. C...PARJ(90) : (D=0.0001 GeV) twice pTmin cutoff for photon emission C... off a lepton. SUBROUTINE PYPTFS(NPART,IPART,PTMAX,PTMIN,PTGEN) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) c INTEGER PYK,PYCHGE,PYCOMP INTEGER PYCOMP C...Parameter statement to help give large particle numbers. PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000, &KEXCIT=4000000,KDIMEN=5000000) C...Parameter statement for maximum size of showers. PARAMETER (MAXNUP=500) C...Commonblocks. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) SAVE /PYJETS/,/PYDAT1/,/PYDAT2/ C...Input array. DIMENSION IPART(MAXNUP) INTEGER NMATCH REAL*8 Q0RES,QXRES,QHARD,PASSKT(20),SCALEQ(20),THORD(20) REAL*8 KTLIMIT COMMON/IPASS/Q0RES,QXRES,QHARD,PASSKT,SCALEQ,THORD,KTLIMIT,NMATCH SAVE /IPASS/ integer itime data itime/0/ save itime C...Local arrays. DIMENSION IPOS(2*MAXNUP),IREC(2*MAXNUP),IFLG(2*MAXNUP), &ISCOL(2*MAXNUP),ISCHG(2*MAXNUP),IMESAV(2*MAXNUP), &PT2SAV(2*MAXNUP),ZSAV(2*MAXNUP),SHTSAV(2*MAXNUP), &MESYS(MAXNUP,0:2),PSUM(5),DPT(5,4) C...Statement functions. SHAT(I,J)=(P(I,4)+P(J,4))**2-(P(I,1)+P(J,1))**2- &(P(I,2)+P(J,2))**2-(P(I,3)+P(J,3))**2 C...Initial values. Check that valid system. PTGEN=0D0 IF(MSTJ(41).NE.1.AND.MSTJ(41).NE.2.AND.MSTJ(41).NE.11.AND. &MSTJ(41).NE.12) RETURN IF(NPART.LT.2) THEN CALL PYERRM(12,'(PYPTFS:) showering system too small') RETURN ENDIF PT2CMX=PTMAX**2 NPART0=NPART C...Mass thresholds and Lambda for QCD evolution. PMB=PMAS(5,1) PMC=PMAS(4,1) ALAM5=PARJ(81) ALAM4=ALAM5*(PMB/ALAM5)**(2D0/25D0) ALAM3=ALAM4*(PMC/ALAM4)**(2D0/27D0) PMBS=PMB**2 PMCS=PMC**2 ALAM5S=ALAM5**2 ALAM4S=ALAM4**2 ALAM3S=ALAM3**2 C...Cutoff scale for QCD evolution. Starting pT2. NFLAV=MAX(0,MIN(5,MSTJ(45))) PT0C=0.5D0*PARJ(82) PT2CMN=MAX(PTMIN,PT0C,1.1D0*ALAM3)**2 C...Parameters for QED evolution. AEM2PI=PARU(101)/PARU(2) PT0EQ=0.5D0*PARJ(83) PT0EL=0.5D0*PARJ(90) C...Begin loop to set up showering partons. Sum four-momenta. NEVOL=0 DO 100 J=1,4 PSUM(J)=0D0 100 CONTINUE DO 180 IP=1,NPART I=IPART(IP) IF(K(I,1).GT.10) GOTO 180 DO 110 J=1,4 PSUM(J)=PSUM(J)+P(I,J) 110 CONTINUE C...Find colour and charge, but skip diquarks. IF(IABS(K(I,2)).GT.1000.AND.IABS(K(I,2)).LT.10000) GOTO 180 KCOL=ISIGN(KCHG(PYCOMP(K(I,2)),2),K(I,2)) KCHA=ISIGN(KCHG(PYCOMP(K(I,2)),1),K(I,2)) C...Either colour or anticolour charge radiates; for gluon both. DO 140 JSGCOL=1,-1,-2 IF(KCOL.EQ.JSGCOL.OR.KCOL.EQ.2) THEN JCOL=4+(1-JSGCOL)/2 JCOLR=9-JCOL C...Basic info about radiating parton. NEVOL=NEVOL+1 IPOS(NEVOL)=I IFLG(NEVOL)=0 ISCOL(NEVOL)=JSGCOL ISCHG(NEVOL)=0 C...Find sister which matching anticolour to the radiating parton. IROLD=I IRNEW=K(IROLD,JCOL)/MSTU(5) 120 IF(K(IRNEW,JCOLR)/MSTU(5).NE.IROLD) THEN C...Step up to mother if radiating parton already branched. IF(K(IRNEW,2).EQ.K(IROLD,2)) THEN IROLD=IRNEW IRNEW=K(IROLD,JCOL)/MSTU(5) GOTO 120 C...Pick sister by history if no anticolour available. ELSE IF(IROLD.GT.1.AND.K(IROLD-1,3).EQ.K(IROLD,3)) THEN IRNEW=IROLD-1 ELSEIF(IROLD.LT.N.AND.K(IROLD+1,3).EQ.K(IROLD,3)) THEN IRNEW=IROLD+1 C...Last resort: pick at random among other primaries. ELSE ISTEP=MAX(1,MIN(NPART-1,INT(1D0+(NPART-1)*PYR(0)))) IRNEW=IPART(1+MOD(IP+ISTEP-1,NPART)) ENDIF ENDIF ENDIF C...Trace down if sister branched. 130 IF(K(IRNEW,1).GT.10) THEN IRNEW=MOD(K(IRNEW,JCOLR),MSTU(5)) GOTO 130 ENDIF IREC(NEVOL)=IRNEW ENDIF 140 CONTINUE C...Also electrical charge may radiate; so far only quarks and leptons. IF((MSTJ(41).EQ.2.OR.MSTJ(41).EQ.12).AND.KCHA.NE.0.AND. & IABS(K(I,2)).LE.18) THEN C...Basic info about radiating parton. NEVOL=NEVOL+1 IPOS(NEVOL)=I IFLG(NEVOL)=0 ISCOL(NEVOL)=0 ISCHG(NEVOL)=KCHA C...Pick sister by history; step up if parton already branched. IROLD=I 150 IF(K(IROLD,3).GT.0.AND.K(K(IROLD,3),2).EQ.K(IROLD,2)) THEN IROLD=K(IROLD,3) GOTO 150 ENDIF IF(IROLD.GT.1.AND.K(IROLD-1,3).EQ.K(IROLD,3)) THEN IRNEW=IROLD-1 ELSEIF(IROLD.LT.N.AND.K(IROLD+1,3).EQ.K(IROLD,3)) THEN IRNEW=IROLD+1 C...Last resort: pick at random among other primaries. ELSE ISTEP=MAX(1,MIN(NPART-1,INT(1D0+(NPART-1)*PYR(0)))) IRNEW=IPART(1+MOD(IP+ISTEP-1,NPART)) ENDIF C...Trace down if sister branched. 160 IF(K(IRNEW,1).GT.10) THEN DO 170 IR=IRNEW+1,N IF(K(IR,3).EQ.IRNEW.AND.K(IR,2).EQ.K(IRNEW,2)) THEN IRNEW=IR GOTO 160 ENDIF 170 CONTINUE ENDIF IREC(NEVOL)=IRNEW ENDIF C...End loop to set up showering partons. System invariant mass. 180 CONTINUE PSUM(5)=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)) C...Check if 3-jet matrix elements to be used. M3JC=0 ALPHA=0.5D0 NMESYS=0 IF(MSTJ(47).GE.1) THEN C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0). KFSRCE=0 IPART1=K(IPART(1),3) IPART2=K(IPART(2),3) 190 IF(IPART1.EQ.IPART2.AND.IPART1.GT.0) THEN KFSRCE=IABS(K(IPART1,2)) ELSEIF(IPART1.GT.IPART2.AND.IPART2.GT.0) THEN IPART1=K(IPART1,3) GOTO 190 ELSEIF(IPART2.GT.IPART1.AND.IPART1.GT.0) THEN IPART2=K(IPART2,3) GOTO 190 ENDIF ITYPES=0 IF(KFSRCE.GE.1.AND.KFSRCE.LE.8) ITYPES=1 IF(KFSRCE.GE.KSUSY1+1.AND.KFSRCE.LE.KSUSY1+8) ITYPES=2 IF(KFSRCE.GE.KSUSY2+1.AND.KFSRCE.LE.KSUSY2+8) ITYPES=2 IF(KFSRCE.GE.21.AND.KFSRCE.LE.24) ITYPES=3 IF(KFSRCE.GE.32.AND.KFSRCE.LE.34) ITYPES=3 IF(KFSRCE.EQ.25.OR.(KFSRCE.GE.35.AND.KFSRCE.LE.37)) ITYPES=4 IF(KFSRCE.GE.KSUSY1+22.AND.KFSRCE.LE.KSUSY1+37) ITYPES=5 IF(KFSRCE.EQ.KSUSY1+21) ITYPES=6 C...Identify two primary showerers. KFLA1=IABS(K(IPART(1),2)) ITYPE1=0 IF(KFLA1.GE.1.AND.KFLA1.LE.8) ITYPE1=1 IF(KFLA1.GE.KSUSY1+1.AND.KFLA1.LE.KSUSY1+8) ITYPE1=2 IF(KFLA1.GE.KSUSY2+1.AND.KFLA1.LE.KSUSY2+8) ITYPE1=2 IF(KFLA1.GE.21.AND.KFLA1.LE.24) ITYPE1=3 IF(KFLA1.GE.32.AND.KFLA1.LE.34) ITYPE1=3 IF(KFLA1.EQ.25.OR.(KFLA1.GE.35.AND.KFLA1.LE.37)) ITYPE1=4 IF(KFLA1.GE.KSUSY1+22.AND.KFLA1.LE.KSUSY1+37) ITYPE1=5 IF(KFLA1.EQ.KSUSY1+21) ITYPE1=6 KFLA2=IABS(K(IPART(2),2)) ITYPE2=0 IF(KFLA2.GE.1.AND.KFLA2.LE.8) ITYPE2=1 IF(KFLA2.GE.KSUSY1+1.AND.KFLA2.LE.KSUSY1+8) ITYPE2=2 IF(KFLA2.GE.KSUSY2+1.AND.KFLA2.LE.KSUSY2+8) ITYPE2=2 IF(KFLA2.GE.21.AND.KFLA2.LE.24) ITYPE2=3 IF(KFLA2.GE.32.AND.KFLA2.LE.34) ITYPE2=3 IF(KFLA2.EQ.25.OR.(KFLA2.GE.35.AND.KFLA2.LE.37)) ITYPE2=4 IF(KFLA2.GE.KSUSY1+22.AND.KFLA2.LE.KSUSY1+37) ITYPE2=5 IF(KFLA2.EQ.KSUSY1+21) ITYPE2=6 C...Order of showerers. Presence of gluino. ITYPMN=MIN(ITYPE1,ITYPE2) ITYPMX=MAX(ITYPE1,ITYPE2) IORD=1 IF(ITYPE1.GT.ITYPE2) IORD=2 IGLUI=0 IF(ITYPE1.EQ.6.OR.ITYPE2.EQ.6) IGLUI=1 C...Require exactly two primary showerers for ME corrections. NPRIM=0 DO 200 I=1,N IF(K(I,3).EQ.IPART1.AND.K(I,2).NE.K(IPART1,2)) NPRIM=NPRIM+1 200 CONTINUE IF(NPRIM.NE.2) THEN C...Predetermined and default matrix element kinds. ELSEIF(MSTJ(38).NE.0) THEN M3JC=MSTJ(38) ALPHA=PARJ(80) MSTJ(38)=0 ELSEIF(MSTJ(47).GE.6) THEN M3JC=MSTJ(47) ELSE ICLASS=1 ICOMBI=4 C...Vector/axial vector -> q + qbar; q -> q + V. IF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.3)) THEN ICLASS=2 IF(KFSRCE.EQ.21.OR.KFSRCE.EQ.22) THEN ICOMBI=1 ELSEIF(KFSRCE.EQ.23.OR.(KFSRCE.EQ.0.AND. & K(IPART(1),2)+K(IPART(2),2).EQ.0)) THEN C...gamma*/Z0: assume e+e- initial state if unknown. EI=-1D0 IF(KFSRCE.EQ.23) THEN IANNFL=K(IPART1,3) IF(K(IANNFL,2).EQ.23) IANNFL=K(IANNFL,3) IF(IANNFL.NE.0) THEN KANNFL=IABS(K(IANNFL,2)) IF(KANNFL.GE.1.AND.KANNFL.LE.18) EI=KCHG(KANNFL,1)/3D0 ENDIF ENDIF AI=SIGN(1D0,EI+0.1D0) VI=AI-4D0*EI*PARU(102) EF=KCHG(KFLA1,1)/3D0 AF=SIGN(1D0,EF+0.1D0) VF=AF-4D0*EF*PARU(102) XWC=1D0/(16D0*PARU(102)*(1D0-PARU(102))) SH=PSUM(5)**2 SQMZ=PMAS(23,1)**2 SQWZ=PSUM(5)*PMAS(23,2) SBWZ=1D0/((SH-SQMZ)**2+SQWZ**2) VECT=EI**2*EF**2+2D0*EI*VI*EF*VF*XWC*SH*(SH-SQMZ)*SBWZ+ & (VI**2+AI**2)*VF**2*XWC**2*SH**2*SBWZ AXIV=(VI**2+AI**2)*AF**2*XWC**2*SH**2*SBWZ ICOMBI=3 ALPHA=VECT/(VECT+AXIV) ELSEIF(KFSRCE.EQ.24.OR.KFSRCE.EQ.0) THEN ICOMBI=4 ENDIF C...For chi -> chi q qbar, use V/A -> q qbar as first approximation. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.5) THEN ICLASS=2 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=3 C...Scalar/pseudoscalar -> q + qbar; q -> q + S. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.4) THEN ICLASS=4 IF(KFSRCE.EQ.25.OR.KFSRCE.EQ.35.OR.KFSRCE.EQ.37) THEN ICOMBI=1 ELSEIF(KFSRCE.EQ.36) THEN ICOMBI=2 ENDIF ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=5 C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S. ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.3)) THEN ICLASS=6 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=7 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.4) THEN ICLASS=8 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=9 C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.5)) THEN ICLASS=10 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=11 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=12 C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g. ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.6) THEN ICLASS=13 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.2)) THEN ICLASS=14 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR. & ITYPES.EQ.1)) THEN ICLASS=15 C...g -> ~g + ~g (eikonal approximation). ELSEIF(ITYPMN.EQ.6.AND.ITYPMX.EQ.6.AND.ITYPES.EQ.0) THEN ICLASS=16 ENDIF M3JC=5*ICLASS+ICOMBI ENDIF C...Store pair that together define matrix element treatment. IF(M3JC.NE.0) THEN NMESYS=1 MESYS(NMESYS,0)=M3JC MESYS(NMESYS,1)=IPART(1) MESYS(NMESYS,2)=IPART(2) ENDIF C...Store qqbar or l+l- pairs for QED radiation. IF(KFLA1.LE.18.AND.KFLA2.LE.18) THEN NMESYS=NMESYS+1 MESYS(NMESYS,0)=101 IF(K(IPART(1),2)+K(IPART(2),2).EQ.0) MESYS(NMESYS,0)=102 MESYS(NMESYS,1)=IPART(1) MESYS(NMESYS,2)=IPART(2) ENDIF C...Store other qqbar/l+l- pairs from g/gamma branchings. DO 240 I1=1,N IF(K(I1,1).GT.10.OR.IABS(K(I1,2)).GT.18) GOTO 240 I1M=K(I1,3) 210 IF(I1M.GT.0.AND.K(I1M,2).EQ.K(I1,2)) THEN I1M=K(I1M,3) GOTO 210 ENDIF IF(K(I1M,2).NE.21.AND.K(I1M,2).NE.22) GOTO 240 DO 230 I2=I1+1,N IF(K(I2,1).GT.10.OR.K(I2,2)+K(I1,2).NE.0) GOTO 230 I2M=K(I2,3) 220 IF(I2M.GT.0.AND.K(I2M,2).EQ.K(I2,2)) THEN I2M=K(I2M,3) GOTO 220 ENDIF IF(I1M.EQ.I2M.AND.I1M.GT.0) THEN NMESYS=NMESYS+1 MESYS(NMESYS,0)=66 MESYS(NMESYS,1)=I1 MESYS(NMESYS,2)=I2 NMESYS=NMESYS+1 MESYS(NMESYS,0)=102 MESYS(NMESYS,1)=I1 MESYS(NMESYS,2)=I2 ENDIF 230 CONTINUE 240 CONTINUE ENDIF C..Loopback point for counting number of emissions. NGEN=0 250 NGEN=NGEN+1 C...Begin loop to evolve all existing partons, if required. 260 IMX=0 PT2MX=0D0 DO 320 IEVOL=1,NEVOL IF(IFLG(IEVOL).EQ.0) THEN C...Basic info on radiator and recoil. I=IPOS(IEVOL) IR=IREC(IEVOL) SHT=SHAT(I,IR) PM2I=P(I,5)**2 PM2R=P(IR,5)**2 C...Invariant mass of "dipole".Starting value for pT evolution. SHTCOR=(SQRT(SHT)-P(IR,5))**2-PM2I C/* PT2CMX depends upon the initiator IF(I.LE.102+NPART0) THEN PT2CMX=SCALEQ(I-100)**2 VETOKT=PASSKT(I-100) ENDIF PT2=MIN(PT2CMX,0.25D0*SHTCOR) c DES if(itime.lt.100) then c DES print*,' starting shower of ',ievol,i,sqrt(pt2cmx), c DES $ sqrt(.25D0*shtcor),vetokt,k(i,2),npart0 c DES endif C...Case of evolution by QCD branching. IF(ISCOL(IEVOL).NE.0) THEN C...If kinematically impossible then do not evolve. IF(PT2.LT.PT2CMN) THEN IFLG(IEVOL)=-1 GOTO 320 ENDIF C...Check if part of system for which ME corrections should be applied. IMESYS=0 DO 270 IME=1,NMESYS IF((I.EQ.MESYS(IME,1).OR.I.EQ.MESYS(IME,2)).AND. & MESYS(IME,0).LT.100) IMESYS=IME 270 CONTINUE C...Special flag for colour octet states. MOCT=0 IF(K(I,2).EQ.21) MOCT=1 IF(K(I,2).EQ.KSUSY1+21) MOCT=2 C...Upper estimate for matrix element weighting and colour factor. C...Note that g->gg and g->qqbar is split on two sides = "dipoles". WTPSGL=2D0 COLFAC=4D0/3D0 IF(MOCT.GE.1) COLFAC=3D0/2D0 IF(IGLUI.EQ.1.AND.IMESYS.EQ.1.AND.MOCT.EQ.0) COLFAC=3D0 WTPSQQ=0.5D0*0.5D0*NFLAV C...Determine overestimated z range: switch at c and b masses. 280 IZRG=1 PT2MNE=PT2CMN B0=27D0/6D0 ALAMS=ALAM3S IF(PT2.GT.1.01D0*PMCS) THEN IZRG=2 PT2MNE=PMCS B0=25D0/6D0 ALAMS=ALAM4S ENDIF IF(PT2.GT.1.01D0*PMBS) THEN IZRG=3 PT2MNE=PMBS B0=23D0/6D0 ALAMS=ALAM5S ENDIF ZMNCUT=0.5D0-SQRT(MAX(0D0,0.25D0-PT2MNE/SHTCOR)) IF(ZMNCUT.LT.1D-8) ZMNCUT=PT2MNE/SHTCOR C...Find evolution coefficients for q->qg/g->gg and g->qqbar. EVEMGL=WTPSGL*COLFAC*LOG(1D0/ZMNCUT-1D0)/B0 EVCOEF=EVEMGL IF(MOCT.EQ.1) THEN EVEMQQ=WTPSQQ*(1D0-2D0*ZMNCUT)/B0 EVCOEF=EVCOEF+EVEMQQ ENDIF C...Pick pT2 (in overestimated z range). 290 PT2=ALAMS*(PT2/ALAMS)**(PYR(0)**(1D0/EVCOEF)) C...Loopback if crossed c/b mass thresholds. IF(IZRG.EQ.3.AND.PT2.LT.PMBS) THEN PT2=PMBS GOTO 280 ENDIF IF(IZRG.EQ.2.AND.PT2.LT.PMCS) THEN PT2=PMCS GOTO 280 ENDIF C...Finish if below lower cutoff. IF(PT2.LT.PT2CMN) THEN IFLG(IEVOL)=-1 GOTO 320 ENDIF C...Pick kind of branching: q->qg/g->gg/X->Xg or g->qqbar. IFLAG=1 IF(MOCT.EQ.1.AND.EVEMGL.LT.PYR(0)*EVCOEF) IFLAG=2 C...Pick z: dz/(1-z) or dz. IF(IFLAG.EQ.1) THEN Z=1D0-ZMNCUT*(1D0/ZMNCUT-1D0)**PYR(0) ELSE Z=ZMNCUT+PYR(0)*(1D0-2D0*ZMNCUT) ENDIF C...Loopback if outside allowed range for given pT2. ZMNNOW=0.5D0-SQRT(MAX(0D0,0.25D0-PT2/SHTCOR)) IF(ZMNNOW.LT.1D-8) ZMNNOW=PT2/SHTCOR IF(Z.LE.ZMNNOW.OR.Z.GE.1D0-ZMNNOW) GOTO 290 PM2=PM2I+PT2/(Z*(1D0-Z)) IF(Z*(1D0-Z).LE.PM2*SHT/(SHT+PM2-PM2R)**2) GOTO 290 C...No weighting for primary partons; to be done later on. IF(IMESYS.GT.0) THEN C...Weighting of q->qg/X->Xg branching. ELSEIF(IFLAG.EQ.1.AND.MOCT.NE.1) THEN IF(1D0+Z**2.LT.WTPSGL*PYR(0)) GOTO 290 C...Weighting of g->gg branching. ELSEIF(IFLAG.EQ.1) THEN IF(1D0+Z**3.LT.WTPSGL*PYR(0)) GOTO 290 C...Flavour choice and weighting of g->qqbar branching. ELSE KFQ=MIN(5,1+INT(NFLAV*PYR(0))) PMQ=PMAS(KFQ,1) ROOTQQ=SQRT(MAX(0D0,1D0-4D0*PMQ**2/PM2)) WTME=ROOTQQ*(Z**2+(1D0-Z)**2) IF(WTME.LT.PYR(0)) GOTO 290 IFLAG=10+KFQ ENDIF C/* Veto itime=itime+1 IF(PT2.GT.VETOKT**2) THEN c DES if(itime.le.100) then c DES print*,' veto ',sqrt(pt2),vetokt c DES endif GOTO 290 ENDIF C...Case of evolution by QED branching. ELSEIF(ISCHG(IEVOL).NE.0) THEN C...If kinematically impossible then do not evolve. PT2EMN=PT0EQ**2 IF(IABS(K(I,2)).GT.10) PT2EMN=PT0EL**2 IF(PT2.LT.PT2EMN) THEN IFLG(IEVOL)=-1 GOTO 320 ENDIF C...Check if part of system for which ME corrections should be applied. IMESYS=0 DO 300 IME=1,NMESYS IF((I.EQ.MESYS(IME,1).OR.I.EQ.MESYS(IME,2)).AND. & MESYS(IME,0).GT.100) IMESYS=IME 300 CONTINUE C...Charge. Matrix element weighting factor. CHG=ISCHG(IEVOL)/3D0 WTPSGA=2D0 C...Determine overestimated z range. Find evolution coefficient. ZMNCUT=0.5D0-SQRT(MAX(0D0,0.25D0-PT2EMN/SHTCOR)) IF(ZMNCUT.LT.1D-8) ZMNCUT=PT2EMN/SHTCOR EVCOEF=AEM2PI*CHG**2*WTPSGA*LOG(1D0/ZMNCUT-1D0) C...Pick pT2 (in overestimated z range). 310 PT2=PT2*PYR(0)**(1D0/EVCOEF) C...Finish if below lower cutoff. IF(PT2.LT.PT2EMN) THEN IFLG(IEVOL)=-1 GOTO 320 ENDIF C...Pick z: dz/(1-z). Z=1D0-ZMNCUT*(1D0/ZMNCUT-1D0)**PYR(0) C...Loopback if outside allowed range for given pT2. ZMNNOW=0.5D0-SQRT(MAX(0D0,0.25D0-PT2/SHTCOR)) IF(ZMNNOW.LT.1D-8) ZMNNOW=PT2/SHTCOR IF(Z.LE.ZMNNOW.OR.Z.GE.1D0-ZMNNOW) GOTO 310 PM2=PM2I+PT2/(Z*(1D0-Z)) IF(Z*(1D0-Z).LE.PM2*SHT/(SHT+PM2-PM2R)**2) GOTO 310 C...Weighting by branching kernel, except if ME weighting later. IF(IMESYS.EQ.0) THEN IF(1D0+Z**2.LT.WTPSGA*PYR(0)) GOTO 310 ENDIF IFLAG=3 ENDIF C...Save acceptable branching. IFLG(IEVOL)=IFLAG IMESAV(IEVOL)=IMESYS PT2SAV(IEVOL)=PT2 ZSAV(IEVOL)=Z SHTSAV(IEVOL)=SHT ENDIF C...Check if branching has highest pT. IF(IFLG(IEVOL).GE.1.AND.PT2SAV(IEVOL).GT.PT2MX) THEN IMX=IEVOL PT2MX=PT2SAV(IEVOL) ENDIF 320 CONTINUE C...Finished if no more branchings to be done. IF(IMX.EQ.0) GOTO 440 C...Restore info on hardest branching to be processed. I=IPOS(IMX) IR=IREC(IMX) KCOL=ISCOL(IMX) KCHA=ISCHG(IMX) IMESYS=IMESAV(IMX) PT2=PT2SAV(IMX) Z=ZSAV(IMX) SHT=SHTSAV(IMX) PM2I=P(I,5)**2 PM2R=P(IR,5)**2 PM2=PM2I+PT2/(Z*(1D0-Z)) C...Special flag for colour octet states. MOCT=0 IF(K(I,2).EQ.21) MOCT=1 IF(K(I,2).EQ.KSUSY1+21) MOCT=2 C...Restore further info for g->qqbar branching. KFQ=0 IF(IFLG(IMX).GT.10) THEN KFQ=IFLG(IMX)-10 PMQ=PMAS(KFQ,1) ROOTQQ=SQRT(MAX(0D0,1D0-4D0*PMQ**2/PM2)) ENDIF C...For branching g include azimuthal asymmetries from polarization. ASYPOL=0D0 IF(MOCT.EQ.1.AND.MOD(MSTJ(46),2).EQ.1) THEN C...Trace grandmother via intermediate recoil copies. KFGM=0 IM=I 330 IF(K(IM,3).NE.K(IM-1,3).AND.K(IM,3).NE.K(IM+1,3).AND. & K(IM,3).GT.0) THEN IM=K(IM,3) GOTO 330 ENDIF IGM=K(IM,3) IF(IGM.GT.0) KFGM=IABS(K(IGM,2)) C...Define approximate energy sharing by identifying aunt. IAU=IM+1 IF(IAU.GT.N-3.OR.K(IAU,3).NE.IGM) IAU=IM-1 IF(KFGM.NE.0.AND.(KFGM.LE.6.OR.KFGM.EQ.21)) THEN ZOLD=P(IM,4)/(P(IM,4)+P(IAU,4)) C...Coefficient from gluon production. IF(KFGM.LE.6) THEN ASYPOL=2D0*(1D0-ZOLD)/(1D0+(1D0-ZOLD)**2) ELSE ASYPOL=((1D0-ZOLD)/(1D0-ZOLD*(1D0-ZOLD)))**2 ENDIF C...Coefficient from gluon decay. IF(KFQ.EQ.0) THEN ASYPOL=ASYPOL*(Z*(1D0-Z)/(1D0-Z*(1D0-Z)))**2 ELSE ASYPOL=-ASYPOL*2D0*Z*(1D0-Z)/(1D0-2D0*Z*(1D0-Z)) ENDIF ENDIF ENDIF C...Create new slots for branching products and recoil. INEW=N+1 IGNEW=N+2 IRNEW=N+3 N=N+3 C...Set status, flavour and mother of new ones. K(INEW,1)=K(I,1) K(IGNEW,1)=3 IF(KCHA.NE.0) K(IGNEW,1)=1 K(IRNEW,1)=K(IR,1) IF(KFQ.EQ.0) THEN K(INEW,2)=K(I,2) K(IGNEW,2)=21 IF(KCHA.NE.0) K(IGNEW,2)=22 ELSE K(INEW,2)=-ISIGN(KFQ,KCOL) K(IGNEW,2)=-K(INEW,2) ENDIF K(IRNEW,2)=K(IR,2) K(INEW,3)=I K(IGNEW,3)=I K(IRNEW,3)=IR C...Find rest frame and angles of branching+recoil. DO 340 J=1,5 P(INEW,J)=P(I,J) P(IGNEW,J)=0D0 P(IRNEW,J)=P(IR,J) 340 CONTINUE BETAX=(P(INEW,1)+P(IRNEW,1))/(P(INEW,4)+P(IRNEW,4)) BETAY=(P(INEW,2)+P(IRNEW,2))/(P(INEW,4)+P(IRNEW,4)) BETAZ=(P(INEW,3)+P(IRNEW,3))/(P(INEW,4)+P(IRNEW,4)) CALL PYROBO(INEW,IRNEW,0D0,0D0,-BETAX,-BETAY,-BETAZ) PHI=PYANGL(P(INEW,1),P(INEW,2)) THETA=PYANGL(P(INEW,3),SQRT(P(INEW,1)**2+P(INEW,2)**2)) C...Derive kinematics of branching: generics (like g->gg). DO 350 J=1,4 P(INEW,J)=0D0 P(IRNEW,J)=0D0 350 CONTINUE PEM=0.5D0*(SHT+PM2-PM2R)/SQRT(SHT) PZM=0.5D0*SQRT(MAX(0D0,(SHT-PM2-PM2R)**2-4D0*PM2*PM2R)/SHT) PT2COR=PM2*(PEM**2*Z*(1D0-Z)-0.25D0*PM2)/PZM**2 PTCOR=SQRT(MAX(0D0,PT2COR)) PZN=(PEM**2*Z-0.5D0*PM2)/PZM PZG=(PEM**2*(1D0-Z)-0.5D0*PM2)/PZM C...Specific kinematics reduction for q->qg with m_q > 0. IF(MOCT.NE.1) THEN PTCOR=(1D0-PM2I/PM2)*PTCOR PZN=PZN+PM2I*PZG/PM2 PZG=(1D0-PM2I/PM2)*PZG C...Specific kinematics reduction for g->qqbar with m_q > 0. ELSEIF(KFQ.NE.0) THEN P(INEW,5)=PMQ P(IGNEW,5)=PMQ PTCOR=ROOTQQ*PTCOR PZN=0.5D0*((1D0+ROOTQQ)*PZN+(1D0-ROOTQQ)*PZG) PZG=PZM-PZN ENDIF C...Pick phi and construct kinematics of branching. 360 PHIROT=PARU(2)*PYR(0) P(INEW,1)=PTCOR*COS(PHIROT) P(INEW,2)=PTCOR*SIN(PHIROT) P(INEW,3)=PZN P(INEW,4)=SQRT(PTCOR**2+P(INEW,3)**2+P(INEW,5)**2) P(IGNEW,1)=-P(INEW,1) P(IGNEW,2)=-P(INEW,2) P(IGNEW,3)=PZG P(IGNEW,4)=SQRT(PTCOR**2+P(IGNEW,3)**2+P(IGNEW,5)**2) P(IRNEW,1)=0D0 P(IRNEW,2)=0D0 P(IRNEW,3)=-PZM P(IRNEW,4)=0.5D0*(SHT+PM2R-PM2)/SQRT(SHT) C...Boost branching system to lab frame. CALL PYROBO(INEW,IRNEW,THETA,PHI,BETAX,BETAY,BETAZ) C...Renew choice of phi angle according to polarization asymmetry. IF(ABS(ASYPOL).GT.1D-3) THEN DO 370 J=1,3 DPT(1,J)=P(I,J) DPT(2,J)=P(IAU,J) DPT(3,J)=P(INEW,J) 370 CONTINUE DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3) DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3) DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2 DO 380 J=1,3 DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/MAX(1D-10,DPMM) DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/MAX(1D-10,DPMM) 380 CONTINUE DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2) DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2) IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1D0*PARJ(82)) THEN CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+ & DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4)) IF(1D0+ASYPOL*(2D0*CAD**2-1D0).LT.PYR(0)*(1D0+ABS(ASYPOL))) & GOTO 360 ENDIF ENDIF C...Matrix element corrections for primary partons when requested. IF(IMESYS.GT.0) THEN M3JC=MESYS(IMESYS,0) C...Identify recoiling partner and set up three-body kinematics. IRP=MESYS(IMESYS,1) IF(IRP.EQ.I) IRP=MESYS(IMESYS,2) IF(IRP.EQ.IR) IRP=IRNEW DO 390 J=1,4 PSUM(J)=P(INEW,J)+P(IRP,J)+P(IGNEW,J) 390 CONTINUE PSUM(5)=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2- & PSUM(3)**2)) X1=2D0*(PSUM(4)*P(INEW,4)-PSUM(1)*P(INEW,1)-PSUM(2)*P(INEW,2)- & PSUM(3)*P(INEW,3))/PSUM(5)**2 X2=2D0*(PSUM(4)*P(IRP,4)-PSUM(1)*P(IRP,1)-PSUM(2)*P(IRP,2)- & PSUM(3)*P(IRP,3))/PSUM(5)**2 R1ME=P(INEW,5)/PSUM(5) R2ME=P(IRP,5)/PSUM(5) C...Matrix elements for gluon emission. IF(M3JC.LT.100) THEN C...Call ME, with right order important for two inequivalent showerers. IF(MESYS(IMESYS,IORD).EQ.I) THEN WME=PYMAEL(M3JC,X1,X2,R1ME,R2ME,ALPHA) ELSE WME=PYMAEL(M3JC,X2,X1,R2ME,R1ME,ALPHA) ENDIF C...Split up total ME when two radiating partons. ISPRAD=1 IF((M3JC.GE.16.AND.M3JC.LE.19).OR.(M3JC.GE.26.AND.M3JC.LE.29) & .OR.(M3JC.GE.36.AND.M3JC.LE.39).OR.(M3JC.GE.46.AND.M3JC.LE.49) & .OR.(M3JC.GE.56.AND.M3JC.LE.64)) ISPRAD=0 IF(ISPRAD.EQ.1) WME=WME*MAX(1D-10,1D0+R1ME**2-R2ME**2-X1)/ & MAX(1D-10,2D0-X1-X2) C...Evaluate shower rate. WPS=2D0/(MAX(1D-10,2D0-X1-X2)* & MAX(1D-10,1D0+R2ME**2-R1ME**2-X2)) IF(IGLUI.EQ.1) WPS=(9D0/4D0)*WPS C...Matrix elements for photon emission: still rather primitive. ELSE C...For generic charge combination currently only massless expression. IF(M3JC.EQ.101) THEN CHG1=KCHG(PYCOMP(K(I,2)),1)*ISIGN(1,K(I,2))/3D0 CHG2=KCHG(PYCOMP(K(IRP,2)),1)*ISIGN(1,K(IRP,2))/3D0 X3=2D0-X1-X2 WME=(CHG1*(1D0-X1)/X3-CHG2*(1D0-X2)/X3)**2*(X1**2+X2**2) WPS=2D0*(CHG1**2*(1D0-X1)/X3+CHG2**2*(1D0-X2)/X3) C...For flavour neutral system assume vector source and include masses. ELSE WME=PYMAEL(11,X1,X2,R1ME,R2ME,0D0)*MAX(1D-10, & 1D0+R1ME**2-R2ME**2-X1)/MAX(1D-10,2D0-X1-X2) WPS=2D0/(MAX(1D-10,2D0-X1-X2)* & MAX(1D-10,1D0+R2ME**2-R1ME**2-X2)) ENDIF ENDIF C...Perform weighting with W_ME/W_PS. IF(WME.LT.PYR(0)*WPS) THEN N=N-3 IFLG(IMX)=0 GOTO 260 ENDIF ENDIF C...Now for sure accepted branching. Save highest pT. IF(NGEN.EQ.1) PTGEN=SQRT(PT2) C...Update status for obsolete ones. Bookkkep the moved original parton C...and new daughter (arbitrary choice for g->gg or g->qqbar). C...Do not bookkeep radiated photon, since it cannot radiate further. K(I,1)=K(I,1)+10 K(IR,1)=K(IR,1)+10 DO 400 IP=1,NPART IF(IPART(IP).EQ.I) IPART(IP)=INEW IF(IPART(IP).EQ.IR) IPART(IP)=IRNEW 400 CONTINUE IF(KCHA.EQ.0) THEN NPART=NPART+1 IPART(NPART)=IGNEW ENDIF C...Initialize colour flow of branching. K(INEW,4)=0 K(IGNEW,4)=0 K(INEW,5)=0 K(IGNEW,5)=0 JCOLP=4+(1-KCOL)/2 JCOLN=9-JCOLP C...Trivial colour flow for l->lgamma and q->qgamma. IF(IABS(KCHA).EQ.3) THEN K(I,4)=INEW K(I,5)=IGNEW ELSEIF(KCHA.NE.0) THEN IF(K(I,4).NE.0) THEN K(I,4)=K(I,4)+INEW K(INEW,4)=MSTU(5)*I ENDIF IF(K(I,5).NE.0) THEN K(I,5)=K(I,5)+INEW K(INEW,5)=MSTU(5)*I ENDIF C...Set colour flow for q->qg and g->gg. ELSEIF(KFQ.EQ.0) THEN K(I,JCOLP)=K(I,JCOLP)+IGNEW K(IGNEW,JCOLP)=MSTU(5)*I K(INEW,JCOLP)=MSTU(5)*IGNEW K(IGNEW,JCOLN)=MSTU(5)*INEW IF(MOCT.GE.1) THEN K(I,JCOLN)=K(I,JCOLN)+INEW K(INEW,JCOLN)=MSTU(5)*I ENDIF C...Set colour flow for g->qqbar. ELSE K(I,JCOLN)=K(I,JCOLN)+INEW K(INEW,JCOLN)=MSTU(5)*I K(I,JCOLP)=K(I,JCOLP)+IGNEW K(IGNEW,JCOLP)=MSTU(5)*I ENDIF C...Colour of recoiling parton sails through unchanged. IF(K(IR,4).NE.0) THEN K(IR,4)=K(IR,4)+IRNEW K(IRNEW,4)=MSTU(5)*IR ENDIF IF(K(IR,5).NE.0) THEN K(IR,5)=K(IR,5)+IRNEW K(IRNEW,5)=MSTU(5)*IR ENDIF C...Vertex information trivial. DO 410 J=1,5 V(INEW,J)=V(I,J) V(IGNEW,J)=V(I,J) V(IRNEW,J)=V(IR,J) 410 CONTINUE C...Update list of old radiators. DO 420 IEVOL=1,NEVOL IF(IPOS(IEVOL).EQ.I.AND.IREC(IEVOL).EQ.IR) THEN IPOS(IEVOL)=INEW IF(KCOL.NE.0.AND.ISCOL(IEVOL).NE.0) IPOS(IEVOL)=IGNEW IREC(IEVOL)=IRNEW IFLG(IEVOL)=0 ELSEIF(IPOS(IEVOL).EQ.I) THEN IPOS(IEVOL)=INEW IFLG(IEVOL)=0 ELSEIF(IPOS(IEVOL).EQ.IR.AND.IREC(IEVOL).EQ.I) THEN IPOS(IEVOL)=IRNEW IREC(IEVOL)=INEW IF(KCOL.NE.0.AND.ISCOL(IEVOL).NE.0) IREC(IEVOL)=IGNEW IFLG(IEVOL)=0 ELSEIF(IPOS(IEVOL).EQ.IR) THEN IPOS(IEVOL)=IRNEW IFLG(IEVOL)=0 ENDIF C...Update links of old connected partons. IF(IREC(IEVOL).EQ.I) THEN IREC(IEVOL)=INEW IFLG(IEVOL)=0 ELSEIF(IREC(IEVOL).EQ.IR) THEN IREC(IEVOL)=IRNEW IFLG(IEVOL)=0 ENDIF 420 CONTINUE C...q->qg or g->gg: create new gluon radiators. IF(KCOL.NE.0.AND.KFQ.EQ.0) THEN NEVOL=NEVOL+1 IPOS(NEVOL)=INEW IREC(NEVOL)=IGNEW IFLG(NEVOL)=0 ISCOL(NEVOL)=KCOL ISCHG(NEVOL)=0 NEVOL=NEVOL+1 IPOS(NEVOL)=IGNEW IREC(NEVOL)=INEW IFLG(NEVOL)=0 ISCOL(NEVOL)=-KCOL ISCHG(NEVOL)=0 ENDIF C...Update matrix elements parton list and add new for g/gamma->qqbar. DO 430 IME=1,NMESYS IF(MESYS(IME,1).EQ.I) MESYS(IME,1)=INEW IF(MESYS(IME,2).EQ.I) MESYS(IME,2)=INEW IF(MESYS(IME,1).EQ.IR) MESYS(IME,1)=IRNEW IF(MESYS(IME,2).EQ.IR) MESYS(IME,2)=IRNEW 430 CONTINUE IF(KFQ.NE.0) THEN NMESYS=NMESYS+1 MESYS(NMESYS,0)=66 MESYS(NMESYS,1)=INEW MESYS(NMESYS,2)=IGNEW NMESYS=NMESYS+1 MESYS(NMESYS,0)=102 MESYS(NMESYS,1)=INEW MESYS(NMESYS,2)=IGNEW ENDIF C...Loopback for more emissions if enough space. PT2CMX=PT2 IF(NPART.LT.MAXNUP-1.AND.NEVOL.LT.2*MAXNUP-2.AND. &NMESYS.LT.MAXNUP-2.AND.N.LT.MSTU(4)-MSTU(32)-5) THEN GOTO 250 ELSE CALL PYERRM(11,'(PYPTFS:) no more memory left for shower') ENDIF C...Done. 440 CONTINUE RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE KTCLUS(IMODE,PP,NN,ECUT,Y,*) IMPLICIT NONE C---DO CLUSTER ANALYSIS OF PARTICLES IN PP C C IMODE = INPUT : DESCRIBED ABOVE C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E C NN = INPUT : NUMBER OF PARTICLES C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED C Y(J) = OUTPUT : VALUE OF Y FOR WHICH EVENT CHANGES FROM BEING C J JET TO J-1 JET C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED (MOST LIKELY DUE TO TOO MANY PARTICLES) C C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER IMODE,NN DOUBLE PRECISION PP(4,*) DOUBLE PRECISION ECUT,Y(*),ONE ONE=1 CALL KTCLUR(IMODE,PP,NN,ONE,ECUT,Y,*999) RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTCLUR(IMODE,PP,NN,R,ECUT,Y,*) IMPLICIT NONE C---DO CLUSTER ANALYSIS OF PARTICLES IN PP C C IMODE = INPUT : DESCRIBED ABOVE C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E C NN = INPUT : NUMBER OF PARTICLES C R = INPUT : ELLIS AND SOPER'S R PARAMETER, SEE ABOVE. C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED C Y(J) = OUTPUT : VALUE OF Y FOR WHICH EVENT CHANGES FROM BEING C J JET TO J-1 JET C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED (MOST LIKELY DUE TO TOO MANY PARTICLES) C C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER NMAX,IM,IMODE,TYPE,ANGL,MONO,RECO,N,I,J,NN, & IMIN,JMIN,KMIN,NUM,HIST,INJET,IABBR,NABBR PARAMETER (NMAX=512,NABBR=7) DOUBLE PRECISION PP(4,*) DOUBLE PRECISION R,ECUT,Y(*),P,KT,ETOT,RSQ,KTP,KTS,KTPAIR,KTSING, & KTMIN,ETSQ,KTLAST,KTMAX,KTTMP LOGICAL FIRST CHARACTER TITLE(4,4)*10 C---KT RECORDS THE KT**2 OF EACH MERGING. C---KTLAST RECORDS FOR EACH MERGING, THE HIGHEST ECUT**2 FOR WHICH THE C RESULT IS NOT MERGED WITH THE BEAM (COULD BE LARGER THAN THE C KT**2 AT WHICH IT WAS MERGED IF THE KT VALUES ARE NOT MONOTONIC). C THIS MAY SOUND POINTLESS, BUT ITS USEFUL FOR DETERMINING WHETHER C SUB-JETS SURVIVED TO SCALE Y=YMAC OR NOT. C---HIST RECORDS MERGING HISTORY: C N=>DELETED TRACK N, M*NMAX+N=>MERGED TRACKS M AND N (M PX,PY,PZ,E C NN = INPUT : NUMBER OF PARTICLES C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED C YCUT = INPUT : Y VALUE AT WHICH TO RECONSTRUCT JET MOMENTA C YMAC = INPUT : Y VALUE USED TO DEFINE MACRO-JETS, TO DETERMINE C WHICH JETS ARE SUB-JETS C PJET(I,J)=OUTPUT : 4-MOMENTUM OF Jth JET AT SCALE YCUT C JET(J) =OUTPUT : THE MACRO-JET WHICH CONTAINS THE Jth JET, C SET TO ZERO IF JET IS NOT A SUB-JET C NJET =OUTPUT : THE NUMBER OF JETS C NSUB =OUTPUT : THE NUMBER OF SUB-JETS (EQUAL TO THE NUMBER OF C NON-ZERO ENTRIES IN JET()) C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER NMAX,RECO,NUM,N,NN,NJET,NSUB,JET(*),HIST,IMIN,JMIN,I,J PARAMETER (NMAX=512) DOUBLE PRECISION PP(4,*),PJET(4,*) DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,ETSQ,YCUT,YMAC,KTLAST, & ROUND PARAMETER (ROUND=0.99999D0) COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM C---CHECK INPUT IF (RECO.LT.1.OR.RECO.GT.3) CALL KTWARN('KTRECO',100,*999) C---COPY PP TO P N=NN IF (NUM.NE.NN) CALL KTWARN('KTRECO',101,*999) CALL KTCOPY(PP,N,P,(RECO.NE.1)) IF (ECUT.EQ.0) THEN ETSQ=1/ETOT**2 ELSE ETSQ=1/ECUT**2 ENDIF C---KEEP MERGING UNTIL YCUT 100 IF (ETSQ*KT(N).LT.ROUND*YCUT) THEN IF (HIST(N).LE.NMAX) THEN CALL KTMOVE(P,KTP,KTS,NMAX,N,HIST(N),0) ELSE IMIN=HIST(N)/NMAX JMIN=HIST(N)-IMIN*NMAX CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO) CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0) ENDIF N=N-1 IF (N.GT.0) GOTO 100 ENDIF C---IF YCUT IS TOO LARGE THERE ARE NO JETS NJET=N NSUB=N IF (N.EQ.0) RETURN C---SET UP OUTPUT MOMENTA DO 210 I=1,NJET IF (RECO.EQ.1) THEN DO 200 J=1,4 PJET(J,I)=P(J,I) 200 CONTINUE ELSE PJET(1,I)=P(6,I)*COS(P(8,I)) PJET(2,I)=P(6,I)*SIN(P(8,I)) PJET(3,I)=P(6,I)*SINH(P(7,I)) PJET(4,I)=P(6,I)*COSH(P(7,I)) ENDIF JET(I)=I 210 CONTINUE C---KEEP MERGING UNTIL YMAC TO FIND THE FATE OF EACH JET 300 IF (ETSQ*KT(N).LT.ROUND*YMAC) THEN IF (HIST(N).LE.NMAX) THEN IMIN=0 JMIN=HIST(N) NSUB=NSUB-1 ELSE IMIN=HIST(N)/NMAX JMIN=HIST(N)-IMIN*NMAX IF (ETSQ*KTLAST(N).LT.ROUND*YMAC) NSUB=NSUB-1 ENDIF DO 310 I=1,NJET IF (JET(I).EQ.JMIN) JET(I)=IMIN IF (JET(I).EQ.N) JET(I)=JMIN 310 CONTINUE N=N-1 IF (N.GT.0) GOTO 300 ENDIF RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTINCL(RECO,PP,NN,PJET,JET,NJET,*) IMPLICIT NONE C---RECONSTRUCT KINEMATICS OF JET SYSTEM, WHICH HAS ALREADY BEEN C ANALYSED BY KTCLUS ACCORDING TO THE INCLUSIVE JET DEFINITION. NOTE C THAT NO CONSISTENCY CHECK IS MADE: USER IS TRUSTED TO USE THE SAME C PP VALUES AS FOR KTCLUS C C RECO = INPUT : RECOMBINATION SCHEME (NEED NOT BE SAME AS KTCLUS) C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E C NN = INPUT : NUMBER OF PARTICLES C PJET(I,J)=OUTPUT : 4-MOMENTUM OF Jth JET AT SCALE YCUT C JET(J) =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE C NJET =OUTPUT : THE NUMBER OF JETS C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER NMAX,RECO,NUM,N,NN,NJET,JET(*),HIST,IMIN,JMIN,I,J PARAMETER (NMAX=512) DOUBLE PRECISION PP(4,*),PJET(4,*) DOUBLE PRECISION P,KT,KTP,KTS,ETOT,RSQ,KTLAST COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM C---CHECK INPUT IF (RECO.LT.1.OR.RECO.GT.3) CALL KTWARN('KTINCL',100,*999) C---COPY PP TO P N=NN IF (NUM.NE.NN) CALL KTWARN('KTINCL',101,*999) CALL KTCOPY(PP,N,P,(RECO.NE.1)) C---INITIALLY EVERY PARTICLE IS IN ITS OWN JET DO 100 I=1,NN JET(I)=I 100 CONTINUE C---KEEP MERGING TO THE BITTER END NJET=0 200 IF (N.GT.0) THEN IF (HIST(N).LE.NMAX) THEN IMIN=0 JMIN=HIST(N) NJET=NJET+1 IF (RECO.EQ.1) THEN DO 300 J=1,4 PJET(J,NJET)=P(J,JMIN) 300 CONTINUE ELSE PJET(1,NJET)=P(6,JMIN)*COS(P(8,JMIN)) PJET(2,NJET)=P(6,JMIN)*SIN(P(8,JMIN)) PJET(3,NJET)=P(6,JMIN)*SINH(P(7,JMIN)) PJET(4,NJET)=P(6,JMIN)*COSH(P(7,JMIN)) ENDIF CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0) ELSE IMIN=HIST(N)/NMAX JMIN=HIST(N)-IMIN*NMAX CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO) CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0) ENDIF DO 400 I=1,NN IF (JET(I).EQ.JMIN) JET(I)=IMIN IF (JET(I).EQ.N) JET(I)=JMIN IF (JET(I).EQ.0) JET(I)=-NJET 400 CONTINUE N=N-1 GOTO 200 ENDIF C---FINALLY EVERY PARTICLE MUST BE IN AN INCLUSIVE JET DO 500 I=1,NN C---IF THERE ARE ANY UNASSIGNED PARTICLES SOMETHING MUST HAVE GONE WRONG IF (JET(I).GE.0) CALL KTWARN('KTINCL',102,*999) JET(I)=-JET(I) 500 CONTINUE RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTISUB(N,NY,YCUT,NSUB,*) IMPLICIT NONE C---COUNT THE NUMBER OF SUB-JETS IN THE Nth INCLUSIVE JET OF AN EVENT C THAT HAS ALREADY BEEN ANALYSED BY KTCLUS. C C N = INPUT : WHICH INCLUSIVE JET TO USE C NY = INPUT : NUMBER OF YCUT VALUES C YCUT(J) = INPUT : Y VALUES AT WHICH NUMBERS OF SUB-JETS ARE COUNTED C NSUB(J) =OUTPUT : NUMBER OF SUB-JETS AT YCUT(J) C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER N,NY,NSUB(NY),NMAX,HIST,I,J,NUM,NM PARAMETER (NMAX=512) DOUBLE PRECISION YCUT(NY),ETOT,RSQ,P,KT,KTP,KTS,KTLAST,ROUND,EPS PARAMETER (ROUND=0.99999D0) COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM DATA EPS/1D-6/ DO 100 I=1,NY NSUB(I)=0 100 CONTINUE C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET NM=0 J=0 DO 110 I=NUM,1,-1 IF (HIST(I).LE.NMAX) J=J+1 IF (J.EQ.N) THEN NM=I GOTO 120 ENDIF 110 CONTINUE 120 CONTINUE C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS IF (NM.EQ.0) CALL KTWARN('KTISUB',100,*999) DO 210 I=NUM,1,-1 DO 200 J=1,NY IF (NSUB(J).EQ.0.AND.RSQ*KT(I).GE.ROUND*YCUT(J)*KT(NM)) & NSUB(J)=I IF (NSUB(J).NE.0.AND.ABS(KTLAST(I)-KTLAST(NM)).GT.EPS) & NSUB(J)=NSUB(J)-1 200 CONTINUE 210 CONTINUE RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTIJOI(N,Y,*) IMPLICIT NONE C---GIVE SAME INFORMATION AS LAST CALL TO KTCLUS EXCEPT THAT ONLY C MERGES OF TWO SUB-JETS INSIDE THE Nth INCLUSIVE JET ARE RECORDED C C N = INPUT : WHICH INCLUSIVE JET TO USE C Y(J) =OUTPUT : Y VALUE WHERE JET CHANGED FROM HAVING C J+1 SUB-JETS TO HAVING J C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER NMAX,HIST,NUM,I,J,N,NM PARAMETER (NMAX=512) DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,Y(*),KTLAST,EPS COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM DATA EPS/1D-6/ C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET NM=0 J=0 DO 100 I=NUM,1,-1 IF (HIST(I).LE.NMAX) J=J+1 IF (J.EQ.N) THEN NM=I GOTO 105 ENDIF 100 CONTINUE 105 CONTINUE C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS IF (NM.EQ.0) CALL KTWARN('KTIJOI',100,*999) J=1 DO 110 I=1,NUM IF (HIST(I).GT.NMAX.AND.ABS(KTLAST(I)-KTLAST(NM)).LT.EPS) THEN Y(J)=RSQ*KT(I)/KT(NM) J=J+1 ENDIF 110 CONTINUE DO 200 I=J,NUM Y(I)=0 200 CONTINUE RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTIREC(RECO,PP,NN,N,YCUT,PSUB,NSUB,*) IMPLICIT NONE C---RECONSTRUCT KINEMATICS OF SUB-JET SYSTEM IN THE Nth INCLUSIVE JET C OF AN EVENT THAT HAS ALREADY BEEN ANALYSED BY KTCLUS C C RECO = INPUT : RECOMBINATION SCHEME (NEED NOT BE SAME AS KTCLUS) C PP(I,J) = INPUT : 4-MOMENTUM OF Jth PARTICLE: I=1,4 => PX,PY,PZ,E C NN = INPUT : NUMBER OF PARTICLES C N = INPUT : WHICH INCLUSIVE JET TO USE C YCUT = INPUT : Y VALUE AT WHICH TO RECONSTRUCT JET MOMENTA C PSUB(I,J)=OUTPUT : 4-MOMENTUM OF Jth SUB-JET AT SCALE YCUT C NSUB =OUTPUT : THE NUMBER OF SUB-JETS C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT THE MOMENTA ARE DECLARED DOUBLE PRECISION, C AND ALL OTHER FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER NMAX,RECO,NUM,NN,NJET,NSUB,JET,HIST,I,J,N,NM PARAMETER (NMAX=512) DOUBLE PRECISION PP(4,*),PSUB(4,*) DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,YCUT,YMAC,KTLAST COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM DIMENSION JET(NMAX) C---FIND WHICH MERGING CORRESPONDS TO THE NTH INCLUSIVE JET NM=0 J=0 DO 100 I=NUM,1,-1 IF (HIST(I).LE.NMAX) J=J+1 IF (J.EQ.N) THEN NM=I GOTO 110 ENDIF 100 CONTINUE 110 CONTINUE C---GIVE UP IF THERE ARE LESS THAN N INCLUSIVE JETS IF (NM.EQ.0) CALL KTWARN('KTIREC',102,*999) C---RECONSTRUCT THE JETS AT THE APPROPRIATE SCALE ECUT=SQRT(KT(NM)/RSQ) YMAC=RSQ CALL KTRECO(RECO,PP,NN,ECUT,YCUT,YMAC,PSUB,JET,NJET,NSUB,*999) C---GET RID OF THE ONES THAT DO NOT END UP IN THE JET WE WANT NSUB=0 DO 210 I=1,NJET IF (JET(I).EQ.HIST(NM)) THEN NSUB=NSUB+1 DO 200 J=1,4 PSUB(J,NSUB)=PSUB(J,I) 200 CONTINUE ENDIF 210 CONTINUE RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTWICH(ECUT,YCUT,JET,NJET,*) IMPLICIT NONE C---GIVE A LIST OF WHICH JET EACH ORIGINAL PARTICLE ENDED UP IN AT SCALE C YCUT, TOGETHER WITH THE NUMBER OF JETS AT THAT SCALE. C C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED C YCUT = INPUT : Y VALUE AT WHICH TO DEFINE JETS C JET(J) =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE, C SET TO ZERO IF IT WAS PUT INTO THE BEAM JETS C NJET =OUTPUT : THE NUMBER OF JETS AT SCALE YCUT (SO JET() C ENTRIES WILL BE IN THE RANGE 0 -> NJET) C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER JET(*),NJET,NTEMP DOUBLE PRECISION ECUT,YCUT CALL KTWCHS(ECUT,YCUT,YCUT,JET,NJET,NTEMP,*999) RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTWCHS(ECUT,YCUT,YMAC,JET,NJET,NSUB,*) IMPLICIT NONE C---GIVE A LIST OF WHICH SUB-JET EACH ORIGINAL PARTICLE ENDED UP IN AT C SCALE YCUT, WITH MACRO-JET SCALE YMAC, TOGETHER WITH THE NUMBER OF C JETS AT SCALE YCUT AND THE NUMBER OF THEM WHICH ARE SUB-JETS. C C ECUT = INPUT : DENOMINATOR OF KT MEASURE. IF ZERO, ETOT IS USED C YCUT = INPUT : Y VALUE AT WHICH TO DEFINE JETS C YMAC = INPUT : Y VALUE AT WHICH TO DEFINE MACRO-JETS C JET(J) =OUTPUT : THE JET WHICH CONTAINS THE Jth PARTICLE, C SET TO ZERO IF IT WAS PUT INTO THE BEAM JETS C NJET =OUTPUT : THE NUMBER OF JETS AT SCALE YCUT (SO JET() C ENTRIES WILL BE IN THE RANGE 0 -> NJET) C NSUB =OUTPUT : THE NUMBER OF SUB-JETS AT SCALE YCUT, WITH C MACRO-JETS DEFINED AT SCALE YMAC (SO ONLY NSUB C OF THE JETS 1 -> NJET WILL APPEAR IN JET()) C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT ALL FLOATING POINT VARIABLES ARE DECLARED DOUBLE PRECISION C INTEGER NMAX,JET(*),NJET,NSUB,HIST,NUM,I,J,JSUB PARAMETER (NMAX=512) DOUBLE PRECISION P1(4,NMAX),P2(4,NMAX) DOUBLE PRECISION ECUT,YCUT,YMAC,ZERO,ETOT,RSQ,P,KTP,KTS,KT,KTLAST COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX), & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM DIMENSION JSUB(NMAX) C---THE MOMENTA HAVE TO BEEN GIVEN LEGAL VALUES, C EVEN THOUGH THEY WILL NEVER BE USED DATA ((P1(J,I),I=1,NMAX),J=1,4),ZERO & /NMAX*1,NMAX*0,NMAX*0,NMAX*1,0/ C---FIRST GET A LIST OF WHICH PARTICLE IS IN WHICH JET AT YCUT CALL KTRECO(1,P1,NUM,ECUT,ZERO,YCUT,P2,JET,NJET,NSUB,*999) C---THEN FIND OUT WHICH JETS ARE SUBJETS CALL KTRECO(1,P1,NUM,ECUT,YCUT,YMAC,P2,JSUB,NJET,NSUB,*999) C---AND MODIFY JET() ACCORDINGLY DO 10 I=1,NUM IF (JET(I).NE.0) THEN IF (JSUB(JET(I)).EQ.0) JET(I)=0 ENDIF 10 CONTINUE RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTFRAM(IOPT,CMF,SIGN,Z,XZ,N,P,Q,*) IMPLICIT NONE C---BOOST PARTICLES IN P TO/FROM FRAME GIVEN BY CMF, Z, XZ. C---IN THIS FRAME CMZ IS STATIONARY, C Z IS ALONG THE (SIGN)Z-AXIS (SIGN=+ OR -) C XZ IS IN THE X-Z PLANE (WITH POSITIVE X COMPONENT) C---IF Z HAS LENGTH ZERO, OR SIGN=0, NO ROTATION IS PERFORMED C---IF XZ HAS ZERO COMPONENT PERPENDICULAR TO Z IN THAT FRAME, C NO AZIMUTHAL ROTATION IS PERFORMED C C IOPT = INPUT : 0=TO FRAME, 1=FROM FRAME C CMF(I) = INPUT : 4-MOMENTUM WHICH IS STATIONARY IN THE FRAME C SIGN = INPUT : DIRECTION OF Z IN THE FRAME, NOTE THAT C ONLY ITS SIGN IS USED, NOT ITS MAGNITUDE C Z(I) = INPUT : 4-MOMENTUM WHICH LIES ON THE (SIGN)Z-AXIS C XZ(I) = INPUT : 4-MOMENTUM WHICH LIES IN THE X-Z PLANE C N = INPUT : NUMBER OF PARTICLES IN P C P(I,J) = INPUT : 4-MOMENTUM OF JTH PARTICLE BEFORE C Q(I,J) = OUTPUT : 4-MOMENTUM OF JTH PARTICLE AFTER C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED C C NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION C C NOTE THAT IT IS SAFE TO CALL WITH P=Q C INTEGER IOPT,I,N DOUBLE PRECISION CMF(4),SIGN,Z(4),XZ(4),P(4,N),Q(4,N), & R(4,4),NEW(4),OLD(4) IF (IOPT.LT.0.OR.IOPT.GT.1) CALL KTWARN('KTFRAM',200,*999) C---FIND BOOST TO GET THERE FROM LAB CALL KTUNIT(R) CALL KTLBST(0,R,CMF,*999) C---FIND ROTATION TO PUT BOOSTED Z ON THE (SIGN)Z AXIS IF (SIGN.NE.0) THEN CALL KTVMUL(R,Z,OLD) IF (OLD(1).NE.0.OR.OLD(2).NE.0.OR.OLD(3).NE.0) THEN NEW(1)=0 NEW(2)=0 NEW(3)=SIGN NEW(4)=ABS(SIGN) CALL KTRROT(R,OLD,NEW,*999) C---FIND ROTATION TO PUT BOOSTED AND ROTATED XZ INTO X-Z PLANE CALL KTVMUL(R,XZ,OLD) IF (OLD(1).NE.0.OR.OLD(2).NE.0) THEN NEW(1)=1 NEW(2)=0 NEW(3)=0 NEW(4)=1 OLD(3)=0 C---NOTE THAT A POTENTIALLY AWKWARD SPECIAL CASE IS AVERTED, BECAUSE IF C OLD AND NEW ARE EXACTLY BACK-TO-BACK, THE ROTATION AXIS IS UNDEFINED C BUT IN THAT CASE KTRROT WILL USE THE Z AXIS, AS REQUIRED CALL KTRROT(R,OLD,NEW,*999) ENDIF ENDIF ENDIF C---INVERT THE TRANSFORMATION IF NECESSARY IF (IOPT.EQ.1) CALL KTINVT(R,R) C---APPLY THE RESULT TO ALL THE VECTORS DO 30 I=1,N CALL KTVMUL(R,P(1,I),Q(1,I)) 30 CONTINUE RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTBREI(IOPT,PLEP,PHAD,POUT,N,P,Q,*) IMPLICIT NONE C---BOOST PARTICLES IN P TO/FROM BREIT FRAME C C IOPT = INPUT : 0/2=TO BREIT FRAME, 1/3=FROM BREIT FRAME C 0/1=NO AZIMUTHAL ROTATION AFTERWARDS C 2/3=LEPTON PLANE ROTATED INTO THE X-Z PLANE C PLEP = INPUT : MOMENTUM OF INCOMING LEPTON IN +Z DIRECTION C PHAD = INPUT : MOMENTUM OF INCOMING HADRON IN +Z DIRECTION C POUT(I) = INPUT : 4-MOMENTUM OF OUTGOING LEPTON C N = INPUT : NUMBER OF PARTICLES IN P C P(I,J) = INPUT : 4-MOMENTUM OF JTH PARTICLE BEFORE C Q(I,J) = OUTPUT : 4-MOMENTUM OF JTH PARTICLE AFTER C LAST ARGUMENT IS LABEL TO JUMP TO IF FOR ANY REASON THE EVENT C COULD NOT BE PROCESSED (MOST LIKELY DUE TO PARTICLES HAVING SMALLER C ENERGY THAN MOMENTUM) C C NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION C C NOTE THAT IT IS SAFE TO CALL WITH P=Q C INTEGER IOPT,N DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N), & CMF(4),Z(4),XZ(4),DOT,QDQ C---CHECK INPUT IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTBREI',200,*999) C---FIND 4-MOMENTUM OF BREIT FRAME (TIMES AN ARBITRARY FACTOR) DOT=ABS(PHAD)*(ABS(PLEP)-POUT(4))-PHAD*(PLEP-POUT(3)) QDQ=(ABS(PLEP)-POUT(4))**2-(PLEP-POUT(3))**2-POUT(2)**2-POUT(1)**2 CMF(1)=DOT*( -POUT(1)) CMF(2)=DOT*( -POUT(2)) CMF(3)=DOT*( PLEP -POUT(3))-QDQ* PHAD CMF(4)=DOT*(ABS(PLEP)-POUT(4))-QDQ*ABS(PHAD) C---FIND ROTATION TO PUT INCOMING HADRON BACK ON Z-AXIS Z(1)=0 Z(2)=0 Z(3)=PHAD Z(4)=ABS(PHAD) XZ(1)=0 XZ(2)=0 XZ(3)=0 XZ(4)=0 C---DO THE BOOST IF (IOPT.LE.1) THEN CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999) ELSE CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999) ENDIF RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTHADR(IOPT,PLEP,PHAD,POUT,N,P,Q,*) IMPLICIT NONE C---BOOST PARTICLES IN P TO/FROM HADRONIC CMF C C ARGUMENTS ARE EXACTLY AS FOR KTBREI C C NOTE THAT ALL MOMENTA ARE DOUBLE PRECISION C C NOTE THAT IT IS SAFE TO CALL WITH P=Q C INTEGER IOPT,N DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N), & CMF(4),Z(4),XZ(4) C---CHECK INPUT IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTHADR',200,*999) C---FIND 4-MOMENTUM OF HADRONIC CMF CMF(1)= -POUT(1) CMF(2)= -POUT(2) CMF(3)= PLEP -POUT(3)+ PHAD CMF(4)=ABS(PLEP)-POUT(4)+ABS(PHAD) C---FIND ROTATION TO PUT INCOMING HADRON BACK ON Z-AXIS Z(1)=0 Z(2)=0 Z(3)=PHAD Z(4)=ABS(PHAD) XZ(1)=0 XZ(2)=0 XZ(3)=0 XZ(4)=0 C---DO THE BOOST IF (IOPT.LE.1) THEN CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999) ELSE CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999) ENDIF RETURN 999 RETURN 1 END C----------------------------------------------------------------------- FUNCTION KTPAIR(ANGL,P,Q,ANGLE) IMPLICIT NONE INTEGER IKTTOPT COMMON/KTSWITCH/IKTTOPT SAVE /KTSWITCH/ C---CALCULATE LOCAL KT OF PAIR, USING ANGULAR SCHEME: C 1=>ANGULAR, 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi) C WHERE f(eta,phi)=2(COSH(eta)-COS(phi)) IS THE QCD EMISSION METRIC C---IF ANGLE<0, IT IS SET TO THE ANGULAR PART OF THE LOCAL KT ON RETURN C IF ANGLE>0, IT IS USED INSTEAD OF THE ANGULAR PART OF THE LOCAL KT INTEGER ANGL DOUBLE PRECISION P(9),Q(9),KTPAIR,R,KTMDPI,ANGLE,ETA,PHI,ESQ C---COMPONENTS OF MOMENTA ARE PX,PY,PZ,E,1/P,PT,ETA,PHI,PT**2 R=ANGLE IF (ANGL.EQ.1) THEN IF (R.LE.0) THEN R=2D0*(1D0-(P(1)*Q(1)+P(2)*Q(2)+P(3)*Q(3))*(P(5)*Q(5))) ENDIF IF(IKTTOPT.EQ.1) THEN ESQ=(P(4)*Q(4)/(P(4)+Q(4)))**2 ELSE ESQ=MIN(P(4),Q(4))**2 ENDIF ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN IF (R.LE.0) THEN ETA=P(7)-Q(7) PHI=KTMDPI(P(8)-Q(8)) IF (ANGL.EQ.2) THEN R=ETA**2+PHI**2 ELSE R=2*(COSH(ETA)-COS(PHI)) ENDIF ENDIF ESQ=MIN(P(9),Q(9)) ELSE CALL KTWARN('KTPAIR',200,*999) STOP ENDIF KTPAIR=ESQ*R IF (ANGLE.LT.0) ANGLE=R 999 END C----------------------------------------------------------------------- FUNCTION KTSING(ANGL,TYPE,P) IMPLICIT NONE C---CALCULATE KT OF PARTICLE, USING ANGULAR SCHEME: C 1=>ANGULAR, 2=>DeltaR, 3=>f(DeltaEta,DeltaPhi) C---TYPE=1 FOR E+E-, 2 FOR EP, 3 FOR PE, 4 FOR PP C FOR EP, PROTON DIRECTION IS DEFINED AS -Z C FOR PE, PROTON DIRECTION IS DEFINED AS +Z INTEGER ANGL,TYPE DOUBLE PRECISION P(9),KTSING,COSTH,R,SMALL DATA SMALL/1D-4/ IF (ANGL.EQ.1) THEN COSTH=P(3)*P(5) IF (TYPE.EQ.2) THEN COSTH=-COSTH ELSEIF (TYPE.EQ.4) THEN COSTH=ABS(COSTH) ELSEIF (TYPE.NE.1.AND.TYPE.NE.3) THEN CALL KTWARN('KTSING',200,*999) STOP ENDIF R=2*(1-COSTH) C---IF CLOSE TO BEAM, USE APPROX 2*(1-COS(THETA))=SIN**2(THETA) IF (R.LT.SMALL) R=(P(1)**2+P(2)**2)*P(5)**2 KTSING=P(4)**2*R ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN KTSING=P(9) ELSE CALL KTWARN('KTSING',201,*999) STOP ENDIF 999 END C----------------------------------------------------------------------- SUBROUTINE KTPMIN(A,NMAX,N,IMIN,JMIN) IMPLICIT NONE C---FIND THE MINIMUM MEMBER OF A(NMAX,NMAX) WITH IMIN < JMIN <= N INTEGER NMAX,N,IMIN,JMIN,KMIN,I,J,K C---REMEMBER THAT A(X+(Y-1)*NMAX)=A(X,Y) C THESE LOOPING VARIABLES ARE J=Y-2, I=X+(Y-1)*NMAX DOUBLE PRECISION A(*),AMIN K=1+NMAX KMIN=K AMIN=A(KMIN) DO 110 J=0,N-2 DO 100 I=K,K+J IF (A(I).LT.AMIN) THEN KMIN=I AMIN=A(KMIN) ENDIF 100 CONTINUE K=K+NMAX 110 CONTINUE JMIN=KMIN/NMAX+1 IMIN=KMIN-(JMIN-1)*NMAX END C----------------------------------------------------------------------- SUBROUTINE KTSMIN(A,NMAX,N,IMIN) IMPLICIT NONE C---FIND THE MINIMUM MEMBER OF A INTEGER N,NMAX,IMIN,I DOUBLE PRECISION A(NMAX) IMIN=1 DO 100 I=1,N IF (A(I).LT.A(IMIN)) IMIN=I 100 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE KTCOPY(A,N,B,ONSHLL) IMPLICIT NONE C---COPY FROM A TO B. 5TH=1/(3-MTM), 6TH=PT, 7TH=ETA, 8TH=PHI, 9TH=PT**2 C IF ONSHLL IS .TRUE. PARTICLE ENTRIES ARE PUT ON-SHELL BY SETTING E=P INTEGER I,N DOUBLE PRECISION A(4,N) LOGICAL ONSHLL DOUBLE PRECISION B(9,N),ETAMAX,SINMIN,EPS DATA ETAMAX,SINMIN,EPS/10,0,1D-6/ C---SINMIN GETS CALCULATED ON FIRST CALL IF (SINMIN.EQ.0) SINMIN=1/COSH(ETAMAX) DO 100 I=1,N B(1,I)=A(1,I) B(2,I)=A(2,I) B(3,I)=A(3,I) B(4,I)=A(4,I) B(5,I)=SQRT(A(1,I)**2+A(2,I)**2+A(3,I)**2) IF (ONSHLL) B(4,I)=B(5,I) IF (B(5,I).EQ.0) B(5,I)=1D-10 B(5,I)=1/B(5,I) B(9,I)=A(1,I)**2+A(2,I)**2 B(6,I)=SQRT(B(9,I)) B(7,I)=B(6,I)*B(5,I) IF (B(7,I).GT.SINMIN) THEN B(7,I)=A(4,I)**2-A(3,I)**2 IF (B(7,I).LE.EPS*B(4,I)**2.OR.ONSHLL) B(7,I)=B(9,I) B(7,I)=LOG((B(4,I)+ABS(B(3,I)))**2/B(7,I))/2 ELSE B(7,I)=ETAMAX+2 ENDIF B(7,I)=SIGN(B(7,I),B(3,I)) IF (A(1,I).EQ.0 .AND. A(2,I).EQ.0) THEN B(8,I)=0 ELSE B(8,I)=ATAN2(A(2,I),A(1,I)) ENDIF 100 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE KTMERG(P,KTP,KTS,NMAX,I,J,N,TYPE,ANGL,MONO,RECO) IMPLICIT NONE C---MERGE THE Jth PARTICLE IN P INTO THE Ith PARTICLE C J IS ASSUMED GREATER THAN I. P CONTAINS N PARTICLES BEFORE MERGING. C---ALSO RECALCULATING THE CORRESPONDING KTP AND KTS VALUES IF MONO.GT.0 C FROM THE RECOMBINED ANGULAR MEASURES IF MONO.GT.1 C---NOTE THAT IF MONO.LE.0, TYPE AND ANGL ARE NOT USED INTEGER ANGL,RECO,TYPE,I,J,K,N,NMAX,MONO DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),PT,PTT, & KTMDPI,KTUP,PI,PJ,ANG,KTPAIR,KTSING,ETAMAX,EPS KTUP(I,J)=KTP(MAX(I,J),MIN(I,J)) DATA ETAMAX,EPS/10,1D-6/ IF (J.LE.I) CALL KTWARN('KTMERG',200,*999) C---COMBINE ANGULAR MEASURES IF NECESSARY IF (MONO.GT.1) THEN DO 100 K=1,N IF (K.NE.I.AND.K.NE.J) THEN IF (RECO.EQ.1) THEN PI=P(4,I) PJ=P(4,J) ELSEIF (RECO.EQ.2) THEN PI=P(6,I) PJ=P(6,J) ELSEIF (RECO.EQ.3) THEN PI=P(9,I) PJ=P(9,J) ELSE CALL KTWARN('KTMERG',201,*999) STOP ENDIF IF (PI.EQ.0.AND.PJ.EQ.0) THEN PI=1 PJ=1 ENDIF KTP(MAX(I,K),MIN(I,K))= & (PI*KTUP(I,K)+PJ*KTUP(J,K))/(PI+PJ) ENDIF 100 CONTINUE ENDIF IF (RECO.EQ.1) THEN C---VECTOR ADDITION P(1,I)=P(1,I)+P(1,J) P(2,I)=P(2,I)+P(2,J) P(3,I)=P(3,I)+P(3,J) P(4,I)=P(4,I)+P(4,J) P(5,I)=SQRT(P(1,I)**2+P(2,I)**2+P(3,I)**2) IF (P(5,I).EQ.0) THEN P(5,I)=1 ELSE P(5,I)=1/P(5,I) ENDIF ELSEIF (RECO.EQ.2) THEN C---PT WEIGHTED ETA-PHI ADDITION PT=P(6,I)+P(6,J) IF (PT.EQ.0) THEN PTT=1 ELSE PTT=1/PT ENDIF P(7,I)=(P(6,I)*P(7,I)+P(6,J)*P(7,J))*PTT P(8,I)=KTMDPI(P(8,I)+P(6,J)*PTT*KTMDPI(P(8,J)-P(8,I))) P(6,I)=PT P(9,I)=PT**2 ELSEIF (RECO.EQ.3) THEN C---PT**2 WEIGHTED ETA-PHI ADDITION PT=P(9,I)+P(9,J) IF (PT.EQ.0) THEN PTT=1 ELSE PTT=1/PT ENDIF P(7,I)=(P(9,I)*P(7,I)+P(9,J)*P(7,J))*PTT P(8,I)=KTMDPI(P(8,I)+P(9,J)*PTT*KTMDPI(P(8,J)-P(8,I))) P(6,I)=P(6,I)+P(6,J) P(9,I)=P(6,I)**2 ELSE CALL KTWARN('KTMERG',202,*999) STOP ENDIF C---IF MONO.GT.0 CALCULATE NEW KT MEASURES. IF MONO.GT.1 USE ANGULAR ONES. IF (MONO.LE.0) RETURN C---CONVERTING BETWEEN 4-MTM AND PT,ETA,PHI IF NECESSARY IF (ANGL.NE.1.AND.RECO.EQ.1) THEN P(9,I)=P(1,I)**2+P(2,I)**2 P(7,I)=P(4,I)**2-P(3,I)**2 IF (P(7,I).LE.EPS*P(4,I)**2) P(7,I)=P(9,I) IF (P(7,I).GT.0) THEN P(7,I)=LOG((P(4,I)+ABS(P(3,I)))**2/P(7,I))/2 IF (P(7,I).GT.ETAMAX) P(7,I)=ETAMAX+2 ELSE P(7,I)=ETAMAX+2 ENDIF P(7,I)=SIGN(P(7,I),P(3,I)) IF (P(1,I).NE.0.AND.P(2,I).NE.0) THEN P(8,I)=ATAN2(P(2,I),P(1,I)) ELSE P(8,I)=0 ENDIF ELSEIF (ANGL.EQ.1.AND.RECO.NE.1) THEN P(1,I)=P(6,I)*COS(P(8,I)) P(2,I)=P(6,I)*SIN(P(8,I)) P(3,I)=P(6,I)*SINH(P(7,I)) P(4,I)=P(6,I)*COSH(P(7,I)) IF (P(4,I).NE.0) THEN P(5,I)=1/P(4,I) ELSE P(5,I)=1 ENDIF ENDIF ANG=0 DO 200 K=1,N IF (K.NE.I.AND.K.NE.J) THEN IF (MONO.GT.1) ANG=KTUP(I,K) KTP(MIN(I,K),MAX(I,K))= & KTPAIR(ANGL,P(1,I),P(1,K),ANG) ENDIF 200 CONTINUE KTS(I)=KTSING(ANGL,TYPE,P(1,I)) 999 END C----------------------------------------------------------------------- SUBROUTINE KTMOVE(P,KTP,KTS,NMAX,N,J,IOPT) IMPLICIT NONE C---MOVE THE Nth PARTICLE IN P TO THE Jth POSITION C---ALSO MOVING KTP AND KTS IF IOPT.GT.0 INTEGER I,J,N,NMAX,IOPT DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX) DO 100 I=1,9 P(I,J)=P(I,N) 100 CONTINUE IF (IOPT.LE.0) RETURN DO 110 I=1,J-1 KTP(I,J)=KTP(I,N) KTP(J,I)=KTP(N,I) 110 CONTINUE DO 120 I=J+1,N-1 KTP(J,I)=KTP(I,N) KTP(I,J)=KTP(N,I) 120 CONTINUE KTS(J)=KTS(N) END C----------------------------------------------------------------------- SUBROUTINE KTUNIT(R) IMPLICIT NONE C SET R EQUAL TO THE 4 BY 4 IDENTITY MATRIX DOUBLE PRECISION R(4,4) INTEGER I,J DO 20 I=1,4 DO 10 J=1,4 R(I,J)=0 IF (I.EQ.J) R(I,J)=1 10 CONTINUE 20 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE KTLBST(IOPT,R,A,*) IMPLICIT NONE C PREMULTIPLY R BY THE 4 BY 4 MATRIX TO C LORENTZ BOOST TO/FROM THE CM FRAME OF A C IOPT=0 => TO C IOPT=1 => FROM C C LAST ARGUMENT IS LABEL TO JUMP TO IF A IS NOT TIME-LIKE C INTEGER IOPT,I,J DOUBLE PRECISION R(4,4),A(4),B(4),C(4,4),M DO 10 I=1,4 B(I)=A(I) 10 CONTINUE M=B(4)**2-B(1)**2-B(2)**2-B(3)**2 IF (M.LE.0) CALL KTWARN('KTLBST',100,*999) M=SQRT(M) B(4)=B(4)+M M=1/(M*B(4)) IF (IOPT.EQ.0) THEN B(4)=-B(4) ELSEIF (IOPT.NE.1) THEN CALL KTWARN('KTLBST',200,*999) STOP ENDIF DO 30 I=1,4 DO 20 J=1,4 C(I,J)=B(I)*B(J)*M IF (I.EQ.J) C(I,J)=C(I,J)+1 20 CONTINUE 30 CONTINUE C(4,4)=C(4,4)-2 CALL KTMMUL(C,R,R) RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTRROT(R,A,B,*) IMPLICIT NONE C PREMULTIPLY R BY THE 4 BY 4 MATRIX TO C ROTATE FROM VECTOR A TO VECTOR B BY THE SHORTEST ROUTE C IF THEY ARE EXACTLY BACK-TO-BACK, THE ROTATION AXIS IS THE VECTOR C WHICH IS PERPENDICULAR TO THEM AND THE X AXIS, UNLESS THEY ARE C PERPENDICULAR TO THE Y AXIS, WHEN IT IS THE VECTOR WHICH IS C PERPENDICULAR TO THEM AND THE Y AXIS. C NOTE THAT THESE CONDITIONS GUARANTEE THAT IF BOTH ARE PERPENDICULAR C TO THE Z AXIS, IT WILL BE USED AS THE ROTATION AXIS. C C LAST ARGUMENT IS LABEL TO JUMP TO IF EITHER HAS LENGTH ZERO C DOUBLE PRECISION R(4,4),M(4,4),A(4),B(4),C(4),D(4),AL,BL,CL,DL,EPS C---SQRT(2*EPS) IS THE ANGLE IN RADIANS OF THE SMALLEST ALLOWED ROTATION C NOTE THAT IF YOU CONVERT THIS PROGRAM TO SINGLE PRECISION, YOU WILL C NEED TO INCREASE EPS TO AROUND 0.5E-4 PARAMETER (EPS=0.5D-6) AL=A(1)**2+A(2)**2+A(3)**2 BL=B(1)**2+B(2)**2+B(3)**2 IF (AL.LE.0.OR.BL.LE.0) CALL KTWARN('KTRROT',100,*999) AL=1/SQRT(AL) BL=1/SQRT(BL) CL=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*AL*BL C---IF THEY ARE COLLINEAR, DON'T NEED TO DO ANYTHING IF (CL.GE.1-EPS) THEN RETURN C---IF THEY ARE BACK-TO-BACK, USE THE AXIS PERP TO THEM AND X AXIS ELSEIF (CL.LE.-1+EPS) THEN IF (ABS(B(2)).GT.EPS) THEN C(1)= 0 C(2)=-B(3) C(3)= B(2) C---UNLESS THEY ARE PERPENDICULAR TO THE Y AXIS, ELSE C(1)= B(3) C(2)= 0 C(3)=-B(1) ENDIF C---OTHERWISE FIND ROTATION AXIS ELSE C(1)=A(2)*B(3)-A(3)*B(2) C(2)=A(3)*B(1)-A(1)*B(3) C(3)=A(1)*B(2)-A(2)*B(1) ENDIF CL=C(1)**2+C(2)**2+C(3)**2 IF (CL.LE.0) CALL KTWARN('KTRROT',101,*999) CL=1/SQRT(CL) C---FIND ROTATION TO INTERMEDIATE AXES FROM A D(1)=A(2)*C(3)-A(3)*C(2) D(2)=A(3)*C(1)-A(1)*C(3) D(3)=A(1)*C(2)-A(2)*C(1) DL=AL*CL M(1,1)=A(1)*AL M(1,2)=A(2)*AL M(1,3)=A(3)*AL M(1,4)=0 M(2,1)=C(1)*CL M(2,2)=C(2)*CL M(2,3)=C(3)*CL M(2,4)=0 M(3,1)=D(1)*DL M(3,2)=D(2)*DL M(3,3)=D(3)*DL M(3,4)=0 M(4,1)=0 M(4,2)=0 M(4,3)=0 M(4,4)=1 CALL KTMMUL(M,R,R) C---AND ROTATION FROM INTERMEDIATE AXES TO B D(1)=B(2)*C(3)-B(3)*C(2) D(2)=B(3)*C(1)-B(1)*C(3) D(3)=B(1)*C(2)-B(2)*C(1) DL=BL*CL M(1,1)=B(1)*BL M(2,1)=B(2)*BL M(3,1)=B(3)*BL M(1,2)=C(1)*CL M(2,2)=C(2)*CL M(3,2)=C(3)*CL M(1,3)=D(1)*DL M(2,3)=D(2)*DL M(3,3)=D(3)*DL CALL KTMMUL(M,R,R) RETURN 999 RETURN 1 END C----------------------------------------------------------------------- SUBROUTINE KTVMUL(M,A,B) IMPLICIT NONE C 4 BY 4 MATRIX TIMES 4 VECTOR: B=M*A. C ALL ARE DOUBLE PRECISION C IT IS SAFE TO CALL WITH B=A C FIRST SUBSCRIPT=ROWS, SECOND=COLUMNS DOUBLE PRECISION M(4,4),A(4),B(4),C(4) INTEGER I,J DO 20 I=1,4 C(I)=0 DO 10 J=1,4 C(I)=C(I)+M(I,J)*A(J) 10 CONTINUE 20 CONTINUE DO 30 I=1,4 B(I)=C(I) 30 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE KTMMUL(A,B,C) IMPLICIT NONE C 4 BY 4 MATRIX MULTIPLICATION: C=A*B. C ALL ARE DOUBLE PRECISION C IT IS SAFE TO CALL WITH C=A OR B. C FIRST SUBSCRIPT=ROWS, SECOND=COLUMNS DOUBLE PRECISION A(4,4),B(4,4),C(4,4),D(4,4) INTEGER I,J,K DO 30 I=1,4 DO 20 J=1,4 D(I,J)=0 DO 10 K=1,4 D(I,J)=D(I,J)+A(I,K)*B(K,J) 10 CONTINUE 20 CONTINUE 30 CONTINUE DO 50 I=1,4 DO 40 J=1,4 C(I,J)=D(I,J) 40 CONTINUE 50 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE KTINVT(A,B) IMPLICIT NONE C---INVERT TRANSFORMATION MATRIX A C C A = INPUT : 4 BY 4 TRANSFORMATION MATRIX C B = OUTPUT : INVERTED TRANSFORMATION MATRIX C C IF A IS NOT A TRANSFORMATION MATRIX YOU WILL GET STRANGE RESULTS C C NOTE THAT IT IS SAFE TO CALL WITH A=B C DOUBLE PRECISION A(4,4),B(4,4),C(4,4) INTEGER I,J C---TRANSPOSE DO 20 I=1,4 DO 10 J=1,4 C(I,J)=A(J,I) 10 CONTINUE 20 CONTINUE C---NEGATE ENERGY-MOMENTUM MIXING TERMS DO 30 I=1,3 C(4,I)=-C(4,I) C(I,4)=-C(I,4) 30 CONTINUE C---OUTPUT DO 50 I=1,4 DO 40 J=1,4 B(I,J)=C(I,J) 40 CONTINUE 50 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE KTWARN(SUBRTN,ICODE,*) C DEALS WITH ERRORS DURING EXECUTION C SUBRTN = NAME OF CALLING SUBROUTINE C ICODE = ERROR CODE: - 99 PRINT WARNING & CONTINUE C 100-199 PRINT WARNING & JUMP C 200- PRINT WARNING & STOP DEAD C----------------------------------------------------------------------- INTEGER ICODE CHARACTER*6 SUBRTN WRITE (6,10) SUBRTN,ICODE 10 FORMAT(/' KTWARN CALLED FROM SUBPROGRAM ',A6,': CODE =',I4/) IF (ICODE.LT.100) RETURN IF (ICODE.LT.200) RETURN 1 STOP END C----------------------------------------------------------------------- C----------------------------------------------------------------------- C-----------------------------------------------------------------------