This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / giface / gmunu.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.38  by  S.Giani
11 *-- Author :
12       SUBROUTINE GMUNU
13 C
14 C *** GENERATION OF MUON-NUCLEUS INTERACTIONS ***
15 C *** NVE 16-MAR-1988 CERN GENEVA ***
16 C
17 C CALLED BY : GTMUON
18 C ORIGIN : H.FESEFELDT (ROUTINE CASMU)
19 C
20 C ***Revised Sep-90 by C. CHIERA, E. LAMANNA:
21 C ***Rebinning of vectors TETAL and XEML
22 C ***to avoid big angle for outgoing muons
23 C
24 C  Rebinning reinstated as original by H-J Trost. The correction of
25 C  the angle of the outgoing muon should take care of the anomalies
26 C  at large angles. 19-JUN-91.
27 C
28 C.    *  This routine is a straigth translation of the Gheisha routine *
29 C.    *  CASMU in the Geant dialect. The inelastic cross section is    *
30 C.    *  taken as 0.0003 mb for E<30 GeV, and is slowly increasing for *
31 C.    *  E>30 GeV. The energy and angle of final state muon is         *
32 C.    *  generated according to the 'free quark parton model'. Instead *
33 C.    *  of the virtual photon a real pion is written on working       *
34 C.    *  common in order to make use of the cascading routines for     *
35 C.    *  hadron production.                                            *
36 C.    ******************************************************************
37 C.
38 #include "geant321/gcbank.inc"
39 #include "geant321/gcphys.inc"
40 #include "geant321/gcjloc.inc"
41 #include "geant321/gckine.inc"
42 #include "geant321/gcking.inc"
43 #include "geant321/gconsp.inc"
44 #include "geant321/gctrak.inc"
45 #include "geant321/gsecti.inc"
46 #include "geant321/gcmulo.inc"
47 #if defined(CERNLIB_USRJMP)
48 #include "geant321/gcjump.inc"
49 #endif
50 C
51 #if !defined(CERNLIB_SINGLE)
52       DOUBLE PRECISION TOTPRO
53 #endif
54       DIMENSION VMUOUT(3),COEF(200),XEML(23),TETAL(35)
55       DIMENSION RNDM(3)
56       LOGICAL FIRST
57       SAVE COEF, FIRST, TETAL, XEML
58       DATA FIRST /.TRUE./
59       DATA TETAL /
60      A   1.0000000,  0.9999995,  0.9999990,  0.9999981,  0.9999962,
61      A   0.9999943,  0.9999905,  0.9999847,  0.9999752,  0.9999599,
62      A   0.9999352,  0.9998951,  0.9998302,  0.9997253,  0.9995556,
63      A   0.9992810,  0.9988368,  0.9981183,  0.9969561,  0.9950773,
64      A   0.9920409,  0.9871377,  0.9792297,  0.9665010,  0.9460785,
65      A   0.9134827,  0.8618938,  0.7813507,  0.6583430,  0.4770452,
66      A   0.2247237, -0.0955139, -0.4461272, -0.7495149, -0.9900000 /
67       DATA XEML /1.,.998,.997,.996,.995,.994,.992,.99,.97,.95,
68      +   .92,.89,.85,.8,.75,.7,.6,.5,.4,.3,.2,.1,.05/
69 C             COEF contains the value of the integral of the
70 C             double differential cross section ds/d(e1) d(cost)
71 C             for the production of the outgoing muon. These
72 C             values are obtained using the function
73 C             GMUSIG and are used to normalize the random value
74 C             used to sample the energy and angle of the outgoing
75 C             muon.
76 C
77       IF(FIRST) THEN
78 *
79 *   Integrate the double differential cross section
80 *
81          DO 8 ICOEF=1, NEKBIN+1
82             EINIT=ELOW(ICOEF)+EMMU
83             TOTPRO=0.0
84             DO 5 ICOST=2,35
85                COST = (TETAL(ICOST)+TETAL(ICOST-1))*0.5
86                DO 3 IEFIN=2,23
87                   EFINAL = (XEML(IEFIN)+XEML(IEFIN-1))*0.5*EINIT
88                   TOTPRO = TOTPRO + GMUSIG(EINIT,EFINAL,COST)*
89      +                     (TETAL(ICOST-1)-TETAL(ICOST))*
90      +                     (XEML(IEFIN-1)-XEML(IEFIN))*EINIT
91    3           CONTINUE
92    5        CONTINUE
93             COEF(ICOEF) = TOTPRO
94    8     CONTINUE
95          FIRST=.FALSE.
96       ENDIF
97 C
98       KCASE  = NAMEC(21)
99       IF(VECT(7).LT.0.01) GO TO 9999
100 C
101 C               Generate 4-momentum of final state muon
102 C
103 C --- IW2TRY loop introduced to avoid W2<0. (HJT/NVE 27-sep-1990) ---
104 C --- In case not successfull within 100 tries ==> No change made ---
105       IW2TRY=0
106  10   CONTINUE
107       IF (IW2TRY .GT. 100) GO TO 9999
108       IW2TRY=IW2TRY+1
109 C
110       TOTPRO = 0.0
111       CALL GRNDM(RNDM,1)
112       RAN    = RNDM(1)
113       HMAX   = (1.-GEKRAT)*COEF(IEKBIN)+GEKRAT*COEF(IEKBIN+1)
114       DO 14 ICOST=2,35
115         COST   = (TETAL(ICOST)+TETAL(ICOST-1))*.5
116         DO 13 IE1=2,23
117           E1     = (XEML(IE1)+XEML(IE1-1))*.5*GETOT
118           TOTPRO = TOTPRO+GMUSIG(GETOT,E1,COST)
119      *       *(TETAL(ICOST-1)-TETAL(ICOST))
120      *       *(XEML(IE1-1)-XEML(IE1))*GETOT
121           IF(RAN*HMAX.LT.TOTPRO) GOTO 15
122   13    CONTINUE
123         IE1    = 23
124   14  CONTINUE
125       ICOST= 35
126   15  CALL GRNDM(RNDM,3)
127       TETA = ACOS(TETAL(ICOST-1))+
128      *     RNDM(1)*(ACOS(TETAL(ICOST))-ACOS(TETAL(ICOST-1)))
129       COST = COS(TETA)
130       E1   = (XEML(IE1)+RNDM(2)*(XEML(IE1-1)-XEML(IE1)))*GETOT
131       IF(E1.LT.EMMU) E1=EMMU+0.0001
132       P1   = SQRT(ABS((E1-EMMU)*(E1+EMMU)))
133       IF (ABS(COST) .GT. 1.) COST=SIGN(1.,COST)
134 C
135 C --- Check value of W2 and in case negative ==> try again ---
136       W2=PMASS*(PMASS+2.*(GETOT-E1))-
137      +    2.*(GETOT*E1-VECT(7)*P1*COST-EMMU**2)
138       IF (W2 .LT. 0.) GO TO 10
139 C
140       SINT   = SQRT(ABS(1.-COST*COST))
141       PHI    = TWOPI*RNDM(3)
142       VMUOUT(1) = P1*SINT*COS(PHI)
143       VMUOUT(2) = P1*SINT*SIN(PHI)
144       VMUOUT(3) = P1*COST
145 C
146 C               Store muon on stack and write virtual photon on
147 C                     result common, rotate muon momenta
148 C
149       CALL GVROT(VECT(4),VMUOUT)
150       IF(IMUNU.EQ.1) THEN
151 C
152 C            Now compute momentum of the outgoing pion
153 C
154         VECT(4) =  VECT(4)*VECT(7)-VMUOUT(1)
155         VECT(5) =  VECT(5)*VECT(7)-VMUOUT(2)
156         VECT(6) =  VECT(6)*VECT(7)-VMUOUT(3)
157         VECT(7) =  SQRT(VECT(4)*VECT(4)+VECT(5)*VECT(5)+VECT(6)*VECT(6))
158         VECT(4) =  VECT(4)/VECT(7)
159         VECT(5) =  VECT(5)/VECT(7)
160         VECT(6) =  VECT(6)/VECT(7)
161 C
162 C            Select pi+ or pi-
163 C
164         IPOLD   = IPART
165         CALL GRNDM(RNDM,1)
166         IPART   = 8.5+RNDM(1)
167         JPA     = LQ(JPART-IPART)
168         DO 16 J=1,5
169   16    NAPART(J) = IQ(JPA+J)
170         ITRTYP    = Q(JPA+6)
171         AMASS     = Q(JPA+7)
172         CHARGE    = Q(JPA+8)
173         TLIFE     = Q(JPA+9)
174         GETOT     = SQRT(VECT(7)*VECT(7)+AMASS*AMASS)
175         GEKIN     = GETOT-AMASS
176         CALL GEKBIN
177 C
178 C             Now force interaction of the pion
179 C
180 C             This piece of code useful only if using the
181 C             Gheisha-Geant interface
182 C
183         STEPHA = BIG
184         IHOLD     = IHADR
185         IF(IHADR.NE.3) IHADR     = 1
186 #if !defined(CERNLIB_USRJMP)
187         CALL GUPHAD
188 #endif
189 #if defined(CERNLIB_USRJMP)
190         CALL JUMPT0(JUPHAD)
191 #endif
192         KK = Q(JMA+11)
193         DO 17 K=1,KK
194 C
195 C             Forbid elastic scattering
196 C
197           ALAM    = ALAM - AIEL(K)
198           AIEL(K) = 0.0
199   17    CONTINUE
200         NMOLD     = NMEC
201 #if !defined(CERNLIB_USRJMP)
202         CALL GUHADR
203 #endif
204 #if defined(CERNLIB_USRJMP)
205         CALL JUMPT0(JUHADR)
206 #endif
207         IHADR     = IHOLD
208         NMEC      = NMOLD
209         ISTOP     = 0
210         IPART     = IPOLD
211         JPA       = LQ(JPART-IPART)
212         DO 26 J=1,5
213   26    NAPART(J) = IQ(JPA+J)
214         ITRTYP    = Q(JPA+6)
215         AMASS     = Q(JPA+7)
216         CHARGE    = Q(JPA+8)
217         TLIFE     = Q(JPA+9)
218       ELSE
219         DESTEP = DESTEP+GETOT-SQRT(P1*P1+AMASS*AMASS)
220       ENDIF
221 C
222 C            Now just put the muon back on the current stack
223 C
224       VECT(7)   =
225      + SQRT(VMUOUT(1)*VMUOUT(1)+VMUOUT(2)*VMUOUT(2)+VMUOUT(3)*VMUOUT(3))
226       VECT(4)   = VMUOUT(1)/VECT(7)
227       VECT(5)   = VMUOUT(2)/VECT(7)
228       VECT(6)   = VMUOUT(3)/VECT(7)
229       GETOT     = SQRT(VECT(7)*VECT(7)+AMASS*AMASS)
230       GEKIN     = GETOT-AMASS
231       CALL GEKBIN
232 C
233 C             Update probabilities
234 C
235   90  SLMUNU=SLENG
236       ZINTMU=GARNDM(DUM)
237       STEPMU=BIG
238 C
239  9999 CONTINUE
240 C
241       RETURN
242       END