* * $Id$ * * $Log$ * Revision 1.1.1.1 1995/10/24 10:21:15 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.38 by S.Giani *-- Author : SUBROUTINE GHEISH C C *** MAIN STEERING FOR HADRON SHOWER DEVELOPMENT *** C *** NVE 15-JUN-1988 CERN GENEVA *** C C CALLED BY : GUHADR (USER ROUTINE) C ORIGIN : F.CARMINATI, H.FESEFELDT C ROUTINES : CALIM 16-SEP-1987 C SETRES 19-AUG-1985 C INTACT 06-OCT-1987 C #include "geant321/gcbank.inc" #include "geant321/gcjloc.inc" #include "geant321/gccuts.inc" #include "geant321/gcflag.inc" #include "geant321/gckine.inc" #include "geant321/gcking.inc" #include "geant321/gcmate.inc" #include "geant321/gcphys.inc" #include "geant321/gctmed.inc" #include "geant321/gctrak.inc" #include "geant321/gsecti.inc" #include "geant321/gcunit.inc" C --- GHEISHA COMMONS --- #include "geant321/mxgkgh.inc" #include "geant321/s_blankp.inc" #include "geant321/s_consts.inc" #include "geant321/s_event.inc" #include "geant321/s_prntfl.inc" C C --- "NEVENT" CHANGED TO "KEVENT" IN COMMON /CURPAR/ DUE TO CLASH --- C --- WITH VARIABLE "NEVENT" IN GEANT COMMON --- C PARAMETER (MXGKCU=MXGKGH) COMMON /CURPAR /WEIGHT(10),DDELTN,IFILE,IRUN,NEVT,KEVENT,SHFLAG, $ ITHST,ITTOT,ITLST,IFRND,TOFCUT,CMOM(5),CENG(5), $ RS,S,ENP(10),NP,NM,NN,NR,NO,NZ,IPA(MXGKCU), $ ATNO2,ZNO2 C C --- "IPART" CHANGED TO "KPART" IN COMMON /RESULT/ DUE TO CLASH --- C --- WITH VARIABLE "IPART" IN GEANT COMMON --- C COMMON /RESULT/ XEND,YEND,ZEND,RCA,RCE,AMAS,NCH,TOF,PX,PY,PZ, $ USERW,INTCT,P,EN,EK,AMASQ,DELTN,ITK,NTK,KPART,IND, $ LCALO,ICEL,SINL,COSL,SINP,COSP, $ XOLD,YOLD,ZOLD,POLD,PXOLD,PYOLD,PZOLD, $ XSCAT,YSCAT,ZSCAT,PSCAT,PXSCAT,PYSCAT,PZSCAT REAL NCH,INTCT C C --- "ABSL(21)" CHANGED TO "ABSLTH(21)" IN COMMON /MAT/ DUE TO CLASH --- C --- WITH VARIABLE "ABSL" IN GEANT COMMON --- C COMMON /MAT/ LMAT, $ DEN(21),RADLTH(21),ATNO(21),ZNO(21),ABSLTH(21), $ CDEN(21),MDEN(21),X0DEN(21),X1DEN(21),RION(21), $ MATID(21),MATID1(21,24),PARMAT(21,10), $ IFRAT,IFRAC(21),FRAC1(21,10),DEN1(21,10), $ ATNO1(21,10),ZNO1(21,10) C DIMENSION IPELOS(35) SAVE IDEOL C C --- TRANSFER GEANT CUT-OFFS INTO GHEISHA VALUES --- DIMENSION CUTS(5) EQUIVALENCE (CUTS(1),CUTGAM) DIMENSION RNDM(1) C #include "geant321/pcodim.inc" #include "geant321/pcodat.inc" C C --- DENOTE STABLE PARTICLES ACCORDING TO GHEISHA CODE --- C --- STABLE : GAMMA, NEUTRINO, ELECTRON, PROTON AND HEAVY FRAGMENTS --- C --- WHEN STOPPING THESE PARTICLES ONLY LOOSE THEIR KINETIC ENERGY --- DATA IPELOS/ $ 1, 1, 0, 1, 0, 0, 0, 0, $ 0, 0, 0, 0, 0, 1, 0, 0, $ 0, 0, 0, 0, 0, 0, 0, 0, $ 0, 0, 0, 0, 0, 1, 1, 1, $ 0, 0, 1/ C C --- LOWERBOUND OF KINETIC ENERGY BIN IN N CROSS-SECTION TABLES --- DATA TEKLOW /0.0001/ C C --- KINETIC ENERGY TO SWITCH FROM "CASN" TO "GNSLWD" FOR N CASCADE --- DATA SWTEKN /0.05/ C DATA IDEOL/0/ C C --- INITIALIZE RELEVANT GHEISHA VARIABLES IN CASE NOT DONE ALREADY --- IF (IFINIT(4) .EQ. 0) CALL GHEINI C C --- SET THE INTERACTION MECHANISM TO "HADR" --- KCASE=NAMEC(12) C C --- SET GHEISHA PRINTING FLAGS ACCORDING TO "DEBUG" STEERING CARD -- IF (IDEOL .EQ. IDEBUG) GO TO 9000 C IF (IDEBUG .NE. 1) GO TO 9001 C C --- SET SELECTED DEBUGGING FLAGS --- DO 9002 LL=1,10 IF ((ISWIT(LL) .LE. 100) .OR. (ISWIT(LL) .GT. 110)) GO TO 9002 JJ=ISWIT(LL)-100 NPRT(JJ)=.TRUE. 9002 CONTINUE GO TO 9000 C C --- NO DEBUGGING SELECTED --- 9001 CONTINUE DO 9003 LL=1,10 NPRT(LL)=.FALSE. 9003 CONTINUE IDEOL=IDEBUG C 9000 CONTINUE C C --- SET THE GHEISHA PARTICLE TYPE TO THE ONE OF GEANT --- IF(IPART.GT.48) THEN IF(ISTOP.EQ.0) GOTO 9999 JPA = LQ(JPART-IPART) AMAS=Q(JPA+7) NCH =Q(JPA+8) KPART=-IPART GOTO 107 ENDIF NETEST=IKPART(KPART) IF ((NETEST .EQ. IPART) .OR. (ISTOP .NE. 0)) GO TO 9004 C PRINT 8881,IPART,KPART,ISTOP 8881 FORMAT(' *GHEISH* IPART,KPART = ',2(I3,1X),' ISTOP = ',I3/ $ ' *GHEISH* ======> PARTICLE TYPES DO NOT MATCH <=======') STOP C 9004 CONTINUE KPART=KIPART(IPART) KKPART=KPART AMAS=RMASS(KPART) NCH=RCHARG(KPART) C C --- TRANSPORT THE TRACK NUMBER TO GHEISHA AND INITIALISE SOME NUMBERS 107 NTK=ITRA INTCT=0.0 NEXT=1 NTOT=0 TOF=0.0 C C --- FILL RESULT COMMON FOR THIS TRACK WITH GEANT VALUES --- C --- CALIM CODE --- XEND=VECT(1) YEND=VECT(2) ZEND=VECT(3) PX=VECT(4) PY=VECT(5) PZ=VECT(6) USERW=UPWGHT C --- SETRES CODE --- P=VECT(7) AMASQ=AMAS*AMAS EN=SQRT(AMASQ+P*P) EK=ABS(EN-ABS(AMAS)) ENOLD=EN C SINL=0.0 COSL=1.0 SINP=0.0 COSP=1.0 C IF (ABS(P) .LE. 1.0E-10) GO TO 1 SINL=PZ COSL=SQRT(ABS(1.0-SINL**2)) C 1 CONTINUE CALL GRNDM(RNDM,1) PHI=RNDM(1)*TWPI IF ((PX .EQ. 0.0) .AND. (PY .EQ. 0.0)) GOTO 3 IF (ABS(PX) .LT. 1.E-10) GOTO 2 PHI=ATAN2(PY,PX) GOTO 3 C 2 CONTINUE IF (PY .GT. 0.0) PHI=PI/2.0 IF (PY .LE. 0.0) PHI=3.0*PI/2.0 C 3 CONTINUE SINP=SIN(PHI) COSP=COS(PHI) C C --- SET GHEISHA INDEX FOR THE CURRENT MEDIUM ALWAYS TO 1 --- IND=1 C C --- TRANSFER GLOBAL MATERIAL CONSTANTS FOR CURRENT MEDIUM --- C --- DETAILED DATA FOR COMPOUNDS IS OBTAINED VIA ROUTINE COMPO --- ATNO(IND+1)=A ZNO(IND+1)=Z DEN(IND+1)=DENS RADLTH(IND+1)=RADL ABSLTH(IND+1)=ABSL C C --- SETUP PARMAT FOR PHYSICS STEERING --- PARMAT(IND+1,5)=0.0 PARMAT(IND+1,8)=IPFIS PARMAT(IND+1,9)=0.0 PARMAT(IND+1,10)=0.0 JTMN=LQ(JTM) IF (JTMN .LE. 0) GO TO 4 PARMAT(IND+1,5)=Q(JTMN+26) 4 CONTINUE C C --- CHECK WHETHER PARTICLE IS STOPPING OR NOT --- IF (ISTOP .EQ. 0) GO TO 5 C IF (NPRT(9)) PRINT 1000,KPART 1000 FORMAT(' *GHEISH* STOPPING GHEISHA PARTICLE ',I3) CALL GHSTOP C --- IN CASE OF DECAY OF PARTICLE OR USER PARTICLE ==> RETURN --- IF (LMEC(NMEC) .EQ. 5 .OR. KPART .LT. 0) GO TO 9999 C --- IN CASE OF HAD. INT. WITH GENERATION OF SEC. ==> GO TO 40 --- IF (IHADR .NE. 2) GO TO 40 C --- ALSO DEPOSIT REST MASS ENERGY FOR IN-STABLE PARTICLES --- IF (IPELOS(KPART) .EQ. 0) DESTEP=DESTEP+ABS(RMASS(KPART)) GO TO 9999 5 CONTINUE C C --- INDICATE LIGHT (<= PI) AND HEAVY PARTICLES (HISTORICALLY) --- C --- CALIM CODE --- J=2 TEST=RMASS(7)-0.001 IF (ABS(AMAS) .LT. TEST) J=1 C C *** DIVISION INTO VARIOUS INTERACTION CHANNELS DENOTED BY "INT" *** C THE CONVENTION FOR "INT" IS THE FOLLOWING C C INT = -1 REACTION CROSS SECTIONS NOT YET TABULATED/PROGRAMMED C = 0 NO INTERACTION C = 1 ELEASTIC SCATTERING C = 2 INELASTIC SCATTERING C = 3 NUCLEAR FISSION WITH INELEASTIC SCATTERING C = 4 NEUTRON CAPTURE C C --- INTACT CODE --- KK=ABS(Q(JMA+11)) ALAM1=0.0 CALL GRNDM(RNDM,1) RAT=RNDM(1)*ALAM NMEC=NMEC+1 ATNO2=A ZNO2 =Z C DO 6 K=1,KK IF (KK .LE. 0) GO TO 6 C IF (KK .EQ. 1) GO TO 7 ATNO2=Q(JMIXT+K) ZNO2 =Q(JMIXT+K+KK) C 7 CONTINUE C C --- TRY FOR ELASTIC SCATTERING --- INT=1 LMEC(NMEC)=13 ALAM1=ALAM1+AIEL(K) IF (RAT .LT. ALAM1) GO TO 8 C C --- TRY FOR INELASTIC SCATTERING --- INT=2 LMEC(NMEC)=20 ALAM1=ALAM1+AIIN(K) IF (RAT .LT. ALAM1) GO TO 8 C C --- TRY FOR NUCLEAR FISSION WITH INELASTIC SCATTERING --- INT=3 LMEC(NMEC)=15 ALAM1=ALAM1+AIFI(K) IF (RAT .LT. ALAM1) GO TO 8 C C --- TRY FOR NEUTRON CAPTURE --- INT=4 LMEC(NMEC)=18 ALAM1=ALAM1+AICA(K) IF (RAT .LT. ALAM1) GO TO 8 C 6 CONTINUE C --- NO REACTION SELECTED ==> ELASTIC SCATTERING --- INT=1 LMEC(NMEC)=13 C C *** TAKE ACTION ACCORDING TO SELECTED REACTION CHANNEL *** C --- FOLLOWING CODE IS A TRANSLATION OF "CALIM" INTO GEANT JARGON --- C 8 CONTINUE IF (NPRT(9)) PRINT 1001,INT 1001 FORMAT(' *GHEISH* INTERACTION TYPE CHOSEN INT = ',I3) C C --- IN CASE OF NO INTERACTION OR UNKNOWN CROSS SECTIONS ==> DONE --- IF (INT .LE. 0) GO TO 40 C C --- IN CASE OF NON-ELASTIC SCATTERING AND NO GENERATION OF SEC. --- C --- PARTICLES DEPOSIT TOTAL PARTICLE ENERGY AND RETURN --- IF ((INT .EQ. 1) .OR. (IHADR .NE. 2)) GO TO 9 ISTOP=2 DESTEP=DESTEP+EN NGKINE=0 GO TO 9999 C 9 CONTINUE IF (INT .NE. 4) GO TO 10 C C --- NEUTRON CAPTURE --- IF (NPRT(9)) PRINT 2000 2000 FORMAT(' *GHEISH* ROUTINE CAPTUR WILL BE CALLED') ISTOP=1 CALL CAPTUR(NOPT) GO TO 40 C 10 CONTINUE IF (INT .NE. 3) GO TO 11 C --- NUCLEAR FISSION --- IF (NPRT(9)) PRINT 2001 2001 FORMAT(' *GHEISH* ROUTINE FISSIO WILL BE CALLED') ISTOP=1 TKIN=FISSIO(EK) GO TO 40 C 11 CONTINUE C C --- ELASTIC AND INELASTIC SCATTERING --- PV( 1,MXGKPV)=P*PX PV( 2,MXGKPV)=P*PY PV( 3,MXGKPV)=P*PZ PV( 4,MXGKPV)=EN PV( 5,MXGKPV)=AMAS PV( 6,MXGKPV)=NCH PV( 7,MXGKPV)=TOF PV( 8,MXGKPV)=KPART PV( 9,MXGKPV)=0. PV(10,MXGKPV)=USERW C C --- ADDITIONAL PARAMETERS TO SIMULATE FERMI MOTION AND EVAPORATION --- DO 111 JENP=1,10 ENP(JENP)=0. 111 CONTINUE ENP(5)=EK ENP(6)=EN ENP(7)=P C IF (INT .NE. 1) GO TO 12 C C *** ELASTIC SCATTERING PROCESSES *** C C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS --- IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35 C C --- NORMAL ELASTIC SCATTERING FOR LIGHT MEDIA --- IF (ATNO2 .LT. 1.5) GO TO 35 C C --- COHERENT ELASTIC SCATTERING FOR HEAVY MEDIA --- IF (NPRT(9)) PRINT 2002 2002 FORMAT(' *GHEISH* ROUTINE COSCAT WILL BE CALLED') CALL COSCAT GO TO 40 C C *** NON-ELASTIC SCATTERING PROCESSES *** 12 CONTINUE C C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS --- IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35 C C *** USE SOMETIMES NUCLEAR REACTION ROUTINE "NUCREC" FOR LOW ENERGY *** C *** PROTON AND NEUTRON SCATTERING *** CALL GRNDM(RNDM,1) TEST1=RNDM(1) TEST2=4.5*(EK-0.01) IF ((KPART .EQ. 14) .AND. (TEST1 .GT. TEST2)) GO TO 85 IF ((KPART .EQ. 16) .AND. (TEST1 .GT. TEST2)) GO TO 86 C C *** FERMI MOTION AND EVAPORATION *** TKIN=CINEMA(EK) PV( 9,MXGKPV)=TKIN ENP(5)=EK+TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS) ENP(7)=SQRT(ABS(ENP(7))) TKIN=FERMI(ENP(5)) ENP(5)=ENP(5)+TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS) ENP(7)=SQRT(ABS(ENP(7))) TKIN=EXNU(ENP(5)) ENP(5)=ENP(5)-TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS) ENP(7)=SQRT(ABS(ENP(7))) C C *** IN CASE OF ENERGY ABOVE CUT-OFF LET THE PARTICLE CASCADE *** TEST=ABS(CHARGE) IF ((TEST .GT. 1.0E-10) .AND. (ENP(5) .GT. CUTHAD)) GO TO 35 IF ((TEST .LE. 1.0E-10) .AND. (ENP(5) .GT. CUTNEU)) GO TO 35 C C --- SECOND CHANCE FOR ANTI-BARYONS DUE TO POSSIBLE ANNIHILATION --- IF ((AMAS .GE. 0.0) .OR. (KPART .LE. 14)) GO TO 13 ANNI=1.3*P IF (ANNI .GT. 0.4) ANNI=0.4 CALL GRNDM(RNDM,1) TEST=RNDM(1) IF (TEST .GT. ANNI) GO TO 35 C C *** PARTICLE WITH ENERGY BELOW CUT-OFF *** C --- ==> ONLY NUCLEAR EVAPORATION AND QUASI-ELASTIC SCATTERING --- 13 CONTINUE C ISTOP=3 C IF (NPRT(9)) PRINT 1002,KPART,EK,EN,P,ENP(5),ENP(6),ENP(7) 1002 FORMAT(' *GHEISH* ENERGY BELOW CUT-OFF FOR GHEISHA PARTICLE ',I3/ $ ' EK,EN,P,ENP(5),ENP(6),ENP(7) = ',6(G12.5,1X)) C IF ((KPART .NE. 14) .AND. (KPART .NE. 16)) GO TO 14 IF (KPART .EQ. 16) GO TO 86 C C --- SLOW PROTON --- 85 CONTINUE IF (NPRT(9)) PRINT 2003,EK,KPART 2003 FORMAT(' *GHEISH* ROUTINE NUCREC WILL BE CALLED', $ ' EK = ',G12.5,' GEV KPART = ',I3) CALL NUCREC(NOPT,2) C IF (NOPT .NE. 0) GO TO 50 C IF (NPRT(9)) PRINT 2004,EK,KPART 2004 FORMAT(' *GHEISH* ROUTINE COSCAT WILL BE CALLED', $ ' EK = ',G12.5,' GEV KPART = ',I3) CALL COSCAT GO TO 40 C C --- SLOW NEUTRON --- 86 CONTINUE IF (NPRT(9)) PRINT 2015 NUCFLG=0 CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW) IF (NUCFLG .NE. 0) GO TO 50 GO TO 40 C C --- OTHER SLOW PARTICLES --- 14 CONTINUE IPA(1)=KPART C --- DECIDE FOR PROTON OR NEUTRON TARGET --- IPA(2)=16 CALL GRNDM(RNDM,1) TEST1=RNDM(1) TEST2=ZNO2/ATNO2 IF (TEST1 .LT. TEST2) IPA(2)=14 AVERN=0.0 NFL=1 IF (IPA(2) .EQ. 16) NFL=2 IPPP=KPART IF (NPRT(9)) PRINT 2005 2005 FORMAT(' *GHEISH* ROUTINE TWOB WILL BE CALLED') CALL TWOB(IPPP,NFL,AVERN) GOTO 40 C C --- INITIALISATION OF CASCADE QUANTITIES --- 35 CONTINUE C C *** CASCADE GENERATION *** C --- CALCULATE FINAL STATE MULTIPLICITY AND LONGITUDINAL AND --- C --- TRANSVERSE MOMENTUM DISTRIBUTIONS --- C C --- FIXED PARTICLE TYPE TO STEER THE CASCADE --- KKPART=KPART C C --- NO CASCADE FOR LEPTONS --- IF (KKPART .LE. 6) GO TO 9999 C C *** WHAT TO DO WITH "NEW PARTICLES" FOR GHEISHA ?????? *** C --- RETURN FOR THE TIME BEING --- IF (KKPART .GE. 35) GO TO 9999 C C --- CASCADE OF HEAVY FRAGMENTS IF ((KKPART .GE. 30) .AND. (KKPART .LE. 32)) GO TO 390 C C --- INITIALIZE THE IPA ARRAY --- CALL VZERO(IPA(1),MXGKCU) C C --- CASCADE OF OMEGA - AND OMEGA - BAR --- IF (KKPART .EQ. 33) GO TO 330 IF (KKPART .EQ. 34) GO TO 331 C NVEPAR=KKPART-17 IF (NVEPAR .LE. 0) GO TO 15 GO TO (318,319,320,321,322,323,324,325,326,327,328,329),NVEPAR C 15 CONTINUE NVEPAR=KKPART-6 GO TO (307,308,309,310,311,312,313,314,315,316,317,318),NVEPAR C C --- PI+ CASCADE --- 307 CONTINUE IF (NPRT(9)) PRINT 2006 2006 FORMAT(' *GHEISH* ROUTINE CASPIP WILL BE CALLED') CALL CASPIP(J,INT,NFL) GO TO 40 C C --- PI0 ==> NO CASCADE --- 308 CONTINUE GO TO 40 C C --- PI- CASCADE --- 309 CONTINUE IF (NPRT(9)) PRINT 2007 2007 FORMAT(' *GHEISH* ROUTINE CASPIM WILL BE CALLED') CALL CASPIM(J,INT,NFL) GO TO 40 C C --- K+ CASCADE --- 310 CONTINUE IF (NPRT(9)) PRINT 2008 2008 FORMAT(' *GHEISH* ROUTINE CASKP WILL BE CALLED') CALL CASKP(J,INT,NFL) GO TO 40 C C --- K0 CASCADE --- 311 CONTINUE IF (NPRT(9)) PRINT 2009 2009 FORMAT(' *GHEISH* ROUTINE CASK0 WILL BE CALLED') CALL CASK0(J,INT,NFL) GO TO 40 C C --- K0 BAR CASCADE --- 312 CONTINUE IF (NPRT(9)) PRINT 2010 2010 FORMAT(' *GHEISH* ROUTINE CASK0B WILL BE CALLED') CALL CASK0B(J,INT,NFL) GO TO 40 C C --- K- CASCADE --- 313 CONTINUE IF (NPRT(9)) PRINT 2011 2011 FORMAT(' *GHEISH* ROUTINE CASKM WILL BE CALLED') CALL CASKM(J,INT,NFL) GO TO 40 C C --- PROTON CASCADE --- 314 CONTINUE IF (NPRT(9)) PRINT 2012 2012 FORMAT(' *GHEISH* ROUTINE CASP WILL BE CALLED') CALL CASP(J,INT,NFL) GO TO 40 C C --- PROTON BAR CASCADE --- 315 CONTINUE IF (NPRT(9)) PRINT 2013 2013 FORMAT(' *GHEISH* ROUTINE CASPB WILL BE CALLED') CALL CASPB(J,INT,NFL) GO TO 40 C C --- NEUTRON CASCADE --- 316 CONTINUE NUCFLG=0 IF (EK .GT. SWTEKN) THEN CALL CASN(J,INT,NFL) IF (NPRT(9)) PRINT 2014 2014 FORMAT(' *GHEISH* ROUTINE CASN WILL BE CALLED') ELSE CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW) IF (NPRT(9)) PRINT 2015 2015 FORMAT(' *GHEISH* ROUTINE GNSLWD WILL BE CALLED') ENDIF IF (NUCFLG .NE. 0) GO TO 50 GO TO 40 C C --- NEUTRON BAR CASCADE --- 317 CONTINUE IF (NPRT(9)) PRINT 2016 2016 FORMAT(' *GHEISH* ROUTINE CASNB WILL BE CALLED') CALL CASNB(J,INT,NFL) GO TO 40 C C --- LAMBDA CASCADE --- 318 CONTINUE IF (NPRT(9)) PRINT 2017 2017 FORMAT(' *GHEISH* ROUTINE CASL0 WILL BE CALLED') CALL CASL0(J,INT,NFL) GO TO 40 C C --- LAMBDA BAR CASCADE --- 319 CONTINUE IF (NPRT(9)) PRINT 2018 2018 FORMAT(' *GHEISH* ROUTINE CASAL0 WILL BE CALLED') CALL CASAL0(J,INT,NFL) GO TO 40 C C --- SIGMA + CASCADE --- 320 CONTINUE IF (NPRT(9)) PRINT 2019 2019 FORMAT(' *GHEISH* ROUTINE CASSP WILL BE CALLED') CALL CASSP(J,INT,NFL) GO TO 40 C C --- SIGMA 0 ==> NO CASCADE --- 321 CONTINUE GO TO 40 C C --- SIGMA - CASCADE --- 322 CONTINUE IF (NPRT(9)) PRINT 2020 2020 FORMAT(' *GHEISH* ROUTINE CASSM WILL BE CALLED') CALL CASSM(J,INT,NFL) GO TO 40 C C --- SIGMA + BAR CASCADE --- 323 CONTINUE IF (NPRT(9)) PRINT 2021 2021 FORMAT(' *GHEISH* ROUTINE CASASP WILL BE CALLED') CALL CASASP(J,INT,NFL) GO TO 40 C C --- SIGMA 0 BAR ==> NO CASCADE --- 324 CONTINUE GO TO 40 C C --- SIGMA - BAR CASCADE --- 325 CONTINUE IF (NPRT(9)) PRINT 2022 2022 FORMAT(' *GHEISH* ROUTINE CASASM WILL BE CALLED') CALL CASASM(J,INT,NFL) GO TO 40 C C --- XI 0 CASCADE --- 326 CONTINUE IF (NPRT(9)) PRINT 2023 2023 FORMAT(' *GHEISH* ROUTINE CASX0 WILL BE CALLED') CALL CASX0(J,INT,NFL) GO TO 40 C C --- XI - CASCADE --- 327 CONTINUE IF (NPRT(9)) PRINT 2024 2024 FORMAT(' *GHEISH* ROUTINE CASXM WILL BE CALLED') CALL CASXM(J,INT,NFL) GO TO 40 C C --- XI 0 BAR CASCADE --- 328 CONTINUE IF (NPRT(9)) PRINT 2025 2025 FORMAT(' *GHEISH* ROUTINE CASAX0 WILL BE CALLED') CALL CASAX0(J,INT,NFL) GO TO 40 C C --- XI - BAR CASCADE --- 329 CONTINUE IF (NPRT(9)) PRINT 2026 2026 FORMAT(' *GHEISH* ROUTINE CASAXM WILL BE CALLED') CALL CASAXM(J,INT,NFL) GO TO 40 C C --- OMEGA - CASCADE --- 330 CONTINUE IF (NPRT(9)) PRINT 2027 2027 FORMAT(' *GHEISH* ROUTINE CASOM WILL BE CALLED') CALL CASOM(J,INT,NFL) GO TO 40 C C --- OMEGA - BAR CASCADE --- 331 CONTINUE IF (NPRT(9)) PRINT 2028 2028 FORMAT(' *GHEISH* ROUTINE CASAOM WILL BE CALLED') CALL CASAOM(J,INT,NFL) GO TO 40 C C --- HEAVY FRAGMENT CASCADE --- 390 CONTINUE IF (NPRT(9)) PRINT 2090 2090 FORMAT(' *GHEISH* ROUTINE CASFRG WILL BE CALLED') NUCFLG=0 CALL CASFRG(NUCFLG,INT,NFL) IF (NUCFLG .NE. 0) GO TO 50 C C *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED *** 40 CONTINUE IF ((NTOT .NE. 0) .OR. (KKPART .NE. KPART)) GO TO 50 C C --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME --- C --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK --- NGKINE=0 TOFG=TOFG+TOF*0.5E-10 C --- In case of crazy momentum value ==> no change to GEANT stack --- IF (P .LT. 0.) GO TO 41 VECT(4)=PX VECT(5)=PY VECT(6)=PZ VECT(7)=P GETOT=EN GEKIN=EK C --- CHECK KINETIC ENERGY --- CALL GEKBIN EDEP=ABS(ENOLD-EN) RMASSI=EN-EK IF (NPRT(9) .AND. (EN .GT. ENOLD)) $ PRINT 8888,EDEP,ENOLD,EN,EK,RMASSI 8888 FORMAT(' *GHEISH* EDEP,ENOLD,EN,EK,M = ',5(G12.5,1X)/ $ ' *GHEISH* =======> EDEP WOULD BE NEGATIVE <========') IF (ISTOP .EQ. 0) DESTEP=DESTEP+EDEP C C --- RE-INITIALIZE THE PROBABILITY FOR HADRONIC INTERACTION --- 41 CONTINUE CALL GRNDM(RNDM,1) IF ((RNDM(1) .LE. 0.) .OR. (RNDM(1) .GE. 1.)) GO TO 41 ZINTHA=-LOG(RNDM(1)) SLHADR=SLENG STEPHA=1.0E10 C NVEDUM=KIPART(IPART) IF (NPRT(9)) PRINT 1003,NTOT,IPART,KPART,KKPART,NVEDUM 1003 FORMAT(' *GHEISH* NO SEC. GEN. NTOT,IPART,KPART,KKPART,KIPART = ', $ 5(I3,1X)/ $ ' CURRENT PARTICLE ON THE STACK AGAIN') GO TO 9999 C C *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND *** C *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED *** 50 CONTINUE C NVEDUM=KIPART(IPART) IF (NPRT(9)) PRINT 1004,NTOT,IPART,KPART,KKPART,NVEDUM 1004 FORMAT(' *GHEISH* SEC. GEN. NTOT,IPART,KPART,KKPART,KIPART = ', $ 5(I3,1X)) C C --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON --- C --- THE GEANT TEMPORARY STACK --- C C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT --- IF ((KPART .NE. 11) .AND. (KPART .NE. 12)) GO TO 52 CALL GRNDM(RNDM,1) KPART=11.5+RNDM(1) C 52 CONTINUE ITY=IKPART(KPART) LNVE=LQ(JPART-ITY) IF (LNVE .LE. 0) PRINT 1234,NTOT,ITY,LNVE 1234 FORMAT('0*GHEISH* 1234 NTOT,ITY,LNVE = ',3(I10,1X)) IF (LNVE .LE. 0) STOP IF (ISTOP .EQ. 0) ISTOP=1 C C --- IN CASE THE NEW PARTICLE IS A NEUTRINO ==> FORGET IT --- IF (KPART .EQ. 2) GO TO 60 C C --- PUT PARTICLE ON THE STACK --- GKIN(1,1)=PX*P GKIN(2,1)=PY*P GKIN(3,1)=PZ*P GKIN(4,1)=SQRT(P*P+RMASS(KPART)**2) GKIN(5,1)=ITY TOFD(1)=TOF*0.5E-10 NGKINE = 1 GPOS(1,1) = VECT(1) GPOS(2,1) = VECT(2) GPOS(3,1) = VECT(3) C IF (NPRT(9)) PRINT 1005,ITY,NGKINE 1005 FORMAT(' *GHEISH* GEANT PART. ',I3,' PUT ONTO STACK AT POS. ',I3) C C *** CHECK WHETHER SECONDARIES HAVE BEEN GENERATED AND COPY THEM *** C *** ALSO ON THE GEANT STACK *** 60 CONTINUE C C --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE --- C --- CONVENTION IS THE FOLLOWING --- C C EVE(INDEX+ 1)= X C EVE(INDEX+ 2)= Y C EVE(INDEX+ 3)= Z C EVE(INDEX+ 4)= NCAL C EVE(INDEX+ 5)= NCELL C EVE(INDEX+ 6)= MASS C EVE(INDEX+ 7)= CHARGE C EVE(INDEX+ 8)= TOF C EVE(INDEX+ 9)= PX C EVE(INDEX+10)= PY C EVE(INDEX+11)= PZ C EVE(INDEX+12)= TYPE C IF (NTOT .LE. 0) GO TO 9999 C C --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED --- DO 61 L=1,NTOT INDEX=(L-1)*12 JND=EVE(INDEX+12) C C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT --- IF ((JND .NE. 11) .AND. (JND .NE. 12)) GO TO 63 CALL GRNDM(RNDM,1) JND=11.5+RNDM(1) C C --- FORGET ABOUT NEUTRINOS --- 63 CONTINUE IF (JND .EQ. 2) GO TO 61 C C --- SWITH TO GEANT QUANTITIES --- ITY=IKPART(JND) JTY=LQ(JPART-ITY) IF (JTY .LE. 0) PRINT 1235,NTOT,ITY,JTY 1235 FORMAT('0*GHEISH* 1235 NTOT,ITY,JTY = ',3(I10,1X)) IF (JTY .LE. 0) STOP * ITRT=Q(JTY+6) PLX=EVE(INDEX+9) PLY=EVE(INDEX+10) PLZ=EVE(INDEX+11) ELT=SQRT(PLX*PLX+PLY*PLY+PLZ*PLZ+Q(JTY+7)**2) C C --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL --- IF (NGKINE .GE. MXGKIN) THEN WRITE(CHMAIL,1236) NTOT, L 1236 FORMAT(' *** GHEISH: ',I9,' particle produced but only ', + I9,' put on the GEANT stack!') CALL GMAIL(1,1) GO TO 9999 ENDIF NGKINE=NGKINE+1 GKIN(1,NGKINE)=PLX GKIN(2,NGKINE)=PLY GKIN(3,NGKINE)=PLZ GKIN(4,NGKINE)=ELT GKIN(5,NGKINE)=ITY TOFD(NGKINE)=EVE(INDEX+8)*0.5E-10 GPOS(1,NGKINE) = VECT(1) GPOS(2,NGKINE) = VECT(2) GPOS(3,NGKINE) = VECT(3) C IF (NPRT(9)) PRINT 1006,ITY,NGKINE,L,(EVE(INDEX+J),J=1,12) 1006 FORMAT(' *GHEISH* GEANT PART. ',I3,' ALSO PUT ONTO STACK AT', $ ' POS. ',I3/ $ ' EVE(',I2,') = '/12(1H ,12X,G12.5/)) C 61 CONTINUE C 9999 CONTINUE C --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW --- NGKINE=MIN(NGKINE,MXGKIN) END