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