PYTHIA/jetset/luboei.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUBOEI(NSAV)
   5
   6 C...Purpose: to modify event so as to approximately take into account
   7 C...Bose-Einstein effects according to a simple phenomenological
   8 C...parametrization.
   9       IMPLICIT DOUBLE PRECISION(D)
  10       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
  11       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  12       SAVE /LUJETS/,/LUDAT1/
  13       DIMENSION DPS(4),KFBE(9),NBE(0:9),BEI(100)
  14       DATA KFBE/211,-211,111,321,-321,130,310,221,331/
  15
  16 C...Boost event to overall CM frame. Calculate CM energy.
  17       IF((MSTJ(51).NE.1.AND.MSTJ(51).NE.2).OR.N-NSAV.LE.1) RETURN
  18       DO 100 J=1,4
  19       DPS(J)=0.
  20   100 CONTINUE
  21       DO 120 I=1,N
  22       KFA=IABS(K(I,2))
  23       IF(K(I,1).LE.10.AND.((KFA.GT.10.AND.KFA.LE.20).OR.KFA.EQ.22).AND.
  24      &K(I,3).GT.0) THEN
  25         KFMA=IABS(K(K(I,3),2))
  26         IF(KFMA.GT.10.AND.KFMA.LE.80) K(I,1)=-K(I,1)
  27       ENDIF
  28       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120
  29       DO 110 J=1,4
  30       DPS(J)=DPS(J)+P(I,J)
  31   110 CONTINUE
  32   120 CONTINUE
  33       CALL LUDBRB(0,0,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4),
  34      &-DPS(3)/DPS(4))
  35       PECM=0.
  36       DO 130 I=1,N
  37       IF(K(I,1).GE.1.AND.K(I,1).LE.10) PECM=PECM+P(I,4)
  38   130 CONTINUE
  39
  40 C...Reserve copy of particles by species at end of record.
  41       NBE(0)=N+MSTU(3)
  42       DO 160 IBE=1,MIN(9,MSTJ(52))
  43       NBE(IBE)=NBE(IBE-1)
  44       DO 150 I=NSAV+1,N
  45       IF(K(I,2).NE.KFBE(IBE)) GOTO 150
  46       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 150
  47       IF(NBE(IBE).GE.MSTU(4)-MSTU(32)-5) THEN
  48         CALL LUERRM(11,'(LUBOEI:) no more memory left in LUJETS')
  49         RETURN
  50       ENDIF
  51       NBE(IBE)=NBE(IBE)+1
  52       K(NBE(IBE),1)=I
  53       DO 140 J=1,3
  54       P(NBE(IBE),J)=0.
  55   140 CONTINUE
  56   150 CONTINUE
  57   160 CONTINUE
  58       IF(NBE(MIN(9,MSTJ(52)))-NBE(0).LE.1) GOTO 280
  59
  60 C...Tabulate integral for subsequent momentum shift.
  61       DO 220 IBE=1,MIN(9,MSTJ(52))
  62       IF(IBE.NE.1.AND.IBE.NE.4.AND.IBE.LE.7) GOTO 180
  63       IF(IBE.EQ.1.AND.MAX(NBE(1)-NBE(0),NBE(2)-NBE(1),NBE(3)-NBE(2))
  64      &.LE.1) GOTO 180
  65       IF(IBE.EQ.4.AND.MAX(NBE(4)-NBE(3),NBE(5)-NBE(4),NBE(6)-NBE(5),
  66      &NBE(7)-NBE(6)).LE.1) GOTO 180
  67       IF(IBE.GE.8.AND.NBE(IBE)-NBE(IBE-1).LE.1) GOTO 180
  68       IF(IBE.EQ.1) PMHQ=2.*ULMASS(211)
  69       IF(IBE.EQ.4) PMHQ=2.*ULMASS(321)
  70       IF(IBE.EQ.8) PMHQ=2.*ULMASS(221)
  71       IF(IBE.EQ.9) PMHQ=2.*ULMASS(331)
  72       QDEL=0.1*MIN(PMHQ,PARJ(93))
  73       IF(MSTJ(51).EQ.1) THEN
  74         NBIN=MIN(100,NINT(9.*PARJ(93)/QDEL))
  75         BEEX=EXP(0.5*QDEL/PARJ(93))
  76         BERT=EXP(-QDEL/PARJ(93))
  77       ELSE
  78         NBIN=MIN(100,NINT(3.*PARJ(93)/QDEL))
  79       ENDIF
  80       DO 170 IBIN=1,NBIN
  81       QBIN=QDEL*(IBIN-0.5)
  82       BEI(IBIN)=QDEL*(QBIN**2+QDEL**2/12.)/SQRT(QBIN**2+PMHQ**2)
  83       IF(MSTJ(51).EQ.1) THEN
  84         BEEX=BEEX*BERT
  85         BEI(IBIN)=BEI(IBIN)*BEEX
  86       ELSE
  87         BEI(IBIN)=BEI(IBIN)*EXP(-(QBIN/PARJ(93))**2)
  88       ENDIF
  89       IF(IBIN.GE.2) BEI(IBIN)=BEI(IBIN)+BEI(IBIN-1)
  90   170 CONTINUE
  91
  92 C...Loop through particle pairs and find old relative momentum.
  93   180 DO 210 I1M=NBE(IBE-1)+1,NBE(IBE)-1
  94       I1=K(I1M,1)
  95       DO 200 I2M=I1M+1,NBE(IBE)
  96       I2=K(I2M,1)
  97       Q2OLD=MAX(0.,(P(I1,4)+P(I2,4))**2-(P(I1,1)+P(I2,1))**2-(P(I1,2)+
  98      &P(I2,2))**2-(P(I1,3)+P(I2,3))**2-(P(I1,5)+P(I2,5))**2)
  99       QOLD=SQRT(Q2OLD)
 100
 101 C...Calculate new relative momentum.
 102       IF(QOLD.LT.1E-3*QDEL) THEN
 103         GOTO 200
 104       ELSEIF(QOLD.LE.QDEL) THEN
 105         QMOV=QOLD/3.
 106       ELSEIF(QOLD.LT.(NBIN-0.1)*QDEL) THEN
 107         RBIN=QOLD/QDEL
 108         IBIN=RBIN
 109         RINP=(RBIN**3-IBIN**3)/(3*IBIN*(IBIN+1)+1)
 110         QMOV=(BEI(IBIN)+RINP*(BEI(IBIN+1)-BEI(IBIN)))*
 111      &  SQRT(Q2OLD+PMHQ**2)/Q2OLD
 112       ELSE
 113         QMOV=BEI(NBIN)*SQRT(Q2OLD+PMHQ**2)/Q2OLD
 114       ENDIF
 115       Q2NEW=Q2OLD*(QOLD/(QOLD+3.*PARJ(92)*QMOV))**(2./3.)
 116
 117 C...Calculate and save shift to be performed on three-momenta.
 118       HC1=(P(I1,4)+P(I2,4))**2-(Q2OLD-Q2NEW)
 119       HC2=(Q2OLD-Q2NEW)*(P(I1,4)-P(I2,4))**2
 120       HA=0.5*(1.-SQRT(HC1*Q2NEW/(HC1*Q2OLD-HC2)))
 121       DO 190 J=1,3
 122       PD=HA*(P(I2,J)-P(I1,J))
 123       P(I1M,J)=P(I1M,J)+PD
 124       P(I2M,J)=P(I2M,J)-PD
 125   190 CONTINUE
 126   200 CONTINUE
 127   210 CONTINUE
 128   220 CONTINUE
 129
 130 C...Shift momenta and recalculate energies.
 131       DO 240 IM=NBE(0)+1,NBE(MIN(9,MSTJ(52)))
 132       I=K(IM,1)
 133       DO 230 J=1,3
 134       P(I,J)=P(I,J)+P(IM,J)
 135   230 CONTINUE
 136       P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
 137   240 CONTINUE
 138
 139 C...Rescale all momenta for energy conservation.
 140       PES=0.
 141       PQS=0.
 142       DO 250 I=1,N
 143       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 250
 144       PES=PES+P(I,4)
 145       PQS=PQS+P(I,5)**2/P(I,4)
 146   250 CONTINUE
 147       FAC=(PECM-PQS)/(PES-PQS)
 148       DO 270 I=1,N
 149       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 270
 150       DO 260 J=1,3
 151       P(I,J)=FAC*P(I,J)
 152   260 CONTINUE
 153       P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
 154   270 CONTINUE
 155
 156 C...Boost back to correct reference frame.
 157   280 CALL LUDBRB(0,0,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),DPS(3)/DPS(4))
 158       DO 290 I=1,N
 159       IF(K(I,1).LT.0) K(I,1)=-K(I,1)
 160   290 CONTINUE
 161
 162       RETURN
 163       END