C...This package contains a replacement routine for the LUBOEI routine C...of the standard JETSET distribution. By using this replacement, C...plus the new LUBOE2 routine (also found below), it is possible to C...gain access to several further options for the treatment of C...Bose-Einstein effects. These options are listed below. C...The LUBOE2 routine is courtesy L. Lonnblad, leif@nordita.dk. C********************************************************************* C...MSTJ(51)=3 => Use new BE algorithm C...MSTJ(51)=4 => Use new BE algorithm but don't treat pairs coming from C... the decay of a single hadron C...MSTJ(53) In e+e- -> W+W-, apply new 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...PARJ(94)(D=0) Epsilon in equation for new BE algorithm - typically C... the fraction of pairs with identical pions w.r.t. C... the total number of pion pairs C...PARJ(95)(R) Set to the energy imbalance after the BE algorithm, C... before rescaling of momenta. 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/ DIMENSION DPS(4),KFBE(9),NBE(0:9),BEI(100) DATA KFBE/211,-211,111,321,-321,130,310,221,331/ C...Boost event to overall CM frame. Calculate CM energy. IF (MSTJ(51).EQ.3.OR.MSTJ(51).EQ.4) THEN CALL LUBOE2(NSAV) RETURN ENDIF 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. NBE(0)=N+MSTU(3) DO 160 IBE=1,MIN(9,MSTJ(52)) NBE(IBE)=NBE(IBE-1) DO 150 I=NSAV+1,N IF(K(I,2).NE.KFBE(IBE)) GOTO 150 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 K(NBE(IBE),1)=I DO 140 J=1,3 P(NBE(IBE),J)=0. 140 CONTINUE 150 CONTINUE 160 CONTINUE IF(NBE(MIN(9,MSTJ(52)))-NBE(0).LE.1) GOTO 280 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)) IF(MSTJ(51).EQ.1) THEN NBIN=MIN(100,NINT(9.*PARJ(93)/QDEL)) BEEX=EXP(0.5*QDEL/PARJ(93)) BERT=EXP(-QDEL/PARJ(93)) ELSE NBIN=MIN(100,NINT(3.*PARJ(93)/QDEL)) 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 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) I2=K(I2M,1) Q2OLD=MAX(0.,(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) QOLD=SQRT(Q2OLD) C...Calculate new relative momentum. IF(QOLD.LT.1E-3*QDEL) THEN GOTO 200 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 Q2NEW=Q2OLD*(QOLD/(QOLD+3.*PARJ(92)*QMOV))**(2./3.) C...Calculate and save shift to be performed on three-momenta. 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 190 J=1,3 PD=HA*(P(I2,J)-P(I1,J)) P(I1M,J)=P(I1M,J)+PD P(I2M,J)=P(I2M,J)-PD 190 CONTINUE 200 CONTINUE 210 CONTINUE 220 CONTINUE C...Shift momenta and recalculate energies. DO 240 IM=NBE(0)+1,NBE(MIN(9,MSTJ(52))) I=K(IM,1) 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) 240 CONTINUE 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 LUBOE2(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 /LUBEDA/ SIEQ(0:1000),SINE(0:1000), $ RAT,XLAM,SIG,SM4,QMAX,QD SAVE /LUBEDA/ DIMENSION DPS(4) C...Part of integrand left for numerical integration for identical (FIP) C...and non-identical (FIM) pairs FIP(Q)=(Q*SQRT(Q**2+SM4)-SM4*LOG(Q+SQRT(Q**2+SM4)))* $ EXP(MIN((Q/SIG)**2,35.0))*(RAT-1.0)*XLAM*Q/ $ ((SIG*(EXP(MIN((Q/SIG)**2,35.0))+RAT*XLAM))**2) FIM(Q)=(Q*SQRT(Q**2+SM4)-SM4*LOG(Q+SQRT(Q**2+SM4)))* $ EXP(MIN((Q/SIG)**2,35.0))*RAT*XLAM*Q/ $ ((SIG*(EXP(MIN((Q/SIG)**2,35.0))+RAT*XLAM))**2) C...Partially integration from 0 to Qprime for identical (FPIP) and C...and non-identical (FPIM) pairs FPIP(Q)=0.5*(Q*SQRT(Q**2+SM4)-SM4*LOG(Q+SQRT(Q**2+SM4)))* $ (1.0+XLAM*EXP(-MIN((Q/SIG)**2,70.0)))/ $ (1.0+RAT*XLAM*EXP(-MIN((Q/SIG)**2,70.0)))+ $ 0.5*SM4*LOG(SQRT(SM4))*(1.0+XLAM)/(1.0+XLAM*RAT) FPIM(Q)=0.5*(Q*SQRT(Q**2+SM4)-SM4*LOG(Q+SQRT(Q**2+SM4)))/ $ (1.0+RAT*XLAM*EXP(-MIN((Q/SIG)**2,70.0)))+ $ 0.5*SM4*LOG(SQRT(SM4))/(1.0+XLAM*RAT) C...Enhancement function f_2 times phase space for identical (FP) and C...non-identical (FM) pairs FP(Q)=((Q**2)/SQRT(Q**2+SM4))* $ (1.0+XLAM*EXP(-MIN((Q/SIG)**2,70.0)))/ $ (1.0+RAT*XLAM*EXP(-MIN((Q/SIG)**2,70.0))) FM(Q)=((Q**2)/SQRT(Q**2+SM4))/ $ (1.0+RAT*XLAM*EXP(-MIN((Q/SIG)**2,70.0))) C...Full integral from 0 to Qprime in the limit of small Qprime for C...identical (FISP) and non-identical (FISM) pairs FISP(Q)=(Q**3)*(1.0+XLAM)/(3.0*SQRT(SM4)*(1.0+RAT*XLAM)) FISM(Q)=(Q**3)/(3.0*SQRT(SM4)*(1.0+RAT*XLAM)) IF((MSTJ(51).NE.3.AND.MSTJ(51).NE.4).OR.N-NSAV.LE.1) RETURN C...Check if parameters have changed, retabulate integral if necessary IF (RAT.NE.PARJ(94).OR.SIG.NE.PARJ(93).OR.XLAM.NE.PARJ(92)) THEN RAT=PARJ(94) SIG=PARJ(93) XLAM=PARJ(92) SM4=4.0*ULMASS(211)**2 SIEQ(0)=0.0 SINE(0)=0.0 QMAX=MAX(SIG,SQRT(SM4))*10.0 QD=QMAX*0.001 DO 100 I=1,1000 QP=(REAL(I)-0.5)*QD SIEQ(I)=SIEQ(I-1)+QD*FIP(QP) SINE(I)=SINE(I-1)+QD*FIM(QP) 100 CONTINUE ENDIF C...Boost event to overall CM frame. Calculate CM energy. Leave out C...some leptons DO 110 J=1,4 DPS(J)=0.0 110 CONTINUE DO 130 I=1,N KFA=IABS(K(I,2)) IF(K(I,1).LE.10.AND.K(I,3).GT.0.AND. $ ((KFA.GT.10.AND.KFA.LE.20).OR.KFA.EQ.22)) 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 130 DO 120 J=1,4 DPS(J)=DPS(J)+P(I,J) 120 CONTINUE 130 CONTINUE CALL LUDBRB(0,0,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4), $ -DPS(3)/DPS(4)) PECM=0. DO 140 I=1,N IF(K(I,1).GE.1.AND.K(I,1).LE.10) PECM=PECM+P(I,4) 140 CONTINUE C...Reserve copy of particles at end of record. And mark entries C...comming from W decays IBE=N+MSTU(3) IBE1=IBE+1 DO 170 I=NSAV+1,N IF (ABS(K(I,2)).NE.211.AND.K(I,2).NE.111) GOTO 170 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 170 IF(IBE.GE.MSTU(4)-MSTU(32)-5) THEN CALL LUERRM(11,'(LUBEC2:) no more memory left in LUJETS') RETURN ENDIF IBE=IBE+1 K(IBE,1)=I IF (MSTJ(53).NE.0) THEN K(IBE,5)=0 IM=I 150 IF (K(IM,3).GT.0) THEN IM=K(IM,3) IF (ABS(K(IM,2)).NE.24) GOTO 150 K(IBE,5)=K(IM,2) ENDIF ENDIF DO 160 J=1,3 P(IBE,J)=0. 160 CONTINUE 170 CONTINUE C...Loop through particle pairs and find old relative momentum. Use C...Newtons method to find new relative momentum. Treat or do not treat C...pairs in e+e- -> W+W- according to MSTJ(53) DO 210 I1M=IBE1,IBE I1=K(I1M,1) IM1=K(I1,3) KFHQ=ABS(KLU(IM1,13)) IF (KFHQ.GT.3) KFHQ=0 DO 200 I2M=I1M+1,IBE I2=K(I2M,1) 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 IF (MSTJ(53).EQ.-1.AND.K(I1M,5).NE.K(I2M,5).AND. $ K(I1,2).NE.K(I2,2)) GOTO 200 IF (MSTJ(51).EQ.4.AND.IM1.EQ.K(I2,3).AND.KFHQ.GT.0) GOTO 200 Q2OLD=MAX(0.0,(P(I1,4)+P(I2,4))**2-(P(I1,3)+P(I2,3))**2- $ (P(I1,2)+P(I2,2))**2-(P(I1,1)+P(I2,1))**2- $ (P(I1,5)+P(I2,5))**2) QOLD=SQRT(Q2OLD) if (k(i1,2).eq.k(i2,2)) then call gfill1(1,QOLD,1.0) else call gfill1(4,QOLD,1.0) endif QNEW=QOLD SINT0=0.5*(QOLD*SQRT(Q2OLD+SM4)+ $ SM4*(LOG(SQRT(SM4)/(QOLD+SQRT(Q2OLD+SM4))))) ILOOP=0 180 ILOOP=ILOOP+1 IQ=MIN(INT(QNEW*1000.0/QMAX)+1,1000) XIQU=MIN((QNEW-REAL(IQ-1)*QD)/QD,1.0) XIQL=MAX((REAL(IQ)*QD-QNEW)/QD,0.0) IF (K(I1,2).EQ.K(I2,2)) THEN IF (QNEW.LT.MIN(SIG,SQRT(SM4))*0.1) THEN QDEL=FISP(QNEW) ELSE QDEL=FPIP(QNEW)-XIQL*SIEQ(IQ-1)-XIQU*SIEQ(IQ) ENDIF QDEL=(QDEL-SINT0)/FP(QNEW) ELSE IF (QNEW.LT.MIN(SIG,SQRT(SM4))*0.1) THEN QDEL=FISM(QNEW) ELSE QDEL=FPIM(QNEW)-XIQL*SINE(IQ-1)-XIQU*SINE(IQ) ENDIF QDEL=(QDEL-SINT0)/FM(QNEW) ENDIF QNEW=QNEW-QDEL IF (ILOOP.GT.100) THEN CALL LUERRM(14,'(LUBOE2:) caught in infinite loop') ELSE IF (ABS(QDEL).GT.0.001) GOTO 180 ENDIF Q2NEW=QNEW**2 C...Calculate and save shift to be performed on three-momenta. 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 190 J=1,3 PD=HA*(P(I2,J)-P(I1,J)) P(I1M,J)=P(I1M,J)+PD P(I2M,J)=P(I2M,J)-PD 190 CONTINUE 200 CONTINUE 210 CONTINUE C...Shift momenta and recalculate energies. DO 230 IM=IBE1,IBE I=K(IM,1) DO 220 J=1,3 P(I,J)=P(I,J)+P(IM,J) 220 CONTINUE P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2) 230 CONTINUE C...Rescale all momenta for energy conservation. PES=0.0 PQS=0.0 DO 240 I=1,N IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 240 PES=PES+P(I,4) PQS=PQS+P(I,5)**2/P(I,4) 240 CONTINUE PARJ(95)=PES-PECM FAC=(PECM-PQS)/(PES-PQS) DO 260 I=1,N IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 260 DO 250 J=1,3 P(I,J)=FAC*P(I,J) 250 CONTINUE P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2) 260 CONTINUE C...Boost back to correct reference frame. CALL LUDBRB(0,0,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),DPS(3)/DPS(4)) DO 270 I=1,N IF(K(I,1).LT.0) K(I,1)=-K(I,1) 270 CONTINUE RETURN END