]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/ghrout/higclu.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / ghrout / higclu.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.2 1996/04/18 15:30:15 ravndal
6* PCM index overflow protection
7*
8* Revision 1.1.1.1 1995/10/24 10:21:13 cernlib
9* Geant
10*
11*
12#include "geant321/pilot.h"
13*CMZ : 3.21/02 29/03/94 15.41.40 by S.Giani
14*FCA : 05/01/99 12:38:13 by Federico Carminati
15* Added protection for uninitialised variables
16*-- Author :
17 SUBROUTINE HIGCLU(IPPP,NFL,AVERN)
18C
19C *** GENERATION OF X- AND PT- VALUES FOR ALL PRODUCED PARTICLES ***
20C *** NVE 01-AUG-1988 CERN GENEVA ***
21C
22C ORIGIN : H.FESEFELDT (11-OCT-1987)
23C
24C A SIMPLE TWO CLUSTER MODEL IS USED
25C THIS SHOULD BE SUFFICIENT FOR LOW ENERGY INTERACTIONS
26C
27C
28#include "geant321/s_defcom.inc"
29#include "geant321/s_genio.inc"
30C
31 REAL NUCSUP
32 DIMENSION SIDE(MXGKCU),C1PAR(5),G1PAR(5),NUCSUP(6)
33 DIMENSION RNDM(3)
34 DATA C1PAR/0.6,0.6,0.35,0.15,0.10/
35 DATA G1PAR/2.6,2.6,1.8,1.30,1.20/
36 DATA NUCSUP/1.0,0.7,0.5,0.4,0.35,0.3/
37 DIMENSION PSUP(6)
38 DATA PSUP/3.,6.,20.,50.,100.,1000./
39C DATA CB/3.0/
40 DATA CB/0.01/
41 BPP(X)=4.000+1.600*LOG(X)
42C
43 MX =MXGKPV-20
44 MX1=MX+1
45 MX2=MX+2
46 MX3=MX+3
47 MX4=MX+4
48 MX5=MX+5
49 MX6=MX+6
50 MX7=MX+7
51 MX8=MX+8
52 EK=ENP(5)
53 EN=ENP(6)
54 P=ENP(7)
55 S=ENP(8)
56 RS=ENP(9)
57 CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
58 SPALL=0
59 IF(P.LT.0.001) GOTO 60
60 NT=0
61C**
62C** CHECK MASS-INDICES FOR ALL PARTICLES
63C**
64 DO 1 I=1,100
65 IF(IPA(I).EQ.0) GOTO 1
66 NT=NT+1
67 IPA(NT)=IPA(I)
68 1 CONTINUE
69 CALL VZERO(IPA(NT+1),MXGKCU-NT)
70C**
71C** SET THE EFFECTICE 4-MOMENTUM-VECTOR FOR INTERACTION
72C**
73 PV( 1,MXGKPV-1)=P*PX
74 PV( 2,MXGKPV-1)=P*PY
75 PV( 3,MXGKPV-1)=P*PZ
76 PV( 4,MXGKPV-1)=EN
77 PV( 5,MXGKPV-1)=AMAS
78 PV( 6,MXGKPV-1)=NCH
79 PV( 7,MXGKPV-1)=TOF
80 PV( 8,MXGKPV-1)=IPART
81 PV( 9,MXGKPV-1)=0.
82 PV(10,MXGKPV-1)=USERW
83 IER(48)=IER(48)+1
84C**
85C** DISTRIBUTE PARTICLES IN FORWARD AND BACKWARD HEMISPHERE OF CMS
86C** OF THE HADRON NUCLEON INTERACTION
87C**
88 SIDE(1)= 1.
89 SIDE(2)=-1.
90 TARG=0.
91 IFOR=1
92 IBACK=1
93 DO 3 I=1,NT
94 IF (I .LE. 2) GO TO 78
95 SIDE(I)=1.
96 CALL GRNDM(RNDM,1)
97 IF (RNDM(1) .LT. 0.5) SIDE(I)=-1.
98 IF (SIDE(I) .LT. 0.) GO TO 76
99C
100C --- PARTICLE IN FORWARD HEMISPHERE ---
101 77 CONTINUE
102 IFOR=IFOR+1
103 IF (IFOR .LE. 18) GO TO 78
104C
105C --- CHANGE IT TO BACKWARD ---
106 SIDE(I)=-1.
107 IFOR=IFOR-1
108 IBACK=IBACK+1
109 GO TO 78
110C
111C --- PARTICLE IN BACKWARD HEMISPHERE ---
112 76 CONTINUE
113 IBACK=IBACK+1
114 IF (IBACK .LE. 18) GO TO 78
115C
116C --- CHANGE IT TO FORWARD ---
117 SIDE(I)=1.
118 IBACK=IBACK-1
119 IFOR=IFOR+1
120C**
121C** SUPPRESSION OF CHARGED PIONS FOR VARIOUS REASONS
122C**
123 78 IF(IPART.EQ.15.OR.IPART.GE.17) GOTO 3
124 IF(ABS(IPA(I)).GE.10) GOTO 3
125 IF(ABS(IPA(I)).EQ. 8) GOTO 3
126 CALL GRNDM(RNDM,1)
127 IF(RNDM(1).GT.(10.-P)/6.) GOTO 3
128 CALL GRNDM(RNDM,1)
129 IF(RNDM(1).GT.ATNO2/300.) GOTO 3
130 IPA(I)=14
131 CALL GRNDM(RNDM,1)
132 IF(RNDM(1).GT.ZNO2/ATNO2) IPA(I)=16
133 TARG=TARG+1.
134 3 CONTINUE
135 TB=2.*IBACK
136 CALL GRNDM(RNDM,1)
137 IF(RS.LT.(2.0+RNDM(1))) TB=(2.*IBACK+NT)/2.
138C**
139C** NUCLEONS + SOME PIONS FROM INTRANUCLEAR CASCADE
140C**
141 AFC=0.312+0.200*LOG(LOG(S))+S**1.5/6000.
142 IF(AFC.GT.0.50) AFC= 0.50
143 XTARG=AFC*(ATNO2**0.33-1.0)*TB
144 IF(XTARG.LE.0.) XTARG=0.01
145 DO 881 IPX=1,6
146 IF(P.LE.PSUP(IPX)) GOTO 882
147 881 CONTINUE
148 IPX=6
149 882 XPNHMF = XTARG*NUCSUP(IPX)
150 XSHHMF = XTARG - XPNHMF
151 IF(XSHHMF.LT.0.01) XSHHMF=0.01
152 IF(XPNHMF.LT.0.01) XPNHMF=0.01
153 SSHHMF=0.5*XSHHMF
154 SPNHMF=0.9*XPNHMF
155 RSHHMF=SSHHMF**2/XSHHMF
156 RPNHMF=SPNHMF**2/XPNHMF
157 IF(RSHHMF.LT.1.1) THEN
158 CALL POISSO(XSHHMF,NSHHMF)
159 GOTO 541
160 ELSE
161 RSHHMF=XSHHMF/(RSHHMF-1.)
162 IF(RSHHMF.LE.20.) THEN
163 CALL SVGAM7(RSHHMF,XHMF)
164 ELSE
165 KRSHMF=IFIX(RSHHMF+0.5)
166 CALL SVERL2(KRSHMF,XHMF)
167 ENDIF
168 XSHHMF=XHMF*XSHHMF/RSHHMF
169 CALL POISSO(XSHHMF,NSHHMF)
170 ENDIF
171 541 IF(RPNHMF.LE.1.1) THEN
172 CALL POISSO(XPNHMF,NPNHMF)
173 GOTO 542
174 ELSE
175 RPNHMF=XPNHMF/(RPNHMF-1.)
176 IF(RPNHMF.LE.20.) THEN
177 CALL SVGAM7(RPNHMF,XHMF)
178 ELSE
179 KRPHMF=IFIX(RPNHMF+0.5)
180 CALL SVERL2(KRPHMF,XHMF)
181 ENDIF
182 XPNHMF=XHMF*XPNHMF/RPNHMF
183 CALL POISSO(XPNHMF,NPNHMF)
184 ENDIF
185 542 NTARG=NSHHMF+NPNHMF
186 NT2=NT+NTARG
187 IF(NT2.LE.MXGKPV-30) GOTO 2
188 NT2=MXGKPV-30
189 NTARG=NT2-NT
190 2 CONTINUE
191 IF(NPRT(4))
192 *WRITE(NEWBCD,3001) NTARG,NT
193 NT1=NT+1
194 IF(NTARG.EQ.0) GOTO 51
195 DO 4 I=NT1,NT2
196 IF(NPNHMF.GT.0) GOTO 52
197 CALL GRNDM(RNDM,1)
198 IPA(I)=-(7+IFIX(RNDM(1)*3.0))
199 SIDE(I)=-1.
200 GOTO 4
201 52 IPA(I)=-16
202 PNRAT=1.-ZNO2/ATNO2
203 CALL GRNDM(RNDM,1)
204 IF(RNDM(1).GT.PNRAT) IPA(I)=-14
205 TARG=TARG+1.
206 SIDE(I)=-2.
207 NPNHMF=NPNHMF-1
208 4 CONTINUE
209 NT=NT2
210C**
211C** CHOOSE MASSES AND CHARGES FOR ALL PARTICLES
212C**
213 51 DO 5 I=1,NT
214 IPA1=ABS(IPA(I))
215 PV(5,I)=RMASS(IPA1)
216 PV(6,I)=RCHARG(IPA1)
217 PV(7,I)=1.
218 IF(PV(5,I).LT.0.) PV(7,I)=-1.
219 PV(5,I)=ABS(PV(5,I))
220 5 CONTINUE
221C**
222C** MARK LEADING STRANGE PARTICLES
223C**
224 LEAD=0
225 IF(IPART.LT.10.OR.IPART.EQ.14.OR.IPART.EQ.16) GOTO 6
226 IPA1=ABS(IPA(1))
227 IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 531
228 LEAD=IPA1
229 GOTO 6
230 531 IPA1=ABS(IPA(2))
231 IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 6
232 LEAD=IPA1
233C**
234C** CHECK AVAILABLE KINETIC ENERGY , CHANGE HEMISPHERE FOR PARTICLES
235C** UNTIL IT FITS
236C**
237 6 IF(NT.LE.1) GOTO 60
238 TAVAI=0.
239 DO 7 I=1,NT
240 IF(SIDE(I).LT.-1.5) GOTO 7
241 TAVAI=TAVAI+ABS(PV(5,I))
242 7 CONTINUE
243 IF(TAVAI.LT.RS) GOTO 12
244 IF(NPRT(4))
245 *WRITE(NEWBCD,3002) (IPA(I),I=1,20),(SIDE(I),I=1,20),TAVAI,RS
246 3002 FORMAT(' *HIGCLU* CHECK AVAILABLE ENERGIES'/
247 * 1H ,20I5/1H ,20F5.0/1H ,'TAVAI,RS ',2F10.3)
248 DO 10 I=1,NT
249 II=NT-I+1
250 IF(SIDE(II).LT.-1.5) GOTO 10
251 IF(II.EQ.NT) GOTO 11
252 NT1=II+1
253 NT2=NT
254 DO 8 J=NT1,NT2
255 IPA(J-1)=IPA(J)
256 SIDE(J-1)=SIDE(J)
257 DO 8 K=1,10
258 8 PV(K,J-1)=PV(K,J)
259 GOTO 11
260 10 CONTINUE
261 11 SIDE(NT)=0.
262 IPA(NT)=0
263 NT=NT-1
264 GOTO 6
265 12 IF(NT.LE.1) GOTO 60
266 B=BPP(P)
267 IF(B.LT.CB) B=CB
268C**
269C** CHOOSE MASSES FOR THE 3 CLUSTER: 1. FORWARD CLUSTER
270C** 2. BACKWARD MESON CLUSTER 3. BACKWARD NUCLEON CLUSTER
271C**
272 RMC0=0.
273 RMD0=0.
274 RME0=0.
275 NTC=0
276 NTD=0
277 NTE=0
278 DO 31 I=1,NT
279 IF(SIDE(I).GT.0.) RMC0=RMC0+ABS(PV(5,I))
280 IF(SIDE(I).GT.0.) NTC =NTC +1
281 IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) RMD0=RMD0+ABS(PV(5,I))
282 IF( SIDE(I).LT.-1.5) RME0=RME0+ABS(PV(5,I))
283 IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) NTD =NTD +1
284 IF( SIDE(I).LT.-1.5) NTE =NTE +1
285 31 CONTINUE
286 32 CALL GRNDM(RNDM,1)
287 RAN=RNDM(1)
288 RMC=RMC0
289 IF(NTC.LE.1) GOTO 33
290 NTC1=NTC
291 IF(NTC1.GT.5) NTC1=5
292 RMC=-LOG(1.-RAN)
293 GPAR=G1PAR(NTC1)
294 CPAR=C1PAR(NTC1)
295 DUMNVE=GPAR
296 IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
297 RMC=RMC0+RMC**CPAR/DUMNVE
298 33 RMD=RMD0
299 IF(NTD.LE.1) GOTO 34
300 NTD1=NTD
301 IF(NTD1.GT.5) NTD1=5
302 CALL GRNDM(RNDM,1)
303 RAN=RNDM(1)
304 RMD=-LOG(1.-RAN)
305 GPAR=G1PAR(NTD1)
306 CPAR=C1PAR(NTD1)
307 DUMNVE=GPAR
308 IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
309 RMD=RMD0+RMD**CPAR/DUMNVE
310 34 IF(RMC+RMD.LE.RS) GOTO 35
311 IF (RMC.LE.RMC0.AND.RMD.LE.RMD0) THEN
312 HNRMDC = 0.999*RS/(RMC+RMD)
313 RMD = RMD*HNRMDC
314 RMC = RMC*HNRMDC
315 ELSE
316 RMC=0.1*RMC0+0.9*RMC
317 RMD=0.1*RMD0+0.9*RMD
318 END IF
319 GOTO 34
320 35 IF(NTE.LE.0) GOTO 38
321 RME=RME0
322 IF(NTE.EQ.1) GOTO 38
323 NTE1=NTE
324 IF(NTE1.GT.5) NTE1=5
325 CALL GRNDM(RNDM,1)
326 RAN=RNDM(1)
327 RME=-LOG(1.-RAN)
328 GPAR=G1PAR(NTE1)
329 CPAR=C1PAR(NTE1)
330 DUMNVE=GPAR
331 IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
332 RME=RME0+RME**CPAR/DUMNVE
333C**
334C** SET BEAM , TARGET OF FIRST INTERACTION IN CMS
335C**
336 38 PV( 1,MX1) =0.
337 PV( 2,MX1) =0.
338 PV( 3,MX1) =P
339 PV( 5,MX1) =ABS(AMAS)
340 PV( 4,MX1) =SQRT(P*P+AMAS*AMAS)
341 PV( 1,MX2) =0.
342 PV( 2,MX2) =0.
343 PV( 3,MX2) =0.
344 PV( 4,MX2) =MP
345 PV( 5,MX2) =MP
346
347C** TRANSFORM INTO CMS.
348
349 CALL ADD(MX1,MX2,MX)
350 CALL LOR(MX1,MX,MX1)
351 CALL LOR(MX2,MX,MX2)
352 PF=(S+RMD*RMD-RMC*RMC)**2 - 4*S*RMD*RMD
353 IF(PF.LT.0.0001) PF=0.0001
354 DUMNVE=2.0*RS
355 IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
356 PF=SQRT(PF)/DUMNVE
357 IF(NPRT(4)) WRITE(6,2002) PF,RMC,RMD,RS
358C**
359C** SET FINAL STATE MASSES AND ENERGIES IN CMS
360C**
361 PV(5,MX3) =RMC
362 PV(5,MX4) =RMD
363 PV(4,MX3) =SQRT(PF*PF+RMC*RMC)
364 PV(4,MX4) =SQRT(PF*PF+RMD*RMD)
365C**
366C** SET |T| AND |TMIN|
367C**
368 T=-1.0E10
369 CALL GRNDM(RNDM,1)
370 IF (B .NE. 0.0) T=LOG(1.-RNDM(1))/B
371 CALL LENGTX(MX1,PIN)
372 TACMIN=(PV(4,MX1) -PV(4,MX3))**2 -(PIN-PF)**2
373C**
374C** CACULATE (SIN(TETA/2.)**2 AND COS(TETA), SET AZIMUTH ANGLE PHI
375C**
376 DUMNVE=4.0*PIN*PF
377 IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
378 CTET=-(T-TACMIN)/DUMNVE
379 CTET=1.0-2.0*CTET
380 IF (CTET .GT. 1.0) CTET=1.0
381 IF (CTET .LT. -1.0) CTET=-1.0
382 DUMNVE=1.0-CTET*CTET
383 IF (DUMNVE .LT. 0.0) DUMNVE=0.0
384 STET=SQRT(DUMNVE)
385 CALL GRNDM(RNDM,1)
386 PHI=RNDM(1)*TWPI
387C**
388C** CALCULATE FINAL STATE MOMENTA IN CMS
389C**
390 PV(1,MX3) =PF*STET*SIN(PHI)
391 PV(2,MX3) =PF*STET*COS(PHI)
392 PV(3,MX3) =PF*CTET
393 PV(1,MX4) =-PV(1,MX3)
394 PV(2,MX4) =-PV(2,MX3)
395 PV(3,MX4) =-PV(3,MX3)
396C**
397C** SIMULATE BACKWARD NUCLEON CLUSTER IN LAB. SYSTEM AND TRANSFORM IN
398C** CMS.
399C**
400 IF(NTE.EQ.0) GOTO 28
401 GA=1.2
402 EKIT1=0.04
403 EKIT2=0.6
404 IF(EK.GT.5.) GOTO 666
405 EKIT1=EKIT1*EK**2/25.
406 EKIT2=EKIT2*EK**2/25.
407 666 A=(1.-GA)/(EKIT2**(1.-GA)-EKIT1**(1.-GA))
408 DO 29 I=1,NT
409 IF(SIDE(I).GT.-1.5) GOTO 29
410 CALL GRNDM(RNDM,3)
411 RAN=RNDM(1)
412 EKIT=(RAN*(1.-GA)/A+EKIT1**(1.-GA))**(1./(1.-GA))
413 PV(4,I)=EKIT+PV(5,I)
414 DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
415 PP=SQRT(DUMNVE)
416 RAN=RNDM(2)
417 COST=LOG(2.23*RAN+0.383)/0.96
418 IF (COST .LT. -1.0) COST=-1.0
419 IF (COST .GT. 1.0) COST=1.0
420 DUMNVE=1.0-COST*COST
421 IF (DUMNVE .LT. 0.0) DUMNVE=0.0
422 SINT=SQRT(DUMNVE)
423 PHI=TWPI*RNDM(3)
424 PV(1,I)=PP*SINT*SIN(PHI)
425 PV(2,I)=PP*SINT*COS(PHI)
426 PV(3,I)=PP*COST
427 CALL LOR(I,MX,I)
428 29 CONTINUE
429C**
430C** FRAGMENTATION OF FORWARD CLUSTER AND BACKWARD MESON CLUSTER
431C**
432 28 PV(1,1)=PV(1,MX3)
433 PV(2,1)=PV(2,MX3)
434 PV(3,1)=PV(3,MX3)
435 PV(4,1)=PV(4,MX3)
436 PV(1,2)=PV(1,MX4)
437 PV(2,2)=PV(2,MX4)
438 PV(3,2)=PV(3,MX4)
439 PV(4,2)=PV(4,MX4)
440 DO 17 I=MX5,MX6
441 DO 16 J=1,3
442 16 PV(J,I)=-PV(J,I-2)
443 DO 17 J=4,5
444 17 PV(J,I)= PV(J,I-2)
445 KGENEV=1
446 IF(NTC.LE.1) GOTO 26
447 TECM=PV(5,MX3)
448 NPG=0
449 DO 18 I=1,NT
450 IF(SIDE(I).LT.0.) GOTO 18
451 IF(NPG.EQ.18) THEN
452 SIDE(I)=-SIDE(I)
453 GOTO 18
454 ENDIF
455 NPG=NPG+1
456 AMASS(NPG)=ABS(PV(5,I))
457 18 CONTINUE
458 IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG)
459 CALL PHASP
460 NPG=0
461 DO 19 I=1,NT
462 IF(SIDE(I).LT.0.OR.NPG.GE.18) GOTO 19
463 NPG=NPG+1
464 PV(1,I)=PCM(1,NPG)
465 PV(2,I)=PCM(2,NPG)
466 PV(3,I)=PCM(3,NPG)
467 PV(4,I)=PCM(4,NPG)
468 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
469 CALL LOR(I,MX5,I)
470 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
471 19 CONTINUE
472 26 IF(NTD.LE.1) GOTO 27
473 TECM=PV(5,MX4)
474 NPG=0
475 DO 20 I=1,NT
476 IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 20
477 IF(NPG.EQ.18) THEN
478 SIDE(I)=-2.
479 PV(4,I)=ABS(PV(5,I))
480 DO 48 J=1,3
481 PV(J,I)=0.
482 48 CONTINUE
483 GOTO 20
484 ENDIF
485 NPG=NPG+1
486 AMASS(NPG)=ABS(PV(5,I))
487 20 CONTINUE
488 IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG)
489 CALL PHASP
490 NPG=0
491 DO 21 I=1,NT
492 IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 21
493 NPG=NPG+1
494 PV(1,I)=PCM(1,NPG)
495 PV(2,I)=PCM(2,NPG)
496 PV(3,I)=PCM(3,NPG)
497 PV(4,I)=PCM(4,NPG)
498 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
499 CALL LOR(I,MX6,I)
500 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
501 21 CONTINUE
502C**
503C** LORENTZ TRANSFORMATION IN LAB SYSTEM
504C**
505 27 TARG=0.
506 DO 36 I=1,NT
507 IF(PV(5,I).GT.0.5) TARG=TARG+1.
508 CALL LOR(I,MX2,I)
509 36 CONTINUE
510 IF(TARG.LT.0.5) TARG=1.
511C**
512C** SOMETIMES THE LEADING STRANGE PARTICLES ARE LOST , SET THEM BACK
513C**
514 IF(LEAD.EQ.0) GOTO 6085
515 DO 6081 I=1,NT
516 IF(ABS(IPA(I)).EQ.LEAD) GOTO 6085
517 6081 CONTINUE
518 I=1
519 IF(LEAD.GE.14.AND.ABS(IPA(2)).GE.14) I=2
520 IF(LEAD.LT.14.AND.ABS(IPA(2)).LT.14) I=2
521 IPA(I)=LEAD
522 EKIN=PV(4,I)-ABS(PV(5,I))
523 PV(5,I)=RMASS(LEAD)
524 PV(7,I)=1.
525 IF(PV(5,I).LT.0.) PV(7,I)=-1.
526 PV(5,I)=ABS(PV(5,I))
527 PV(6,I)=RCHARG(LEAD)
528 PV(4,I)=PV(5,I)+EKIN
529 CALL LENGTX(I,PP)
530 DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
531 PP1=SQRT(DUMNVE)
532C
533 IF (PP .GE. 1.0E-6) GO TO 8000
534 CALL GRNDM(RNDM,2)
535 RTHNVE=PI*RNDM(1)
536 PHINVE=TWPI*RNDM(2)
537 PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
538 PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
539 PV(3,I)=PP1*COS(RTHNVE)
540 GO TO 8001
541 8000 CONTINUE
542 PV(1,I)=PV(1,I)*PP1/PP
543 PV(2,I)=PV(2,I)*PP1/PP
544 PV(3,I)=PV(3,I)*PP1/PP
545 8001 CONTINUE
546C
547C** FOR VARIOUS REASONS, THE ENERGY BALANCE IS NOT SUFFICIENT,
548C** CHECK THAT, ENERGY BALANCE, ANGLE OF FINAL SYSTEM E.T.C.
549 6085 KGENEV=1
550 PV(1,MX4) =0.
551 PV(2,MX4) =0.
552 PV(3,MX4) =P
553 PV(4,MX4) =SQRT(P*P+AMAS*AMAS)
554 PV(5,MX4) =ABS(AMAS)
555 EKIN0=PV(4,MX4) -PV(5,MX4)
556 PV(1,MX5) =0.
557 PV(2,MX5) =0.
558 PV(3,MX5) =0.
559 PV(4,MX5) =MP*TARG
560 PV(5,MX5) =PV(4,MX5)
561 EKIN=PV(4,MX4) +PV(4,MX5)
562 I=MX4
563 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
564 I=MX5
565 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
566 CALL ADD(MX4,MX5,MX6)
567 CALL LOR(MX4,MX6,MX4)
568 CALL LOR(MX5,MX6,MX5)
569 TECM=PV(4,MX4) +PV(4,MX5)
570 NPG=NT
571 PV(1,MX8) =0.
572 PV(2,MX8) =0.
573 PV(3,MX8) =0.
574 PV(4,MX8) =0.
575 PV(5,MX8) =0.
576 EKIN1=0.
577 DO 598 I=1,NPG
578 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
579 CALL ADD(MX8,I,MX8)
580 EKIN1=EKIN1+PV(4,I)-PV(5,I)
581 EKIN=EKIN-PV(5,I)
582 IF(I.GT.18) GOTO 598
583 AMASS(I)=PV(5,I)
584 598 CONTINUE
585 IF(NPG.GT.18) GOTO 597
586 CALL PHASP
587 EKIN=0.
588 DO 599 I=1,NPG
589 PV(1,MX7)=PCM(1,I)
590 PV(2,MX7)=PCM(2,I)
591 PV(3,MX7)=PCM(3,I)
592 PV(4,MX7)=PCM(4,I)
593 PV(5,MX7)=AMASS(I)
594 CALL LOR(MX7,MX5,MX7)
595 599 EKIN=EKIN+PV(4,MX7)-PV(5,MX7)
596 CALL ANG(MX8,MX4,COST,TETA)
597 IF(NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1,EKIN
598C**
599C** MAKE SHURE, THAT KINETIC ENERGIES ARE CORRECT
600C** THE 3. CLUSTER IS NOT PRODUCED WITHIN PROPER KINEMATICS!!!
601C** EKIN= KINETIC ENERGY THEORETICALLY
602C** EKIN1= KINETIC ENERGY SIMULATED
603C**
604 597 IF(EKIN1.EQ.0.) GOTO 600
605 PV(1,MX7) =0.
606 PV(2,MX7) =0.
607 PV(3,MX7) =0.
608 PV(4,MX7) =0.
609 PV(5,MX7) =0.
610 WGT=EKIN/EKIN1
611 EKIN1=0.
612 DO 602 I=1,NT
613 EKIN=PV(4,I)-PV(5,I)
614 EKIN=EKIN*WGT
615 PV(4,I)=EKIN+PV(5,I)
616 DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
617 PP=SQRT(DUMNVE)
618 CALL LENGTX(I,PP1)
619C
620 IF (PP1 .GE. 1.0E-6) GO TO 8002
621 CALL GRNDM(RNDM,2)
622 RTHNVE=PI*RNDM(1)
623 PHINVE=TWPI*RNDM(2)
624 PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE)
625 PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE)
626 PV(3,I)=PP*COS(RTHNVE)
627 GO TO 8003
628 8002 CONTINUE
629 PV(1,I)=PV(1,I)*PP/PP1
630 PV(2,I)=PV(2,I)*PP/PP1
631 PV(3,I)=PV(3,I)*PP/PP1
632 8003 CONTINUE
633C
634 EKIN1=EKIN1+EKIN
635 CALL ADD(MX7,I,MX7)
636 602 CONTINUE
637 CALL ANG(MX7,MX4,COST,TETA)
638 IF(NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1
639C**
640C** ROTATE IN DIRECTION OF Z-AXIS, SEE COMMENTS IN 'GENXPT'
641C**
642 600 PV(1,MX7)=0.
643 PV(2,MX7)=0.
644 PV(3,MX7)=0.
645 PV(4,MX7)=0.
646 PV(5,MX7)=0.
647 DO 596 I=1,NT
648 596 CALL ADD(MX7,I,MX7)
649* call rannor(ran1,ran2)
650 CALL GRNDM(RNDM,2)
651 RY=RNDM(1)
652 RZ=RNDM(2)
653 RX=6.283185*RZ
654 A1=SQRT(-2.*LOG(RY))
655 RAN1=A1*SIN(RX)
656 RAN2=A1*COS(RX)
657 PV(1,MX7)=PV(1,MX7)+RAN1*0.020*TARG
658 PV(2,MX7)=PV(2,MX7)+RAN2*0.020*TARG
659 CALL DEFS(MX4,MX7,MX8)
660 PV(1,MX7)=0.
661 PV(2,MX7)=0.
662 PV(3,MX7)=0.
663 PV(4,MX7)=0.
664 PV(5,MX7)=0.
665 DO 595 I=1,NT
666 CALL TRAC(I,MX8,I)
667 595 CALL ADD(MX7,I,MX7)
668 CALL ANG(MX7,MX4,COST,TETA)
669 IF(NPRT(4)) WRITE(NEWBCD,2003) TETA
670C**
671C** ROTATE IN DIRECTION OF PRIMARY PARTICLE
672C**
673 DEKIN=0.
674 NPIONS=0
675 EK1=0.
676 DO 25 I=1,NT
677 CALL DEFS1(I,MXGKPV-1,I)
678 IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
679 IF(ATNO2.LT.1.5) GOTO 25
680 CALL LENGTX(I,PP)
681 EKIN=PV(4,I)-ABS(PV(5,I))
682 CALL NORMAL(RAN)
683 EKIN=EKIN-CFA*(1.+0.5*RAN)
684 IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
685 CALL STEEQ(XXH,I)
686 DEKIN=DEKIN+EKIN*(1.-XXH)
687 EKIN=EKIN*XXH
688 IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) NPIONS=NPIONS+1
689 IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) EK1=EK1+EKIN
690 PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
691 PV(4,I)=EKIN+ABS(PV(5,I))
692C
693 IF (PP .GE. 1.0E-6) GO TO 8004
694 CALL GRNDM(RNDM,2)
695 RTHNVE=PI*RNDM(1)
696 PHINVE=TWPI*RNDM(2)
697 PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
698 PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
699 PV(3,I)=PP1*COS(RTHNVE)
700 GO TO 8005
701 8004 CONTINUE
702 PV(1,I)=PV(1,I)*PP1/PP
703 PV(2,I)=PV(2,I)*PP1/PP
704 PV(3,I)=PV(3,I)*PP1/PP
705 8005 CONTINUE
706C
707 25 CONTINUE
708 IF(EK1.EQ.0.) GOTO 23
709 IF(NPIONS.LE.0) GOTO 23
710 DEKIN=1.+DEKIN/EK1
711 DO 22 I=1,NT
712 IF(ABS(IPA(I)).LT.7.OR.ABS(IPA(I)).GT.9) GOTO 22
713 CALL LENGTX(I,PP)
714 EKIN=PV(4,I)-ABS(PV(5,I))
715 EKIN=EKIN*DEKIN
716 IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
717 PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
718 PV(4,I)=EKIN+ABS(PV(5,I))
719C
720 IF (PP .GE. 1.0E-6) GO TO 8006
721 CALL GRNDM(RNDM,2)
722 RTHNVE=PI*RNDM(1)
723 PHINVE=TWPI*RNDM(2)
724 PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
725 PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
726 PV(3,I)=PP1*COS(RTHNVE)
727 GO TO 8007
728 8006 CONTINUE
729 PV(1,I)=PV(1,I)*PP1/PP
730 PV(2,I)=PV(2,I)*PP1/PP
731 PV(3,I)=PV(3,I)*PP1/PP
732 8007 CONTINUE
733C
734 22 CONTINUE
735 23 IF(ATNO2.LT.1.5) GOTO 40
736C**
737C** ADD BLACK TRACK PARTICLES
738C**
739 CALL HIGHAB(SPROB)
740 CALL GRNDM(RNDM,1)
741 IF(RNDM(1).LT.SPROB) GOTO 40
742 TEX=ENP(1)
743 SPALL=TARG
744 IF(TEX.LT.0.001) GOTO 445
745 BLACK=(1.5+1.25*TARG)*ENP(1)/(ENP(1)+ENP(3))
746 CALL POISSO(BLACK,NBL)
747 IF(NPRT(4))
748 *WRITE(NEWBCD,3003) NBL,TEX
749 IF(IFIX(TARG)+NBL.GT.ATNO2) NBL=ATNO2-TARG
750 IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
751 IF(NBL.LE.0) GOTO 445
752 EKIN=TEX/NBL
753 EKIN2=0.
754 CALL STEEP(XX)
755 DO 441 I=1,NBL
756 CALL GRNDM(RNDM,1)
757 IF(RNDM(1).LT.SPROB) GOTO 441
758 IF(NT.EQ.MXGKPV-2) GOTO 441
759 IF(EKIN2.GT.TEX) GOTO 443
760 CALL GRNDM(RNDM,1)
761 RAN1=RNDM(1)
762 CALL NORMAL(RAN2)
763 EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
764 IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
765 EKIN1=EKIN1*XX
766 EKIN2=EKIN2+EKIN1
767 IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
768 IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
769 IPA1=16
770 PNRAT=1.-ZNO2/ATNO2
771 CALL GRNDM(RNDM,3)
772 IF(RNDM(1).GT.PNRAT) IPA1=14
773 NT=NT+1
774 SPALL=SPALL+1.
775 COST=-1.0+RNDM(2)*2.0
776 DUMNVE=1.0-COST*COST
777 IF (DUMNVE .LT. 0.0) DUMNVE=0.0
778 SINT=SQRT(DUMNVE)
779 PHI=TWPI*RNDM(3)
780 IPA(NT)=-IPA1
781 SIDE(NT)=-4.
782 PV(5,NT)=ABS(RMASS(IPA1))
783 PV(6,NT)=RCHARG(IPA1)
784 PV(7,NT)=1.
785 PV(4,NT)=EKIN1+PV(5,NT)
786 DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
787 PP=SQRT(DUMNVE)
788 PV(1,NT)=PP*SINT*SIN(PHI)
789 PV(2,NT)=PP*SINT*COS(PHI)
790 PV(3,NT)=PP*COST
791 441 CONTINUE
792 443 IF(ATNO2.LT.10.) GOTO 445
793 IF(EK.GT.2.0) GOTO 445
794 II=NT+1
795 KK=0
796 EKA=EK
797 IF(EKA.GT.1.) EKA=EKA*EKA
798 IF(EKA.LT.0.1) EKA=0.1
799 IKA=3.6*EXP((ZNO2**2/ATNO2-35.56)/6.45)/EKA
800 IF(IKA.LE.0) GO TO 445
801 DO 444 I=1,NT
802 II=II-1
803 IF(IPA(II).NE.-14) GOTO 444
804 IPA(II)=-16
805 IPA1 = 16
806 PV(5,II)=ABS(RMASS(IPA1))
807 PV(6,II)=RCHARG(IPA1)
808 KK=KK+1
809 IF(KK.GT.IKA) GOTO 445
810 444 CONTINUE
811 445 TEX=ENP(3)
812 IF(TEX.LT.0.001) GOTO 40
813 BLACK=(1.5+1.25*TARG)*ENP(3)/(ENP(1)+ENP(3))
814 CALL POISSO(BLACK,NBL)
815 IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
816 IF(NBL.LE.0) GOTO 40
817 EKIN=TEX/NBL
818 EKIN2=0.
819 CALL STEEP(XX)
820 IF(NPRT(4))
821 *WRITE(NEWBCD,3004) NBL,TEX
822 DO 442 I=1,NBL
823 CALL GRNDM(RNDM,1)
824 IF(RNDM(1).LT.SPROB) GOTO 442
825 IF(NT.EQ.MXGKPV-2) GOTO 442
826 IF(EKIN2.GT.TEX) GOTO 40
827 CALL GRNDM(RNDM,1)
828 RAN1=RNDM(1)
829 CALL NORMAL(RAN2)
830 EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
831 IF(EKIN1.LT.0.0) EKIN1=-0.005*LOG(RAN1)
832 EKIN1=EKIN1*XX
833 EKIN2=EKIN2+EKIN1
834 IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
835 IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
836 CALL GRNDM(RNDM,3)
837 COST=-1.0+RNDM(1)*2.0
838 DUMNVE=1.0-COST*COST
839 IF (DUMNVE .LT. 0.0) DUMNVE=0.0
840 SINT=SQRT(DUMNVE)
841 PHI=TWPI*RNDM(2)
842 RAN=RNDM(3)
843 IPA(NT+1)=-30
844 IF(RAN.GT.0.60) IPA(NT+1)=-31
845 IF(RAN.GT.0.90) IPA(NT+1)=-32
846 SIDE(NT+1)=-4.
847 PV(5,NT+1)=(ABS(IPA(NT+1))-28)*MP
848 SPALL=SPALL+PV(5,NT+1)*1.066
849 IF(SPALL.GT.ATNO2) GOTO 40
850 NT=NT+1
851 PV(6,NT)=1.
852 IF(IPA(NT).EQ.-32) PV(6,NT)=2.
853 PV(7,NT)=1.
854 PV(4,NT)=PV(5,NT)+EKIN1
855 DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
856 PP=SQRT(DUMNVE)
857 PV(1,NT)=PP*SINT*SIN(PHI)
858 PV(2,NT)=PP*SINT*COS(PHI)
859 PV(3,NT)=PP*COST
860 442 CONTINUE
861C**
862C** STORE ON EVENT COMMON
863C**
864 40 CALL GRNDM(RNDM,1)
865 IF(RS.GT.(4.+RNDM(1)*1.)) GOTO 42
866 DO 41 I=1,NT
867 CALL LENGTX(I,ETB)
868 IF(ETB.LT.P) GOTO 41
869 ETF=P
870 PV(4,I)=SQRT(PV(5,I)**2+ETF**2)
871 DUMNVE=ETB
872 IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
873 ETF=ETF/DUMNVE
874 PV(1,I)=PV(1,I)*ETF
875 PV(2,I)=PV(2,I)*ETF
876 PV(3,I)=PV(3,I)*ETF
877 41 CONTINUE
878 42 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
879 EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
880 EKIN2=0.
881 CALL TDELAY(TOF1)
882 CALL GRNDM(RNDM,1)
883 RAN=RNDM(1)
884 TOF=TOF-TOF1*LOG(RAN)
885 DO 44 I=1,NT
886 IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I)
887 PV(7,I)=TOF
888 PV(8,I)=ABS(IPA(I))
889 PV(9,I)=0.
890 44 PV(10,I)=0.
891 CALL GHETUN(NT)
892 DO 45 I=1,NT
893 EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
894 45 CONTINUE
895 EKIN2=(EKIN2-EKIN)/EKIN
896 IF(NPRT(4))
897 $ WRITE(NEWBCD,2006) NT,EKIN,ENP(1),ENP(3),EKIN1,EKIN2
898 IF(EKIN2.GT.0.2) GOTO 60
899 INTCT=INTCT+1.
900 NMODE=3
901 IF(SPALL.LT.0.5.AND.ATNO2.GT.1.5) NMODE=14
902 CALL SETCUR(NT)
903 NTK=NTK+1
904 IF(NT.EQ.1) GOTO 300
905 DO 50 II=2,NT
906 I=II-1
907 IF(NTOT.LT.NSIZE/12) GOTO 43
908 GO TO 9999
909 43 CALL SETTRK(I)
910 50 CONTINUE
911 300 CONTINUE
912 GO TO 9999
913C**
914C** IT IS NOT POSSIBLE TO PRODUCE A PROPER TWO CLUSTER FINAL STATE.
915C** CONTINUE WITH QUASI ELASTIC SCATTERING
916C**
917 60 IF(NPRT(4)) WRITE(NEWBCD,2005)
918 DO 61 I=3,MXGKCU
919 61 IPA(I)=0
920 IPA(1)=IPART
921 IPA(2)=14
922 IF(NFL.EQ.2) IPA(2)=16
923 CALL TWOB(IPPP,NFL,AVERN)
924 GO TO 9999
925C
926 2000 FORMAT(' *HIGCLU* CMS PARAMETERS OF FINAL STATE PARTICLES',
927 $ ' AFTER ',I3,' TRIALS')
928 2001 FORMAT(' *HIGCLU* TRACK',2X,I3,2X,10F8.2,2X,I3,2X,F3.0)
929 2002 FORMAT(' *HIGCLU* MOMENTUM ',F8.3,' MASSES ',2F8.4,' RS ',F8.4)
930 2003 FORMAT(' *HIGCLU* TETA,EKIN0,EKIN1,EKIN ',4F10.4)
931 2004 FORMAT(' *HIGCLU* TECM,NPB,MASSES: ',F10.4,1X,I3,1X,8F10.4/
932 $ 1H ,26X,15X,8F10.4)
933 2005 FORMAT(' *HIGCLU* NUMBER OF FINAL STATE PARTICLES',
934 $ ' LESS THAN 2 ==> CONTINUE WITH 2-BODY SCATTERING')
935 2006 FORMAT(' *HIGCLU* COMP.',1X,I5,1X,5F7.2)
936 3001 FORMAT(' *HIGCLU* NUCLEAR EXCITATION ',I5,' PARTICLES PRODUCED',
937 $ ' IN ADDITION TO',I5,' NORMAL PARTICLES')
938 3003 FORMAT(' *HIGCLU* ',I3,' BLACK TRACK PARTICLES PRODUCED',
939 $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV')
940 3004 FORMAT(' *HIGCLU* ',I5,' HEAVY FRAGMENTS WITH TOTAL ENERGY OF ',
941 $ F8.4,' GEV')
942C
943 9999 CONTINUE
944 END