This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / fiface / flufin.F.ori
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.3  1997/06/20 18:32:55  japost
6 * A summary of the problem buf fix and its original report:
7 *
8 * ------------------------------------------------------------------------------
9 * From: Shawn McKee - (313) 764-4395 <mckee@pooh.physics.lsa.umich.edu>
10 * Subject: RE: Problem with FLUKA interface
11 * Date: Thu, 12 Jun 97 18:09:31 METDST
12 *
13 *
14 * I have seen this error with FLUKA.  I traced it to
15 * the following statement:
16 *
17 *      IF (RNDEVT .GE. SINE/FSIG) THEN
18 *
19 * Apparently FSIG (for this event) was not initialized
20 * properly (it is 0.0).  My fix was to insert a check
21 * on FSIG in the FLUFIN routine source.
22 *
23 * ------------------------------------------------------------------------------
24 * From: pniemine@estwm0.wm.estec.esa.nl
25 * Subject: Problem with FLUKA interface - cern.heplib #5667
26 * Date: Thu, 12 Jun 97 17:21:14 METDST
27 *
28 * In trying to run a GEANT program with FLUKA interface in our VMS-AXP (Alpha)
29 * cluster, I get the following error message (originating from FLUFIN) once
30 * a hadronic interaction is to be simulated:
31 *
32 * arithmetic trap, floating/decimal divide by zero
33 * ------------------------------------------------------------------------------
34 *
35 * Revision 1.2  1996/02/27 14:12:57  ravndal
36 * Correction for overstopped antiprotons
37 *
38 * Revision 1.1.1.1  1995/10/24 10:19:53  cernlib
39 * Geant
40 *
41 *
42 #include "geant321/pilot.h"
43 *CMZ :  3.21/02 19/05/94  13.35.12  by  S.Ravndal
44 *-- Author :
45       SUBROUTINE FLUFIN
46 #include "geant321/gcbank.inc"
47 #include "geant321/gccuts.inc"
48 #include "geant321/gcjloc.inc"
49 #include "geant321/gcflag.inc"
50 #include "geant321/gckine.inc"
51 #include "geant321/gcking.inc"
52 #include "geant321/gcmate.inc"
53 #include "geant321/gcphys.inc"
54 #include "geant321/gctrak.inc"
55 #include "geant321/gsecti.inc"
56 #include "geant321/gctmed.inc"
57 #include "geant321/gcunit.inc"
58 #include "geant321/dimpar.inc"
59  
60 #if !defined(CERNLIB_SINGLE)
61 #include "geant321/finuct.inc"
62 #endif
63 #include "geant321/finuc.inc"
64       REAL RNDM(1)
65 #if !defined(CERNLIB_SINGLE)
66       DOUBLE PRECISION AOCMBM, AMSS , ZTAR, RHO , ZLIN, ZLEL, ZLRAD,
67      +ZUL
68 #endif
69       COMMON / FKMAPA / AOCMBM (MXXMDF), AMSS (MXXMDF), ZTAR   (MXXMDF),
70      +                  RHO    (MXXMDF), ZLIN (MXXMDF), ZLEL   (MXXMDF),
71      +                  ZLRAD  (MXXMDF), ZUL  (MXXMDF), MEDIUM (MXXRGN),
72      +                  MULFLG (MXXMDF),IFCOMP(MXXMDF), MSSNUM (MXXMDF),
73      +                  NREGS, NMATF, MTBSNM
74 #if !defined(CERNLIB_SINGLE)
75 #include "geant321/part2t.inc"
76 #endif
77 #include "geant321/part2.inc"
78 #if !defined(CERNLIB_SINGLE)
79 #include "geant321/comcont.inc"
80 #endif
81 #include "geant321/comcon.inc"
82 #if !defined(CERNLIB_SINGLE)
83 #include "geant321/fheavyt.inc"
84 #endif
85 #include "geant321/fheavy.inc"
86 #include "geant321/paprop.inc"
87 #if !defined(CERNLIB_SINGLE)
88 #include "geant321/papropt.inc"
89 #endif
90 #include "geant321/gfkdis.inc"
91 #if !defined(CERNLIB_SINGLE)
92       DOUBLE PRECISION POO,EKE,TXI,TYI,TZI,AMM,WE,ONE,PGEANT,DMOD
93 #endif
94       PARAMETER (ONE=1)
95       DIMENSION IGTOFL(49),IFLTOG(39),IHVTOG(6),ZSAMP(50)
96       DATA IGTOFL / 0, 0, 0, 0, 0, 0,23,13,14,12, 15,16, 8, 1, 2,19, 0,
97      +17,21,22, 20, 34, 36, 38, 9,18, 31, 32, 33, 35, 37, 39, 17*0/
98  
99       DATA IFLTOG /14,15, 3, 2, 4, 4, 1,13,25, 5, 6,10, 8, 9,11,12,18,
100      +26,16,21, 19,20, 7, 7*0, 27, 28, 29, 22, 30, 23, 31, 24, 32/
101       DATA IHVTOG /13,14,45,46,49,47/
102 *
103       NP = 0
104       NPHEAV = 0
105 *
106 *    Stopped particles:
107 *    o Neutral particles are sent to GHSTOP
108 *    o pi+ and K+/K- are forced to decay
109 *    o pi-, antiprotons and antineutrons are sent to FLUKA
110 *      for annihilation (not here but later in this routine)
111       IF ((IGF.EQ.1).OR.
112      +      (GEKIN.EQ.0..AND.ITRTYP.EQ.3.AND.IPART.NE.25)) THEN
113          CALL GHEISH
114          IGF = 0
115          GOTO 999
116       ELSE IF (GEKIN.EQ.0..AND.
117      +        (IPART.EQ.8.OR.IPART.EQ.12.OR.IPART.EQ.11)) THEN
118          CALL GDECAY
119          NMEC=NMEC+1
120          LMEC(NMEC)=5
121          ISTOP=1
122          GOTO 999
123       ENDIF
124 *
125       IF (IFINIT(5) .EQ. 0) CALL FLINIT
126       INT=0
127       IJ=IGTOFL(IPART)
128       IF(IJ.EQ.0) GOTO 110
129       NMEC = NMEC + 1
130       EKE = GEKIN
131       TXI = VECT(4)
132       TYI = VECT(5)
133       TZI = VECT(6)
134       DMOD = ONE/SQRT(TXI**2+TYI**2+TZI**2)
135       TXI = TXI*DMOD
136       TYI = TYI*DMOD
137       TZI = TZI*DMOD
138       WE  = 1.
139       JMA = LQ(JMATE-NMAT)
140       NCOMP = ABS (Q(JMA+11))
141       AMM = Q(JMA+6)
142       JMIXT = LQ(JMA-5)
143  
144 *    Antiprotons, antineutrons and pi- are sent to
145 *    eventv for annihilation
146       IF (GEKIN.EQ.0..AND.
147      +         (IPART.EQ.15.OR.IPART.EQ.9.OR.IPART.EQ.25)) THEN
148          IF(NCOMP.LE.1) THEN
149             AMSS(1) = Q(JMA+6)
150             ZTAR(1) = Q(JMA+7)
151             MSSNUM(1) = 0
152             RHO(1) = Q(JMA+8)
153          ELSE
154             ZSAMP(1) = 0.
155             DO 10 I=1,NCOMP
156                ZSAMP(I+1) = ZSAMP(I) + Q(JMIXT+NCOMP+I)
157    10       CONTINUE
158             CALL GRNDM(RNDM,1)
159             ZCONT=ZSAMP(NCOMP+1)*RNDM(1)
160             DO 20 I=1,NCOMP
161                IF(ZCONT.LE.ZSAMP(I+1)) GO TO 30
162    20       CONTINUE
163             I = NCOMP
164    30       CONTINUE
165             AMSS(1)   = Q(JMIXT+I)
166             MSSNUM(1) = 0
167             ZTAR(1)   = Q(JMIXT+NCOMP+I)
168             RHO(1)    = Q(JMIXT+2*NCOMP+I)*DENS
169          END IF
170          EKE = 1E-9
171          POO=SQRT(EKE*(EKE+2*AM(IJ)))
172          CALL EVENTV(IJ,POO,EKE,TXI,TYI,TZI,WE,1)
173          NMEC=NMEC+1
174          LMEC(NMEC)=17
175          GOTO 80
176       ELSE IF (GEKIN.LE.CUTHAD .AND. ITRTYP.EQ.4) THEN
177          DESTEP = DESTEP + GEKIN
178          GEKIN  = 0.
179          GETOT  = AMASS
180          VECT(7) = 0.
181          ISTOP = 1
182          LMEC(NMEC)=30
183          GO TO 110
184       ENDIF
185 *
186       CALL GRNDM(RNDM,1)
187       RNDEVT=RNDM(1)
188
189       IF (FSIG.LE.0.) THEN
190          SRAT = 1.E7
191       ELSE
192          SRAT = SINE/FSIG
193       ENDIF
194       IF ( RNDEVT .GE. SRAT) THEN
195          IF (GEKIN .GT. 0.02) THEN
196             POO=SQRT(EKE*(EKE+2*AM(IJ)))
197          ELSE
198             GO TO 110
199          END IF
200          INT=1
201          LMEC(NMEC)=13
202          IF(NCOMP.LE.1) THEN
203             CALL NUCREL(IJ,POO,EKE,TXI,TYI,TZI,AMM,WE)
204          ELSE
205             CALL GRNDM(RNDM,1)
206             RCONT=ELXNOR*RNDM(1)
207             DO 40  I=1,NCOMP
208                IF(RCONT.LE.CABELX(I)) GO TO 50
209    40       CONTINUE
210             I=NCOMP
211    50       CONTINUE
212             CALL NUCREL(IJ,POO,EKE,TXI,TYI,TZI,ONE*Q(JMIXT+I),WE)
213          END IF
214       ELSE
215          LMEC(NMEC)=20
216          IF (IHADR.EQ.2) THEN
217             ISTOP = 2
218             DESTEP = DESTEP + GETOT
219             GO TO 110
220          ENDIF
221          IF (GEKIN .GT. 0.02) THEN
222             POO=SQRT(EKE*(EKE+2*AM(IJ)))
223          ELSE
224             IF ((IJ.EQ.2 .OR. IJ.EQ.9 .OR. IJ.EQ.14 .OR. IJ.EQ.16)
225      +            .AND. GEKIN .GT. 0.0) THEN
226                POO=SQRT(EKE*(EKE+2*AM(IJ)))
227             ELSE
228                NMEC=NMEC-1
229                GO TO 110
230             END IF
231          END IF
232          INT=2
233          IF(NCOMP.LE.1) THEN
234             AMSS(1) = Q(JMA+6)
235             ZTAR(1) = Q(JMA+7)
236             MSSNUM(1) = 0
237             RHO(1) = Q(JMA+8)
238          ELSE
239             CALL GRNDM(RNDM,1)
240             RCONT=ANXNOR*RNDM(1)
241             DO 60  I=1,NCOMP
242                IF(RCONT.LE.CABINX(I)) GO TO 70
243    60       CONTINUE
244             I=NCOMP
245    70       CONTINUE
246             AMSS(1)   = Q(JMIXT+I)
247             MSSNUM(1) = 0
248             ZTAR(1)   = Q(JMIXT+NCOMP+I)
249             RHO(1)    = Q(JMIXT+2*NCOMP+I)*DENS
250          END IF
251          CALL EVENTV(IJ,POO,EKE,TXI,TYI,TZI,WE,1)
252       END IF
253 *
254    80 IF(NP.EQ.1.AND.NPHEAV.EQ.0.AND.KPART(1).EQ.IJ) THEN
255          VECT(4)=CXR(1)
256          VECT(5)=CYR(1)
257          VECT(6)=CZR(1)
258          VECT(7)=SQRT(TKI(1)*(TKI(1)+2*AMASS))
259          GETOT=TKI(1)+AMASS
260          GEKIN=TKI(1)
261       ELSE
262          ISTOP=1
263          NSTAK1 = MIN(NP,MXGKIN-NGKINE)
264          IF(NP.GT.NSTAK1) THEN
265             WRITE(CHMAIL,10000) NP-NSTAK1
266             CALL GMAIL(0,0)
267          ENDIF
268          DO 90  K=1,NSTAK1
269             NGKINE = NGKINE + 1
270             IF (KPART(K) .EQ. 24 .OR. KPART(K) .EQ. 25) THEN
271                KPART(K) = 19
272                CALL GRNDM(RNDM,1)
273                IF (RNDM(1) .GT. 0.5) KPART(K) = 12
274             END IF
275             IGEPAR = IFLTOG(KPART(K))
276             JPA = LQ(JPART-IGEPAR)
277             AGEMAS = Q(JPA+7)
278             PGEANT = SQRT(TKI(K)*(TKI(K)+2*AGEMAS))
279             GKIN(1,NGKINE)=CXR(K)*PGEANT
280             GKIN(2,NGKINE)=CYR(K)*PGEANT
281             GKIN(3,NGKINE)=CZR(K)*PGEANT
282             GKIN(4,NGKINE)=TKI(K)+AGEMAS
283             GKIN(5,NGKINE)=IGEPAR
284             TOFD(NGKINE)=0.0
285             GPOS(1,NGKINE) = VECT(1)
286             GPOS(2,NGKINE) = VECT(2)
287             GPOS(3,NGKINE) = VECT(3)
288    90    CONTINUE
289 *
290          NSTAK2 = MIN(NPHEAV,MXGKIN-NGKINE)
291          IF(NPHEAV.GT.NSTAK2) THEN
292             WRITE(CHMAIL,10100) NPHEAV-NSTAK2
293             CALL GMAIL(0,0)
294          ENDIF
295          DO 100 K=1,NSTAK2
296             NGKINE = NGKINE + 1
297             IGEPAR = IHVTOG(KHEAVY(K))
298             JPA = LQ(JPART-IGEPAR)
299             AGEMAS = Q(JPA+7)
300             PGEANT = SQRT(TKHEAV(K)*(TKHEAV(K)+2*AGEMAS))
301             GKIN(1,NGKINE)=CXHEAV(K)*PGEANT
302             GKIN(2,NGKINE)=CYHEAV(K)*PGEANT
303             GKIN(3,NGKINE)=CZHEAV(K)*PGEANT
304             GKIN(4,NGKINE)=TKHEAV(K)+AGEMAS
305             GKIN(5,NGKINE)=IGEPAR
306             TOFD(NGKINE)=0.0
307             GPOS(1,NGKINE) = VECT(1)
308             GPOS(2,NGKINE) = VECT(2)
309             GPOS(3,NGKINE) = VECT(3)
310   100    CONTINUE
311 *
312          KCASE=NAMEC(12)
313       END IF
314   110 CONTINUE
315       ZINTHA = GARNDM(DUMMY)
316       SLHADR = SLENG
317       STEPHA = 1.0E10
318 10000 FORMAT(' **** FLUFIN: Stack overflow, ',I6,' particles lost')
319 10100 FORMAT(' **** FLUFIN: Stack overflow, ',I6,
320      +' heavy particles lost')
321   999 END