C...This file collects software to study string recoupling in WW events. C...The main program is just one example how to study some consequences C...for the reconstructed W mass, while the subrouties have more lasting C...value. C...Copyright Torbjorn Sjostrand, CERN, Geneva, 1993. COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5) COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200) COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 COMMON/LUDATR/MRLU(6),RRLU(100) DIMENSION NEVT(15),IQQ(4),PMWR(4,4),IJJ(4),DISTSV(4,4),NACC(5,5), &PDIFF(5),NACC2(5,5) DOUBLE PRECISION WACC(5,5),WACC2(5,5) DATA NEVT/15*0/,NACC/25*0/,WACC/25*0D0/,NACC2/25*0/,WACC2/25*0D0/ C...Initialize. Number of events. MRLU(1)=19930926 ECM=170. NEV=40000 NRE=5 MDCY(111,1)=0 CALL PYWWIN(ECM) C***With current Jetset versions you need to set switch below to C***reproduce correctly the behaviour of some of the toy models. C***Corrected 1997-12-17. MSTJ(50)=0 C...Set up cluster finding for 4 clusters. MSTU(41)=1 MSTU(47)=4 PARU(44)=8. C...Histograms. DO 100 IRE=1,NRE CALL GBOOK1(IRE,'Prob(number of jets)',20,-0.5,19.5) CALL GBOOK1(IRE+5,'Prob(reconnect all)',5,-0.5,4.5) CALL GBOOK1(IRE+10,'Prob(reconnect 4-jet)',5,-0.5,4.5) CALL GBOOK1(IRE+15,'Prob(reconnect survive)',5,-0.5,4.5) CALL GBOOK1(IRE+20,'(true)',100,70.,90.) CALL GBOOK1(IRE+25,'(LUCLUS)',100,70.,90.) CALL GBOOK1(IRE+30,'(LUCLUS-true)',100,-10.,10.) CALL GBOOK1(IRE+35,'(LUCLUS min dist)',100,70.,90.) CALL GBOOK1(IRE+40,'(LUCLUS min dist-true)',100,-10.,10.) CALL GBOOK1(IRE+45,'(LUCLUS min dist2)',100,70.,90.) CALL GBOOK1(IRE+50,'(LUCLUS min dist2-true)',100,-10.,10.) CALL GBOOK1(IRE+55,'(LUCLUS min dist3)',100,70.,90.) CALL GBOOK1(IRE+60,'(LUCLUS min dist3-true)',100,-10.,10.) CALL GBOOK1(IRE+65,'(LUCLUS min dist4)',100,70.,90.) CALL GBOOK1(IRE+70,'(LUCLUS min dist4-true)',100,-10.,10.) 100 CONTINUE C...Loop through the cases of interest. DO 300 IRE=1,NRE MODE=0 PARA=0. IF(IRE.EQ.1) THEN IALT=5 MODE=1 ELSEIF(IRE.EQ.2) THEN IALT=4 MODE=2 PARA=0.6 ELSEIF(IRE.EQ.3) THEN IALT=0 ELSEIF(IRE.EQ.4) THEN IALT=2 MODE=1 ELSE IALT=2 MODE=2 ENDIF C...Generate events for each case. DO 290 IEV=1,NEV CALL PYWWEV(IALT,MODE,PARA) c IF(IEV.EQ.1) CALL LULIST(1) C...Do cluster search. IACC=0 CALL LUCLUS(NJET) CALL GFILL1(IRE,FLOAT(NJET),1.) CALL GFILL1(IRE+5,FLOAT(IRR),1.) IF(NJET.NE.4) GOTO 290 NEVT(IRE)=NEVT(IRE)+1 CALL GFILL1(IRE+10,FLOAT(IRR),1.) C...Check if all clusters well separated and have good energy. DSPMIN=0.5 EMIN=20. ISEP=1 DO 120 I1=N+1,N+3 DO 120 I2=I1+1,N+4 PA1=SQRT(P(I1,1)**2+P(I1,2)**2+P(I1,3)**2) PA2=SQRT(P(I2,1)**2+P(I2,2)**2+P(I2,3)**2) P12=P(I1,1)*P(I2,1)+P(I1,2)*P(I2,2)+P(I1,3)*P(I2,3) DSEP=ACOS(MAX(-1.,MIN(1.,P12/(PA1*PA2)))) DISTSV(I1-N,I2-N)=DSEP DISTSV(I2-N,I1-N)=DSEP 120 IF(DSEP.LT.DSPMIN) ISEP=0 DO 130 I1=N+1,N+4 130 IF(P(I1,4).LT.EMIN) ISEP=0 IF(ISEP.EQ.0) GOTO 290 NEVT(IRE+5)=NEVT(IRE+5)+1 CALL GFILL1(IRE+15,FLOAT(IRR),1.) C...Fill true average W mass. PMAVG=0.5*(P(IW1,5)+P(IW2,5)) CALL GFILL1(IRE+20,PMAVG,1.) C...Relate clusters to partons. IQQ(1)=IQ1 IQQ(2)=IQ2 IQQ(3)=IQ3 IQQ(4)=IQ4 DO 150 I1=N+1,N+4 DO 150 I2=1,4 PMWRS=(P(I1,4)+P(IQQ(I2),4))**2-(P(I1,1)+P(IQQ(I2),1))**2- &(P(I1,2)+P(IQQ(I2),2))**2-(P(I1,3)+P(IQQ(I2),3))**2 150 PMWR(I1-N,I2)=SQRT(MAX(0.,PMWRS)) PROMIN=1E10 DO 160 IJ1=1,4 DO 160 IJ2=1,4 DO 160 IJ3=1,4 DO 160 IJ4=1,4 IF(IJ1.EQ.IJ2.OR.IJ1.EQ.IJ3.OR.IJ1.EQ.IJ4.OR.IJ2.EQ.IJ3.OR. &IJ2.EQ.IJ4.OR.IJ3.EQ.IJ4) GOTO 160 PRONOW=PMWR(IJ1,1)*PMWR(IJ2,2)*PMWR(IJ3,3)*PMWR(IJ4,4) IF(PRONOW.LT.PROMIN) THEN PROMIN=PRONOW IJJ(1)=N+IJ1 IJJ(2)=N+IJ2 IJJ(3)=N+IJ3 IJJ(4)=N+IJ4 ENDIF 160 CONTINUE C...Fill reconstructed average W masses. PMWRS=(P(IJJ(1),4)+P(IJJ(2),4))**2-(P(IJJ(1),1)+P(IJJ(2),1))**2- &(P(IJJ(1),2)+P(IJJ(2),2))**2-(P(IJJ(1),3)+P(IJJ(2),3))**2 PMWW1=SQRT(MAX(0.,PMWRS)) PMWRS=(P(IJJ(3),4)+P(IJJ(4),4))**2-(P(IJJ(3),1)+P(IJJ(4),1))**2- &(P(IJJ(3),2)+P(IJJ(4),2))**2-(P(IJJ(3),3)+P(IJJ(4),3))**2 PMWW2=SQRT(MAX(0.,PMWRS)) PMWWA=0.5*(PMWW1+PMWW2) CALL GFILL1(IRE+25,PMWWA,1.) CALL GFILL1(IRE+30,PMWWA-PMAVG,1.) C...Pick jet combination which gives masses closest to nominal ones. PMBEST=0. PMDEVI=80. DMAXI=0. PMDEVJ=80. DO 180 ICA=1,3 IF(ICA.EQ.1) THEN I1=N+1 I2=N+2 I3=N+3 I4=N+4 ELSEIF(ICA.EQ.2) THEN I1=N+1 I2=N+3 I3=N+2 I4=N+4 ELSE I1=N+1 I2=N+4 I3=N+2 I4=N+3 ENDIF PMWRS1=(P(I1,4)+P(I2,4))**2-(P(I1,1)+P(I2,1))**2- &(P(I1,2)+P(I2,2))**2-(P(I1,3)+P(I2,3))**2 PMWR1=SQRT(MAX(0.,PMWRS1)) PMWRS2=(P(I3,4)+P(I4,4))**2-(P(I3,1)+P(I4,1))**2- &(P(I3,2)+P(I4,2))**2-(P(I3,3)+P(I4,3))**2 PMWR2=SQRT(MAX(0.,PMWRS2)) PMWRA=0.5*(PMWR1+PMWR2) IF(ABS(PMWRA-80.).LT.ABS(PMBEST-80.)) PMBEST=PMWRA PMDEWW=ABS(PMWR1-80.)+ABS(PMWR2-80.) IF(PMDEWW.LT.PMDEVI) THEN PMDEVI=PMDEWW PMBESA=PMWRA ENDIF DTHE=DISTSV(I1-N,I2-N)+DISTSV(I3-N,I4-N) IF(DTHE.GT.DMAXI) THEN DMAXI=DTHE PMBESB=PMWRA ENDIF PMDEW=ABS(MAX(PMWR1,PMWR2)-80.8) IF(PMDEW.LT.PMDEVJ) THEN PMDEVJ=PMDEW PMBESC=PMWRA ENDIF 180 CONTINUE C...Fill average masses for different best choices. CALL GFILL1(IRE+35,PMBEST,1.) CALL GFILL1(IRE+40,PMBEST-PMAVG,1.) CALL GFILL1(IRE+45,PMBESA,1.) CALL GFILL1(IRE+50,PMBESA-PMAVG,1.) CALL GFILL1(IRE+55,PMBESB,1.) CALL GFILL1(IRE+60,PMBESB-PMAVG,1.) CALL GFILL1(IRE+65,PMBESC,1.) CALL GFILL1(IRE+70,PMBESC-PMAVG,1.) C...Numerical statistics; not histograms. PDIFF(1)=PMWWA-PMAVG PDIFF(2)=PMBEST-PMAVG PDIFF(3)=PMBESA-PMAVG PDIFF(4)=PMBESB-PMAVG PDIFF(5)=PMBESC-PMAVG DO 200 J=1,5 IF(ABS(PDIFF(J)).LT.10.) THEN NACC(IRE,J)=NACC(IRE,J)+1 WACC(IRE,J)=WACC(IRE,J)+PDIFF(J) ENDIF IF(ABS(PDIFF(J)).LT.4.) THEN NACC2(IRE,J)=NACC2(IRE,J)+1 WACC2(IRE,J)=WACC2(IRE,J)+PDIFF(J) ENDIF 200 CONTINUE C...End event and case loop. 290 CONTINUE 300 CONTINUE C...Numerical statistics. DO 320 IRE=1,5 DO 320 J=1,5 WACC(IRE,J)=WACC(IRE,J)/NACC(IRE,J) WACC2(IRE,J)=WACC2(IRE,J)/NACC2(IRE,J) 320 CONTINUE DO 330 J=1,5 WRITE(6,'(I5,5I6,5F10.3)') J,(NACC(IRE,J),IRE=1,5), &(WACC(IRE,J),IRE=1,5) WRITE(6,'(I5,5I6,5F10.3)') J,(NACC2(IRE,J),IRE=1,5), &(WACC2(IRE,J),IRE=1,5) 330 CONTINUE C...Produce output. C CALL PYSTAT(1) WRITE(6,*) 'NUMBER OF EVENTS:',NEV,(NEVT(IRE+5),IRE=1,NRE) DO 400 IRE=1,NRE FAC=1./NEV CALL GSCALE(IRE,FAC) CALL GSCALE(IRE+5,FAC) FAC=1./MAX(1,NEVT(IRE)) CALL GSCALE(IRE+10,FAC) FAC=1./MAX(1,NEVT(IRE+5)) CALL GSCALE(IRE+15,FAC) CALL GSCALE(IRE+20,5.*FAC) CALL GSCALE(IRE+25,5.*FAC) CALL GSCALE(IRE+30,5.*FAC) CALL GSCALE(IRE+35,5.*FAC) CALL GSCALE(IRE+40,5.*FAC) CALL GSCALE(IRE+45,5.*FAC) CALL GSCALE(IRE+50,5.*FAC) CALL GSCALE(IRE+55,5.*FAC) CALL GSCALE(IRE+60,5.*FAC) CALL GSCALE(IRE+65,5.*FAC) CALL GSCALE(IRE+70,5.*FAC) 400 CONTINUE c CALL GDUMP('gout.85') CALL GCLEAR END C********************************************************************* SUBROUTINE PYWWIN(ECM) C...Subroutine to set up e+ e- -> W+ W- -> q1 q2bar q3 q4bar C...in mode suitable for event generation with string recoupling. COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5) COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200) COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) COMMON/PYWWCM/ECMWW,PMW,PGW,HBAR,TFRAG,RHAD C...Process. MSEL=0 MSUB(25)=1 C...Set W to decay hadronically and not to top C...(modified 19/6 1995). DO 100 IDCY=MDCY(24,2),MDCY(24,2)+MDCY(24,3)-1 KFD1=IABS(KFDP(IDCY,1)) KFD2=IABS(KFDP(IDCY,2)) IF(KFD1.EQ.6.OR.KFD2.EQ.6) MDME(IDCY,1)=-1 IF(KFD1.GE.11.AND.KFD1.LE.18) MDME(IDCY,1)=MIN(0,MDME(IDCY,1)) 100 CONTINUE C...Switch off components to be done after string reconnection. MSTP(111)=0 MSTP(128)=1 C...Initialize. CALL PYINIT('CMS','E+','E-',ECM) C...Save ECM; set up some default values (modified 19/6 1995). ECMWW=ECM PMW=PMAS(24,1) PGW=PMAS(24,2) HBAR=0.2 TFRAG=1.5 RHAD=0.5 END C********************************************************************* SUBROUTINE PYWWEV(IALT,MODE,PARA) C...Main driver routine to generate events with string recoupling. C IALT = 0 : standard event generation without any ado. C MODE, PARA : dummy. C = 1 : two strings at rest crossing at an angle; u+d quarks. C MODE = 0 : no reconnection. C = 1 : reconnect after shower but before fragmentation. C = 2 : reconnect before shower and fragmentation. C PARA : angle in radians between q1 and qbar4. C = 2 : reconnect at origin of event. C MODE = 0 : no reconnection. C = 1 : reconnect after shower but before fragmentation. C = 2 : reconnect before shower and fragmentation. C PARA : dummy. C = 3 : reconnect at origin of event based on spherical overlap. C MODE = 0 : no reconnection. C = 1 : reconnect, instantaneous sphere. C = 2 : reconnect, time-retarded sphere. C PARA : strength of reconnection probability. C = 4 : reconnect where strings overlap based on cylinder geometry. C MODE = 0 : no reconnection. C = 1 : reconnect, instantaneous string pieces. C = 2 : reconnect, time-retarded string pieces. C PARA : strength of reconnection probability. C = 5 : reconnect when strings cross, a la type II superconductor. C MODE = 0 : no recoupling. C = 1 : recoupling at first crossing. C = 2 : recoupling at first crossing which decreases C total string length. C PARA : dummy. COMMON/PYWWCM/ECMWW,PMW,PGW,HBAR,TFRAG,RHAD COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 C...Simple alternative: never recoupling. IRR=0 DCUT1=0. DCUT2=0. IF(IALT.EQ.0) THEN CALL PYWWA0 C...Two strings at rest - toy model. ELSEIF(IALT.EQ.1) THEN CALL PYWWA1(MODE,PARA,0.5*ECMWW) C...Reconnect two strings at origin of event. ELSEIF(IALT.EQ.2) THEN CALL PYWWA2(MODE) C...Reconnect strings after shower according to spherical overlap. ELSEIF(IALT.EQ.3) THEN CALL PYWWA3(MODE,PARA) C...Reconnect strings after shower according to cylindrical overlap. ELSEIF(IALT.EQ.4) THEN CALL PYWWA4(MODE,PARA) C...Type II superconductor - reconnect when cross. ELSEIF(IALT.EQ.5) THEN CALL PYWWA5(MODE) ENDIF RETURN END C********************************************************************* SUBROUTINE PYWWA0 C...Subroutine to generate simpleminded WW events without reconnection. COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 C...Generate event. CALL PYEVNT C...Find position of W's and quarks. IW1=0 IW2=0 IQ1=0 IQ2=0 IQ3=0 IQ4=0 DO 100 I=7,12 IF(K(I,2).EQ.24) IW1=I IF(K(I,2).EQ.-24) IW2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.24) IQ1=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.24) IQ2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.-24) IQ3=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.-24) IQ4=I 100 CONTINUE C...Do fragmentation. CALL LUEXEC CALL PYLAMB(RLEN) ALBEF=RLEN ALAFT=RLEN RETURN END C********************************************************************* SUBROUTINE PYWWA1(MODE,THETA,EW) C...Subroutine to generate simpleminded WW event and reconnect strings. C...Assumed to be at threshold, i.e. mW = Ecm/2, and into u+dbar+d+ubar. C...MODE = 0 : no string reconnection. C... = 1 : string reconnection after shower but before fragmentation. C... = 2 : string reconnection before shower and fragmentation. COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 DIMENSION IJOIN(2) C...Set up two W's for documentation. IW1=1 IW2=2 MSTU(10)=1 P(IW1,5)=EW CALL LU1ENT(IW1,24,EW,0.,0.) MSTU(10)=1 P(IW2,5)=EW CALL LU1ENT(IW2,-24,EW,0.,0.) MSTU(10)=2 K(IW1,1)=21 K(IW2,1)=21 C...Set up two jet systems. IQ1=3 IQ2=4 IQ3=5 IQ4=6 CALL LU2ENT(IQ1,2,-1,EW) CALL LU2ENT(IQ3,1,-2,EW) CALL LUDBRB(IQ1,IQ2,0.5*THETA,0.,0D0,0D0,0D0) CALL LUDBRB(IQ3,IQ4,3.1416-0.5*THETA,0.,0D0,0D0,0D0) K(IQ1,3)=IW1 K(IQ2,3)=IW1 K(IQ3,3)=IW2 K(IQ4,3)=IW2 C...(Re)connect strings. IF(MODE.EQ.0) THEN IJOIN(1)=IQ1 IJOIN(2)=IQ2 CALL LUJOIN(2,IJOIN) IJOIN(1)=IQ3 IJOIN(2)=IQ4 CALL LUJOIN(2,IJOIN) ELSE IJOIN(1)=IQ1 IJOIN(2)=IQ4 CALL LUJOIN(2,IJOIN) IJOIN(1)=IQ2 IJOIN(2)=IQ3 CALL LUJOIN(2,IJOIN) ENDIF C...Do parton shower. IF(MODE.LE.1) THEN CALL LUSHOW(IQ1,IQ2,EW) CALL LUSHOW(IQ3,IQ4,EW) ELSE PPM2=(P(IQ1,4)+P(IQ4,4))**2-(P(IQ1,1)+P(IQ4,1))**2- & (P(IQ1,2)+P(IQ4,2))**2-(P(IQ1,3)+P(IQ4,3))**2 PPM=SQRT(MAX(0.,PPM2)) CALL LUSHOW(IQ1,IQ4,PPM) PPM2=(P(IQ3,4)+P(IQ2,4))**2-(P(IQ3,1)+P(IQ2,1))**2- & (P(IQ3,2)+P(IQ2,2))**2-(P(IQ3,3)+P(IQ2,3))**2 PPM=SQRT(MAX(0.,PPM2)) CALL LUSHOW(IQ3,IQ2,PPM) ENDIF C...Do fragmentation. CALL LUEXEC IF(MODE.GE.1) IRR=1 CALL PYLAMB(RLEN) ALBEF=RLEN ALAFT=RLEN RETURN END C********************************************************************* SUBROUTINE PYWWA2(MODE) C...Subroutine to generate WW event and reconnect strings. C...MODE = 0 : no string reconnection. C... = 1 : string reconnection after shower but before fragmentation. C... = 2 : string reconnection before shower and fragmentation. COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 DIMENSION IJOIN(2) C...Generate event. MSTP71=MSTP(71) MSTP(71)=0 CALL PYEVNT C...Find position of W's and quarks. IW1=0 IW2=0 IQ1=0 IQ2=0 IQ3=0 IQ4=0 DO 100 I=9,N IF(K(I,2).EQ.24) IW1=I IF(K(I,2).EQ.-24) IW2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.24) IQ1=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.24) IQ2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.-24) IQ3=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.-24) IQ4=I 100 CONTINUE C...(Re)connect strings. IF(MODE.EQ.0) THEN IJOIN(1)=IQ1 IJOIN(2)=IQ2 CALL LUJOIN(2,IJOIN) IJOIN(1)=IQ3 IJOIN(2)=IQ4 CALL LUJOIN(2,IJOIN) ELSE IJOIN(1)=IQ1 IJOIN(2)=IQ4 CALL LUJOIN(2,IJOIN) IJOIN(1)=IQ2 IJOIN(2)=IQ3 CALL LUJOIN(2,IJOIN) ENDIF C...Do parton shower. IF(MODE.LE.1) THEN NPS1=N CALL LUSHOW(IQ1,IQ2,P(IW1,5)) NPS2=N CALL LUSHOW(IQ3,IQ4,P(IW2,5)) C...Amount of energy to either side of recoupling. FRACP=(P(NPS1+1,5)**2+P(NPS1+2,5)**2-P(NPS1+3,5)**2)/ & (2.*P(NPS1+1,5)**2) DCUT1=2.*MAX(FRACP,1.-FRACP)-1. FRACM=(P(NPS2+1,5)**2+P(NPS2+2,5)**2-P(NPS2+3,5)**2)/ & (2.*P(NPS2+1,5)**2) DCUT2=2.*MAX(FRACM,1.-FRACM)-1. ELSE PPM2=(P(IQ1,4)+P(IQ4,4))**2-(P(IQ1,1)+P(IQ4,1))**2- & (P(IQ1,2)+P(IQ4,2))**2-(P(IQ1,3)+P(IQ4,3))**2 PPM=SQRT(MAX(0.,PPM2)) CALL LUSHOW(IQ1,IQ4,PPM) PPM2=(P(IQ3,4)+P(IQ2,4))**2-(P(IQ3,1)+P(IQ2,1))**2- & (P(IQ3,2)+P(IQ2,2))**2-(P(IQ3,3)+P(IQ2,3))**2 PPM=SQRT(MAX(0.,PPM2)) CALL LUSHOW(IQ3,IQ2,PPM) ENDIF C...Do fragmentation. CALL LUEXEC MSTP(71)=MSTP71 IF(MODE.GE.1) IRR=1 CALL PYLAMB(RLEN) ALBEF=RLEN ALAFT=RLEN RETURN END C********************************************************************* SUBROUTINE PYWWA3(MODE,PARA) C...Subroutine to reconnect strings at origin of event, after shower, C...based on spherical overlap. C...Input: MODE = 0 : no reconnection. C = 1 : reconnect, instantaneous sphere. C = 2 : reconnect, time-retarded sphere. C PARA : strength of reconnection probability. PARAMETER (PI=3.141592) PARAMETER (NPT=100) COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) COMMON/PYWWCM/ECMWW,PMW,PGW,HBAR,TFRAG,RHAD COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 DIMENSION IJOIN(2) C...Generate event. MREC=0 MSTP71=MSTP(71) MSTP(71)=0 CALL PYEVNT C...Find position of W's and quarks. IW1=0 IW2=0 IQ1=0 IQ2=0 IQ3=0 IQ4=0 DO 100 I=9,N IF(K(I,2).EQ.24) IW1=I IF(K(I,2).EQ.-24) IW2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.24) IQ1=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.24) IQ2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.-24) IQ3=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.-24) IQ4=I 100 CONTINUE IF(MODE.EQ.0) GOTO 200 C...Kinematics: boost factors etc. PM1=P(IW1,5) PM2=P(IW2,5) PE1=P(IW1,4) PE2=P(IW2,4) PCOM=SQRT(P(IW1,1)**2+P(IW1,2)**2+P(IW1,3)**2) BE1=PCOM/PE1 GA1=PE1/PM1 BE2=PCOM/PE2 GA2=PE2/PM2 C...Select W decay times. TAU1=HBAR*(-LOG(RLU(0)))*PM1/ &SQRT((PM1**2-PMW**2)**2+(PM1**2*PGW/PMW)**2) TAU2=HBAR*(-LOG(RLU(0)))*PM2/ &SQRT((PM2**2-PMW**2)**2+(PM2**2*PGW/PMW)**2) GTMAX=MAX(GA1*TAU1,GA2*TAU2) C...Loop over number of space-time points. SUM=0. DO 110 IPT=1,NPT C...Pick x,y,z,t Gaussian (width RHAD and TFRAG, respectively). R=SQRT(-LOG(RLU(0))) PHI=2.*PI*RLU(0) X=RHAD*R*COS(PHI) Y=RHAD*R*SIN(PHI) R=SQRT(-LOG(RLU(0))) PHI=2.*PI*RLU(0) Z=RHAD*R*COS(PHI) T=GTMAX+SQRT(0.5)*TFRAG*R*ABS(SIN(PHI)) C...Weight for sample distribution and W+ and W- distributions. WTSMP=EXP(-(X**2+Y**2+Z**2)/RHAD**2)* &EXP(-2.*(T-GTMAX)**2/TFRAG**2) R2=X**2+Y**2+GA1**2*(Z-BE1*T)**2 TR=GA1*(T-BE1*Z)-TAU1 WT1=EXP(-R2/(2.*RHAD**2))*EXP(-TR**2/TFRAG**2) IF(MODE.EQ.2.AND.TR-SQRT(R2).LT.0.) WT1=0. R2=X**2+Y**2+GA2**2*(Z+BE2*T)**2 TR=GA2*(T+BE2*Z)-TAU2 WT2=EXP(-R2/(2.*RHAD**2))*EXP(-TR**2/TFRAG**2) IF(MODE.EQ.2.AND.TR-SQRT(R2).LT.0.) WT2=0. C...Result of integration. WT=WT1*WT2/WTSMP SUM=SUM+WT 110 CONTINUE RES=SUM/NPT C...Decide whether to reconnect. PREC=0.5*(1.-EXP(-PARA*RES)) IF(PREC.GT.RLU(0)) MREC=1 C...(Re)connect strings. 200 IF(MREC.EQ.0) THEN IJOIN(1)=IQ1 IJOIN(2)=IQ2 CALL LUJOIN(2,IJOIN) IJOIN(1)=IQ3 IJOIN(2)=IQ4 CALL LUJOIN(2,IJOIN) ELSE IJOIN(1)=IQ1 IJOIN(2)=IQ4 CALL LUJOIN(2,IJOIN) IJOIN(1)=IQ2 IJOIN(2)=IQ3 CALL LUJOIN(2,IJOIN) ENDIF C...Do parton shower. NPS1=N CALL LUSHOW(IQ1,IQ2,P(IW1,5)) NPS2=N CALL LUSHOW(IQ3,IQ4,P(IW2,5)) C...Do fragmentation. CALL LUEXEC MSTP(71)=MSTP71 IF(MREC.EQ.1) IRR=1 CALL PYLAMB(RLEN) ALBEF=RLEN ALAFT=RLEN C...Amount of energy to either side of recoupling. FRACP=(P(NPS1+1,5)**2+P(NPS1+2,5)**2-P(NPS1+3,5)**2)/ &(2.*P(NPS1+1,5)**2) DCUT1=2.*MAX(FRACP,1.-FRACP)-1. FRACM=(P(NPS2+1,5)**2+P(NPS2+2,5)**2-P(NPS2+3,5)**2)/ &(2.*P(NPS2+1,5)**2) DCUT2=2.*MAX(FRACM,1.-FRACM)-1. RETURN END C********************************************************************* SUBROUTINE PYWWA4(MODE,PARA) C...Subroutine to generate WW event and reconnect strings C...according to overlap of cylindrical string piece volumes. C...Input: MODE = 0 : no reconnection. C = 1 : reconnect, instantaneous cylinders. C = 2 : reconnect, time-retarded cylinders. C PARA : strength of reconnection probability. PARAMETER (PI=3.141592) PARAMETER (NPT=100) PARAMETER (BLOWR=2.5) PARAMETER (BLOWT=2.) COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/PYWWCM/ECMWW,PMW,PGW,HBAR,TFRAG,RHAD COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 DIMENSION XP(3),XM(3),INP(50),INM(50),IJOIN(100), &V1(3),V2(3),BETP(50,4),DIRP(50,3),BETM(50,4),DIRM(50,3), &XD(4),XB(4),PSUMP(4),PSUMM(4) DIMENSION IAP(NPT),IAM(NPT),WTA(NPT) C...Generate event. IACC=0 CALL PYEVNT C...Find position of W's and quarks. IW1=0 IW2=0 IQ1=0 IQ2=0 IQ3=0 IQ4=0 DO 100 I=7,14 IF(K(I,2).EQ.24) IW1=I IF(K(I,2).EQ.-24) IW2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.24) IQ1=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.24) IQ2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.-24) IQ3=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.-24) IQ4=I 100 CONTINUE IF(MODE.EQ.0) GOTO 320 C...Select decay vertices of W+ and W-. TP=HBAR*(-LOG(RLU(0)))*P(IW1,4)/ &SQRT((P(IW1,5)**2-PMW**2)**2+(P(IW1,5)**2*PGW/PMW)**2) TM=HBAR*(-LOG(RLU(0)))*P(IW2,4)/ &SQRT((P(IW2,5)**2-PMW**2)**2+(P(IW2,5)**2*PGW/PMW)**2) DO 110 J=1,3 XP(J)=TP*P(IW1,J)/P(IW1,4) XM(J)=TM*P(IW2,J)/P(IW2,4) 110 CONTINUE GTMAX=MAX(TP,TM) C...Find partons pointing back to W+ and W-; store them with quark C...end of string first. NNP=0 NNM=0 ISGP=0 ISGM=0 DO 140 I=9,N IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 140 C...Do not count photons (added 19/6 1995). IF(K(I,2).EQ.22) GOTO 140 IF(K(I,3).EQ.IW1) THEN IF(ISGP.EQ.0) ISGP=ISIGN(1,K(I,2)) NNP=NNP+1 IF(ISGP.EQ.1) THEN INP(NNP)=I ELSE DO 120 I1=NNP,2,-1 120 INP(I1)=INP(I1-1) INP(1)=I ENDIF IF(K(I,1).EQ.1) ISGP=0 ELSEIF(K(I,3).EQ.IW2) THEN IF(ISGM.EQ.0) ISGM=ISIGN(1,K(I,2)) NNM=NNM+1 IF(ISGM.EQ.1) THEN INM(NNM)=I ELSE DO 130 I1=NNM,2,-1 130 INM(I1)=INM(I1-1) INM(1)=I ENDIF IF(K(I,1).EQ.1) ISGM=0 ENDIF 140 CONTINUE C...Reconstruct velocity and direction of W+ string pieces. DO 180 IIP=1,NNP-1 IF(K(INP(IIP),2).LT.0) GOTO 180 I1=INP(IIP) I2=INP(IIP+1) P1A=SQRT(P(I1,1)**2+P(I1,2)**2+P(I1,3)**2) P2A=SQRT(P(I2,1)**2+P(I2,2)**2+P(I2,3)**2) DO 150 J=1,3 V1(J)=P(I1,J)/P1A 150 V2(J)=P(I2,J)/P2A DO 160 J=1,3 BETP(IIP,J)=0.5*(V1(J)+V2(J)) 160 DIRP(IIP,J)=V1(J)-V2(J) BETP(IIP,4)=1./SQRT(1.-BETP(IIP,1)**2-BETP(IIP,2)**2- &BETP(IIP,3)**2) DIRL=SQRT(DIRP(IIP,1)**2+DIRP(IIP,2)**2+DIRP(IIP,3)**2) DO 170 J=1,3 170 DIRP(IIP,J)=DIRP(IIP,J)/DIRL 180 CONTINUE C...Reconstruct velocity and direction of W- string pieces. DO 230 IIM=1,NNM-1 IF(K(INM(IIM),2).LT.0) GOTO 230 I1=INM(IIM) I2=INM(IIM+1) P1A=SQRT(P(I1,1)**2+P(I1,2)**2+P(I1,3)**2) P2A=SQRT(P(I2,1)**2+P(I2,2)**2+P(I2,3)**2) DO 190 J=1,3 V1(J)=P(I1,J)/P1A 190 V2(J)=P(I2,J)/P2A DO 200 J=1,3 BETM(IIM,J)=0.5*(V1(J)+V2(J)) 200 DIRM(IIM,J)=V1(J)-V2(J) BETM(IIM,4)=1./SQRT(1.-BETM(IIM,1)**2-BETM(IIM,2)**2- &BETM(IIM,3)**2) DIRL=SQRT(DIRM(IIM,1)**2+DIRM(IIM,2)**2+DIRM(IIM,3)**2) DO 220 J=1,3 220 DIRM(IIM,J)=DIRM(IIM,J)/DIRL 230 CONTINUE C...Loop over number of space-time points. NACC=0 SUM=0. DO 280 IPT=1,NPT C...Pick x,y,z,t Gaussian (width RHAD and TFRAG, respectively). R=SQRT(-LOG(RLU(0))) PHI=2.*PI*RLU(0) X=BLOWR*RHAD*R*COS(PHI) Y=BLOWR*RHAD*R*SIN(PHI) R=SQRT(-LOG(RLU(0))) PHI=2.*PI*RLU(0) Z=BLOWR*RHAD*R*COS(PHI) T=GTMAX+BLOWT*SQRT(0.5)*TFRAG*R*ABS(SIN(PHI)) C...Weight for sample distribution. WTSMP=EXP(-(X**2+Y**2+Z**2)/(BLOWR*RHAD)**2)* &EXP(-2.*(T-GTMAX)**2/(BLOWT*TFRAG)**2) C...Loop over W+ string pieces and find one with largest weight. IMAXP=0 WTMAXP=1E-10 XD(1)=X-XP(1) XD(2)=Y-XP(2) XD(3)=Z-XP(3) XD(4)=T-TP DO 250 IIP=1,NNP-1 IF(K(INP(IIP),2).LT.0) GOTO 250 BED=BETP(IIP,1)*XD(1)+BETP(IIP,2)*XD(2)+BETP(IIP,3)*XD(3) BEDG=BETP(IIP,4)*(BETP(IIP,4)*BED/(1.+BETP(IIP,4))-XD(4)) DO 240 J=1,3 240 XB(J)=XD(J)+BEDG*BETP(IIP,J) XB(4)=BETP(IIP,4)*(XD(4)-BED) SR2=XB(1)**2+XB(2)**2+XB(3)**2 SZ2=(DIRP(IIP,1)*XB(1)+DIRP(IIP,2)*XB(2)+DIRP(IIP,3)*XB(3))**2 WTP=EXP(-(SR2-SZ2)/(2.*RHAD**2))*EXP(-(XB(4)**2-SZ2)/TFRAG**2) IF(MODE.EQ.1.AND.XB(4)-SQRT(SZ2).LT.0.) WTP=0. IF(MODE.EQ.2.AND.XB(4)-SQRT(SR2).LT.0.) WTP=0. IF(WTP.GT.WTMAXP) THEN IMAXP=IIP WTMAXP=WTP ENDIF 250 CONTINUE C...Loop over W- string pieces and find one with largest weight. IMAXM=0 WTMAXM=1E-10 XD(1)=X-XM(1) XD(2)=Y-XM(2) XD(3)=Z-XM(3) XD(4)=T-TM DO 270 IIM=1,NNM-1 IF(K(INM(IIM),2).LT.0) GOTO 270 BED=BETM(IIM,1)*XD(1)+BETM(IIM,2)*XD(2)+BETM(IIM,3)*XD(3) BEDG=BETM(IIM,4)*(BETM(IIM,4)*BED/(1.+BETM(IIM,4))-XD(4)) DO 260 J=1,3 260 XB(J)=XD(J)+BEDG*BETM(IIM,J) XB(4)=BETM(IIM,4)*(XD(4)-BED) SR2=XB(1)**2+XB(2)**2+XB(3)**2 SZ2=(DIRM(IIM,1)*XB(1)+DIRM(IIM,2)*XB(2)+DIRM(IIM,3)*XB(3))**2 WTM=EXP(-(SR2-SZ2)/(2.*RHAD**2))*EXP(-(XB(4)**2-SZ2)/TFRAG**2) IF(MODE.EQ.1.AND.XB(4)-SQRT(SZ2).LT.0.) WTM=0. IF(MODE.EQ.2.AND.XB(4)-SQRT(SR2).LT.0.) WTM=0. IF(WTM.GT.WTMAXM) THEN IMAXM=IIM WTMAXM=WTM ENDIF 270 CONTINUE C...Result of integration. WT=0. IF(IMAXP.NE.0.AND.IMAXM.NE.0) THEN WT=WTMAXP*WTMAXM/WTSMP SUM=SUM+WT NACC=NACC+1 IAP(NACC)=IMAXP IAM(NACC)=IMAXM WTA(NACC)=WT ENDIF 280 CONTINUE RES=BLOWR**3*BLOWT*SUM/NPT C...Decide whether to reconnect. PREC=1.-EXP(-PARA*RES) MREC=0 IF(PREC.GT.RLU(0)) MREC=1 IF(MREC.EQ.0) GOTO 320 C...Decide which pair of strings to reconnect. RSUM=RLU(0)*SUM DO 290 IA=1,NACC IACC=IA RSUM=RSUM-WTA(IA) IF(RSUM.LE.0.) GOTO 300 290 CONTINUE 300 IIP=IAP(IACC) IIM=IAM(IACC) CALL PYLAMB(RLEN) ALBEF=RLEN C...Recouple strings. NJOIN=0 DO 310 IS=1,NNP+NNM NJOIN=NJOIN+1 IF(IS.LE.IIP) THEN I=INP(IS) ELSEIF(IS.LE.IIP+NNM-IIM) THEN I=INM(IS-IIP+IIM) ELSEIF(IS.LE.IIP+NNM) THEN I=INM(IS-IIP-NNM+IIM) ELSE I=INP(IS-NNM) ENDIF IJOIN(NJOIN)=I IF(K(I,2).LT.0) THEN CALL LUJOIN(NJOIN,IJOIN) NJOIN=0 ENDIF 310 CONTINUE C...Amount of energy to either side of recoupling. DO 330 J=1,4 PSUMP(J)=0. 330 PSUMM(J)=0. DO 340 IS=1,IIP I=INP(IS) DO 340 J=1,4 340 PSUMP(J)=PSUMP(J)+P(I,J) DO 350 IS=1,IIM I=INM(IS) DO 350 J=1,4 350 PSUMM(J)=PSUMM(J)+P(I,J) FRACP=(P(IW1,4)*PSUMP(4)-P(IW1,1)*PSUMP(1)-P(IW1,2)*PSUMP(2)- &P(IW1,3)*PSUMP(3))/P(IW1,5)**2 DCUT1=2.*MAX(FRACP,1.-FRACP)-1. FRACM=(P(IW2,4)*PSUMM(4)-P(IW2,1)*PSUMM(1)-P(IW2,2)*PSUMM(2)- &P(IW2,3)*PSUMM(3))/P(IW2,5)**2 DCUT2=2.*MAX(FRACM,1.-FRACM)-1. C...Do fragmentation. 320 CONTINUE CALL LUEXEC IF(IACC.NE.0) IRR=1 CALL PYLAMB(RLEN) ALAFT=RLEN RETURN END C********************************************************************* SUBROUTINE PYWWA5(MODE) C...Subroutine to generate WW event and reconnect strings C...when they cross (type II superconductor). C...Input: MODE = 0 : no string reconnection. C = 1 : string reconnection. C = 2 : string reconnection if string length decreased. COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/PYWWCM/ECMWW,PMW,PGW,HBAR,TFRAG,RHAD COMMON/PYWWII/IW1,IW2,IQ1,IQ2,IQ3,IQ4,IRR,ALBEF,ALAFT,DCUT1,DCUT2 DIMENSION XP(3),XM(3),V1P(3),V2P(3),V1M(3),V2M(3) DIMENSION INP(50),INM(50),IPC(20),IMC(20),TC(0:20),TPC(20), &TMC(20),IJOIN(100),PSUMP(4),PSUMM(4) C...Function to give four product. FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3) C...Generate event. IACC=0 CALL PYEVNT C...Find position of W's and quarks. IW1=0 IW2=0 IQ1=0 IQ2=0 IQ3=0 IQ4=0 DO 100 I=7,14 IF(K(I,2).EQ.24) IW1=I IF(K(I,2).EQ.-24) IW2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.24) IQ1=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.24) IQ2=I IF(K(I,2).GT.0.AND.K(I,2).LT.10.AND.K(K(I,3),2).EQ.-24) IQ3=I IF(K(I,2).GT.-10.AND.K(I,2).LT.0.AND.K(K(I,3),2).EQ.-24) IQ4=I 100 CONTINUE IF(MODE.EQ.0) GOTO 250 C...Select decay vertices of W+ and W-. TP=HBAR*(-LOG(RLU(0)))*P(IW1,4)/ &SQRT((P(IW1,5)**2-PMW**2)**2+(P(IW1,5)**2*PGW/PMW)**2) TM=HBAR*(-LOG(RLU(0)))*P(IW2,4)/ &SQRT((P(IW2,5)**2-PMW**2)**2+(P(IW2,5)**2*PGW/PMW)**2) DO 110 J=1,3 XP(J)=TP*P(IW1,J)/P(IW1,4) XM(J)=TM*P(IW2,J)/P(IW2,4) 110 CONTINUE C...Find partons pointing back to W+ and W-; store them with quark C...end of string first. NNP=0 NNM=0 ISGP=0 ISGM=0 DO 140 I=9,N IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 140 C...Do not count photons (added 19/6 1995). IF(K(I,2).EQ.22) GOTO 140 IF(K(I,3).EQ.IW1) THEN IF(ISGP.EQ.0) ISGP=ISIGN(1,K(I,2)) NNP=NNP+1 IF(ISGP.EQ.1) THEN INP(NNP)=I ELSE DO 120 I1=NNP,2,-1 120 INP(I1)=INP(I1-1) INP(1)=I ENDIF IF(K(I,1).EQ.1) ISGP=0 ELSEIF(K(I,3).EQ.IW2) THEN IF(ISGM.EQ.0) ISGM=ISIGN(1,K(I,2)) NNM=NNM+1 IF(ISGM.EQ.1) THEN INM(NNM)=I ELSE DO 130 I1=NNM,2,-1 130 INM(I1)=INM(I1-1) INM(1)=I ENDIF IF(K(I,1).EQ.1) ISGM=0 ENDIF 140 CONTINUE C...Loop through all string pieces, one from W+ and one from W-. NPAIR=0 NCROSS=0 TC(0)=0. DO 200 IIP=1,NNP-1 IF(K(INP(IIP),2).LT.0) GOTO 200 I1P=INP(IIP) I2P=INP(IIP+1) DO 190 IIM=1,NNM-1 IF(K(INM(IIM),2).LT.0) GOTO 190 I1M=INM(IIM) I2M=INM(IIM+1) NPAIR=NPAIR+1 C...Find endpoint velocity vectors and determine crossing. DO 150 J=1,3 V1P(J)=P(I1P,J)/P(I1P,4) V2P(J)=P(I2P,J)/P(I2P,4) V1M(J)=P(I1M,J)/P(I1M,4) V2M(J)=P(I2M,J)/P(I2M,4) 150 CONTINUE CALL PYWWXX(TP,TM,XP,XM,V1P,V2P,V1M,V2M,IANSW,T,TAUP,TAUM) C...Order crossings by time. End loop over crossings. IF(IANSW.EQ.1.AND.NCROSS.LT.20) THEN NCROSS=NCROSS+1 DO 160 I1=NCROSS,1,-1 IF(T.GT.TC(I1-1).OR.I1.EQ.1) THEN IPC(I1)=IIP IMC(I1)=IIM TC(I1)=T TPC(I1)=TAUP TMC(I1)=TAUM GOTO 170 ELSE IPC(I1)=IPC(I1-1) IMC(I1)=IMC(I1-1) TC(I1)=TC(I1-1) TPC(I1)=TPC(I1-1) TMC(I1)=TMC(I1-1) ENDIF 160 CONTINUE 170 CONTINUE ENDIF 190 CONTINUE 200 CONTINUE C...Loop over crossings; find first (if any) acceptable one. IF(NCROSS.GE.1) THEN DO 210 IC=1,NCROSS PNFRAG=EXP(-(TPC(IC)**2+TMC(IC)**2)/TFRAG**2) IF(PNFRAG.GT.RLU(0)) THEN IF(MODE.EQ.1) THEN IACC=IC GOTO 220 ELSE IIP=IPC(IC) IIM=IMC(IC) I1P=INP(IIP) I2P=INP(IIP+1) I1M=INM(IIM) I2M=INM(IIM+1) ELOLD=FOUR(I1P,I2P)*FOUR(I1M,I2M) ELNEW=FOUR(I1P,I2M)*FOUR(I1M,I2P) IF(ELNEW.LT.ELOLD) THEN IACC=IC GOTO 220 ENDIF ENDIF ENDIF 210 CONTINUE 220 CONTINUE ENDIF CALL PYLAMB(RLEN) ALBEF=RLEN C...Recouple strings. IF(IACC.NE.0) THEN IIP=IPC(IACC) IIM=IMC(IACC) NJOIN=0 DO 230 IS=1,NNP+NNM NJOIN=NJOIN+1 IF(IS.LE.IIP) THEN I=INP(IS) ELSEIF(IS.LE.IIP+NNM-IIM) THEN I=INM(IS-IIP+IIM) ELSEIF(IS.LE.IIP+NNM) THEN I=INM(IS-IIP-NNM+IIM) ELSE I=INP(IS-NNM) ENDIF IJOIN(NJOIN)=I IF(K(I,2).LT.0) THEN CALL LUJOIN(NJOIN,IJOIN) NJOIN=0 ENDIF 230 CONTINUE ENDIF C...Amount of energy to either side of recoupling. DO 330 J=1,4 PSUMP(J)=0. 330 PSUMM(J)=0. DO 340 IS=1,IIP I=INP(IS) DO 340 J=1,4 340 PSUMP(J)=PSUMP(J)+P(I,J) DO 350 IS=1,IIM I=INM(IS) DO 350 J=1,4 350 PSUMM(J)=PSUMM(J)+P(I,J) FRACP=(P(IW1,4)*PSUMP(4)-P(IW1,1)*PSUMP(1)-P(IW1,2)*PSUMP(2)- &P(IW1,3)*PSUMP(3))/P(IW1,5)**2 DCUT1=2.*MAX(FRACP,1.-FRACP)-1. FRACM=(P(IW2,4)*PSUMM(4)-P(IW2,1)*PSUMM(1)-P(IW2,2)*PSUMM(2)- &P(IW2,3)*PSUMM(3))/P(IW2,5)**2 DCUT2=2.*MAX(FRACM,1.-FRACM)-1. C...Do fragmentation. 250 CALL LUEXEC IF(IACC.NE.0) IRR=1 CALL PYLAMB(RLEN) ALAFT=RLEN RETURN END C********************************************************************* SUBROUTINE PYWWXX(TP,TM,XP,XM,V1P,V2P,V1M,V2M,IANSW,T,TAUP,TAUM) C...Subroutine to determine if two string pieces cross. C...P stands for W+ related quantities, M for W- related ones. C...Input: TP, TM times of creation of string pieces. C XP, XM vertices of creation of string pieces. C V1P, V2P, V1M, V2M velocity vectors of endpoint partons. C...Output: IANSW = 0 if not cross; = 1 if do cross; -1 if error. C T time of crossing (when they do cross, else 0). C TAUP, TAUM invariant times of strings at crossing point. DIMENSION XP(3),XM(3),V1P(3),V2P(3),V1M(3),V2M(3) DIMENSION Q(4,3),XPP(3),XMM(3) C...Function to do determinant of three rows of q matrix. DETER(I,J,K)=Q(I,1)*Q(J,2)*Q(K,3)-Q(I,1)*Q(K,2)*Q(J,3)+ &Q(J,1)*Q(K,2)*Q(I,3)-Q(J,1)*Q(I,2)*Q(K,3)+ &Q(K,1)*Q(I,2)*Q(J,3)-Q(K,1)*Q(J,2)*Q(I,3) C...Define q matrix and find t. DO 100 J=1,3 Q(1,J)=V2P(J)-V1P(J) Q(2,J)=-(V2M(J)-V1M(J)) Q(3,J)=XP(J)-XM(J)-TP*V1P(J)+TM*V1M(J) Q(4,J)=V1P(J)-V1M(J) 100 CONTINUE T=-DETER(1,2,3)/DETER(1,2,4) C...Find alpha and beta. S11=Q(1,1)*(T-TP) S12=Q(2,1)*(T-TM) S13=Q(3,1)+Q(4,1)*T S21=Q(1,2)*(T-TP) S22=Q(2,2)*(T-TM) S23=Q(3,2)+Q(4,2)*T DEN=S11*S22-S12*S21 ALP=(S12*S23-S22*S13)/DEN BET=(S21*S13-S11*S23)/DEN C...Check if solution acceptable. IANSW=1 IF(T.LT.TP.OR.T.LT.TM) IANSW=0 IF(ALP.LT.0..OR.ALP.GT.1.) IANSW=0 IF(BET.LT.0..OR.BET.GT.1.) IANSW=0 C...Find point of crossing and check that not inconsistent. DO 110 J=1,3 XPP(J)=XP(J)+(V1P(J)+ALP*(V2P(J)-V1P(J)))*(T-TP) XMM(J)=XM(J)+(V1M(J)+BET*(V2M(J)-V1M(J)))*(T-TM) 110 CONTINUE D2PM=(XPP(1)-XMM(1))**2+(XPP(2)-XMM(2))**2+(XPP(3)-XMM(3))**2 D2P=XPP(1)**2+XPP(2)**2+XPP(3)**2 D2M=XMM(1)**2+XMM(2)**2+XMM(3)**2 IF(D2PM.GT.1E-4*(D2P+D2M)) IANSW=-1 C...Find string eigentimes at crossing. IF(IANSW.EQ.1) THEN TAUP=SQRT(MAX(0.,(T-TP)**2-(XPP(1)-XP(1))**2- & (XPP(2)-XP(2))**2-(XPP(3)-XP(3))**2)) TAUM=SQRT(MAX(0.,(T-TM)**2-(XMM(1)-XM(1))**2- & (XMM(2)-XM(2))**2-(XMM(3)-XM(3))**2)) ELSE TAUP=0. TAUM=0. ENDIF RETURN END C********************************************************************* SUBROUTINE PYLAMB(RLEN) C...Find lambda measure length of string. COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) PARAMETER (PMRHO=0.8) C...Loop through partons; find where there is a string piece. RLEN=0. DO 100 I=1,N IF(K(I,1).NE.2.AND.K(I,1).NE.12) GOTO 100 C...Calculate string piece mass and add. PM2=(P(I,4)+P(I+1,4))**2-(P(I,1)+P(I+1,1))**2- &(P(I,2)+P(I+1,2))**2-(P(I,3)+P(I+1,3))**2 DLEN=LOG(1.+SQRT(MAX(0.,PM2))/PMRHO) RLEN=RLEN+DLEN 100 CONTINUE RETURN END