This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / giface / ghesig.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:15  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.37  by  S.Giani
11 *-- Author :
12       FUNCTION GHESIG(PPART,EKIN,AVER,A,Z,W,KK,DENS,QCOR,LPART)
13 C
14 C *** CALCULATION OF THE PROBABILITIES FOR (IN)ELASTIC INTERACTIONS ***
15 C *** OF STABLE PARTICLES ON PROTON AND NEUTRON                     ***
16 C *** NVE 07-APR-1988 ***
17 C
18 C CALLED BY : GPGHEI
19 C ORIGIN : F.CARMINATI, H.FESEFELDT (ROUTINE INTACT 06-OCT-1987)
20 C
21 C *** IPART DENOTES THE GHEISHA PARTICLE INDEX ***
22 C
23 C CONVENTION :
24 C
25 C   PARTICLE                 IPART
26 C   ------------------------------
27 C   GAMMA                    1
28 C   NEUTRINO                 2
29 C   POSITRON                 3
30 C   ELECTRON                 4
31 C   MUON +                   5
32 C   MUON -                   6
33 C   PION +                   7
34 C   PION 0                   8
35 C   PION -                   9
36 C   KAON +                  10
37 C   KAON 0 S                11
38 C   KAON 0 L                12
39 C   KAON -                  13
40 C   PROTON                  14
41 C   PROTON BAR              15
42 C   NEUTRON                 16
43 C   NEUTRON BAR             17
44 C   LAMBDA                  18
45 C   LAMBDA BAR              19
46 C   SIGMA +                 20
47 C   SIGMA 0                 21
48 C   SIGMA -                 22
49 C   SIGMA + BAR             23
50 C   SIGMA 0 BAR             24
51 C   SIGMA - BAR             25
52 C   XSI 0                   26
53 C   XSI -                   27
54 C   XSI 0 BAR               28
55 C   XSI - BAR               29
56 C   DEUTERON                30
57 C   TRITON                  31
58 C   ALPHA                   32
59 C   OMEGA -                 33
60 C   OMEGA - BAR             34
61 C   NEW PARTICLES           35
62 C
63 #include "geant321/gcunit.inc"
64 #include "geant321/gconsp.inc"
65 C
66 C --- FOR CROSS-SECTION INFORMATION SEE "PCSDATA" ---
67 C
68 #include "geant321/gsecti.inc"
69 C --- GHEISHA COMMONS ---
70 #include "geant321/s_result.inc"
71 #include "geant321/s_prntfl.inc"
72 C
73       DIMENSION A(KK),Z(KK),W(KK)
74 C
75       DIMENSION ALPHA(35),ALPHAC(41),IPART2(7),CSA(4)
76       DIMENSION PARTEL(35),PARTIN(35),ICORR(35),INTRC(35)
77 C
78       PARAMETER (ONETHR=1./3.)
79       SAVE ALPHA,ALPHAC,PARTEL,PARTIN,CSA,IPART2,ICORR,INTRC
80 C
81 #include "geant321/s_csdim.inc"
82 #include "geant321/pcodim.inc"
83 #include "geant321/s_csdat.inc"
84 #include "geant321/pcodat.inc"
85 C
86       DATA ALPHA    / 6*0.7,
87      +                0.75 ,0.75 ,0.75 ,
88      +                0.76,0.76 ,0.76 ,0.76 ,
89      +                0.685,0.63 ,0.685,0.63,0.685,0.63,
90      +                3*0.685,3*0.63,2*0.685,2*0.63,
91      +                3*0.7,0.685,0.63,0.7/
92       DATA ALPHAC    /1.2,1.2,1.2,1.15,0.90,0.91,0.98,1.06,1.10,1.11,
93      +                1.10,1.08,1.05,1.01,0.985,0.962,0.945,0.932,
94      +                0.925,0.920,0.920,0.921,0.922,0.923,0.928,0.931,
95      +                0.940,0.945,0.950,0.955,0.958,0.962,0.965,0.976,
96      +                0.982,0.988,0.992,1.010,1.020,1.030,1.040/
97       DATA PARTEL/6*0.,29*1./
98       DATA PARTIN/6*0.,1.00,0.00,1.05,1.20,1.35,1.30,1.20,1.00,1.30,
99      +            1.00,1.30,1.00,1.30,1.00,1.00,1.00,1.30,1.30,1.30,
100      +            1.00,1.00,1.30,1.30,1.00,1.,1.,1.,1.3,1./
101       DATA ICORR /14*1, 0, 1, 0, 1, 0, 3*1, 3*0, 2*1, 2*0, 4*1, 2*0/
102       DATA INTRC /6*0, 1, 0, 12*1, 0, 2*2, 0, 10*1, 0/
103 C
104 C CROSS SECTIONS ON NUCLEUS ARE KNOWN ONLY FOR PIONS AND PROTONS.
105 C THE GENERAL LAW SIGMA(A)=1.25*SIGMA(TOT,PROTON)*A**ALPHA IS VALID
106 C ONLY FOR MOMENTA > 2 GEV.THE PARAMETRIZATION DONE HERE GIVES ONLY
107 C A BEHAVIOUR AVERAGED OVER MOMENTA AND PARTICLE TYPES.
108 C FOR A DETECTOR WITH ONLY A FEW MATERIALS IT'S OF COURSE MUCHBETTER
109 C TO USE TABLES OF THE MEASURED CROSS SECTIONS .
110 C FOR ELEMENTS WITH THE FOLLOWING ATOMIC NUMBERS MEASURED CROSS
111 C SECTIONS ARE AVAILABLE (SEE "PCSDATA").
112 C
113 C                 H   AL     CU     PB
114       DATA  CSA  /1. ,27.00 ,63.54 ,207.19 /
115       DATA IPART2/9,8,7,11,10,13,12/
116 C
117 C
118 C --- Initialise GHESIG and switch to GHEISHA particle code ---
119       GHESIG=0.0
120       IF(LPART.GT.48)GO TO 160
121       IPART=KIPART(LPART)
122 C
123 C --- No interaction for gammas, neutrinos, electrons, positrons, muons,
124 C --- neutral pions, neutral sigmas and antisigmas and new particles.
125       IF(INTRC(IPART).EQ. 0) GO TO 160
126       P=PPART
127       EK=EKIN
128 C
129 C --- Initialise the cross-sections with 0.0 ---
130       DO 10  K=1,KK
131          AIEL(K)=0.0
132          AIIN(K)=0.0
133          AIFI(K)=0.0
134          AICA(K)=0.0
135    10 CONTINUE
136 C
137       IF ((IPART .GE. 30) .AND. (IPART .LE. 32)) THEN
138 C
139 C --- Take geometrical cross sections for inelastic scattering ---
140 C --- of deuterons, tritons and alphas ---
141          IF (IPART .EQ. 30) THEN
142             APART=2.0**ONETHR
143          ELSEIF (IPART .EQ. 31) THEN
144             APART=3.0**ONETHR
145          ELSEIF (IPART .EQ. 32) THEN
146             APART=4.0**ONETHR
147          END IF
148          DO 20 K=1,KK
149             AIIN(K)=49.0*(APART+A(K)**ONETHR)**2
150    20    CONTINUE
151          IF (NPRT(9)) THEN
152             WRITE(CHMAIL,10000)
153             CALL GMAIL(0,0)
154          END IF
155 C
156       ELSE IF ((IPART .EQ. 16) .AND. (EK .LE. 0.0327)) THEN
157 C
158 C --- Use tables for low energy neutrons ---
159 C --- get energy bin ---
160          JE2=17
161          DO 30 J=2,17
162             IF (EK .LT. ELAB(J)) THEN
163                JE2=J
164                GO TO 40
165             END IF
166    30    CONTINUE
167 C
168    40    JE1=JE2-1
169          EKX=MAX(EK,1.0E-9)
170          DELAB=ELAB(JE2)-ELAB(JE1)
171          DO 70 K=1,KK
172 C
173 C --- Get A bin ---
174             JA2=15
175             DO 50 J=2,15
176                IF (A(K) .LT. CNLWAT(J)) THEN
177                   JA2=J
178                   GO TO 60
179                END IF
180    50       CONTINUE
181 C
182    60       JA1=JA2-1
183             DNLWAT=CNLWAT(JA2)-CNLWAT(JA1)
184 C
185 C --- Use linear interpolation or extrapolation by Y=RCE*X+RCA*X+B ---
186 C
187 C --- Elastic cross section ---
188 C --- E interpolation or extrapolation at JA1 ---
189             DY=CNLWEL(JA1,JE2)-CNLWEL(JA1,JE1)
190             RCE=DY/DELAB
191 C --- A interpolation or extrapolation at JE1 ---
192             DY=CNLWEL(JA2,JE1)-CNLWEL(JA1,JE1)
193             RCA=DY/DNLWAT
194             B=CNLWEL(JA1,JE1)-RCE*ELAB(JE1)-RCA*CNLWAT(JA1)
195             AIEL(K)=RCE*EK+RCA*A(K)+B
196 C
197 C --- Inelastic cross section ---
198 C --- E interpolation or extrapolation at JA1 ---
199             DY=CNLWIN(JA1,JE2)-CNLWIN(JA1,JE1)
200             RCE=DY/DELAB
201 C --- A interpolation or extrapolation at JE1 ---
202             DY=CNLWIN(JA2,JE1)-CNLWIN(JA1,JE1)
203             RCA=DY/DNLWAT
204             B=CNLWIN(JA1,JE1)-RCE*ELAB(JE1)-RCA*CNLWAT(JA1)
205             AIIN(K)=RCE*EK+RCA*A(K)+B
206 C
207             IZNO=Z(K)+0.01
208             AICA(K)=11.12*CSCAP(IZNO)/(EKX*1.0E6)**0.577
209    70    CONTINUE
210          IF (NPRT(9)) THEN
211             WRITE(CHMAIL,10100)
212             CALL GMAIL(0,0)
213          END IF
214       ELSE
215 C
216 C --- Use parametrization of cross section data for all other cases ---
217 C
218          IF (NPRT(9)) THEN
219             WRITE(CHMAIL,10200)
220             CALL GMAIL(0,0)
221          END IF
222 C
223 C --- Get momentum bin ---
224          J=40
225          DO 80 I=2,41
226             IF (P .LT. PLAB(I)) THEN
227                J=I-1
228                GO TO 90
229             END IF
230    80    CONTINUE
231 C
232 C --- Start with  cross sections for scattering on free protons ---
233 C --- use linear interpolation or extrapolation by Y=RC*X+B     ---
234    90    DX=PLAB(J+1)-PLAB(J)
235 C --- Elastic cross section ---
236          DY=CSEL(IPART,J+1)-CSEL(IPART,J)
237          RC=DY/DX
238          B=CSEL(IPART,J)-RC*PLAB(J)
239          AIELIN=RC*P+B
240 C --- Inelastic cross section ---
241          DY=CSIN(IPART,J+1)-CSIN(IPART,J)
242          RC=DY/DX
243          B=CSIN(IPART,J)-RC*PLAB(J)
244          AIININ=RC*P+B
245          ALPH=ALPHA(IPART)
246          IF(IPART.LT.14) THEN
247             DY=ALPHAC(J+1)-ALPHAC(J)
248             RC=DY/DX
249             B=ALPHAC(J)-RC*PLAB(J)
250             CORFAC=RC*P+B
251             ALPH=ALPH*CORFAC
252 C
253             IPART3=IPART2(IPART-6)
254 C
255 C --- Elastic cross section ---
256             DY=CSEL(IPART3,J+1)-CSEL(IPART3,J)
257             RC=DY/DX
258             B=CSEL(IPART3,J)-RC*PLAB(J)
259             XSECEL=RC*P+B
260 C --- Inelastic cross section ---
261             DY=CSIN(IPART3,J+1)-CSIN(IPART3,J)
262             RC=DY/DX
263             B=CSIN(IPART3,J)-RC*PLAB(J)
264             XSECIN=RC*P+B
265 C
266          END IF
267 C
268          DO 100 K=1,KK
269             AIEL(K)=AIELIN
270             AIIN(K)=AIININ
271 C
272             IF (A(K) .GE. 1.5) THEN
273 C
274 C --- A-dependence from parametrization ---
275                CREL=1.0
276                CRIN=1.0
277 C --- Get medium bin  1=Hydr.  2=Al  3=Cu  4=Pb ---
278                I=3
279                IF (A(K) .LT. 50.0) I=2
280                IF (A(K) .GT. 100.0) I=4
281                IF ((IPART .EQ. 14) .OR. (IPART .EQ. 16)) THEN
282 C
283 C --- Protons and neutrons ---
284 C
285 C --- Elastic cross section ---
286                   DY=CSPNEL(I-1,J+1)-CSPNEL(I-1,J)
287                   RC=DY/DX
288                   B=CSPNEL(I-1,J)-RC*PLAB(J)
289                   XSECEL=RC*P+B
290 C --- Inelastic cross section ---
291                   DY=CSPNIN(I-1,J+1)-CSPNIN(I-1,J)
292                   RC=DY/DX
293                   B=CSPNIN(I-1,J)-RC*PLAB(J)
294                   XSECIN=RC*P+B
295                   IF (AIEL(K) .GE. 0.001) CREL=XSECEL/(0.36*AIEL(K)*
296      +            CSA(I)**1.17)
297                   AITOT=AIEL(K)+AIIN(K)
298                   IF (AITOT .GE. 0.001) CRIN=XSECIN/(AITOT*CSA(I)**
299      +            ALPH)
300 C
301                ELSEIF (IPART .LT. 15) THEN
302 C
303 C --- Calculate correction factors from values on Al,Cu,Pb for all ---
304 C --- mesons use linear interpolation or extrapolation by Y=RC*X+B ---
305 C --- Note that data is only available for pions and protons
306                   WGCH=0.5
307                   IF (A(K) .LT. 20.0) WGCH=0.5+0.5*EXP(-(A(K)-1.0))
308                   AIEL(K)=WGCH*AIEL(K)+(1.0-WGCH)*XSECEL
309                   AIIN(K)=WGCH*AIIN(K)+(1.0-WGCH)*XSECIN
310 C
311 C --- This section not for kaons ---
312                   IF (IPART .LT. 10) THEN
313 C
314 C --- Elastic cross section ---
315                      DY=CSPIEL(I-1,J+1)-CSPIEL(I-1,J)
316                      RC=DY/DX
317                      B=CSPIEL(I-1,J)-RC*PLAB(J)
318                      XSPIEL=RC*P+B
319 C --- Inelastic cross section ---
320                      DY=CSPIIN(I-1,J+1)-CSPIIN(I-1,J)
321                      RC=DY/DX
322                      B=CSPIIN(I-1,J)-RC*PLAB(J)
323                      XSPIIN=RC*P+B
324 C
325                      IF (AIEL(K) .GE. 0.001) CREL=XSPIEL/(0.36* AIEL(K)
326      +               *CSA(I)**1.17)
327                      AITOT=AIEL(K)+AIIN(K)
328                      IF (AITOT .GE. 0.001) CRIN=XSPIIN/(AITOT*CSA(I)
329      +               **ALPH)
330                   END IF
331                END IF
332                AIIN(K)=CRIN*(AIIN(K)+AIEL(K))*A(K)**ALPH
333                AIEL(K)=CREL*0.36*AIEL(K)*A(K)**1.17
334                AIEL(K)=AIEL(K)*PARTEL(IPART)
335                AIIN(K)=AIIN(K)*PARTIN(IPART)
336             END IF
337   100    CONTINUE
338 C
339 C --- Fission cross sections ---
340 C --- A-dependence given by  sigma(3 MeV)=-67.0+38.7*Z**(4/3)/A ---
341       END IF
342       I=21
343       DO 110 II=1,21
344          IF (EK .LT. EKFISS(II)) THEN
345             I=II
346             GO TO 120
347          END IF
348   110 CONTINUE
349 C
350   120 CONTINUE
351       DO 130 K=1,KK
352 C
353 C --- No fission for materials with A < 230 ---
354          IF (A(K) .GE. 230.) THEN
355 C
356 C --- Only data for U(233), U(235) and Pu(239) for Ek .le. 0.01 GeV ---
357             J=4
358             IF (EK .LE. 0.01) THEN
359 C
360 C --- Distinguish U(233), U(235), Pu(239) and U(238) and rest by J=1,4
361                IF ((Z(K) .EQ. 92.0).AND.(ABS(A(K)-233.0).LT.0.5))THEN
362                   J=1
363                ELSEIF((Z(K).EQ.92.0).AND.(ABS(A(K)-235.0).LT.0.5))THEN
364                   J=2
365                ELSEIF((Z(K).EQ.94.0).AND.(ABS(A(K)-239.0).LT.0.5))THEN
366                   J=3
367                END IF
368             END IF
369             IF(J.EQ.4) THEN
370 C
371                Z43BA=Z(K)**(4.0*ONETHR)/A(K)
372                Z43BA=MAX(-67.0+38.7*Z43BA,0.)
373             ELSE
374                Z43BA=1.
375             ENDIF
376 C
377 C --- Energy dependence taken from U(238) ---
378 C --- approximated as step-function ---
379             AIFI(K)=CSFISS(J,I)*Z43BA
380          END IF
381   130 CONTINUE
382 C
383 C --- Corrections for compounds ---
384 C --- These corrections should only be applied to inorganic scintill. ---
385 C --- apply the correction only if user selected it within GEANT ---
386       IF (QCOR .GT. 0.0.AND.KK.GT.1) THEN
387 C --- Do not apply corrections for anti-baryons ---
388          IF (ICORR(IPART).EQ.1) THEN
389 C
390 C --- ACC40 between 0.3 and 0.5 for Pi+ and Pi0, 0.2 for other mesons ---
391             ACC40=0.325
392 C           IF (IPART .GE. 9) ACC40=0.20
393 C --- ACC40 = 0.15 for baryons ---
394             IF (IPART .GE. 14) ACC40=0.15
395 C --- ACCA=0.08 for all pions, 0.02 for all other particles ---
396             ACCA= 0.08
397             IF (IPART .GE. 10) ACCA=0.02
398 C
399             ACCB=0.32*(ACC40-ACCA)
400             ACC=ACCA-ACCB*LOG(EK)
401             IF (ACC .GT. 0.5) ACC=0.5
402             IF (ACC .GT. 0.0) THEN
403 C
404                CAVER=AVER**ACC
405                CAVER=1.+(CAVER-1.)*QCOR
406                DO 140 K=1,KK
407                   AIEL(K)=AIEL(K)*CAVER
408                   AIIN(K)=AIIN(K)*CAVER
409                   AIFI(K)=AIFI(K)*CAVER
410                   AICA(K)=AICA(K)*CAVER
411   140          CONTINUE
412             END IF
413          END IF
414       END IF
415 C
416 C --- Calculate interaction probability ---
417 C
418 C --- Correction factor for high (P > 100 GeV/C) energies ---
419       CORH=1.0
420 C --- Assume a LOG(P) dependence of the correction factor with values ---
421 C P = 100 GeV/C  ==> CORH = 1.
422 C P =   1 TeV/C  ==> CORH = 1.25
423 C     DX=LOG(1000.)-LOG(100.)
424 C     DY=1.25-1.
425 C     RC=DY/DX
426 C     B=1.-RC*LOG(100.)
427 C     CORH=RC*LOG(P)+B
428       IF (P .GT. 100.) CORH=0.1085736156*LOG(P)+0.5
429       ALAM=0.0
430       DO 150 K=1,KK
431          AFACT=AVO*1E-3*DENS*W(K)/A(K)
432          AIEL(K)=MAX(CORH*AIEL(K)*AFACT,0.)
433          AIIN(K)=MAX(CORH*AIIN(K)*AFACT,0.)
434          AIFI(K)=MAX(CORH*AIFI(K)*AFACT,0.)
435          AICA(K)=MAX(CORH*AICA(K)*AFACT,0.)
436 C
437          ALAM=ALAM+AIEL(K)+AIIN(K)+AIFI(K)+AICA(K)
438   150 CONTINUE
439 C
440 C --- Pass the interaction probability to GEANT ---
441       GHESIG=ALAM
442 C
443       GO TO 999
444 C
445 C --- Printout of skipped particles in case of interface debug ---
446   160 IF (NPRT(9)) THEN
447          WRITE(CHMAIL,10300) IPART
448          CALL GMAIL(0,0)
449       END IF
450 10000 FORMAT(' *GHESIG* GEOM X-SECT. FOR INEL. SCAT. OF D,T AND ALPHA')
451 10100 FORMAT(' *GHESIG* X-SECT. FROM LOW ENERGY NEUTRON TABLES')
452 10200 FORMAT(' *GHESIG* X-SECT. FROM PARAMETRIZATION OF DATA')
453 10300 FORMAT('0*GHESIG* GHEISHA PARTICLE ',I3,' SKIPPED')
454   999 END