[u/mrichter/AliRoot.git] / GEANT321 / fluka / bamjev.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:19:53  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.41  by  S.Giani
*-- Author :
*$ CREATE BAMJEV.FOR
*COPY BAMJEV
*
*=== bamjev ===========================================================*
*
      SUBROUTINE BAMJEV ( IHAD, KFA1, KFA2, KFA3, KFA4, AE0, IOPT )
 
#include "geant321/dblprc.inc"
#include "geant321/dimpar.inc"
#include "geant321/iounit.inc"
*
*----------------------------------------------------------------------*
*  Bamjet90: slight revision by A. Ferrari                             *
*----------------------------------------------------------------------*
*
*----------------------------------------------------------------------*
*     Ihad   = number of final hadrons and hadron resonances           *
*     Ae0    = initial energy in GeV                                   *
*     Kfai   = initial quark flavours (u=1,d=2,s=3,c=4,ubar=7,dbar=8,  *
*                                      sbar=9, cbar=10)                *
*     Iopt   = 1,2,3,4,5 means:                                        *
*              1: single (anti)quark jet,   (kfa1)                     *
*              2: single (anti)diquark jet, (kfa1-kfa2)                *
*              3: complete quark antiquark twojet event, (kfa1,kfa2)   *
*              4: complete (anti)quark-(anti)diquark two jet event,    *
*                 (kfa1,kfa2-kfa3)                                     *
*              5: complete diquark-(anti)diquark two jet event,        *
*                 (kfa1-kfa2,kfa3-kfa4)                                *
*     Common/finpar/ contains the momenta,energies and quantum numbers *
*     of the created hadrons                                           *
*     Iv     = actual vertex,iv=1,4,5,6,9,10 are meson verteces        *
*              iv=2,3,7,8 are baryon verteces                          *
*     La     = 1 means cut-off                                         *
*     Ll     = 0,1 means quark jet, antiquark jet, (diquark jet, anti- *
*              diquark jet)                                            *
*     Common/remain/ contains rest jet energy,momenta and quantumnum-  *
*                    bers                                              *
*----------------------------------------------------------------------*
*
#include "geant321/bamjcm.inc"
#include "geant321/finpar2.inc"
#include "geant321/part.inc"
#include "geant321/inpdat.inc"
 
      COMMON/FKREMA/ RPXR,RPYR,RPZR,RER,KR1R,KR2R
      SAVE NCOU, ITMX
*
      DATA NCOU/0/
      DATA ITMX/0/
      NCOU=NCOU+1
      A1SAVE = A1
      B1SAVE = B1
      B2SAVE = B2
      B3SAVE = B3
      IF ( AE0 .LE. 4.D+00 )THEN
         A1 = 0.88D+00
      ELSE IF ( AE0 .LE. 8.0D+00 ) THEN
         A1 = 0.88D+00
      ELSE IF ( AE0 .LT. 30.D+00 ) THEN
         A1 = 0.88D+00
      ELSE
         A1 = 0.88D+00
      END IF
         FRB12 = 0.5D+00
         B1 = 4.D+00 + 1.D+00 / ( FRB12 * AE0 )**2
         B2 = 4.D+00 + 1.D+00 / ( FRB12 * AE0 )**2
      E00 = 40.D+00
* The following is consistent with B3=6
      B3 = B3 * LOG10 ( E00 ) /
     &   ( LOG10 ( 1.D+00 + ( AE0 / E00 )**2 ) + LOG10 ( E00 ) )
C     IF (NCOU.EQ.4701)LT=1
*or   IF (LT.EQ.1)WRITE(LUNOUT,3399)IOPT,NCOU
 3399 FORMAT('  BAMJET',2I10 )
