[u/mrichter/AliRoot.git] / GEANT321 / giface / gmunu.F

*
* $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 GMUNU
C
C *** GENERATION OF MUON-NUCLEUS INTERACTIONS ***
C *** NVE 16-MAR-1988 CERN GENEVA ***
C
C CALLED BY : GTMUON
C ORIGIN : H.FESEFELDT (ROUTINE CASMU)
C
C ***Revised Sep-90 by C. CHIERA, E. LAMANNA:
C ***Rebinning of vectors TETAL and XEML
C ***to avoid big angle for outgoing muons
C
C  Rebinning reinstated as original by H-J Trost. The correction of
C  the angle of the outgoing muon should take care of the anomalies
C  at large angles. 19-JUN-91.
C
C.    *  This routine is a straigth translation of the Gheisha routine *
C.    *  CASMU in the Geant dialect. The inelastic cross section is    *
C.    *  taken as 0.0003 mb for E<30 GeV, and is slowly increasing for *
C.    *  E>30 GeV. The energy and angle of final state muon is         *
C.    *  generated according to the 'free quark parton model'. Instead *
C.    *  of the virtual photon a real pion is written on working       *
C.    *  common in order to make use of the cascading routines for     *
C.    *  hadron production.                                            *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gcphys.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gckine.inc"
#include "geant321/gcking.inc"
#include "geant321/gconsp.inc"
#include "geant321/gctrak.inc"
#include "geant321/gsecti.inc"
#include "geant321/gcmulo.inc"
#if defined(CERNLIB_USRJMP)
#include "geant321/gcjump.inc"
#endif
C
#if !defined(CERNLIB_SINGLE)
      DOUBLE PRECISION TOTPRO
#endif
      DIMENSION VMUOUT(3),COEF(200),XEML(23),TETAL(35)
      DIMENSION RNDM(3)
      LOGICAL FIRST
      SAVE COEF, FIRST, TETAL, XEML
      DATA FIRST /.TRUE./
      DATA TETAL /
     A   1.0000000,  0.9999995,  0.9999990,  0.9999981,  0.9999962,
     A   0.9999943,  0.9999905,  0.9999847,  0.9999752,  0.9999599,
     A   0.9999352,  0.9998951,  0.9998302,  0.9997253,  0.9995556,
     A   0.9992810,  0.9988368,  0.9981183,  0.9969561,  0.9950773,
     A   0.9920409,  0.9871377,  0.9792297,  0.9665010,  0.9460785,
     A   0.9134827,  0.8618938,  0.7813507,  0.6583430,  0.4770452,
     A   0.2247237, -0.0955139, -0.4461272, -0.7495149, -0.9900000 /
      DATA XEML /1.,.998,.997,.996,.995,.994,.992,.99,.97,.95,
     +   .92,.89,.85,.8,.75,.7,.6,.5,.4,.3,.2,.1,.05/
C             COEF contains the value of the integral of the
C             double differential cross section ds/d(e1) d(cost)
C             for the production of the outgoing muon. These
C             values are obtained using the function
C             GMUSIG and are used to normalize the random value
C             used to sample the energy and angle of the outgoing
C             muon.
C
      IF(FIRST) THEN
*
*   Integrate the double differential cross section
*
         DO 8 ICOEF=1, NEKBIN+1
            EINIT=ELOW(ICOEF)+EMMU
            TOTPRO=0.0
            DO 5 ICOST=2,35
               COST = (TETAL(ICOST)+TETAL(ICOST-1))*0.5
               DO 3 IEFIN=2,23
                  EFINAL = (XEML(IEFIN)+XEML(IEFIN-1))*0.5*EINIT
                  TOTPRO = TOTPRO + GMUSIG(EINIT,EFINAL,COST)*
     +                     (TETAL(ICOST-1)-TETAL(ICOST))*
     +                     (XEML(IEFIN-1)-XEML(IEFIN))*EINIT
   3           CONTINUE
   5        CONTINUE
            COEF(ICOEF) = TOTPRO
   8     CONTINUE
         FIRST=.FALSE.
      ENDIF
C
      KCASE  = NAMEC(21)
      IF(VECT(7).LT.0.01) GO TO 9999
C
C               Generate 4-momentum of final state muon
C
C --- IW2TRY loop introduced to avoid W2<0. (HJT/NVE 27-sep-1990) ---
C --- In case not successfull within 100 tries ==> No change made ---
      IW2TRY=0
 10   CONTINUE
      IF (IW2TRY .GT. 100) GO TO 9999
      IW2TRY=IW2TRY+1
C
      TOTPRO = 0.0
      CALL GRNDM(RNDM,1)
      RAN    = RNDM(1)
      HMAX   = (1.-GEKRAT)*COEF(IEKBIN)+GEKRAT*COEF(IEKBIN+1)
      DO 14 ICOST=2,35
        COST   = (TETAL(ICOST)+TETAL(ICOST-1))*.5
        DO 13 IE1=2,23
          E1     = (XEML(IE1)+XEML(IE1-1))*.5*GETOT
          TOTPRO = TOTPRO+GMUSIG(GETOT,E1,COST)
     *       *(TETAL(ICOST-1)-TETAL(ICOST))
     *       *(XEML(IE1-1)-XEML(IE1))*GETOT
          IF(RAN*HMAX.LT.TOTPRO) GOTO 15
  13    CONTINUE
        IE1    = 23
  14  CONTINUE
      ICOST= 35
  15  CALL GRNDM(RNDM,3)
      TETA = ACOS(TETAL(ICOST-1))+
     *     RNDM(1)*(ACOS(TETAL(ICOST))-ACOS(TETAL(ICOST-1)))
      COST = COS(TETA)
      E1   = (XEML(IE1)+RNDM(2)*(XEML(IE1-1)-XEML(IE1)))*GETOT
      IF(E1.LT.EMMU) E1=EMMU+0.0001
      P1   = SQRT(ABS((E1-EMMU)*(E1+EMMU)))
      IF (ABS(COST) .GT. 1.) COST=SIGN(1.,COST)
