HIJING/hipyset1_35/luboei_hijing.F

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