*or   IF (LT.EQ.1)WRITE(LUNOUT,288)IHAD,KFA1,KFA2,KFA3,KFA4,AE0,IOPT
  288 FORMAT (5I5,2E12.4,' BAMJET,IHAD,KFA1,KFA2,KFA3,KFA4,AE0,IOPT')
   63 CONTINUE
      DO 62  I=1,ITMX
         KFR1(I) = 0
         KFR2(I) = 0
   62 CONTINUE
      ITMX=0
   60 CONTINUE
      IYY=0
      IHAD=0
      IT=0
      E0=AE0*.5D0
      IF(IOPT.EQ.1.OR.IOPT.EQ.2) E0=AE0
      LL=0
      IF(KFA1.GT.6.AND.IOPT.EQ.1) LL=1
      IF(KFA1.LE.6.AND.IOPT.EQ.2) LL=1
      IF(KFA1.GT.6.AND.IOPT.EQ.4) LL=1
      PGX = 0.D0
      PGY = 0.D0
      PGZ = 0.D0
      RPX0 = 0.D0
      RPY0 = 0.D0
* The following 6 initializations might be useless, but they make the
* code much clearer
      KR1R = 0
      KR2R = 0
      KR1L = 0
      KR2L = 0
      RER  = 0.D+00
      REL  = 0.D+00
      DO 10 I=1,KMXJCM
         KFR1(I) = 0
         KFR2(I) = 0
         LA = 0
         IT = IT+1
         J  = IT-1
* The following line seems useless
*        K  = IT+1
   40    CONTINUE
C*****CUT OFF TASK
*  | Abbrch is called to cut the chain
         CALL ABBRCH(IT,LL,LA,LT,E0,PGX,PGY,PGZ,KFR1,KFR2,RE,
     &   KR1R,KR2R,KR1L,KR2L,RPX,RPY,RPZ,RPXR,RPYR,RPZR,RPXL,RPYL,RPZL,
     &   RER,REL,IV,B1,B2,KFA1,KFA2,KFA3,KFA4,IOPT,IYY)
         ITMX=MAX(ITMX,IT)
         IF(LA .EQ. 0) GO TO 20
         IT=IT-1
         IF(IOPT.EQ.3.AND.LL.EQ.0) GO TO 70
         IF(IOPT.EQ.4.AND.KFA1.LE.6.AND.LL.EQ.0) GO TO 70
         IF(IOPT.EQ.4.AND.KFA1.GT.6.AND.LL.EQ.1) GO TO 70
         IF(IOPT.EQ.5.AND.LL.EQ.0) GO TO 70
         GO TO 50
   70    CONTINUE
         IYY = 1
         LL  = 1
         IF(IOPT.EQ.4.AND.KFA1.GT.6) LL = 0
         IAR=IT
         GO TO 30
   20    CONTINUE
C*****CHOICE OF THE VERTEX
         CALL FKVERT(IT,LT,LL,KFA1,E0,IV,RE,KFR1,KFR2,AME,IOPT)
C*****CHOICE OF THE FLAVOUR
         CALL FKFLAV(IT,LT,LL,E0,IV,RE,KFR1,KFR2,ISU,BET,KFA1,KFA2,
     &               KFA3,KFA4,IOPT)
C*****CLASSIFICATION OF THE PARTICLES
         CALL HKLASS(IT,LT,LA,LL,KFR1,KFR2,KR1R,KR2R,KR1L,KR2L,IV,IMPS,
     &         IMVE,IB08,IA08,IB10,IA10,AS,B8,KFA1,KFA2,KFA3,KFA4,IOPT)
         ITMX=MAX(ITMX,IT)
         IF (IT .EQ. 1) RX = E0
         IF (IT .GT. 1) RX = RE(J)
         IF(AMF(IT) .GT. RX) GO TO 63
         IF(AMF(IT) .LE. RX) GO TO 19
         LA = 1
         GO TO 40
   19    CONTINUE
         IHAD = IHAD + 1
*or         IF(LT .EQ. 0) GO TO 31
*or         WRITE(LUNOUT,32)IHAD
*or   31    CONTINUE
C*****CHOICE OF THE ENERGY
         CALL ENERGI(IT,LL,LT,IV,RE,HMA,HE,E0,A1)