C
C --- Check value of W2 and in case negative ==> try again ---
      W2=PMASS*(PMASS+2.*(GETOT-E1))-
     +    2.*(GETOT*E1-VECT(7)*P1*COST-EMMU**2)
      IF (W2 .LT. 0.) GO TO 10
C
      SINT   = SQRT(ABS(1.-COST*COST))
      PHI    = TWOPI*RNDM(3)
      VMUOUT(1) = P1*SINT*COS(PHI)
      VMUOUT(2) = P1*SINT*SIN(PHI)
      VMUOUT(3) = P1*COST
C
C               Store muon on stack and write virtual photon on
C                     result common, rotate muon momenta
C
      CALL GVROT(VECT(4),VMUOUT)
      IF(IMUNU.EQ.1) THEN
C
C            Now compute momentum of the outgoing pion
C
        VECT(4) =  VECT(4)*VECT(7)-VMUOUT(1)
        VECT(5) =  VECT(5)*VECT(7)-VMUOUT(2)
        VECT(6) =  VECT(6)*VECT(7)-VMUOUT(3)
        VECT(7) =  SQRT(VECT(4)*VECT(4)+VECT(5)*VECT(5)+VECT(6)*VECT(6))
        VECT(4) =  VECT(4)/VECT(7)
        VECT(5) =  VECT(5)/VECT(7)
        VECT(6) =  VECT(6)/VECT(7)
C
C            Select pi+ or pi-
C
        IPOLD   = IPART
        CALL GRNDM(RNDM,1)
        IPART   = 8.5+RNDM(1)
        JPA     = LQ(JPART-IPART)
        DO 16 J=1,5
  16    NAPART(J) = IQ(JPA+J)
        ITRTYP    = Q(JPA+6)
        AMASS     = Q(JPA+7)
        CHARGE    = Q(JPA+8)
        TLIFE     = Q(JPA+9)
        GETOT     = SQRT(VECT(7)*VECT(7)+AMASS*AMASS)
        GEKIN     = GETOT-AMASS
        CALL GEKBIN
C
C             Now force interaction of the pion
C
C             This piece of code useful only if using the
C             Gheisha-Geant interface
C
        STEPHA = BIG
        IHOLD     = IHADR
        IF(IHADR.NE.3) IHADR     = 1
#if !defined(CERNLIB_USRJMP)
        CALL GUPHAD
#endif
#if defined(CERNLIB_USRJMP)
        CALL JUMPT0(JUPHAD)
#endif
        KK = Q(JMA+11)
        DO 17 K=1,KK
C
C             Forbid elastic scattering
C
          ALAM    = ALAM - AIEL(K)
          AIEL(K) = 0.0
  17    CONTINUE
        NMOLD     = NMEC
#if !defined(CERNLIB_USRJMP)
        CALL GUHADR
#endif
#if defined(CERNLIB_USRJMP)
        CALL JUMPT0(JUHADR)
#endif
        IHADR     = IHOLD
        NMEC      = NMOLD
        ISTOP     = 0
        IPART     = IPOLD
        JPA       = LQ(JPART-IPART)
        DO 26 J=1,5
  26    NAPART(J) = IQ(JPA+J)
        ITRTYP    = Q(JPA+6)
        AMASS     = Q(JPA+7)
        CHARGE    = Q(JPA+8)
        TLIFE     = Q(JPA+9)
      ELSE
        DESTEP = DESTEP+GETOT-SQRT(P1*P1+AMASS*AMASS)
      ENDIF
C
C            Now just put the muon back on the current stack
C
      VECT(7)   =
     + SQRT(VMUOUT(1)*VMUOUT(1)+VMUOUT(2)*VMUOUT(2)+VMUOUT(3)*VMUOUT(3))
      VECT(4)   = VMUOUT(1)/VECT(7)
      VECT(5)   = VMUOUT(2)/VECT(7)
      VECT(6)   = VMUOUT(3)/VECT(7)
      GETOT     = SQRT(VECT(7)*VECT(7)+AMASS*AMASS)
      GEKIN     = GETOT-AMASS
      CALL GEKBIN
C
C             Update probabilities
C
  90  SLMUNU=SLENG
      ZINTMU=GARNDM(DUM)
      STEPMU=BIG
C
 9999 CONTINUE
C
      RETURN
      END
Commit	Line	Data
fe4da5cc	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.5EINIT
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)+GEKRATCOEF(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)).5GETOT
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.(GETOTE1-VECT(7)P1COST-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) = P1SINTCOS(PHI)
143	VMUOUT(2) = P1SINTSIN(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)+AMASSAMASS)
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(P1P1+AMASSAMASS)
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)+AMASSAMASS)
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