5 * Revision 1.1.1.1 1995/10/24 10:20:14 cernlib
9 #include "geant321/pilot.h"
10 *CMZ : 3.21/04 11/01/95 08.57.12 by S.Ravndal
12 SUBROUTINE GFTMAT(IMATE,IPATT,CHMECA,KDIM,TKIN,VALUE,PCUT,IXST)
14 C. ******************************************************************
16 C. * FETCH and INTERPOLATE the DE/DX and Cross sections *
17 C. * tabulated in JMATE banks coresponding to : *
18 C. * material IMATE, particle IPATT, mecanism name CHMECA, *
19 C. * kinetic energies TKIN *
21 C. * The CHMECAnism name can be : *
22 C. * 'HADF' 'INEF' 'ELAF' 'FISF' 'CAPF' *
23 C. * 'HADG' 'INEG' 'ELAG' 'FISG' 'CAPG' *
24 C. * 'LOSS' 'PHOT' 'ANNI' 'COMP' 'BREM' *
25 C. * 'PAIR' 'DRAY' 'PFIS' 'RAYL' 'HADG' *
26 C. * 'MUNU' 'RANG' 'STEP' *
28 C. * For Hadronic particles it also computes the *
29 C. * hadronic cross section from FLUKA ( '***F' ) or *
30 C. * GHEISHA ( '***G' ) programs: *
31 C. * HADF or HADG -- total *
32 C. * INEF or INEG -- inelastic *
33 C. * ELAF or ELAG -- elastic *
34 C. * FISF or FISG -- fission (0.0 for FLUKA) *
35 C. * CAPF or CAPG -- neutron capture (0.0 for FLUKA) *
37 C. * Input parameters *
38 C. * IMATE Geant material number *
39 C. * IPATT Geant particle number *
40 C. * CHMECA mechanism name of the bank to be fetched *
41 C. * KDIM dimension of the arrays TKIN , VALUE *
42 C. * TKIN array of kinetic energy of incident particle (in Gev) *
44 C. * Output parameters *
45 C. * VALUE array of energy loss (in Mev/cm) , *
46 C. * or stopping range (cm) ,or continuous step (cm) *
47 C. * or macroscopic cross section (in 1/cm) *
48 C. * PCUT(5) array of the physical cuts in material IMATE (Gev) *
49 C. * IXST flag = 1 if the array VALUE is filled , =0 otherwise *
51 C. * ==>Called by : <USER> GPLMAT GRPMAT *
52 C. * Authors R.Brun, M.Maire ********* *
54 C. ******************************************************************
56 #include "geant321/gcbank.inc"
57 #include "geant321/gcnum.inc"
58 #include "geant321/gconst.inc"
59 #include "geant321/gcunit.inc"
60 #include "geant321/gcmulo.inc"
61 #include "geant321/gsecti.inc"
62 #include "geant321/gcflag.inc"
63 #include "geant321/gctrak.inc"
64 #include "geant321/gcmate.inc"
65 #include "geant321/gctmed.inc"
66 #include "geant321/gfkdis.inc"
67 #include "geant321/gcking.inc"
68 #include "geant321/gckine.inc"
72 DIMENSION TKIN(KDIM), VALUE(KDIM), PCUT(5)
73 #include "geant321/gcnmec.inc"
75 * ------------------------------------------------------------------
78 IF(KDIM.LE.0) GO TO 999
81 IF(CHMECA.EQ.CHNMEC(KMECA)) THEN
86 WRITE(CHMAIL,'('' *** GFTMAT: Mechanism '',A, '
87 + //' ''not implemented'')') CHMECA
98 IF(JMATE.LE.0) GO TO 999
99 IF(IMATE.LE.0) GO TO 999
100 IF(IMATE.GT.NMATE) GO TO 110
101 JMA = LQ(JMATE-IMATE)
102 IF(JMA.LE.0) GO TO 110
105 IF(Z.LT.1.) GO TO 999
118 IF(JTMED.LE.0) GO TO 999
119 IF(NTMED.LE.0) GO TO 999
123 IF(JTM.LE.0) GO TO 40
126 IF(IMAT.EQ.IMATE) THEN
128 IF(JTMN.NE.0) JBANK = JTMN
132 50 CALL UCOPY( Q(JBANK+6),PCUT(1),5)
143 IF(JPART.LE.0) GO TO 999
144 IF(IPATT.LE.0) GO TO 999
145 IF(IPATT.GT.NPART) GO TO 110
146 JPA = LQ(JPART-IPATT)
147 IF(JPA.LE.0) GO TO 110
152 * *** Find the correct pointer
161 IF (CHMECA.EQ.'PHOT') JBANK = LQ(JMA- 6)
162 IF (CHMECA.EQ.'COMP') JBANK = LQ(JMA- 8)
163 IF (CHMECA.EQ.'PAIR') JBANK = LQ(JMA-10)
164 IF (CHMECA.EQ.'PFIS') JBANK = LQ(JMA-12)
165 IF (CHMECA.EQ.'RAYL') JBANK = LQ(JMA-13)
167 * *** Electrons / positons
169 ELSE IF (ITYPE.EQ.2) THEN
170 IF (CHMECA.EQ.'LOSS') THEN
172 IF (CHARGE.GT.0.) ISHIF = NEK1
173 ELSE IF (CHMECA.EQ.'RANG') THEN
175 IF (CHARGE.GT.0.) ISHIF = NEK1
176 ELSE IF (CHMECA.EQ.'STEP') THEN
178 ELSE IF ((CHMECA.EQ.'ANNI').AND.(CHARGE.GT.0.)) THEN
180 ELSE IF (CHMECA.EQ.'BREM') THEN
182 ELSE IF (CHMECA.EQ.'DRAY') THEN
184 IF (CHARGE.GT.0.) ISHIF = NEK1
187 * *** Neutral hadrons (***F for FLUKA cross sections
188 * ***G for GHEISHA cross sections
189 * LOWN and N*** for MICAP cross sections)
191 ELSE IF (ITYPE.EQ.3) THEN
192 IF((CHMECA.EQ.'HADF').OR.(CHMECA.EQ.'INEF')
193 + .OR.(CHMECA.EQ.'ELAF') .OR.(CHMECA.EQ.'FISF')
194 + .OR.(CHMECA.EQ.'CAPF')) THEN
196 IF (IFINIT(5) .EQ. 0) CALL FLINIT
197 ELSE IF((CHMECA.EQ.'HADG').OR.(CHMECA.EQ.'INEG')
198 + .OR.(CHMECA.EQ. 'ELAG').OR.(CHMECA.EQ.'FISG')
199 + .OR.(CHMECA.EQ.'CAPG')) THEN
202 ELSE IF(IMECA.GE.IBLOWN.AND.IPATT.EQ.13) THEN
203 IF (IFINIT(7) .EQ. 0) CALL GMORIN
208 * *** Charged hadrons (***F for FLUKA cross sections
209 * ***G for GHEISHA cross sections)
212 ELSE IF (ITYPE.EQ.4.OR.ITYPE.EQ.8) THEN
214 IF (CHMECA.EQ.'LOSS') THEN
216 ELSE IF (CHMECA.EQ.'RANG') THEN
217 JBANK = LQ(JMA- 16) + NEK1
218 ELSE IF (CHMECA.EQ.'STEP') THEN
220 JRANG = LQ(JMA -16) + NEK1
221 CUTPRO = CUTHAD*RMASS
222 CUTPRO = MAX(ELOW(1), MIN( CUTPRO, ELOW(NEK1)*0.99))
223 IKCUT = GEKA*LOG10(CUTPRO) + GEKB
224 GKC = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT))
225 STOPC = (1.-GKC)*Q(JRANG+IKCUT) + GKC*Q(JRANG+IKCUT+1)
226 ELSE IF (CHMECA.EQ.'DRAY') THEN
231 ELSE IF(((CHMECA.EQ.'HADF').OR.(CHMECA.EQ.'INEF')
232 + .OR.(CHMECA.EQ. 'ELAF').OR.(CHMECA.EQ.'FISF')
233 + .OR.(CHMECA.EQ.'CAPF')).AND.ITYPE.NE.8) THEN
235 IF (IFINIT(5) .EQ. 0) CALL FLINIT
236 ELSE IF(((CHMECA.EQ.'HADG').OR.(CHMECA.EQ.'INEG')
237 + .OR.(CHMECA.EQ. 'ELAG').OR.(CHMECA.EQ.'FISG')
238 + .OR.(CHMECA.EQ.'CAPG')).AND.ITYPE.NE.8) THEN
246 ELSE IF (ITYPE.EQ.5) THEN
247 IF (CHMECA.EQ.'LOSS') THEN
249 ELSE IF (CHMECA.EQ.'RANG') THEN
251 ELSE IF (CHMECA.EQ.'STEP') THEN
253 ELSE IF (CHMECA.EQ.'MUNU') THEN
255 ELSE IF (CHMECA.EQ.'BREM') THEN
258 ELSE IF (CHMECA.EQ.'PAIR') THEN
261 ELSE IF (CHMECA.EQ.'DRAY') THEN
268 ELSEIF (ITYPE.EQ.6) THEN
275 ELSEIF (ITYPE.EQ.7) THEN
276 IF (CHMECA.EQ.'LABS') THEN
278 * *** Not implemented yet!
284 IF(CHARGE.NE.0.AND.ITCKOV.NE.0) THEN
285 IF(IQ(JTM-2).GE.3) THEN
286 IF(LQ(JTM-3).NE.0.AND.LQ(LQ(JTM-3)-3).NE.0) THEN
288 * *** In this tracking medium Cerenkov photons are generated and
289 * *** tracked. Set to 1 the corresponding flag and calculate the
290 * *** relevant pointers.
294 JABSCO = LQ(JTCKOV-1)
295 JEFFIC = LQ(JTCKOV-2)
296 JINDEX = LQ(JTCKOV-3)
297 JCURIN = LQ(JTCKOV-4)
303 IF(JBANK.EQ.0) GO TO 999
307 JBANK = JBANK + ISHIF
309 * Find bin number in table JMATE
310 EKP = TKIN(IKB)*RMASS
311 EKP = MAX(ELOW(1), MIN( EKP, ELOW(NEK1)*0.99))
312 IKP=GEKA*LOG10(EKP)+GEKB +0.001
313 GKRA=(EKP-ELOW(IKP))/(ELOW(IKP+1)-ELOW(IKP))
316 * Retieve value from bank JMATE
317 VALUE(IKB) = (1.-GKRA)*Q(JBANK+IKP) + GKRA*Q(JBANK+IKP+1)
318 IF ((CHMECA.EQ.'PHOT').AND.(EKP.GE.0.05 )) VALUE(IKB) =
320 IF ((CHMECA.EQ.'MUNU').AND.(EKP.LT.0.05 )) VALUE(IKB) =
322 IF ( CHMECA.EQ.'LOSS') THEN
323 VALUE(IKB) = VALUE(IKB)*CHARGE**2*1.E+3
324 ELSE IF (CHMECA.EQ.'RANG') THEN
325 VALUE(IKB) = VALUE(IKB)/(RMASS*CHARGE*CHARGE)
326 ELSE IF (CHMECA.NE.'STEP') THEN
327 IF (VALUE(IKB).GT.0.) THEN
328 VALUE(IKB) = 1./VALUE(IKB)
334 ELSEIF (JBANK.EQ.-1) THEN
335 * Compute step due to muls + loss + field
338 PMOM =SQRT(GEKIN*(GETOT+AMASS))
343 IF (IFIELD*FIELDM.NE.0.)
344 + SFIELD = 3333.*DEGRAD*TMAXFD*PMOM/ABS(FIELDM*CHARGE)
346 + SMULS = MIN (2232.*RADL*((PMOM**2)/(GETOT*CHARGE))**2 ,
348 IF (ILOSS*DEEMAX.GT.0.) THEN
349 STOPP = (1.-GKRA)*Q(JRANG+IKP) + GKRA*Q(JRANG+IKP+1)
350 STOPMX = (STOPP - STOPC)/(RMASS*CHARGE*CHARGE)
351 IF (STOPMX.LT.0.) STOPMX = 0.
352 EKF = MAX ( ELOW(1) , (1.-DEEMAX)*EKP )
353 IKF = GEKA*LOG10(EKF) + GEKB
354 GKF = (EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF))
355 SLOSP= STOPP - (1.-GKF)*Q(JRANG+IKF) - GKF*Q(JRANG+IKF+1)
356 SLOSS = SLOSP/(RMASS*CHARGE*CHARGE)
359 VECT(7)=SQRT(TKIN(IKB)*(TKIN(IKB)+2*AMASS))
361 STCKOV = MXPHOT/MAX(3.*DNDL,1E-10)
366 IF (STOPMX.LE.STMIN) THEN
369 VALUE(IKB) = MAX(STMIN,MIN(STCKOV,SLOSS,SFIELD,SMULS,
373 ELSEIF (JBANK.EQ.-2) THEN
374 * Compute delta ray cross section for hadrons
378 BET2=GEKIN*GAMASS/(GETOT*GETOT)
379 TMAX=EMASS*GEKIN*GAMASS/(0.5*AMASS*AMASS+EMASS*GETOT)
380 IF(TMAX.GT.DCUTM)THEN
382 SIG=(1.-Y+BET2*Y*LOG(Y))/DCUTM
383 IF(AMASS.GT.0.9)SIG=SIG+0.5*(TMAX-DCUTM)/(GETOT*GETOT)
384 VALUE(IKB)=SIG*AZRO*CHARGE*CHARGE*EMASS/BET2
389 * compute hadronic cross section from FLUKA code
391 ELSEIF (JBANK.EQ.-3) THEN
393 PMOM = SQRT(GEKIN*(GEKIN+2*AMASS))
396 IF (IPATT.NE.IPART) THEN
404 IF (CHMECA.EQ.'HADF') VALUE(IKB)= FSIG
405 IF (CHMECA.EQ.'INEF') VALUE(IKB)= SINE
406 IF (CHMECA.EQ.'ELAF') VALUE(IKB)= SELA
407 IF (CHMECA.EQ.'FISF') VALUE(IKB)= 0.0
408 IF (CHMECA.EQ.'CAPF') VALUE(IKB)= 0.0
410 * compute hadronic cross section from GHEISHA code
412 ELSEIF (JBANK.EQ.-4) THEN
414 PMOM = SQRT(GEKIN*(GEKIN+2*AMASS))
419 IF (JTMN.GT.0) HHHH=Q(JTMN+26)
420 VALUE(IKB)=GHESIG(PMOM,GEKIN,A,Q(JMIXT+1), Q(JMIXT+NLM+
421 + 1),Q(JMIXT+2*NLM+1),NLM,DENS,HHHH,IPATT)
422 IF (CHMECA .EQ. 'INEG') THEN
425 VALUE(IKB) = AIIN(K) + VALUE(IKB)
427 ELSE IF (CHMECA .EQ. 'ELAG') THEN
430 VALUE(IKB) = AIEL(K) + VALUE(IKB)
432 ELSE IF (CHMECA .EQ. 'FISG') THEN
435 VALUE(IKB) = AIFI(K) + VALUE(IKB)
437 ELSE IF (CHMECA .EQ. 'CAPG') THEN
440 VALUE(IKB) = AICA(K) + VALUE(IKB)
445 VALUE(IKB)=GHESIG(PMOM,GEKIN,A,A,Z,1.,1,DENS,0.,IPATT)
446 IF (CHMECA .EQ. 'INEG') VALUE(IKB) = AIIN(1)
447 IF (CHMECA .EQ. 'ELAG') VALUE(IKB) = AIEL(1)
448 IF (CHMECA .EQ. 'FISG') VALUE(IKB) = AIFI(1)
449 IF (CHMECA .EQ. 'CAPG') VALUE(IKB) = AICA(1)
453 * compute the cross-section for low-energy neutrons
456 ELSEIF (JBANK.EQ.-5) THEN
458 IF (GEKIN.LE..02) THEN
459 IF (CHMECA .EQ. 'LOWN') THEN
460 VALUE(IKB) = SIGMOR(GEKIN*1.E+9,IMATE)
463 CALL GMXSEC (IMECA,VALUE(IKB))
473 110 WRITE(CHMAIL,10100) IMATE ,IPATT
476 10000 FORMAT(' ***** GFTMAT: No processes active for geantinos')
477 10100 FORMAT(' ***** GFTMAT error : material',I4,
478 + ' or particle',I4,' not defined' )