C*****CHOICE OF THE MOMENTUM
*  |
*  |  He is the total energy, hma the mass one (input) hpx, hpy, hpz
*  |  the momentum components (output values), hps the transversal
*  |  momentum (output)
*  |
         CALL FKIMPU(HE,HMA,HPS,HPX,HPY,HPZ,LT,LL,B3)
         IF (IT .GT. 1) GO TO 13
         RPX(IT)=RPX0-HPX
         RPY(IT)=RPY0-HPY
         GO TO 14
   13    RPX(IT)=RPX(J)-HPX
         RPY(IT)=RPY(J)-HPY
   14    CONTINUE
         IF (IOPT.EQ.1.AND.LL.EQ.1)HPZ=-HPZ
         IF(IOPT.EQ.2.AND.LL.EQ.1) HPZ=-HPZ
         IF(IOPT.EQ.4.AND.KFA1.GT.6) HPZ=-HPZ
         IF(IOPT.EQ.5) HPZ=-HPZ
         PGX=PGX+HPX
         PGY=PGY+HPY
         PGZ=PGZ+HPZ
         PXF(IT)=HPX
         PYF(IT)=HPY
         PZF(IT)=HPZ
*or         IF (LT .EQ. 0) GO TO 15
*or            WRITE(LUNOUT,16)PGX,PGY,PGZ
*or   16       FORMAT(1H0,12HPGX,PGY,PGZ=,3F8.4)
*or   15    CONTINUE
   30    CONTINUE
 
   10 CONTINUE
*
* we suppose that exiting from loop must be achieved via " go to 50
*
      WRITE (LUNERR,*)' BAMJEV: EXITING FROM LOOP ABNORMALLY!!!! '
      WRITE (LUNOUT,*)' BAMJEV: EXITING FROM LOOP ABNORMALLY!!!! '
   50 CONTINUE
      ITMX=MAX(ITMX,IT)
      IF(IOPT.EQ.1.OR.IOPT.EQ.2) GO TO 51
C*****PUT THE RIGHT AND LEFT JET TOGETHER
      CALL VEREIN(IT,LA,LT,RER,REL,RPXR,RPYR,RPZR,RPXL,RPYL,RPZL,
     &KR1R,KR2R,KR1L,KR2L,IHAD,LL,KFR1,KFR2,IMPS,IMVE,IB08,IA08,
     &IB10,IA10,B3,AS,B8,IAR,KFA1,KFA2,KFA3,KFA4,IOPT)
      IF(LA.EQ.3) GO TO 63
      IF(LA.EQ.2) GO TO 63
   51 CONTINUE
      IF(IOPT.EQ.3.OR.IOPT.EQ.4.OR.IOPT.EQ.5) GO TO 52
      IF(LL.EQ.0) GO TO 52
      RPXR=RPXL
      RPYR=RPYL
      RPZR=RPZL
      RER=REL
      KR1R=KR1L
      KR2R=KR2L
   52 CONTINUE
      IF(LE.EQ.0) GO TO 1000
      WRITE(LUNOUT,92)
   92 FORMAT(2X,'NF,NAME,MASS,IQ,IB,PX,PY,PZ,E')
      DO 91 J=1,IT
         WRITE(LUNOUT,90)NREF(J),ANF(J),AMF(J),ICHF(J),IBARF(J),PXF(J),
     &   PYF(J),PZF(J),HEF(J)
   90    FORMAT(2X,I3,A6,F6.3,2I4,4F8.4)
   91 CONTINUE
 1000 CONTINUE
   64 FORMAT(1H0,38HNUMBER OF EVENTS WITH PREST GT. EREST=,I4,
     */,21HNUMBER OF ALL EVENTS=,I4)
 2000 FORMAT(1H0,'NUMBER OF EVENTS WITH ONLY ONE PARTICLE=',I4)
