source: trunk/MagicSoft/Simulation/Corsika/Mmcs/cghei.f@ 9872

Last change on this file since 9872 was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 24.2 KB
Line 
1 SUBROUTINE CGHEI
2
3C-----------------------------------------------------------------------
4C C(ORSIKA) GHE(ISHA) I(NTERFACE)
5C
6C MAIN STEERING ROUTINE FOR HADRON PACKAGE GHEISHA ***
7C THIS SUBROUTINE IS CALLED FROM NUCINT
8C
9C ORIGIN : F.CARMINATI, H.FESEFELDT (ROUTINE GHESIG)
10C REDESIGN: P. GABRIEL IK1 FZK KARLSRUHE
11C-----------------------------------------------------------------------
12
13*KEEP,CGCOMP.
14 PARAMETER (KK=3)
15 COMMON/CGCOMP/ ACOMP,ZCOMP,WCOMP
16 REAL ACOMP(KK),ZCOMP(KK),WCOMP(KK)
17*KEEP,ELABCT.
18 COMMON /ELABCT/ ELCUT
19 DOUBLE PRECISION ELCUT(4)
20*KEEP,ELADPM.
21 COMMON /ELADPM/ ELMEAN,ELMEAA,IELDPM,IELDPA
22 DOUBLE PRECISION ELMEAN(37),ELMEAA(37)
23 INTEGER IELDPM(37,13),IELDPA(37,13)
24*KEEP,ELASTY.
25 COMMON /ELASTY/ ELAST,IELIS,IELHM,IELNU,IELPI
26 DOUBLE PRECISION ELAST
27 INTEGER IELIS(20),IELHM(20),IELNU(20),IELPI(20)
28*KEEP,GENER.
29 COMMON /GENER/ GEN,ALEVEL
30 DOUBLE PRECISION GEN,ALEVEL
31*KEEP,MULT.
32 COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
33 DOUBLE PRECISION EKINL
34 INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
35*KEEP,PAM.
36 COMMON /PAM/ PAMA,SIGNUM
37 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
38*KEEP,PARPAR.
39 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
40 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
41 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
42 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
43 INTEGER ITYPE,LEVL
44*KEEP,RUNPAR.
45 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
46 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
47 * MONIOU,MDEBUG,NUCNUC,
48 * CETAPE,
49 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
50 * N1STTR,MDBASE,
51 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
52 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
53 * ,GHEISH,GHESIG
54 COMMON /RUNPAC/ DSN,HOST,USER
55 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
56 REAL STEPFC
57 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
58 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
59 * N1STTR,MDBASE
60 INTEGER CETAPE
61 CHARACTER*79 DSN
62 CHARACTER*20 HOST,USER
63
64 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
65 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
66 * ,GHEISH,GHESIG
67*KEND.
68
69 DOUBLE PRECISION ELASTI,ELABOR,PLX,PLY,PLZ,PLSQ,PLTOT,RMASSK
70
71 COMMON/GSECTI/ AIEL(20),AIIN(20),AIFI(20),AICA(20),ALAM,K0FLAG
72 INTEGER K0FLAG
73 REAL AIEL,AIIN,AIFI,AICA,ALAM
74
75C --- GHEISHA COMMONS ---
76 PARAMETER (MXGKGH=100)
77 PARAMETER (MXGKPV=MXGKGH)
78 COMMON /VECUTY/ PV(10,MXGKPV)
79
80 COMMON/CONSTS/ PI,TWPI,PIBTW,MP,MPI,MMU,MEL,MKCH,MK0,SMP,SMPI,
81 $ SMU,CT,CTKCH,CTK0,
82 $ ML0,MSP,MS0,MSM,MX0,MXM,CTL0,CTSP,CTSM,CTX0,CTXM,
83 $ RMASS(35),RCHARG(35)
84
85 REAL MP,MPI,MMU,MEL,MKCH,MK0,
86 * ML0,MSP,MS0,MSM,MX0,MXM
87
88 PARAMETER (MXEVEN=12*MXGKGH)
89 COMMON/EVENT / NSIZE,NCUR,NEXT,NTOT,EVE(MXEVEN)
90
91 COMMON/PRNTFL/INBCD,NEWBCD,INBIN,NEWBIN,NPEVT,NEVTP,LPRT,NPRT(10)
92 LOGICAL LPRT,NPRT
93
94
95C --- "NEVENT" CHANGED TO "KEVENT" IN COMMON /CURPAR/ DUE TO CLASH ---
96C --- WITH VARIABLE "NEVENT" IN GEANT COMMON ---
97
98 PARAMETER (MXGKCU=MXGKGH)
99 COMMON /CURPAR/ WEIGHT(10),DDELTN,IFILE,IRUN,NEVT,KEVENT,SHFLAG,
100 $ ITHST,ITTOT,ITLST,IFRND,TOFCUT,CMOM(5),CENG(5),
101 $ RS,S,ENP(10),NP,NM,NN,NR,NO,NZ,IPA(MXGKCU),
102 $ ATNO2,ZNO2
103
104C --- "IPART" CHANGED TO "KPART" IN COMMON /RESULT/ DUE TO CLASH ---
105C --- WITH VARIABLE "IPART" IN GEANT COMMON ---
106
107 COMMON /RESULT/ XEND,YEND,ZEND,RCA,RCE,AMAS,NCH,TOF,PX,PY,PZ,
108 $ USERW,INTCT,P,EN,EK,AMASQ,DELTN,ITK,NTK,KPART,IND,
109 $ LCALO,ICEL,SINL,COSL,SINP,COSP,
110 $ XOLD,YOLD,ZOLD,POLD,PXOLD,PYOLD,PZOLD,
111 $ XSCAT,YSCAT,ZSCAT,PSCAT,PXSCAT,PYSCAT,PZSCAT
112 REAL NCH,INTCT
113
114C --- "ABSL(21)" CHANGED TO "ABSLTH(21)" IN COMMON /MAT/ DUE TO CLASH ---
115C --- WITH VARIABLE "ABSL" IN GEANT COMMON ---
116
117 COMMON /MAT/ LMAT,
118 $ DEN(21),RADLTH(21),ATNO(21),ZNO(21),ABSLTH(21),
119 $ CDEN(21),MDEN(21),X0DEN(21),X1DEN(21),RION(21),
120 $ MATID(21),MATID1(21,24),PARMAT(21,10),
121 $ IFRAT,IFRAC(21),FRAC1(21,10),DEN1(21,10),
122 $ ATNO1(21,10),ZNO1(21,10)
123
124 DIMENSION IPELOS(35)
125 REAL EMAX,EEESQ
126 SAVE IDEOL
127
128 DIMENSION RNDM(1)
129
130C --- DIMENSION STMTS. FOR GEANT/GHEISHA PARTICLE CODE CONVERSIONS ---
131C --- KIPART(I)=GHEISHA CODE CORRESPONDING TO GEANT CODE I ---
132C --- IKPART(I)=GEANT CODE CORRESPONDING TO GHEISHA CODE I ---
133
134 DIMENSION KIPART(48),IKPART(35)
135
136C --- DATA STMTS. FOR GEANT/GHEISHA PARTICLE CODE CONVERSIONS ---
137C --- KIPART(I)=GHEISHA CODE CORRESPONDING TO GEANT CODE I ---
138C --- IKPART(I)=GEANT CODE CORRESPONDING TO GHEISHA CODE I ---
139
140 DATA KIPART/
141 $ 1, 3, 4, 2, 5, 6, 8, 7,
142 $ 9, 12, 10, 13, 16, 14, 15, 11,
143 $ 35, 18, 20, 21, 22, 26, 27, 33,
144 $ 17, 19, 23, 24, 25, 28, 29, 34,
145 $ 35, 35, 35, 35, 35, 35, 35, 35,
146 $ 35, 35, 35, 35, 30, 31, 32, 35/
147
148 DATA IKPART/
149 $ 1, 4, 2, 3, 5, 6, 8, 7,
150 $ 9, 11, 16, 10, 12, 14, 15, 13,
151 $ 25, 18, 26, 19, 20, 21, 27, 28,
152 $ 29, 22, 23, 30, 31, 45, 46, 47,
153 $ 24, 32, 48/
154
155
156C --- DENOTE STABLE PARTICLES ACCORDING TO GHEISHA CODE ---
157C --- STABLE : GAMMA, NEUTRINO, ELECTRON, PROTON AND HEAVY FRAGMENTS ---
158C --- WHEN STOPPING THESE PARTICLES ONLY LOOSE THEIR KINETIC ENERGY ---
159 DATA IPELOS/
160 $ 1, 1, 0, 1, 0, 0, 0, 0,
161 $ 0, 0, 0, 0, 0, 1, 0, 0,
162 $ 0, 0, 0, 0, 0, 0, 0, 0,
163 $ 0, 0, 0, 0, 0, 1, 1, 1,
164 $ 0, 0, 1/
165
166C --- LOWERBOUND OF KINETIC ENERGY BIN IN N CROSS-SECTION TABLES ---
167 DATA TEKLOW /0.0001/
168
169C --- KINETIC ENERGY TO SWITCH FROM "CASN" TO "GNSLWD" FOR N CASCADE ---
170 DATA SWTEKN /0.05/
171
172 DATA IDEOL/0/
173C-----------------------------------------------------------------------
174
175 IF ( DEBUG ) WRITE(MDEBUG,*) 'CGHEI :'
176
177C --- DEFINE PARTICLE TYPE
178 IF ( ITYPE .LE. 48 ) THEN
179 IPART = ITYPE
180 ELSEIF ( ITYPE .EQ. 201 ) THEN
181 IPART = 45
182 ELSEIF ( ITYPE .EQ. 301 ) THEN
183 IPART = 46
184 ELSEIF ( ITYPE .EQ. 402 ) THEN
185 IPART = 47
186 ELSE
187 WRITE (MONIOU,7795) ITYPE
188 7795 FORMAT (//,' *CGHEI* ILLEGAL PARTICLE TYPE OCCURS =',I5)
189 IPART = 48
190 ENDIF
191
192 NETEST=IKPART(KPART)
193 IF ( NETEST .EQ. IPART ) GO TO 9004
194
195 WRITE(MONIOU,8881) IPART,KPART
196 8881 FORMAT(' *CGHEI* IPART,KPART = ',2(I3,1X)/
197 $ ' *CGHEI* ======> PARTICLE TYPES DO NOT MATCH <=======')
198 STOP
199
200 9004 CONTINUE
201 KPART=KIPART(IPART)
202 KKPART=KPART
203
204C --- TRANSPORT THE TRACK NUMBER TO GHEISHA AND INITIALISE SOME NUMBERS
205C --- NTK=ITRA ITRA = CURRENT TRACK NUMBER IN GEANT (GCKINE)
206 NTK=0
207 INTCT=0.0
208 NEXT=1
209 NTOT=0
210 INT=0
211 TOF=0.0
212
213C --- STORE COORDINATES FOR SECONDARIES AND RESET ITYPE
214 SECPAR(1) = 0.
215 DO 7001 LK = 5, 8
216 SECPAR(LK) = CURPAR(LK)
217 7001 CONTINUE
218
219
220C --- FILL RESULT COMMON FOR THIS TRACK WITH CORSIKA VALUES ---
221
222 AMAS=RMASS(KPART)
223 NCH=RCHARG(KPART)
224 107 XEND = CURPAR(7)
225 YEND = CURPAR(8)
226 ZEND = CURPAR(5)
227 SINL = -CURPAR(3)
228 PHI = CURPAR(4)
229 USERW=0.0
230
231 AMASQ=AMAS*AMAS
232 EN = CURPAR(2) * ABS(AMAS)
233 EK = ABS ( EN - ABS(AMAS) )
234 ENOLD=EN
235 EMAX = 0.
236 P = SQRT ( EN*EN - AMASQ )
237 ELABOR = EN
238
239 SINP = SIN(PHI)
240 COSP = COS(PHI)
241 COSL = SQRT ( ABS(1.-SINL**2) )
242
243 PX = COSL * COSP
244 PY = COSL * SINP
245 PZ = SINL
246
247C --- SET GHEISHA INDEX FOR THE CURRENT MEDIUM ALWAYS TO 1 ---
248 IND=1
249
250C --- TRANSFER GLOBAL MATERIAL CONSTANTS FOR CURRENT MEDIUM ---
251C --- DETAILED DATA FOR COMPOUNDS IS OBTAINED VIA ROUTINE COMPO ---
252 ATNO(IND+1) = 14.56
253 ZNO(IND+1) = 7.265
254 DEN(IND+1) = 0.0
255 RADLTH(IND+1)= 0.0
256 ABSLTH(IND+1)= 0.0
257
258C --- SETUP PARMAT FOR PHYSICS STEERING ---
259 PARMAT(IND+1,10)=0.0
260
261 5 CONTINUE
262
263C --- INDICATE LIGHT (<= PI) AND HEAVY PARTICLES (HISTORICALLY) ---
264C --- CALIM CODE ---
265 J=2
266 TEST=RMASS(7)-0.001
267 IF (ABS(AMAS) .LT. TEST) J=1
268
269C *** DIVISION INTO VARIOUS INTERACTION CHANNELS DENOTED BY "INT" ***
270C THE CONVENTION FOR "INT" IS THE FOLLOWING
271
272C INT = -1 REACTION CROSS SECTIONS NOT YET TABULATED/PROGRAMMED
273C = 0 NO INTERACTION
274C = 1 ELEASTIC SCATTERING
275C = 2 INELASTIC SCATTERING
276C = 3 NUCLEAR FISSION WITH INELEASTIC SCATTERING
277C = 4 NEUTRON CAPTURE
278C INT = 3, 4 SHOULD BE DELETED FOR AIR TARGET
279
280C --- INTACT CODE ---
281 ALAM1=0.0
282 CALL GRNDM(RNDM,1)
283 RAT=RNDM(1)*ALAM
284 ATNO2 = 14.56
285 ZNO2 = 7.265
286
287 DO 6 K=1,KK
288 ATNO2 = ACOMP(K)
289 ZNO2 = ZCOMP(K)
290
291C --- TRY FOR ELASTIC SCATTERING ---
292 INT=1
293 ALAM1=ALAM1+AIEL(K)
294 IF (RAT .LT. ALAM1) GO TO 8
295
296C --- TRY FOR INELASTIC SCATTERING ---
297 INT=2
298 ALAM1=ALAM1+AIIN(K)
299 IF (RAT .LT. ALAM1) GO TO 8
300
301C --- TRY FOR NEUTRON CAPTURE ---
302 INT=4
303 ALAM1=ALAM1+AICA(K)
304 IF (RAT .LT. ALAM1) GO TO 8
305
306 6 CONTINUE
307
308C --- NO REACTION SELECTED ==> ELASTIC SCATTERING ---
309 INT=1
310
311C *** TAKE ACTION ACCORDING TO SELECTED REACTION CHANNEL ***
312C --- FOLLOWING CODE IS A TRANSLATION OF "CALIM" INTO GEANT JARGON ---
313
314 8 CONTINUE
315 IF (NPRT(9)) WRITE(MDEBUG,1001) INT
316 1001 FORMAT(' *CGHEI* INTERACTION TYPE CHOSEN INT = ',I3)
317
318 IF (INT .NE. 4) GO TO 10
319
320C --- NEUTRON CAPTURE ---
321 IF (NPRT(9)) WRITE(MDEBUG,2000)
322 2000 FORMAT(' *CGHEI* ROUTINE CAPTUR WILL BE CALLED')
323 CALL CAPTUR(NOPT)
324 GO TO 40
325
326 10 CONTINUE
327
328C --- ELASTIC AND INELASTIC SCATTERING ---
329 PV(1,MXGKPV)=P*PX
330 PV(2,MXGKPV)=P*PY
331 PV(3,MXGKPV)=P*PZ
332 PV(4,MXGKPV)=EN
333 PV(5,MXGKPV)=AMAS
334 PV(6,MXGKPV)=NCH
335 PV(7,MXGKPV)=TOF
336 PV(8,MXGKPV)=KPART
337 PV(9,MXGKPV)=0.
338 PV(10,MXGKPV)=USERW
339
340C --- ADDITIONAL PARAMETERS TO SIMULATE FERMI MOTION AND EVAPORATION ---
341 DO 111 JENP=1,10
342 ENP(JENP)=0.
343 111 CONTINUE
344 ENP(5)=EK
345 ENP(6)=EN
346 ENP(7)=P
347
348 IF (INT .NE. 1) GO TO 12
349
350C *** ELASTIC SCATTERING PROCESSES ***
351
352C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS ---
353 IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35
354
355C --- NORMAL ELASTIC SCATTERING FOR LIGHT MEDIA ---
356 IF (ATNO2 .LT. 1.5) GO TO 35
357
358C --- COHERENT ELASTIC SCATTERING FOR HEAVY MEDIA ---
359 IF (NPRT(9)) WRITE(MDEBUG,2002)
360 2002 FORMAT(' *CGHEI* ROUTINE COSCAT WILL BE CALLED')
361 CALL COSCAT
362 GO TO 40
363
364C *** NON-ELASTIC SCATTERING PROCESSES ***
365 12 CONTINUE
366
367C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS ---
368 IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35
369
370C *** USE SOMETIMES NUCLEAR REACTION ROUTINE "NUCREC" FOR LOW ENERGY ***
371C *** PROTON AND NEUTRON SCATTERING ***
372 CALL GRNDM(RNDM,1)
373 TEST1=RNDM(1)
374 TEST2=4.5*(EK-0.01)
375 IF ((KPART .EQ. 14) .AND. (TEST1 .GT. TEST2)) GO TO 85
376 IF ((KPART .EQ. 16) .AND. (TEST1 .GT. TEST2)) GO TO 86
377
378C *** FERMI MOTION AND EVAPORATION ***
379 TKIN=CINEMA(EK)
380 PV(9,MXGKPV) = TKIN
381 ENP(5)=EK+TKIN
382C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
383 IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
384 ENP(6)=ENP(5)+ABS(AMAS)
385 ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
386 ENP(7)=SQRT(ABS(ENP(7)))
387 TKIN=FERMI(ENP(5))
388 ENP(5)=ENP(5)+TKIN
389C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
390 IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
391 ENP(6)=ENP(5)+ABS(AMAS)
392 ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
393 ENP(7)=SQRT(ABS(ENP(7)))
394 TKIN=EXNU(ENP(5))
395 ENP(5)=ENP(5)-TKIN
396C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
397 IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
398 ENP(6)=ENP(5)+ABS(AMAS)
399 ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
400 ENP(7)=SQRT(ABS(ENP(7)))
401
402C *** IN CASE OF ENERGY ABOVE CUT-OFF LET THE PARTICLE CASCADE ***
403 IF ( ENP(5) .GT. ELCUT(1)) GO TO 35
404
405C --- SECOND CHANCE FOR ANTI-BARYONS DUE TO POSSIBLE ANNIHILATION ---
406 IF ((AMAS .GE. 0.0) .OR. (KPART .LE. 14)) GO TO 13
407 ANNI=1.3*P
408 IF (ANNI .GT. 0.4) ANNI=0.4
409 CALL GRNDM(RNDM,1)
410 TEST=RNDM(1)
411 IF (TEST .GT. ANNI) GO TO 35
412
413C *** PARTICLE WITH ENERGY BELOW CUT-OFF ***
414C --- ==> ONLY NUCLEAR EVAPORATION AND QUASI-ELASTIC SCATTERING ---
415 13 CONTINUE
416
417 IF (NPRT(9)) WRITE(MDEBUG,1002) KPART,EK,EN,P,ENP(5),ENP(6),ENP(7)
418 1002 FORMAT(' *CGHEI* ENERGY BELOW CUT-OFF FOR GHEISHA PARTICLE ',I3/
419 $ ' EK,EN,P,ENP(5),ENP(6),ENP(7) = ',6(G12.5,1X))
420
421 IF ((KPART .NE. 14) .AND. (KPART .NE. 16)) GO TO 14
422 IF (KPART .EQ. 16) GO TO 86
423
424C --- SLOW PROTON ---
425 85 CONTINUE
426 IF (NPRT(9)) WRITE(MDEBUG,2003) EK,KPART
427 2003 FORMAT(' *CGHEI* ROUTINE NUCREC WILL BE CALLED',
428 $ ' EK = ',G12.5,' GEV KPART = ',I3)
429 CALL NUCREC(NOPT,2)
430
431 IF (NOPT .NE. 0) GO TO 50
432
433 IF (NPRT(9)) WRITE(MDEBUG,2004)EK,KPART
434 2004 FORMAT(' *CGHEI* ROUTINE COSCAT WILL BE CALLED',
435 $ ' EK = ',G12.5,' GEV KPART = ',I3)
436 CALL COSCAT
437 GO TO 40
438
439C --- SLOW NEUTRON ---
440 86 CONTINUE
441 IF (NPRT(9)) WRITE(MDEBUG,2015)
442 NUCFLG=0
443 CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW)
444 IF (NUCFLG .NE. 0) GO TO 50
445 GO TO 40
446
447C --- OTHER SLOW PARTICLES ---
448 14 CONTINUE
449 IPA(1)=KPART
450C --- DECIDE FOR PROTON OR NEUTRON TARGET ---
451 IPA(2)=16
452 CALL GRNDM(RNDM,1)
453 TEST1=RNDM(1)
454 TEST2=ZNO2/ATNO2
455 IF (TEST1 .LT. TEST2) IPA(2)=14
456 AVERN=0.0
457 NFL=1
458 IF (IPA(2) .EQ. 16) NFL=2
459 IPPP=KPART
460 IF (NPRT(9)) WRITE(MDEBUG,2005)
461 2005 FORMAT(' *CGHEI* ROUTINE TWOB WILL BE CALLED')
462 CALL TWOB(IPPP,NFL,AVERN)
463 GOTO 40
464
465C --- INITIALISATION OF CASCADE QUANTITIES ---
466 35 CONTINUE
467
468C *** CASCADE GENERATION ***
469C --- CALCULATE FINAL STATE MULTIPLICITY AND LONGITUDINAL AND ---
470C --- TRANSVERSE MOMENTUM DISTRIBUTIONS ---
471
472C --- FIXED PARTICLE TYPE TO STEER THE CASCADE ---
473 KKPART=KPART
474
475C --- NO CASCADE FOR LEPTONS ---
476 IF (KKPART .LE. 6) GO TO 9999
477
478C *** WHAT TO DO WITH "NEW PARTICLES" FOR GHEISHA ?????? ***
479C --- RETURN FOR THE TIME BEING ---
480 IF (KKPART .GE. 35) GO TO 9999
481
482C --- CASCADE OF HEAVY FRAGMENTS
483 IF ((KKPART .GE. 30) .AND. (KKPART .LE. 32)) GO TO 390
484
485C --- INITIALIZE THE IPA ARRAY ---
486 CALL VZERO(IPA(1),MXGKCU)
487
488C --- CASCADE OF OMEGA - AND OMEGA - BAR ---
489 IF (KKPART .EQ. 33) GO TO 330
490 IF (KKPART .EQ. 34) GO TO 331
491
492 NVEPAR=KKPART-17
493 IF (NVEPAR .LE. 0) GO TO 15
494 GO TO (318,319,320,321,322,323,324,325,326,327,328,329),NVEPAR
495
496 15 CONTINUE
497 NVEPAR=KKPART-6
498 GO TO (307,308,309,310,311,312,313,314,315,316,317,318),NVEPAR
499
500C --- PI+ CASCADE ---
501 307 CONTINUE
502 IF (NPRT(9)) WRITE(MDEBUG,2006)
503 2006 FORMAT(' *CGHEI* ROUTINE CASPIP WILL BE CALLED')
504 CALL CASPIP(J,INT,NFL)
505 GO TO 40
506
507C --- PI0 ==> NO CASCADE ---
508 308 CONTINUE
509 GO TO 40
510
511C --- PI- CASCADE ---
512 309 CONTINUE
513 IF (NPRT(9)) WRITE(MDEBUG,2007)
514 2007 FORMAT(' *CGHEI* ROUTINE CASPIM WILL BE CALLED')
515 CALL CASPIM(J,INT,NFL)
516 GO TO 40
517
518C --- K+ CASCADE ---
519 310 CONTINUE
520 IF (NPRT(9)) WRITE(MDEBUG,2008)
521 2008 FORMAT(' *CGHEI* ROUTINE CASKP WILL BE CALLED')
522 CALL CASKP(J,INT,NFL)
523 GO TO 40
524
525C --- K0 CASCADE ---
526 311 CONTINUE
527 IF (NPRT(9)) WRITE(MDEBUG,2009)
528 2009 FORMAT(' *CGHEI* ROUTINE CASK0 WILL BE CALLED')
529 CALL CASK0(J,INT,NFL)
530 GO TO 40
531
532C --- K0 BAR CASCADE ---
533 312 CONTINUE
534 IF (NPRT(9)) WRITE(MDEBUG,2010)
535 2010 FORMAT(' *CGHEI* ROUTINE CASK0B WILL BE CALLED')
536 CALL CASK0B(J,INT,NFL)
537 GO TO 40
538
539C --- K- CASCADE ---
540 313 CONTINUE
541 IF (NPRT(9)) WRITE(MDEBUG,2011)
542 2011 FORMAT(' *CGHEI* ROUTINE CASKM WILL BE CALLED')
543 CALL CASKM(J,INT,NFL)
544 GO TO 40
545
546C --- PROTON CASCADE ---
547 314 CONTINUE
548 IF (NPRT(9)) WRITE(MDEBUG,2012)
549 2012 FORMAT(' *CGHEI* ROUTINE CASP WILL BE CALLED')
550 CALL CASP(J,INT,NFL)
551 GO TO 40
552
553C --- PROTON BAR CASCADE ---
554 315 CONTINUE
555 IF (NPRT(9)) WRITE(MDEBUG,2013)
556 2013 FORMAT(' *CGHEI* ROUTINE CASPB WILL BE CALLED')
557 CALL CASPB(J,INT,NFL)
558 GO TO 40
559
560C --- NEUTRON CASCADE ---
561 316 CONTINUE
562 NUCFLG=0
563 IF (EK .GT. SWTEKN) THEN
564 CALL CASN(J,INT,NFL)
565 IF (NPRT(9)) WRITE(MDEBUG,2014)
566 2014 FORMAT(' *CGHEI* ROUTINE CASN WILL BE CALLED')
567 ELSE
568 CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW)
569 IF (NPRT(9)) WRITE(MDEBUG,2015)
570 2015 FORMAT(' *CGHEI* ROUTINE GNSLWD WILL BE CALLED')
571 ENDIF
572 IF (NUCFLG .NE. 0) GO TO 50
573 GO TO 40
574
575C --- NEUTRON BAR CASCADE ---
576 317 CONTINUE
577 IF (NPRT(9)) WRITE(MDEBUG,2016)
578 2016 FORMAT(' *CGHEI* ROUTINE CASNB WILL BE CALLED')
579 CALL CASNB(J,INT,NFL)
580 GO TO 40
581
582C --- LAMBDA CASCADE ---
583 318 CONTINUE
584 IF (NPRT(9)) WRITE(MDEBUG,2017)
585 2017 FORMAT(' *CGHEI* ROUTINE CASL0 WILL BE CALLED')
586 CALL CASL0(J,INT,NFL)
587 GO TO 40
588
589C --- LAMBDA BAR CASCADE ---
590 319 CONTINUE
591 IF (NPRT(9)) WRITE(MDEBUG,2018)
592 2018 FORMAT(' *CGHEI* ROUTINE CASAL0 WILL BE CALLED')
593 CALL CASAL0(J,INT,NFL)
594 GO TO 40
595
596C --- SIGMA + CASCADE ---
597 320 CONTINUE
598 IF (NPRT(9)) WRITE(MDEBUG,2019)
599 2019 FORMAT(' *CGHEI* ROUTINE CASSP WILL BE CALLED')
600 CALL CASSP(J,INT,NFL)
601 GO TO 40
602
603C --- SIGMA 0 ==> NO CASCADE ---
604 321 CONTINUE
605 GO TO 40
606
607C --- SIGMA - CASCADE ---
608 322 CONTINUE
609 IF (NPRT(9)) WRITE(MDEBUG,2020)
610 2020 FORMAT(' *CGHEI* ROUTINE CASSM WILL BE CALLED')
611 CALL CASSM(J,INT,NFL)
612 GO TO 40
613
614C --- SIGMA + BAR CASCADE ---
615 323 CONTINUE
616 IF (NPRT(9)) WRITE(MDEBUG,2021)
617 2021 FORMAT(' *CGHEI* ROUTINE CASASP WILL BE CALLED')
618 CALL CASASP(J,INT,NFL)
619 GO TO 40
620
621C --- SIGMA 0 BAR ==> NO CASCADE ---
622 324 CONTINUE
623 GO TO 40
624
625C --- SIGMA - BAR CASCADE ---
626 325 CONTINUE
627 IF (NPRT(9)) WRITE(MDEBUG,2022)
628 2022 FORMAT(' *CGHEI* ROUTINE CASASM WILL BE CALLED')
629 CALL CASASM(J,INT,NFL)
630 GO TO 40
631
632C --- XI 0 CASCADE ---
633 326 CONTINUE
634 IF (NPRT(9)) PRINT 2023
635 2023 FORMAT(' *CGHEI* ROUTINE CASX0 WILL BE CALLED')
636 CALL CASX0(J,INT,NFL)
637 GO TO 40
638
639C --- XI - CASCADE ---
640 327 CONTINUE
641 IF (NPRT(9)) PRINT 2024
642 2024 FORMAT(' *CGHEI* ROUTINE CASXM WILL BE CALLED')
643 CALL CASXM(J,INT,NFL)
644 GO TO 40
645
646C --- XI 0 BAR CASCADE ---
647 328 CONTINUE
648 IF (NPRT(9)) PRINT 2025
649 2025 FORMAT(' *CGHEI* ROUTINE CASAX0 WILL BE CALLED')
650 CALL CASAX0(J,INT,NFL)
651 GO TO 40
652
653C --- XI - BAR CASCADE ---
654 329 CONTINUE
655 IF (NPRT(9)) PRINT 2026
656 2026 FORMAT(' *CGHEI* ROUTINE CASAXM WILL BE CALLED')
657 CALL CASAXM(J,INT,NFL)
658 GO TO 40
659
660C --- OMEGA - CASCADE ---
661 330 CONTINUE
662 IF (NPRT(9)) PRINT 2027
663 2027 FORMAT(' *CGHEI* ROUTINE CASOM WILL BE CALLED')
664 CALL CASOM(J,INT,NFL)
665 GO TO 40
666
667C --- OMEGA - BAR CASCADE ---
668 331 CONTINUE
669 IF (NPRT(9)) PRINT 2028
670 2028 FORMAT(' *CGHEI* ROUTINE CASAOM WILL BE CALLED')
671 CALL CASAOM(J,INT,NFL)
672 GO TO 40
673
674C --- HEAVY FRAGMENT CASCADE ---
675 390 CONTINUE
676 IF (NPRT(9)) WRITE(MDEBUG,2090)
677 2090 FORMAT(' *CGHEI* ROUTINE CASFRG WILL BE CALLED')
678 NUCFLG=0
679 CALL CASFRG(NUCFLG,INT,NFL)
680 IF (NUCFLG .NE. 0) GO TO 50
681
682C *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
683 40 CONTINUE
684 IF ((NTOT .NE. 0) .OR. (KKPART .NE. KPART)) GO TO 50
685
686 50 CONTINUE
687
688 NVEDUM=KIPART(IPART)
689 IF (NPRT(9)) WRITE(MDEBUG,1004)NTOT,IPART,KPART,KKPART,NVEDUM
690 1004 FORMAT(' *CGHEI* SEC. GEN. NTOT,IPART,KPART,KKPART,KIPART = ',
691 $ 5(I3,1X))
692
693C --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
694C --- THE TEMPORARY STACK ---
695
696C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT ---
697 IF ((KPART .NE. 11) .AND. (KPART .NE. 12)) GO TO 52
698 CALL GRNDM(RNDM,1)
699 KPART=11.5+RNDM(1)
700
701 52 CONTINUE
702
703C --- IN CASE THE NEW PARTICLE IS A NEUTRINO ==> FORGET IT ---
704 IF (KPART .EQ. 2) GO TO 60
705
706C --- PUT CURRENT GHEISHA PARTICLE ON THE CORSIKA STACK
707C --- ( IF SURVIVING ANGLE CUT ! )
708 NGKINE = 1
709 SECPAR(3) = -PZ
710
711C --- CALCULATE ELASTICITY
712 IF ( EN .GT. EMAX ) THEN
713 EMAX = EN
714 ENDIF
715
716 IF ( SECPAR(3) .GT. C(29) ) THEN
717
718 ITY=IKPART(KPART)
719 IF ( ITY .LT. 45 ) THEN
720 SECPAR(1) = DBLE(ITY)
721 ELSEIF ( ITY .EQ. 45 ) THEN
722 SECPAR(1) = 201.D0
723 ELSEIF ( ITY .EQ. 46 ) THEN
724 SECPAR(1) = 301.D0
725 ELSEIF ( ITY .EQ. 47 ) THEN
726 SECPAR(1) = 402.D0
727 ENDIF
728 IF ( ABS(AMAS) .LT. 1.E-9 ) THEN
729 SECPAR(2) = EN
730 ELSE
731 SECPAR(2) = DBLE(EN) / DBLE(ABS(AMAS))
732 ENDIF
733 IF ( PX .NE. 0. .OR. PY .NE. 0. ) THEN
734 SECPAR(4) = ATAN2 ( DBLE(PY) , DBLE(PX) )
735 ELSE
736 SECPAR(4) = 0.D0
737 ENDIF
738
739 CALL TSTACK
740
741 ENDIF
742
743C *** CHECK WHETHER SECONDARIES HAVE BEEN GENERATED AND COPY THEM ***
744C *** ALSO ON THE GEANT STACK ***
745 60 CONTINUE
746
747C --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
748C --- CONVENTION IS THE FOLLOWING ---
749C
750C EVE(INDEX+ 1)= X
751C EVE(INDEX+ 2)= Y
752C EVE(INDEX+ 3)= Z
753C EVE(INDEX+ 4)= NCAL
754C EVE(INDEX+ 5)= NCELL
755C EVE(INDEX+ 6)= MASS
756C EVE(INDEX+ 7)= CHARGE
757C EVE(INDEX+ 8)= TOF
758C EVE(INDEX+ 9)= PX
759C EVE(INDEX+10)= PY
760C EVE(INDEX+11)= PZ
761C EVE(INDEX+12)= TYPE
762
763 IF (NTOT .LE. 0) GO TO 9999
764
765C --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
766 DO 61 L=1,NTOT
767 INDEX=(L-1)*12
768 JND=EVE(INDEX+12)
769
770C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT ---
771 IF ((JND .NE. 11) .AND. (JND .NE. 12)) GO TO 63
772 CALL GRNDM(RNDM,1)
773 JND=11.5+RNDM(1)
774
775C --- FORGET ABOUT NEUTRINOS ---
776 63 CONTINUE
777 IF (JND .EQ. 2) GO TO 61
778
779C --- SWITCH TO CORSIKA QUANTITIES ---
780 ITY=IKPART(JND)
781 IF (NPRT(9)) WRITE(MDEBUG,1006) ITY,NGKINE,L,(EVE(INDEX+J),J=1,12)
782 1006 FORMAT(' *CGHEI* GEANT PART. ',I3,' ALSO PUT ONTO STACK AT',
783 $ ' POS. ',I3/
784 $ ' EVE(',I2,') = ',(' ',10G12.5))
785
786 PLX = EVE(INDEX+9)
787 PLY = EVE(INDEX+10)
788 PLZ = EVE(INDEX+11)
789 PLSQ = PLX**2 + PLY**2 + PLZ**2
790 PLTOT = SQRT (PLSQ)
791 RMASSK = ABS(RMASS(JND))
792
793C FIND HIGHEST ENERGY PARTICLE FOR ELASTICITY
794 EEESQ = PLSQ + RMASSK**2
795 IF ( EEESQ .GT. EMAX**2 ) THEN
796 EMAX = SQRT(EEESQ)
797 ENDIF
798
799C --- APPLY ANGLE CUT AND
800C --- ADD PARTICLE TO THE CORSIKA STACK (RESTRICTED TO 100) ---
801 IF ( PLTOT .LE. 1.D-10 ) GOTO 61
802 SECPAR(3) = -PLZ / PLTOT
803 IF ( SECPAR(3) .LE. C(29) ) GOTO 61
804
805 IF (NGKINE .GE. MXGKGH) GO TO 9999
806 NGKINE=NGKINE+1
807 IF ( ITY .LT. 45 ) THEN
808 SECPAR(1) = DBLE(ITY)
809 ELSEIF ( ITY .EQ. 45 ) THEN
810 SECPAR(1) = 201.D0
811 ELSEIF ( ITY .EQ. 46 ) THEN
812 SECPAR(1) = 301.D0
813 ELSEIF ( ITY .EQ. 47 ) THEN
814 SECPAR(1) = 402.D0
815 ELSE
816 SECPAR(1) = 0.D0
817 WRITE(MONIOU,*) '*CGHEI* ILLEGAL PARTICLE TYPE'
818 ENDIF
819 IF ( RMASSK .LT. 1.D-9 ) THEN
820 SECPAR(2) = PLTOT
821 ELSE
822 SECPAR(2) = SQRT (PLSQ+RMASSK**2) / RMASSK
823 ENDIF
824 IF ( PLX .NE. 0.D0 .OR. PLY .NE. 0.D0 ) THEN
825 SECPAR(4) = ATAN2 ( PLY,PLX )
826 ELSE
827 SECPAR(4) = 0.D0
828 ENDIF
829
830 CALL TSTACK
831
832 61 CONTINUE
833
834C --- COUNTER FOR ENERGY-MULTIPLICITY MATRIX
835 MSMM = MSMM + NTOT
836
837C --- FILL ELASTICITY IN MATRICES
838 ELASTI = EMAX/ENOLD
839 MELL = MIN ( 1.D0+10.D0* MAX( 0.D0, ELASTI ) , 11.D0 )
840 MEN = MIN ( 4.D0+ 3.D0*LOG10(MAX( .1D0, EKINL )), 37.D0 )
841 IELDPM(MEN,MELL) = IELDPM(MEN,MELL) + 1
842 IELDPA(MEN,MELL) = IELDPA(MEN,MELL) + 1
843 IF ( ELASTI .LT. 1. ) THEN
844 ELMEAN(MEN) = ELMEAN(MEN) + ELASTI
845 ELMEAA(MEN) = ELMEAA(MEN) + ELASTI
846 ENDIF
847
848 9999 CONTINUE
849C --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
850 NGKINE=MIN(NGKINE,MXGKGH)
851
852 RETURN
853 END
Note: See TracBrowser for help on using the repository browser.