C...This file contains routines for Bose-Einstein studies. C...The routines are intended to replace the standard LUBOEI routine C...found in JETSET 7.4. C...The physics of the routines is described in C...Leif Lonnblad and Torbjorn Sjostrand, C..."Modelling Bose-Einstein correlations at LEP 2" C...preprint no. LU TP 97-30, NORDITA-97/75 P and hep-ph/9711460, C...to appear in European Physical Journal C. C...Here you also find further information underlying the meaning C...of the switches below and how to combine them sensibly. C********************************************************************* C...MSTJ(53) In e+e- -> W+W-, apply BE algorithm C... 0 => on all pion pairs C... 1 => only on pairs were both pions come from the same W C... 2 => only on pairs were the pions come from different Ws C... -1 => on all pairs except unequal pions coming from C... different Ws C... -2 => When calculating balancing shifts for pions from same C... W, only consider pairs from this W C...MSTJ(54) Alternative local energy compensation C... 0 => Global energy compensation C... 1 => Compensate with identical pairs by negative BE C... enhancement with a third of the radius C... -5 => Compansate with pair giving the smallest string length C... -2 => Compensate with pair giving the smallest invariant mass C...MSTJ(55) Calculation of difference vector C... 0 => In the lab frame C... 1 => In the CMS of the given pair. C...MSTJ(56) In e+e- -> W+W-, include distance between W's C... 0 => Radius is the same for all pairs C... 1 => Radius for pairs from different W's is R+deltaR_WW C...MSTJ(57) Penalty for shifting particles with close by identical C... neighbors in local energy compensation. C... 0 => No penalty C... 1 => Penalty C...PARJ(95)(R) Set to the energy imbalance after the BE algorithm, C... before rescaling of momenta. C...PARJ(96)(R) Set to the alpha needed to retain energy-momentum C... conservation in each event for relevant models. C...PARJ(97)(D=0) If larger than 0 used to get alternative shape for C... MSTJ(54)=1 C********************************************************************* SUBROUTINE LUBOEI(NSAV) C...Purpose: to modify event so as to approximately take into account C...Bose-Einstein effects according to a simple phenomenological C...parametrization. IMPLICIT DOUBLE PRECISION(D) COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) SAVE /LUJETS/,/LUDAT1/ COMMON /LUDAT2/ KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) SAVE /LUDAT2/ DIMENSION DPS(4),KFBE(9),NBE(0:10),BEI(100),BEI3(100), $ BEIW(100),BEI3W(100) DATA KFBE/211,-211,111,321,-321,130,310,221,331/ SDIP(I,J)=((P(I,4)+P(J,4))**2-(P(I,3)+P(J,3))**2- $ (P(I,2)+P(J,2))**2-(P(I,1)+P(J,1))**2) C...Boost event to overall CM frame. Calculate CM energy. IF((MSTJ(51).NE.1.AND.MSTJ(51).NE.2).OR.N-NSAV.LE.1) RETURN DO 100 J=1,4 DPS(J)=0. 100 CONTINUE DO 120 I=1,N KFA=IABS(K(I,2)) IF(K(I,1).LE.10.AND.((KFA.GT.10.AND.KFA.LE.20).OR.KFA.EQ.22).AND. &K(I,3).GT.0) THEN KFMA=IABS(K(K(I,3),2)) IF(KFMA.GT.10.AND.KFMA.LE.80) K(I,1)=-K(I,1) ENDIF IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120 DO 110 J=1,4 DPS(J)=DPS(J)+P(I,J) 110 CONTINUE 120 CONTINUE CALL LUDBRB(0,0,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4), &-DPS(3)/DPS(4)) PECM=0. DO 130 I=1,N IF(K(I,1).GE.1.AND.K(I,1).LE.10) PECM=PECM+P(I,4) 130 CONTINUE C...Reserve copy of particles by species at end of record. IWP=0 IWN=0 NBE(0)=N+MSTU(3) NMAX=NBE(0) SMMIN=PECM DO 160 IBE=1,MIN(10,MSTJ(52)+1) NBE(IBE)=NBE(IBE-1) DO 150 I=NSAV+1,N IF (IBE.EQ.MIN(10,MSTJ(52)+1)) THEN DO 152 IIBE=1,IBE-1 IF (K(I,2).EQ.KFBE(IIBE)) GOTO 150 152 CONTINUE ELSE IF(K(I,2).NE.KFBE(IBE)) GOTO 150 ENDIF IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 150 IF(NBE(IBE).GE.MSTU(4)-MSTU(32)-5) THEN CALL LUERRM(11,'(LUBOEI:) no more memory left in LUJETS') RETURN ENDIF NBE(IBE)=NBE(IBE)+1 NMAX=NBE(IBE) K(NBE(IBE),1)=I K(NBE(IBE),5)=0 SMMIN=MIN(SMMIN,P(I,5)) IF (MSTJ(53).NE.0.OR.MSTJ(56).GT.0) THEN IM=I 151 IF (K(IM,3).GT.0) THEN IM=K(IM,3) IF (ABS(K(IM,2)).NE.24) GOTO 151 K(NBE(IBE),5)=K(IM,2) IF (IWP.EQ.0.AND.K(IM,2).EQ.24) IWP=IM IF (IWN.EQ.0.AND.K(IM,2).EQ.-24) IWN=IM ENDIF ENDIF DO 140 J=1,3 P(NBE(IBE),J)=0. V(NBE(IBE),J)=0. 140 CONTINUE P(NBE(IBE),5)=-1.0 150 CONTINUE 160 CONTINUE IF(NBE(MIN(9,MSTJ(52)))-NBE(0).LE.1) GOTO 280 C...Calculate separation between W+ and W- SIGW=PARJ(93) IF (IWP.GT.0.AND.IWN.GT.0.AND.MSTJ(56).GT.0) THEN DMW=PMAS(24,1) DGW=PMAS(24,2) DMP=P(IWP,5) DMN=P(IWN,5) TAUPD=DMP/SQRT((DMP**2-DMW**2)**2+(DGW*(DMP**2)/DMW)**2) TAUND=DMN/SQRT((DMN**2-DMW**2)**2+(DGW*(DMN**2)/DMW)**2) TAUP=-TAUPD*LOG(RLU(IDUM)) TAUN=-TAUND*LOG(RLU(IDUM)) DXP=TAUP*PLU(IWP,8)/DMP DXN=TAUN*PLU(IWN,8)/DMN DX=DXP+DXN SIGW=1.0/(1.0/PARJ(93)+REAL(MSTJ(56))*DX) ELSE SIGW=PARJ(93) ENDIF IF (MSTJ(57).EQ.1) THEN DO 221 IBE=1,MIN(9,MSTJ(52)) DO 222 I1M=NBE(IBE-1)+1,NBE(IBE)-1 Q2MIN=PECM**2 I1=K(I1M,1) DO 223 I2M=NBE(IBE-1)+1,NBE(IBE)-1 IF (I2M.EQ.I1M) GOTO 223 I2=K(I2M,1) Q2=(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- $ (P(I1,5)+P(I2,5))**2 IF (Q2.GT.0.0.AND.Q2.LT.Q2MIN) THEN Q2MIN=Q2 ENDIF 223 CONTINUE P(I1M,5)=Q2MIN 222 CONTINUE 221 CONTINUE ENDIF C...Tabulate integral for subsequent momentum shift. DO 220 IBE=1,MIN(9,MSTJ(52)) IF(IBE.NE.1.AND.IBE.NE.4.AND.IBE.LE.7) GOTO 180 IF(IBE.EQ.1.AND.MAX(NBE(1)-NBE(0),NBE(2)-NBE(1),NBE(3)-NBE(2)) &.LE.1) GOTO 180 IF(IBE.EQ.4.AND.MAX(NBE(4)-NBE(3),NBE(5)-NBE(4),NBE(6)-NBE(5), &NBE(7)-NBE(6)).LE.1) GOTO 180 IF(IBE.GE.8.AND.NBE(IBE)-NBE(IBE-1).LE.1) GOTO 180 IF(IBE.EQ.1) PMHQ=2.*ULMASS(211) IF(IBE.EQ.4) PMHQ=2.*ULMASS(321) IF(IBE.EQ.8) PMHQ=2.*ULMASS(221) IF(IBE.EQ.9) PMHQ=2.*ULMASS(331) QDEL=0.1*MIN(PMHQ,PARJ(93)) QDEL3=0.1*MIN(PMHQ,PARJ(93)*3.0) QDELW=0.1*MIN(PMHQ,SIGW) QDEL3W=0.1*MIN(PMHQ,SIGW*3.0) IF(MSTJ(51).EQ.1) THEN NBIN=MIN(100,NINT(9.*PARJ(93)/QDEL)) NBIN3=MIN(100,NINT(27.*PARJ(93)/QDEL3)) NBINW=MIN(100,NINT(9.*SIGW/QDELW)) NBIN3W=MIN(100,NINT(27.*SIGW/QDEL3W)) BEEX=EXP(0.5*QDEL/PARJ(93)) BEEX3=EXP(0.5*QDEL3/(3.0*PARJ(93))) BEEXW=EXP(0.5*QDELW/SIGW) BEEX3W=EXP(0.5*QDEL3W/(3.0*SIGW)) BERT=EXP(-QDEL/PARJ(93)) BERT3=EXP(-QDEL3/(3.0*PARJ(93))) BERTW=EXP(-QDELW/SIGW) BERT3W=EXP(-QDEL3W/(3.0*SIGW)) ELSE NBIN=MIN(100,NINT(3.*PARJ(93)/QDEL)) NBIN3=MIN(100,NINT(9.*PARJ(93)/QDEL3)) NBINW=MIN(100,NINT(3.*SIGW/QDELW)) NBIN3W=MIN(100,NINT(9.*SIGW/QDEL3W)) ENDIF DO 170 IBIN=1,NBIN QBIN=QDEL*(IBIN-0.5) BEI(IBIN)=QDEL*(QBIN**2+QDEL**2/12.)/SQRT(QBIN**2+PMHQ**2) IF(MSTJ(51).EQ.1) THEN BEEX=BEEX*BERT BEI(IBIN)=BEI(IBIN)*BEEX ELSE BEI(IBIN)=BEI(IBIN)*EXP(-(QBIN/PARJ(93))**2) ENDIF IF(IBIN.GE.2) BEI(IBIN)=BEI(IBIN)+BEI(IBIN-1) 170 CONTINUE DO 173 IBIN=1,NBIN3 QBIN=QDEL3*(IBIN-0.5) BEI3(IBIN)=QDEL3*(QBIN**2+QDEL3**2/12.)/SQRT(QBIN**2+PMHQ**2) IF(MSTJ(51).EQ.1) THEN BEEX3=BEEX3*BERT3 BEI3(IBIN)=BEI3(IBIN)*BEEX3 ELSE BEI3(IBIN)=BEI3(IBIN)*EXP(-(QBIN/(3.0*PARJ(93)))**2) ENDIF IF(IBIN.GE.2) BEI3(IBIN)=BEI3(IBIN)+BEI3(IBIN-1) 173 CONTINUE DO 171 IBIN=1,NBINW QBIN=QDELW*(IBIN-0.5) BEIW(IBIN)=QDELW*(QBIN**2+QDELW**2/12.)/SQRT(QBIN**2+PMHQ**2) IF(MSTJ(51).EQ.1) THEN BEEXW=BEEXW*BERTW BEIW(IBIN)=BEIW(IBIN)*BEEXW ELSE BEIW(IBIN)=BEIW(IBIN)*EXP(-(QBIN/SIGW)**2) ENDIF IF(IBIN.GE.2) BEIW(IBIN)=BEIW(IBIN)+BEIW(IBIN-1) 171 CONTINUE DO 174 IBIN=1,NBIN3W QBIN=QDEL3W*(IBIN-0.5) BEI3W(IBIN)=QDEL3W*(QBIN**2+QDEL3W**2/12.)/SQRT(QBIN**2+PMHQ**2) IF(MSTJ(51).EQ.1) THEN BEEX3W=BEEX3W*BERT3W BEI3W(IBIN)=BEI3W(IBIN)*BEEX3W ELSE BEI3W(IBIN)=BEI3W(IBIN)*EXP(-(QBIN/(3.0*SIGW))**2) ENDIF IF(IBIN.GE.2) BEI3W(IBIN)=BEI3W(IBIN)+BEI3W(IBIN-1) 174 CONTINUE C...Loop through particle pairs and find old relative momentum. 180 DO 210 I1M=NBE(IBE-1)+1,NBE(IBE)-1 I1=K(I1M,1) DO 200 I2M=I1M+1,NBE(IBE) IF (MSTJ(53).EQ.1.AND.K(I1M,5).NE.K(I2M,5)) GOTO 200 IF (MSTJ(53).EQ.2.AND.K(I1M,5).EQ.K(I2M,5)) GOTO 200 I2=K(I2M,1) Q2OLD=(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-(P(I1,5)+P(I2,5))**2 IF (Q2OLD.LE.0.0) GOTO 200 QOLD=SQRT(Q2OLD) C...Calculate new relative momentum. QMOV=0.0 QMOV3=0.0 QMOVW=0.0 QMOV3W=0.0 IF(QOLD.LT.1E-3*QDEL) THEN GOTO 203 ELSEIF(QOLD.LE.QDEL) THEN QMOV=QOLD/3. ELSEIF(QOLD.LT.(NBIN-0.1)*QDEL) THEN RBIN=QOLD/QDEL IBIN=RBIN RINP=(RBIN**3-IBIN**3)/(3*IBIN*(IBIN+1)+1) QMOV=(BEI(IBIN)+RINP*(BEI(IBIN+1)-BEI(IBIN)))* & SQRT(Q2OLD+PMHQ**2)/Q2OLD ELSE QMOV=BEI(NBIN)*SQRT(Q2OLD+PMHQ**2)/Q2OLD ENDIF 203 Q2NEW=Q2OLD*(QOLD/(QOLD+3.*PARJ(92)*QMOV))**(2./3.) IF(QOLD.LT.1E-3*QDEL3) THEN GOTO 204 ELSEIF(QOLD.LE.QDEL3) THEN QMOV3=QOLD/3. ELSEIF(QOLD.LT.(NBIN3-0.1)*QDEL3) THEN RBIN3=QOLD/QDEL3 IBIN3=RBIN3 RINP3=(RBIN3**3-IBIN3**3)/(3*IBIN3*(IBIN3+1)+1) QMOV3=(BEI3(IBIN3)+RINP3*(BEI3(IBIN3+1)-BEI3(IBIN3)))* & SQRT(Q2OLD+PMHQ**2)/Q2OLD ELSE QMOV3=BEI3(NBIN3)*SQRT(Q2OLD+PMHQ**2)/Q2OLD ENDIF 204 Q2NEW3=Q2OLD*(QOLD/(QOLD+3.*PARJ(92)*QMOV3))**(2./3.) RSCALE=1.0 IF (PARJ(97).GT.0.0) $ RSCALE=1.0-EXP(-(QOLD/(PARJ(93)*PARJ(97)))**2) IF (MSTJ(56).LE.0.OR.IWP.EQ.0.OR.IWN.EQ.0.OR. $ K(I1M,5).EQ.K(I2M,5)) GOTO 207 IF(QOLD.LT.1E-3*QDELW) THEN GOTO 205 ELSEIF(QOLD.LE.QDELW) THEN QMOVW=QOLD/3. ELSEIF(QOLD.LT.(NBINW-0.1)*QDELW) THEN RBINW=QOLD/QDELW IBINW=RBINW RINPW=(RBINW**3-IBINW**3)/(3*IBINW*(IBINW+1)+1) QMOVW=(BEIW(IBINW)+RINPW*(BEIW(IBINW+1)-BEIW(IBINW)))* & SQRT(Q2OLD+PMHQ**2)/Q2OLD ELSE QMOVW=BEIW(NBINW)*SQRT(Q2OLD+PMHQ**2)/Q2OLD ENDIF 205 Q2NEW=Q2OLD*(QOLD/(QOLD+3.*PARJ(92)*QMOVW))**(2./3.) IF(QOLD.LT.1E-3*QDEL3W) THEN GOTO 206 ELSEIF(QOLD.LE.QDEL3W) THEN QMOV3W=QOLD/3. ELSEIF(QOLD.LT.(NBIN3W-0.1)*QDEL3W) THEN RBIN3W=QOLD/QDEL3W IBIN3W=RBIN3W RINP3W=(RBIN3W**3-IBIN3W**3)/(3*IBIN3W*(IBIN3W+1)+1) QMOV3W=(BEI3W(IBIN3W)+RINP3W*(BEI3W(IBIN3W+1)-BEI3W(IBIN3W)))* & SQRT(Q2OLD+PMHQ**2)/Q2OLD ELSE QMOV3W=BEI3W(NBIN3W)*SQRT(Q2OLD+PMHQ**2)/Q2OLD ENDIF 206 Q2NEW3=Q2OLD*(QOLD/(QOLD+3.*PARJ(92)*QMOV3W))**(2./3.) IF (PARJ(97).GT.0.0) $ RSCALE=1.0-EXP(-(QOLD/(SIGW*PARJ(97)))**2) 207 CALL LUBESQ(I1,I2,NMAX,Q2OLD,Q2NEW) DO 190 J=1,3 P(I1M,J)=P(I1M,J)+P(NMAX+1,J) P(I2M,J)=P(I2M,J)+P(NMAX+2,J) 190 CONTINUE IF(MSTJ(54).EQ.1) THEN CALL LUBESQ(I1,I2,NMAX,Q2OLD,Q2NEW3) DO 193 J=1,3 V(I1M,J)=V(I1M,J)+P(NMAX+1,J)*RSCALE V(I2M,J)=V(I2M,J)+P(NMAX+2,J)*RSCALE 193 CONTINUE ELSEIF(MSTJ(54).LE.-1) THEN EDEL=P(I1,4)+P(I2,4)- $ SQRT(MAX(Q2NEW-Q2OLD+(P(I1,4)+P(I2,4))**2,0.0)) A2=(P(I1,1)-P(I2,1))**2+(P(I1,2)-P(I2,2))**2+ $ (P(I1,3)-P(I2,3))**2 WMAX=-1.0E20 MI3=0 MI4=0 IF (MSTJ(54).EQ.-10) THEN SMIN=PECM**2 II1=I1 II2=I2 IF (RLU(IDUM).GT.0.5) THEN II1=I2 II2=I1 ENDIF DO 1095 I3M=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1)) IF (I3M.EQ.I1M.OR.I3M.EQ.I2M) GOTO 1095 IF (MSTJ(53).EQ.1.AND.K(I3M,5).NE.K(I1M,5)) GOTO 1095 IF (MSTJ(53).EQ.-2.AND.K(I1M,5).EQ.K(I2M,5).AND. $ K(I3M,5).NE.K(I1M,5)) GOTO 1095 I3=K(I3M,1) IF(K(I3,2).EQ.K(I1,2)) GOTO 1095 S13=(P(II1,4)+P(I3,4))**2- $ (P(II1,3)+P(I3,3))**2- $ (P(II1,2)+P(I3,2))**2- $ (P(II1,1)+P(I3,1))**2 IF (S13.GE.SMIN) GOTO 1095 MI3=I3M SMIN=S13 1095 CONTINUE IF (MI3.GT.0) THEN I3=K(MI3,1) SMIN=PECM**2 DO 1096 I4M=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1)) IF (I4M.EQ.I1M.OR.I4M.EQ.I2M.OR.I4M.EQ.MI3) GOTO 1096 IF (MSTJ(53).EQ.1.AND.K(I4M,5).NE.K(I1M,5)) GOTO 1096 IF (MSTJ(53).EQ.-2.AND.K(I1M,5).EQ.K(I2M,5).AND. $ K(I4M,5).NE.K(I1M,5)) GOTO 1096 I4=K(I4M,1) IF(K(I4,2).EQ.K(I1,2).OR.K(I4,2).EQ.K(I3,2)) GOTO 1096 S24=(P(II2,4)+P(I4,4))**2- $ (P(II2,3)+P(I4,3))**2- $ (P(II2,2)+P(I4,2))**2- $ (P(II2,1)+P(I4,1))**2 IF (S24.GE.SMIN) GOTO 1096 IF ((P(I3,4)+P(I4,4)+EDEL)**2.LT.(P(I3,1)+P(I4,1))**2+ $ (P(I3,2)+P(I4,2))**2+(P(I3,3)+P(I4,3))**2+ $ (P(I3,5)+P(I4,5))**2) GOTO 1096 MI4=I4M SMIN=S24 1096 CONTINUE ENDIF ELSE S12=SDIP(I1,I2) SM1=(P(I1,5)+SMMIN)**2 DO 195 I3M=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1)) IF (I3M.EQ.I1M.OR.I3M.EQ.I2M) GOTO 195 IF (MSTJ(53).EQ.1.AND.K(I3M,5).NE.K(I1M,5)) GOTO 195 IF (MSTJ(53).EQ.-2.AND.K(I1M,5).EQ.K(I2M,5).AND. $ K(I3M,5).NE.K(I1M,5)) GOTO 195 I3=K(I3M,1) IF(K(I3,2).EQ.K(I1,2))GOTO 195 S13=SDIP(I1,I3) S23=SDIP(I2,I3) SM3=(P(I3,5)+SMMIN)**2 IF (MSTJ(54).EQ.-5) THEN WI=1.0/(MIN(S12*SM3,S13*MIN(SM1,SM3), $ S23*MIN(SM1,SM3))*SM1) ELSE WI=1.0/((P(I1,4)+P(I2,4)+P(I3,4))**2- $ (P(I1,3)+P(I2,3)+P(I3,3))**2- $ (P(I1,2)+P(I2,2)+P(I3,2))**2- $ (P(I1,1)+P(I2,1)+P(I3,1))**2) ENDIF IF (MSTJ(57).EQ.1.AND.P(I3M,5).GT.0) $ WI=WI*(1.0-EXP(P(I3M,5)/(PARJ(93)**2))) IF (WI.LE.WMAX) GOTO 195 DO 196 I4M=I3M+1,NBE(MIN(10,MSTJ(52)+1)) IF (I4M.EQ.I1M.OR.I4M.EQ.I2M) GOTO 196 IF (MSTJ(53).EQ.1.AND.K(I4M,5).NE.K(I1M,5)) GOTO 196 IF (MSTJ(53).EQ.-2.AND.K(I1M,5).EQ.K(I2M,5).AND. $ K(I4M,5).NE.K(I1M,5)) GOTO 196 I4=K(I4M,1) IF (K(I3,2).EQ.K(I4,2).OR.K(I4,2).EQ.K(I1,2)) GOTO 196 IF ((P(I3,4)+P(I4,4)+EDEL)**2.LT.(P(I3,1)+P(I4,1))**2+ $ (P(I3,2)+P(I4,2))**2+(P(I3,3)+P(I4,3))**2+ $ (P(I3,5)+P(I4,5))**2) GOTO 196 IF (MSTJ(54).EQ.-5) THEN S14=SDIP(I1,I4) S24=SDIP(I2,I4) S34=SDIP(I3,I4) W=S12*MIN(MIN(S23,S24),MIN(S13,S14))*S34 W=MIN(W,S13*MIN(MIN(S23,S34),S12)*S24) W=MIN(W,S14*MIN(MIN(S24,S34),S12)*S23) W=MIN(W,MIN(S23,S24)*S13*S14) W=1.0/W ELSE C...weight=1-cos(theta)/mtot2 S1234=(P(I1,4)+P(I2,4)+P(I3,4)+P(I4,4))**2- $ (P(I1,3)+P(I2,3)+P(I3,3)+P(I4,3))**2- $ (P(I1,2)+P(I2,2)+P(I3,2)+P(I4,2))**2- $ (P(I1,1)+P(I2,1)+P(I3,1)+P(I4,1))**2 W=1.0/S1234 IF (W.LE.WMAX) GOTO 196 IF (MSTJ(54).LE.-2) THEN AB2=(P(I1,1)-P(I2,1)+P(I3,1)-P(I4,1))**2+ $ (P(I1,2)-P(I2,2)+P(I3,2)-P(I4,2))**2+ $ (P(I1,3)-P(I2,3)+P(I3,3)-P(I4,3))**2 B2=(P(I3,1)-P(I4,1))**2+(P(I3,2)-P(I4,2))**2+ $ (P(I3,3)-P(I4,3))**2 CTH=ABS(0.5*(AB2-A2-B2)/SQRT(A2*B2)) W=(1.0-CTH)/S1234 IF (MSTJ(54).EQ.-3) W=W*(1.0-CTH) IF (MSTJ(54).EQ.-4) W=W*(1.0-CTH)**3 ENDIF ENDIF IF (MSTJ(57).EQ.1.AND.P(I3M,5).GT.0) $ W=W*(1.0-EXP(P(I3M,5)/(PARJ(93)**2))) IF (MSTJ(57).EQ.1.AND.P(I4M,5).GT.0) $ W=W*(1.0-EXP(P(I4M,5)/(PARJ(93)**2))) IF (W.LE.WMAX) GOTO 196 MI3=I3M MI4=I4M WMAX=W 196 CONTINUE 195 CONTINUE ENDIF IF (MI4.EQ.0) GOTO 200 I3=K(MI3,1) I4=K(MI4,1) EOLD=P(I3,4)+P(I4,4) ENEW=EOLD+EDEL P2=(P(I3,1)+P(I4,1))**2+(P(I3,2)+P(I4,2))**2+ $ (P(I3,3)+P(I4,3))**2 Q2NEWP=MAX(0.0,ENEW**2-P2-(P(I3,5)+P(I4,5))**2) Q2OLDP=MAX(0.0,EOLD**2-P2-(P(I3,5)+P(I4,5))**2) CALL LUBESQ(I3,I4,NMAX,Q2OLDP,Q2NEWP) DO 197 J=1,3 V(MI3,J)=V(MI3,J)+P(NMAX+1,J) V(MI4,J)=V(MI4,J)+P(NMAX+2,J) 197 CONTINUE ENDIF 200 CONTINUE 210 CONTINUE 220 CONTINUE C...Shift momenta and recalculate energies. ESUMP=0.0 ESUM=0.0 PROD=0.0 DO 240 IM=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1)) I=K(IM,1) ESUMP=ESUMP+P(I,4) DO 230 J=1,3 P(I,J)=P(I,J)+P(IM,J) 230 CONTINUE P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2) ESUM=ESUM+P(I,4) DO 233 J=1,3 PROD=PROD+V(IM,J)*P(I,J)/P(I,4) 233 CONTINUE 240 CONTINUE PARJ(96)=0.0 IF (MSTJ(54).NE.0.AND.PROD.NE.0.0) THEN 333 ALPHA=(ESUMP-ESUM)/PROD PARJ(96)=PARJ(96)+ALPHA PROD=0.0 ESUM=0.0 DO 300 IM=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1)) I=K(IM,1) DO 310 J=1,3 P(I,J)=P(I,J)+ALPHA*V(IM,J) 310 CONTINUE P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2) ESUM=ESUM+P(I,4) DO 313 J=1,3 PROD=PROD+V(IM,J)*P(I,J)/P(I,4) 313 CONTINUE 300 CONTINUE IF (PROD.NE.0.0.AND.ABS(ESUMP-ESUM)/PECM.GT.0.00001) GOTO 333 ENDIF C...Rescale all momenta for energy conservation. PES=0. PQS=0. DO 250 I=1,N IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 250 PES=PES+P(I,4) PQS=PQS+P(I,5)**2/P(I,4) 250 CONTINUE PARJ(95)=PES-PECM FAC=(PECM-PQS)/(PES-PQS) DO 270 I=1,N IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 270 DO 260 J=1,3 P(I,J)=FAC*P(I,J) 260 CONTINUE P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2) 270 CONTINUE C...Boost back to correct reference frame. 280 CALL LUDBRB(0,0,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),DPS(3)/DPS(4)) DO 290 I=1,N IF(K(I,1).LT.0) K(I,1)=-K(I,1) 290 CONTINUE RETURN END C********************************************************************* SUBROUTINE LUBESQ(I1,I2,NI,Q2OLD,Q2NEW) C...Purpose: To calculate the momentum shift i a system of two C...particles assuming the relative momentum squared C...should be shifted to Q2NEW. NI is the last position occupied in C.../LUJETS/ IMPLICIT DOUBLE PRECISION(D) COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) SAVE /LUJETS/,/LUDAT1/ DIMENSION DP(5) SAVE HC1 double precision dq2,dp2,dp12,se,de IF (MSTJ(55).EQ.-1) THEN HC1=(P(I1,4)+P(I2,4))**2-(Q2OLD-Q2NEW) HC2=(Q2OLD-Q2NEW)*(P(I1,4)-P(I2,4))**2 HA=0.5*(1.-SQRT(HC1*Q2NEW/(HC1*Q2OLD-HC2))) DO 100 J=1,3 PD=HA*(P(I2,J)-P(I1,J)) P(NI+1,J)=PD P(NI+2,J)=-PD 100 CONTINUE RETURN ENDIF IF (MSTJ(55).EQ.0) THEN DQ2=Q2NEW-Q2OLD DP2=(P(I1,1)-P(I2,1))**2+(P(I1,2)-P(I2,2))**2+ $ (P(I1,3)-P(I2,3))**2 DP12=P(I1,1)**2+P(I1,2)**2+P(I1,3)**2 $ -P(I2,1)**2-P(I2,2)**2-P(I2,3)**2 SE=P(I1,4)+P(I2,4) DE=P(I1,4)-P(I2,4) DQ2SE=DQ2+SE**2 DA=SE*DE*DP12-DP2*DQ2SE DB=DP2*DQ2SE-DP12**2 HA=(DA+SQRT(MAX(DA**2+DQ2*(DQ2+SE**2-DE**2)*DB,0D0)))/(2D0*DB) DO 101 J=1,3 PD=HA*(P(I1,J)-P(I2,J)) P(NI+1,J)=PD P(NI+2,J)=-PD 101 CONTINUE RETURN ENDIF K(NI+1,1)=1 K(NI+2,1)=1 DO 200 J=1,5 P(NI+1,J)=P(I1,J) P(NI+2,J)=P(I2,J) DP(J)=DBLE(P(I1,J))+DBLE(P(I2,J)) 200 CONTINUE C...Boost to cms and rotate first particle to z-axis CALL LUDBRB(NI+1,NI+2,0.0,0.0, $ -DP(1)/DP(4),-DP(2)/DP(4),-DP(3)/DP(4)) PHI=ULANGL(P(NI+1,1),P(NI+1,2)) THE=ULANGL(P(NI+1,3),SQRT(P(NI+1,1)**2+P(NI+1,2)**2)) S=Q2NEW+(P(I1,5)+P(I2,5))**2 PZ=0.5*SQRT(Q2NEW*(S-(P(I1,5)-P(I2,5))**2)/S) P(NI+1,1)=0.0 P(NI+1,2)=0.0 P(NI+1,3)=PZ P(NI+1,4)=SQRT(PZ**2+P(I1,5)**2) P(NI+2,1)=0.0 P(NI+2,2)=0.0 P(NI+2,3)=-PZ P(NI+2,4)=SQRT(PZ**2+P(I2,5)**2) DP(4)=SQRT(DP(1)**2+DP(2)**2+DP(3)**2+DBLE(S)) CALL LUDBRB(NI+1,NI+2,THE,PHI, $ DP(1)/DP(4),DP(2)/DP(4),DP(3)/DP(4)) DO 300 J=1,3 P(NI+1,J)=P(NI+1,J)-P(I1,J) P(NI+2,J)=P(NI+2,J)-P(I2,J) 300 CONTINUE RETURN END