*or      IF(LT.EQ.0) GO TO 17
*or      WRITE(LUNOUT,18)IHAD
*or   18 FORMAT(1H0,15HMULTIPLIZITAET=,I3)
*or   32 FORMAT(1H0,13HHADRONANZAHL=,I3)
*or   17 CONTINUE
      A1 = A1SAVE
      B1 = B1SAVE
      B2 = B2SAVE
      B3 = B3SAVE
      RETURN
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:19:53 cernlib
	6	* Geant
	7	*
	8	*
	9	#include "geant321/pilot.h"
	10	*CMZ : 3.21/02 29/03/94 15.41.41 by S.Giani
	11	*-- Author :
	12	*$ CREATE BAMJEV.FOR
	13	*COPY BAMJEV
	14	*
	15	=== bamjev ===========================================================
	16	*
	17	SUBROUTINE BAMJEV ( IHAD, KFA1, KFA2, KFA3, KFA4, AE0, IOPT )
	18
	19	#include "geant321/dblprc.inc"
	20	#include "geant321/dimpar.inc"
	21	#include "geant321/iounit.inc"
	22	*
	23	----------------------------------------------------------------------
	24	* Bamjet90: slight revision by A. Ferrari *
	25	----------------------------------------------------------------------
	26	*
	27	----------------------------------------------------------------------
	28	* Ihad = number of final hadrons and hadron resonances *
	29	* Ae0 = initial energy in GeV *
	30	* Kfai = initial quark flavours (u=1,d=2,s=3,c=4,ubar=7,dbar=8, *
	31	* sbar=9, cbar=10) *
	32	* Iopt = 1,2,3,4,5 means: *
	33	* 1: single (anti)quark jet, (kfa1) *
	34	* 2: single (anti)diquark jet, (kfa1-kfa2) *
	35	* 3: complete quark antiquark twojet event, (kfa1,kfa2) *
	36	* 4: complete (anti)quark-(anti)diquark two jet event, *
	37	* (kfa1,kfa2-kfa3) *
	38	* 5: complete diquark-(anti)diquark two jet event, *
	39	* (kfa1-kfa2,kfa3-kfa4) *
	40	* Common/finpar/ contains the momenta,energies and quantum numbers *
	41	* of the created hadrons *
	42	* Iv = actual vertex,iv=1,4,5,6,9,10 are meson verteces *
	43	* iv=2,3,7,8 are baryon verteces *
	44	* La = 1 means cut-off *
	45	* Ll = 0,1 means quark jet, antiquark jet, (diquark jet, anti- *
	46	* diquark jet) *
	47	* Common/remain/ contains rest jet energy,momenta and quantumnum- *
	48	* bers *
	49	----------------------------------------------------------------------
	50	*
	51	#include "geant321/bamjcm.inc"
	52	#include "geant321/finpar2.inc"
	53	#include "geant321/part.inc"
	54	#include "geant321/inpdat.inc"
	55
	56	COMMON/FKREMA/ RPXR,RPYR,RPZR,RER,KR1R,KR2R
	57	SAVE NCOU, ITMX
	58	*
	59	DATA NCOU/0/
	60	DATA ITMX/0/
	61	NCOU=NCOU+1
	62	A1SAVE = A1
	63	B1SAVE = B1
	64	B2SAVE = B2
