GEANT321/fluka/verein.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:03  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.44  by  S.Giani
  11 *-- Author :
  12 *$ CREATE VEREIN.FOR
  13 *COPY VEREIN
  14 *
  15 *=== verein ===========================================================*
  16 *
  17       SUBROUTINE VEREIN(IT,LA,LT,RER,REL,RPXR,RPYR,RPZR,RPXL,RPYL,RPZL,
  18      *KR1R,KR2R,KR1L,KR2L,IHAD,LL,KFR1,KFR2,IMPS,IMVE,IB08,IA08,
  19      *IB10,IA10,B3,AS,B8,IAR,KFA1,KFA2,KFA3,KFA4,IOPT)
  20
  21 #include "geant321/dblprc.inc"
  22 #include "geant321/dimpar.inc"
  23 #include "geant321/iounit.inc"
  24 *
  25 *----------------------------------------------------------------------*
  26 *  Verein89: slight revision by A. Ferrari                             *
  27 *----------------------------------------------------------------------*
  28 *
  29 #include "geant321/finpar2.inc"
  30 #include "geant321/part.inc"
  31       DIMENSION KFR1(*),KFR2(*),IMPS(6,6)
  32 *
  33       PARAMETER (KMXJCM = 100)
  34       DIMENSION IV(KMXJCM)
  35       DIMENSION IMVE(6,6),IB08(6,21),IA08(6,21),IB10(6,21),IA10(6,21)
  36       REAL RNDM(1)
  37 C*****VEREIN COMBINES THE TWO JETS INTO A COMPLETE E+E- EVENT OR INTO
  38 C*****A COMPLETE QQQ (AQ,AQ,AQ) EVENT
  39 C*****RER,REL,RPXR,RPXL,RPYR,RPYL,RPZR,RPZL ARE REST JET ENERGIES AND
  40 C*****MOMENTA , KR1R,KR1L,KR2R,KR2L ARE REST JET FLAVOURS
  41 *
  42 *     R is for right, L for left!
  43 *
  44 *  Initialize Iv
  45 *
  46       DATA IV /KMXJCM*0/
  47       IF(LT.EQ.0) GO TO 300
  48       WRITE(LUNOUT,202)IT
  49   202 FORMAT(1H0,I5,3H=IT)
  50   300 CONTINUE
  51 C*****ONLY ONE RESONANCE WILL BE CREATED IF IT=0
  52       IF(IT.EQ.0) GO TO 100
  53       IF(KR2R.EQ.0.OR.KR1R.EQ.0)GO TO 10
  54       IF(KR1L.EQ.0.OR.KR2L.EQ.0)GO TO 10
  55 C*****RESTJET CONTAINS FOUR QUARKS,THEREFORE TWO MESONS WILL BE FORMED
  56 *
  57 *  Reg is et equal to the sum of the left and right jet energy rest,
  58 *  the same for momenta
  59 *
  60       REG=RER+REL
  61       RPXG=RPXR+RPXL
  62       RPYG=RPYR+RPYL
  63       RPZG=RPZR+RPZL
  64       J=IT+1
  65       I=IT+2
  66       IT=J
  67 *
  68 *  Hklass uses always only Iv(i) with i >= it, so from the previous
  69 *  instructions, it is always > than the original it , so the used
  70 *  Iv's should be always 0 also in the bamjet original array ?????
  71 *
  72       CALL HKLASS(IT,LT,LA,LL,KFR1,KFR2,KR1R,KR2R,KR1L,KR2L,IV,IMPS,
  73      &IMVE,IB08,IA08,IB10,IA10,AS,B8,KFA1,KFA2,KFA3,KFA4,IOPT)
  74       IHAD=IHAD+2
  75       IF(LT.EQ.1) WRITE(LUNOUT,301)IHAD
  76   301 FORMAT(1H0,13HHADRONANZAHL=,I5)
  77       GO TO 20
  78 C*****RESTJET CONTAINS TWO OR THREE QUARKS,THEREFORE ONE MESON OR ONE
  79 C*****BARYON OR ONE ANTIBARYON  WILL BE FORMED
  80    10 CONTINUE
  81       IAK=IT
  82       CALL GRNDM(RNDM,1)
  83       X=RNDM(1)
  84       IF(X.LE.0.5D0.AND.IAR.NE.0) IAK=IAR
  85       REG=RER+REL+HEF(IAK)
  86       RPXG=RPXR+RPXL+PXF(IAK)
  87       RPYG=RPYR+RPYL+PYF(IAK)
  88       RPZG=RPZR+RPZL+PZF(IAK)
  89       I=IT+1
  90       J=IAK
  91       IT=I
  92 *
  93 *  Hklass uses always only Iv(i) with i >= it, so from the previous
  94 *  instructions, it is always > than the original it , so the used
  95 *  Iv's should be always 0 also in the bamjet original array ?????
  96 *
  97       CALL HKLASS(IT,LT,LA,LL,KFR1,KFR2,KR1R,KR2R,KR1L,KR2L,IV,IMPS,
  98      &IMVE,IB08,IA08,IB10,IA10,AS,B8,KFA1,KFA2,KFA3,KFA4,IOPT)
  99 *  Hp1 is not used, also Hp2 !!!
 100 *     HP2=HP1
 101       IHAD=IHAD+1
 102       IF(LT.EQ.1) WRITE(LUNOUT,301)IHAD
 103    20 CONTINUE
 104       RPG=SQRT(RPXG**2+RPYG**2+RPZG**2)
 105       IF(RPG.GT.REG) LA=3
 106       IF(LA.EQ.3.AND.LT.GT.0) WRITE(LUNOUT,71)
 107    71 FORMAT(2X,'REST JET MOMENTUM IS GREATER THAN REST JET ENERGY')
 108       IF(LA.EQ.3) RETURN
 109       RMG=SQRT((REG-RPG)*(REG+RPG))
 110       HM1=AMF(J)
 111       HM2=AMF(I)
 112       IF(RMG.LT.HM1+HM2) LA=3
 113       IF(LA.EQ.3.AND.LT.GT.0) WRITE(LUNOUT,71)
 114       IF(LA.EQ.3) RETURN
 115       IF(LT.EQ.0) GO TO 30
 116       WRITE(LUNOUT,60)REG,RMG,RPG,HM1,HM2,J,I
 117    60 FORMAT(1H0,19HEG,MG,PG,M1,M2,J,I=,5F8.4,2I3)
 118 C     LORENZ PARAMETERS
 119    30 CONTINUE
 120       GAA=REG/RMG
 121       GABE=RPG/RMG
 122 C     DECAY INTO 2 HADRONS IN THE CMS OF THE REMAINDER
 123       HE1=(RMG**2+HM1**2-HM2**2)/(2.D0*RMG)
 124       HE2=RMG-HE1
 125 *   Hp1 is not used !
 126 *     HP1=SQRT(HE1**2-HM1**2)
 127 C     SAMPLES THE MOMENTUM DIRECTIONS OF THE LAST PARTICLES
 128       HE=HE1
 129       HMA=HM1
 130       CALL FKIMPU(HE,HMA,HPS,HPX,HPY,HPZ,LT,LL,B3)
 131       HP1X=HPX
 132       HP1Y=HPY
 133       HP1Z=HPZ
 134       HP2X=-HP1X
 135       HP2Y=-HP1Y
 136       IF (IOPT.EQ.4.AND.HM2.GT.HM1)  GO TO 999
 137       HP1Z=-HPZ
 138       HP2Z=HPZ
 139   999 CONTINUE
 140       HP2Z=-HP1Z
 141       IF(RPG.EQ.0.D0)GO TO 80
 142 C     ROTATION BACK TO THE CMS*
 143       HEX=HE1*GAA+HP1Z*GABE
 144       HEY=HE2*GAA+HP2Z*GABE
 145       HP1Z=HE1*GABE+HP1Z*GAA
 146       HP2Z=HE2*GABE+HP2Z*GAA
 147       HE1=HEX
 148       HE2=HEY
 149 C     ROTATION INTO THE CMS OF THE E+ E-  COLLISION
 150       X=HP1X
 151       Y=HP1Y
 152       Z=HP1Z
 153       COTE=RPZG/RPG
 154       SITE=SQRT((1.D0-COTE)*(1.D0+COTE))
 155       IF(SITE.EQ.0.D0) GO TO 11
 156       COV=RPXG/(RPG*SITE)
 157       SIV=RPYG/(RPG*SITE)
 158       GO TO 12
 159    11 CONTINUE
 160       SIV=1.D0
 161       COV=0.D0
 162    12 CONTINUE
 163       COPS=-SIV
 164       SIPS=COV
 165       CALL DRELAB(X,Y,Z,COTE,SITE,COPS,SIPS)
 166       HP1X=X
 167       HP1Y=Y
 168       HP1Z=Z
 169       X=HP2X
 170       Y=HP2Y
 171       Z=HP2Z
 172       CALL DRELAB(X,Y,Z,COTE,SITE,COPS,SIPS)
 173       HP2X=X
 174       HP2Y=Y
 175       HP2Z=Z
 176    80 CONTINUE
 177       HEF(J)=HE1
 178       PXF(J)=HP1X
 179       PYF(J)=HP1Y
 180       PZF(J)=HP1Z
 181       HEF(I)=HE2
 182       PXF(I)=HP2X
 183       PYF(I)=HP2Y
 184       PZF(I)=HP2Z
 185       IF(LT.EQ.0)GO TO 13
 186       WRITE(LUNOUT,50)HEF(J),PXF(J),PYF(J),PZF(J)
 187       WRITE(LUNOUT,50)HEF(I),PXF(I),PYF(I),PZF(I)
 188    50 FORMAT(1H0,11HE,PX,PY,PZ=,4F8.4)
 189    13 CONTINUE
 190       GO TO 200
 191   100 CONTINUE
 192 C*****ONLY ONE RESONANCE WILL BE CREATED IF IT=0
 193       IF(LT.EQ.0)GO TO 14
 194       WRITE(LUNOUT,70)
 195    70 FORMAT(1H0,'CUT OFF OF THE LEFT AND RIGHT JET BEFOR THE FIRST
 196      *DECAY STEP',/,'IF ALLOWED ONLY ONE RESONANCE COULD BE CREATED')
 197    14 CONTINUE
 198       IHAD=1
 199 *
 200 *  La = 2 means resonance creation
 201 *
 202       LA=2
 203   200 CONTINUE
 204       RETURN
 205       END