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