65	B3SAVE = B3
66	IF ( AE0 .LE. 4.D+00 )THEN
67	A1 = 0.88D+00
68	ELSE IF ( AE0 .LE. 8.0D+00 ) THEN
69	A1 = 0.88D+00
70	ELSE IF ( AE0 .LT. 30.D+00 ) THEN
71	A1 = 0.88D+00
72	ELSE
73	A1 = 0.88D+00
74	END IF
75	FRB12 = 0.5D+00
76	B1 = 4.D+00 + 1.D+00 / ( FRB12 * AE0 )**2
77	B2 = 4.D+00 + 1.D+00 / ( FRB12 * AE0 )**2
78	E00 = 40.D+00
79	* The following is consistent with B3=6
80	B3 = B3 * LOG10 ( E00 ) /
81	& ( LOG10 ( 1.D+00 + ( AE0 / E00 )**2 ) + LOG10 ( E00 ) )
82	C IF (NCOU.EQ.4701)LT=1
83	*or IF (LT.EQ.1)WRITE(LUNOUT,3399)IOPT,NCOU
84	3399 FORMAT(' BAMJET',2I10 )
85	*or IF (LT.EQ.1)WRITE(LUNOUT,288)IHAD,KFA1,KFA2,KFA3,KFA4,AE0,IOPT
86	288 FORMAT (5I5,2E12.4,' BAMJET,IHAD,KFA1,KFA2,KFA3,KFA4,AE0,IOPT')
87	63 CONTINUE
88	DO 62 I=1,ITMX
89	KFR1(I) = 0
90	KFR2(I) = 0
91	62 CONTINUE
92	ITMX=0
93	60 CONTINUE
94	IYY=0
95	IHAD=0
96	IT=0
97	E0=AE0*.5D0
98	IF(IOPT.EQ.1.OR.IOPT.EQ.2) E0=AE0
99	LL=0
100	IF(KFA1.GT.6.AND.IOPT.EQ.1) LL=1
101	IF(KFA1.LE.6.AND.IOPT.EQ.2) LL=1
102	IF(KFA1.GT.6.AND.IOPT.EQ.4) LL=1
103	PGX = 0.D0
104	PGY = 0.D0
105	PGZ = 0.D0
106	RPX0 = 0.D0
107	RPY0 = 0.D0
108	* The following 6 initializations might be useless, but they make the
109	* code much clearer
110	KR1R = 0
111	KR2R = 0
112	KR1L = 0
113	KR2L = 0
114	RER = 0.D+00
115	REL = 0.D+00
116	DO 10 I=1,KMXJCM
117	KFR1(I) = 0
118	KFR2(I) = 0
119	LA = 0
120	IT = IT+1
121	J = IT-1
122	* The following line seems useless
123	* K = IT+1
124	40 CONTINUE
125	C*****CUT OFF TASK
126	* \| Abbrch is called to cut the chain
127	CALL ABBRCH(IT,LL,LA,LT,E0,PGX,PGY,PGZ,KFR1,KFR2,RE,
128	& KR1R,KR2R,KR1L,KR2L,RPX,RPY,RPZ,RPXR,RPYR,RPZR,RPXL,RPYL,RPZL,
129	& RER,REL,IV,B1,B2,KFA1,KFA2,KFA3,KFA4,IOPT,IYY)
130	ITMX=MAX(ITMX,IT)
131	IF(LA .EQ. 0) GO TO 20
132	IT=IT-1
133	IF(IOPT.EQ.3.AND.LL.EQ.0) GO TO 70
134	IF(IOPT.EQ.4.AND.KFA1.LE.6.AND.LL.EQ.0) GO TO 70
135	IF(IOPT.EQ.4.AND.KFA1.GT.6.AND.LL.EQ.1) GO TO 70
136	IF(IOPT.EQ.5.AND.LL.EQ.0) GO TO 70
137	GO TO 50
138	70 CONTINUE
139	IYY = 1
140	LL = 1
141	IF(IOPT.EQ.4.AND.KFA1.GT.6) LL = 0
142	IAR=IT
143	GO TO 30
144	20 CONTINUE
145	C*****CHOICE OF THE VERTEX
146	CALL FKVERT(IT,LT,LL,KFA1,E0,IV,RE,KFR1,KFR2,AME,IOPT)
147	C*****CHOICE OF THE FLAVOUR
148	CALL FKFLAV(IT,LT,LL,E0,IV,RE,KFR1,KFR2,ISU,BET,KFA1,KFA2,
149	& KFA3,KFA4,IOPT)
150	C*****CLASSIFICATION OF THE PARTICLES
151	CALL HKLASS(IT,LT,LA,LL,KFR1,KFR2,KR1R,KR2R,KR1L,KR2L,IV,IMPS,
152	& IMVE,IB08,IA08,IB10,IA10,AS,B8,KFA1,KFA2,KFA3,KFA4,IOPT)
153	ITMX=MAX(ITMX,IT)
154	IF (IT .EQ. 1) RX = E0
155	IF (IT .GT. 1) RX = RE(J)
156	IF(AMF(IT) .GT. RX) GO TO 63
157	IF(AMF(IT) .LE. RX) GO TO 19
158	LA = 1
159	GO TO 40
160	19 CONTINUE
161	IHAD = IHAD + 1
162	*or IF(LT .EQ. 0) GO TO 31
163	*or WRITE(LUNOUT,32)IHAD
164	*or 31 CONTINUE
165	C*****CHOICE OF THE ENERGY
166	CALL ENERGI(IT,LL,LT,IV,RE,HMA,HE,E0,A1)
167	C*****CHOICE OF THE MOMENTUM
168	* \|
169	* \| He is the total energy, hma the mass one (input) hpx, hpy, hpz
170	* \| the momentum components (output values), hps the transversal
171	* \| momentum (output)
172	* \|
173	CALL FKIMPU(HE,HMA,HPS,HPX,HPY,HPZ,LT,LL,B3)
174	IF (IT .GT. 1) GO TO 13
175	RPX(IT)=RPX0-HPX
176	RPY(IT)=RPY0-HPY
177	GO TO 14
178	13 RPX(IT)=RPX(J)-HPX
179	RPY(IT)=RPY(J)-HPY
180	14 CONTINUE
181	IF (IOPT.EQ.1.AND.LL.EQ.1)HPZ=-HPZ
182	IF(IOPT.EQ.2.AND.LL.EQ.1) HPZ=-HPZ
183	IF(IOPT.EQ.4.AND.KFA1.GT.6) HPZ=-HPZ
184	IF(IOPT.EQ.5) HPZ=-HPZ
185	PGX=PGX+HPX
186	PGY=PGY+HPY
187	PGZ=PGZ+HPZ
188	PXF(IT)=HPX
189	PYF(IT)=HPY
190	PZF(IT)=HPZ
191	*or IF (LT .EQ. 0) GO TO 15
192	*or WRITE(LUNOUT,16)PGX,PGY,PGZ
193	*or 16 FORMAT(1H0,12HPGX,PGY,PGZ=,3F8.4)
194	*or 15 CONTINUE
195	30 CONTINUE
196
197	10 CONTINUE
198	*
199	* we suppose that exiting from loop must be achieved via " go to 50
200	*
201	WRITE (LUNERR,*)' BAMJEV: EXITING FROM LOOP ABNORMALLY!!!! '
202	WRITE (LUNOUT,*)' BAMJEV: EXITING FROM LOOP ABNORMALLY!!!! '
203	50 CONTINUE
204	ITMX=MAX(ITMX,IT)
205	IF(IOPT.EQ.1.OR.IOPT.EQ.2) GO TO 51
206	C*****PUT THE RIGHT AND LEFT JET TOGETHER
207	CALL VEREIN(IT,LA,LT,RER,REL,RPXR,RPYR,RPZR,RPXL,RPYL,RPZL,
208	&KR1R,KR2R,KR1L,KR2L,IHAD,LL,KFR1,KFR2,IMPS,IMVE,IB08,IA08,
209	&IB10,IA10,B3,AS,B8,IAR,KFA1,KFA2,KFA3,KFA4,IOPT)
210	IF(LA.EQ.3) GO TO 63
211	IF(LA.EQ.2) GO TO 63
212	51 CONTINUE
213	IF(IOPT.EQ.3.OR.IOPT.EQ.4.OR.IOPT.EQ.5) GO TO 52
214	IF(LL.EQ.0) GO TO 52
215	RPXR=RPXL
216	RPYR=RPYL
217	RPZR=RPZL
218	RER=REL
219	KR1R=KR1L
220	KR2R=KR2L
221	52 CONTINUE
222	IF(LE.EQ.0) GO TO 1000
223	WRITE(LUNOUT,92)
224	92 FORMAT(2X,'NF,NAME,MASS,IQ,IB,PX,PY,PZ,E')
225	DO 91 J=1,IT
226	WRITE(LUNOUT,90)NREF(J),ANF(J),AMF(J),ICHF(J),IBARF(J),PXF(J),
227	& PYF(J),PZF(J),HEF(J)
228	90 FORMAT(2X,I3,A6,F6.3,2I4,4F8.4)
229	91 CONTINUE
230	1000 CONTINUE
231	64 FORMAT(1H0,38HNUMBER OF EVENTS WITH PREST GT. EREST=,I4,
232	*/,21HNUMBER OF ALL EVENTS=,I4)
233	2000 FORMAT(1H0,'NUMBER OF EVENTS WITH ONLY ONE PARTICLE=',I4)
234	*or IF(LT.EQ.0) GO TO 17
235	*or WRITE(LUNOUT,18)IHAD
236	*or 18 FORMAT(1H0,15HMULTIPLIZITAET=,I3)
237	*or 32 FORMAT(1H0,13HHADRONANZAHL=,I3)
238	*or 17 CONTINUE
239	A1 = A1SAVE
240	B1 = B1SAVE
241	B2 = B2SAVE
242	B3 = B3SAVE
243	RETURN
244	END