Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / gheisha / casp.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:21:00 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10*CMZ : 3.21/02 29/03/94 15.41.39 by S.Giani
11*-- Author :
12 SUBROUTINE CASP(K,INT,NFL)
13C
14C *** CASCADE OF PROTON ***
15C *** NVE 04-MAY-1988 CERN GENEVA ***
16C
17C ORIGIN : H.FESEFELDT (13-SEP-1987)
18C
19C P UNDERGOES INTERACTION WITH NUCLEON WITHIN NUCLEUS.
20C CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE PIONS/KAONS.
21C IF NOT ASSUME NUCLEAR EXCITATION OCCURS AND INPUT PARTICLE
22C IS DEGRADED IN ENERGY. NO OTHER PARTICLES PRODUCED.
23C IF REACTION IS POSSIBLE FIND CORRECT NUMBER OF PIONS/PROTONS/
24C NEUTRONS PRODUCED USING AN INTERPOLATION TO MULTIPLICITY DATA.
25C REPLACE SOME PIONS OR PROTONS/NEUTRONS BY KAONS OR STRANGE BARYONS
26C ACCORDING TO AVERAGE MULTIPLICITY PER INELASTIC REACTIONS.
27C
28#include "geant321/mxgkgh.inc"
29#include "geant321/s_consts.inc"
30#include "geant321/s_curpar.inc"
31#include "geant321/s_result.inc"
32#include "geant321/s_prntfl.inc"
33#include "geant321/s_kginit.inc"
34#include "geant321/limits.inc"
35C
36 REAL N
37 DIMENSION PMUL(2,1200),ANORM(2,60),SUPP(10),CECH(10),B(2)
38 DIMENSION RNDM(1)
39 SAVE PMUL,ANORM
40 DATA SUPP/0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98/
41 DATA CECH/0.50,0.45,0.40,0.35,0.30,0.25,0.06,0.04,0.005,0./
42 DATA B/0.70,0.35/,C/1.25/
43C
44C --- INITIALIZATION INDICATED BY KGINIT(10) ---
45 IF (KGINIT(10) .NE. 0) GO TO 10
46 KGINIT(10)=1
47C
48C --- INITIALIZE PMUL AND ANORM ARRAYS ---
49 DO 9000 J=1,1200
50 DO 9001 I=1,2
51 PMUL(I,J)=0.0
52 IF (J .LE. 60) ANORM(I,J)=0.0
53 9001 CONTINUE
54 9000 CONTINUE
55C
56C** COMPUTE NORMALIZATION CONSTANTS
57C** FOR P AS TARGET
58C
59 L=0
60 DO 1 NP1=1,20
61 NP=NP1-1
62 NMM1=NP1-2
63 IF(NMM1.LE.1) NMM1=1
64 DO 1 NM1=NMM1,NP1
65 NM=NM1-1
66 DO 1 NZ1=1,20
67 NZ=NZ1-1
68 L=L+1
69 IF(L.GT.1200) GOTO 1
70 NPROT=2-NP+NM
71 NNEUT=2-NPROT
72 NT=NP+NM+NZ
73 IF(NT.LE.0.OR.NT.GT.60) GOTO 1
74 PMUL(1,L)=PMLTPC(NP,NM,NZ,NT,B(1),C)
75 NPROTF=NFAC(NPROT)
76 NNEUTF=NFAC(NNEUT)
77 PMUL(1,L)=PMUL(1,L)/(NPROTF*NNEUTF)
78 ANORM(1,NT)=ANORM(1,NT)+PMUL(1,L)
79 1 CONTINUE
80C** FOR N AS TARGET
81 L=0
82 DO 2 NP1=1,20
83 NP=NP1-1
84 NMM1=NP1-1
85 IF(NMM1.LE.1) NMM1=1
86 NPP1=NP1+1
87 DO 2 NM1=NMM1,NPP1
88 NM=NM1-1
89 DO 2 NZ1=1,20
90 NZ=NZ1-1
91 L=L+1
92 IF(L.GT.1200) GOTO 2
93 NPROT=1-NP+NM
94 NNEUT=2-NPROT
95 NT=NP+NM+NZ
96 IF(NT.LE.0.OR.NT.GT.60) GOTO 2
97 PMUL(2,L)=PMLTPC(NP,NM,NZ,NT,B(2),C)
98 NPROTF=NFAC(NPROT)
99 NNEUTF=NFAC(NNEUT)
100 PMUL(2,L)=PMUL(2,L)/(NPROTF*NNEUTF)
101 ANORM(2,NT)=ANORM(2,NT)+PMUL(2,L)
102 2 CONTINUE
103 DO 3 I=1,60
104 IF(ANORM(1,I).GT.0.) ANORM(1,I)=1./ANORM(1,I)
105 IF(ANORM(2,I).GT.0.) ANORM(2,I)=1./ANORM(2,I)
106 3 CONTINUE
107 IF(.NOT.NPRT(10)) GOTO 10
108 WRITE(NEWBCD,2001)
109 DO 4 NFL=1,2
110 WRITE(NEWBCD,2002) NFL
111 WRITE(NEWBCD,2003) (ANORM(NFL,I),I=1,60)
112 WRITE(NEWBCD,2003) (PMUL(NFL,I),I=1,1200)
113 4 CONTINUE
114C** CHOOSE PROTON OR NEUTRON AS TARGET
115 10 NFL=2
116 CALL GRNDM(RNDM,1)
117 IF(RNDM(1).LT.ZNO2/ATNO2) NFL=1
118 TARMAS=RMASS(14)
119 IF (NFL .EQ. 2) TARMAS=RMASS(16)
120 S=AMASQ+TARMAS**2+2.0*TARMAS*EN
121 RS=SQRT(S)
122 ENP(8)=AMASQ+TARMAS**2+2.0*TARMAS*ENP(6)
123 ENP(9)=SQRT(ENP(8))
124 EAB=RS-TARMAS-RMASS(14)
125C** ELASTIC SCATTERING
126 NCECH=0
127 NP=0
128 NM=0
129 NZ=0
130 N=0.
131 IF(INT.EQ.2) GOTO 20
132C** INTRODUCE CHARGE EXCHANGE REACTION PN --> NP
133 IF(NFL.EQ.1) GOTO 100
134 IPLAB=IFIX(P*2.5)+1
135 IF(IPLAB.GT.10) IPLAB=10
136 CALL GRNDM(RNDM,1)
137 IF(RNDM(1).GT.CECH(IPLAB)/ATNO2**0.42) GOTO 100
138 NCECH=1
139 GOTO 100
140C** CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION IN REACT.
141 20 IF (EAB .LE. RMASS(7)) GOTO 55
142C** SUPPRESSION OF HIGH MULTIPLICITY EVENTS AT LOW MOMENTUM
143 IEAB=IFIX(EAB*5.)+1
144 IF(IEAB.GT.10) GOTO 22
145 CALL GRNDM(RNDM,1)
146 IF(RNDM(1).LT.SUPP(IEAB)) GOTO 22
147 N=1.
148 GOTO (23,24),NFL
149 23 CONTINUE
150 TEST=-(1+B(1))**2/(2.0*C**2)
151 IF (TEST .LT. EXPXL) TEST=EXPXL
152 IF (TEST .GT. EXPXU) TEST=EXPXU
153 W0=EXP(TEST)/2.0
154 WP=EXP(TEST)
155 CALL GRNDM(RNDM,1)
156 RAN=RNDM(1)
157 NP=0
158 NM=0
159 NZ=1
160 IF(RAN.LT.W0/(W0+WP)) GOTO 100
161 NP=1
162 NM=0
163 NZ=0
164 GOTO 100
165 24 CONTINUE
166 TEST=-(1+B(2))**2/(2.0*C**2)
167 IF (TEST .LT. EXPXL) TEST=EXPXL
168 IF (TEST .GT. EXPXU) TEST=EXPXU
169 W0=EXP(TEST)
170 WP=EXP(TEST)/2.0
171 TEST=-(-1+B(2))**2/(2.0*C**2)
172 IF (TEST .LT. EXPXL) TEST=EXPXL
173 IF (TEST .GT. EXPXU) TEST=EXPXU
174 WM=EXP(TEST)/2.0
175 WT=W0+WP+WM
176 WP=W0+WP
177 CALL GRNDM(RNDM,1)
178 RAN=RNDM(1)
179 NP=0
180 NM=0
181 NZ=1
182 IF(RAN.LT.W0/WT) GOTO 100
183 NP=1
184 NM=0
185 NZ=0
186 IF(RAN.LT.WP/WT) GOTO 100
187 NP=0
188 NM=1
189 NZ=0
190 GOTO 100
191C
192 22 ALEAB=LOG(EAB)
193C** NO. OF TOTAL PARTICLES VS SQRT(S)-2*MP
194 N=3.62567+0.665843*ALEAB+0.336514*ALEAB*ALEAB
195 * +0.117712*ALEAB*ALEAB*ALEAB+0.0136912*ALEAB*ALEAB*ALEAB*ALEAB
196 N=N-2.
197C** NORMALIZATION CONSTANT FOR KNO-DISTRIBUTION
198 ANPN=0.
199 DO 21 NT=1,60
200 TEST=-(PI/4.0)*(NT/N)**2
201 IF (TEST .LT. EXPXL) TEST=EXPXL
202 IF (TEST .GT. EXPXU) TEST=EXPXU
203 DUM1=PI*NT/(2.0*N*N)
204 DUM2=ABS(DUM1)
205 DUM3=EXP(TEST)
206 ADDNVE=0.0
207 IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
208 IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
209 ANPN=ANPN+ADDNVE
210 21 CONTINUE
211 ANPN=1./ANPN
212C** P OR N AS TARGET
213 CALL GRNDM(RNDM,1)
214 RAN=RNDM(1)
215 EXCS=0.
216 GOTO (30,40),NFL
217C** FOR P AS TARGET
218 30 L=0
219 DO 31 NP1=1,20
220 NP=NP1-1
221 NMM1=NP1-2
222 IF(NMM1.LE.1) NMM1=1
223 DO 31 NM1=NMM1,NP1
224 NM=NM1-1
225 DO 31 NZ1=1,20
226 NZ=NZ1-1
227 L=L+1
228 IF(L.GT.1200) GOTO 31
229 NT=NP+NM+NZ
230 IF(NT.LE.0.OR.NT.GT.60) GOTO 31
231 TEST=-(PI/4.0)*(NT/N)**2
232 IF (TEST .LT. EXPXL) TEST=EXPXL
233 IF (TEST .GT. EXPXU) TEST=EXPXU
234 DUM1=ANPN*PI*NT*PMUL(1,L)*ANORM(1,NT)/(2.0*N*N)
235 DUM2=ABS(DUM1)
236 DUM3=EXP(TEST)
237 ADDNVE=0.0
238 IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
239 IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
240 EXCS=EXCS+ADDNVE
241 IF(RAN.LT.EXCS) GOTO 100
242 31 CONTINUE
243 GOTO 80
244C** FOR N AS TARGET
245 40 L=0
246 DO 41 NP1=1,20
247 NP=NP1-1
248 NMM1=NP1-1
249 IF(NMM1.LE.1) NMM1=1
250 NPP1=NP1+1
251 DO 41 NM1=NMM1,NPP1
252 NM=NM1-1
253 DO 41 NZ1=1,20
254 NZ=NZ1-1
255 L=L+1
256 IF(L.GT.1200) GOTO 41
257 NT=NP+NM+NZ
258 IF(NT.LE.0.OR.NT.GT.60) GOTO 41
259 TEST=-(PI/4.0)*(NT/N)**2
260 IF (TEST .LT. EXPXL) TEST=EXPXL
261 IF (TEST .GT. EXPXU) TEST=EXPXU
262 DUM1=ANPN*PI*NT*PMUL(2,L)*ANORM(2,NT)/(2.0*N*N)
263 DUM2=ABS(DUM1)
264 DUM3=EXP(TEST)
265 ADDNVE=0.0
266 IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
267 IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
268 EXCS=EXCS+ADDNVE
269 IF(RAN.LT.EXCS) GOTO 100
270 41 CONTINUE
271 GOTO 80
272 50 IF(NPRT(4))
273 *WRITE(NEWBCD,1003) EAB,N,NFL,NP,NM,NZ
274 CALL STPAIR
275 IF(INT.EQ.1) CALL TWOB(14,NFL,N)
276 IF(INT.EQ.2) CALL GENXPT(14,NFL,N)
277 GO TO 9999
278 55 IF(NPRT(4))
279 *WRITE(NEWBCD,1001)
280 GOTO 53
281C** EXCLUSIVE REACTION NOT FOUND,ASSUME ELASTIC SCATTERING
282 80 IF(NPRT(4))
283 *WRITE(NEWBCD,1004) RS,N
284 53 INT=1
285 NP=0
286 NM=0
287 NZ=0
288 100 DO 101 I=1,60
289 101 IPA(I)=0
290 IF(INT.LE.0) GOTO 131
291 NPROT=2-NP+NM+(1-NFL)
292 NNEUT=2-NPROT
293 GOTO (102,112),NFL
294 102 GOTO (103,104),INT
295 103 IPA(1)=14
296 IPA(2)=14
297 NT=2
298 GOTO 130
299 104 IF(NNEUT.EQ.1) GOTO 105
300 IF(NNEUT.EQ.2) GOTO 106
301 IPA(1)=14
302 IPA(2)=14
303 GOTO 120
304 105 IPA(1)=14
305 IPA(2)=16
306 CALL GRNDM(RNDM,1)
307 IF(RNDM(1).LT.0.5) GOTO 120
308 IPA(1)=16
309 IPA(2)=14
310 GOTO 120
311 106 IPA(1)=16
312 IPA(2)=16
313 GOTO 120
314 112 GOTO (113,114),INT
315 113 IPA(1)=14
316 IPA(2)=16
317 NT=2
318 IF(NCECH.EQ.0) GOTO 130
319 IPA(1)=16
320 IPA(2)=14
321 GOTO 130
322 114 IF(NNEUT.EQ.1) GOTO 115
323 IF(NNEUT.EQ.2) GOTO 116
324 IPA(1)=14
325 IPA(2)=14
326 GOTO 120
327 115 IPA(1)=16
328 IPA(2)=14
329 CALL GRNDM(RNDM,1)
330 IF(RNDM(1).LT.0.33) GOTO 120
331 IPA(1)=14
332 IPA(2)=16
333 GOTO 120
334 116 IPA(1)=16
335 IPA(2)=16
336 120 NT=2
337 IF(NP.EQ.0) GOTO 122
338 DO 121 I=1,NP
339 NT=NT+1
340 121 IPA(NT)=7
341 122 IF(NM.EQ.0) GOTO 124
342 DO 123 I=1,NM
343 NT=NT+1
344 123 IPA(NT)=9
345 124 IF(NZ.EQ.0) GOTO 130
346 DO 125 I=1,NZ
347 NT=NT+1
348 125 IPA(NT)=8
349 130 IF(NPRT(4))
350 *WRITE(NEWBCD,2004) NT,(IPA(I),I=1,20)
351 GOTO 50
352 131 IF(NPRT(4))
353 *WRITE(NEWBCD,2005)
354C
3551001 FORMAT('0*CASP* CASCADE ENERGETICALLY NOT POSSIBLE',
356 $ ' CONTINUE WITH QUASI-ELASTIC SCATTERING')
3571003 FORMAT(' *CASP* PROTON-INDUCED CASCADE,',
358 $ ' AVAIL. ENERGY',2X,F8.4,
359 $ 2X,'<NTOT>',2X,F8.4,2X,'FROM',4(2X,I3),2X,'PARTICLES')
3601004 FORMAT(' *CASP* PROTON-INDUCED CASCADE,',
361 $ ' EXCLUSIVE REACTION NOT FOUND',
362 $ ' TRY ELASTIC SCATTERING AVAIL. ENERGY',2X,F8.4,2X,
363 $ ' <NTOT>',2X,F8.4)
3642001 FORMAT('0*CASP* TABLES FOR MULT. DATA PROTON INDUCED REACTION',
365 $ ' FOR DEFINITION OF NUMBERS SEE FORTRAN CODING')
3662002 FORMAT(' *CASP* TARGET PARTICLE FLAG',2X,I5)
3672003 FORMAT(1H ,10E12.4)
3682004 FORMAT(' *CASP* ',I3,2X,'PARTICLES , MASS INDEX ARRAY',2X,20I4)
3692005 FORMAT(' *CASP* NO PARTICLES PRODUCED')
370C
371 9999 CONTINUE
372 END