This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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