TGeant3/GeantPatch.F

   1 *CMZ :          02/02/99  18.09.13  by  Federico Carminati
   2 *-- Author :    Federico Carminati   02/02/99
   3       SUBROUTINE GFLUCT(DEMEAN,DE)
   4 C.
   5 C.    ******************************************************************
   6 C.    *                                                                *
   7 C.    *   Subroutine to decide which method is used to simulate        *
   8 C.    *   the straggling around the mean energy loss.                  *
   9 C.    *                                                                *
  10 C.    *                                                                *
  11 C.    *   DNMIN:  <---------->1<-------->30<--------->50<--------->    *
  12 C.    *                                                                *
  13 C.    *   LOSS=2  :                                                    *
  14 C.    *   STRA=0  <----------GLANDZ-------------------><--GLANDO-->    *
  15 C.    *   LOSS=1,3:                                                    *
  16 C.    *   STRA=0  <---------------------GLANDZ-------------------->    *
  17 C.    *                                                                *
  18 C.    *   STRA=1  <-----------PAI---------------------><--GLANDZ-->    *
  19 C.    *                                                                *
  20 C.    *   DNMIN :  an estimation of the number of collisions           *
  21 C.    *            with energy close to the ionization energy          *
  22 C.    *            (see PHYS333)                                       *
  23 C.    *                                                                *
  24 C.    *   Input  : DEMEAN (mean energy loss)                           *
  25 C.    *   Output : DE   (energy loss in the current step)              *
  26 C.    *                                                                *
  27 C.    *    ==>Called by : GTELEC,GTMUON,GTHADR                         *
  28 C.    *                                                                *
  29 C.    ******************************************************************
  30 C.
  31 #include "geant321/pilot.h"
  32 #undef CERNLIB_GEANT321_GCBANK_INC
  33 #undef CERNLIB_GEANT321_GCLINK_INC
  34 #include "geant321/gcbank.inc"
  35 #undef CERNLIB_GEANT321_GCJLOC_INC
  36 #include "geant321/gcjloc.inc"
  37 #undef CERNLIB_GEANT321_GCONSP_INC
  38 #include "geant321/gconsp.inc"
  39 #undef CERNLIB_GEANT321_GCMATE_INC
  40 #include "geant321/gcmate.inc"
  41 #undef CERNLIB_GEANT321_GCCUTS_INC
  42 #include "geant321/gccuts.inc"
  43 #undef CERNLIB_GEANT321_GCKINE_INC
  44 #include "geant321/gckine.inc"
  45 #undef CERNLIB_GEANT321_GCMULO_INC
  46 #include "geant321/gcmulo.inc"
  47 #undef CERNLIB_GEANT321_GCPHYS_INC
  48 #include "geant321/gcphys.inc"
  49 #undef CERNLIB_GEANT321_GCTRAK_INC
  50 #include "geant321/gctrak.inc"
  51 *KEND.
  52       PARAMETER (EULER=0.577215,GAM1=EULER-1)
  53       PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
  54      +           P5=.74442,P6=1.1934)
  55       PARAMETER (DGEV=0.153536 E-3, DNLIM=50)
  56 * These parameters are needed by M.Kowalski's fluctuation algorithm
  57       PARAMETER (FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2)
  58       PARAMETER (XEXPO=-EEXPO+1, YEXPO=1/XEXPO)
  59 * These parameters are needed by M.Kowalski's fluctuation algorithm
  60       DIMENSION RNDM(2)
  61       DE2(DPOT,RAN)=(DPOT**XEXPO*(1-RAN)+EEND**XEXPO*RAN)**YEXPO
  62       FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5)
  63       IF(STEP.LE.0) THEN
  64          DE=DEMEAN
  65       ELSE
  66          DEDX = DEMEAN/STEP
  67          POTI=Q(JPROB+9)
  68          IF(ISTRA.EQ.0.AND.(ILOSS.EQ.1.OR.ILOSS.EQ.3)) THEN
  69             CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10))
  70          ELSEIF (ILOSS.EQ.5) THEN
  71 * This is Marek Kowalski's fluctuation algorithm, it works only when
  72 * the step size has been limited to one ionisation on average
  73             CALL GRNDM(RNDM,1)
  74             DE=DE2(FPOT,RNDM(1))
  75 *
  76          ELSE
  77 * *** mean ionization potential (GeV)
  78 *        POTI=16E-9*Z**0.9
  79             GAMMA = GETOT/AMASS
  80             BETA = VECT(7)/GETOT
  81             BET2 = BETA**2
  82 * ***    low energy transfer
  83             XI = DGEV*CHARGE**2*STEP*DENS*Z/(A*BET2)
  84 * ***    regime
  85 * ***    ISTRA = 1 --> PAI + URBAN
  86             DNMIN = MIN(XI,DEMEAN)/POTI
  87             IF (ISTRA.EQ.0) THEN
  88                IF(DNMIN.GE.DNLIM) THEN
  89 *  Energy straggling using Gaussian, Landau & Vavilov theories.
  90 *  STEP   =  current step-length (cm)
  91 *  DELAND =  DE/DX - <DE/DX>     (GeV)
  92 *  Author      : G.N. Patrick
  93                   IF(STEP.LT.1.E-7)THEN
  94                      DELAND=0.
  95                   ELSE
  96 *     Maximum energy transfer to atomic electron (GeV).
  97                      ETA = BETA*GAMMA
  98                      RATIO = EMASS/AMASS
  99 *     Calculate Kappa significance ratio.
 100 *                 EMAX=(2*EMASS*ETA**2)/(1+2*RATIO*GAMMA+RATIO**2)
 101 *                 CAPPA = XI/EMAX
 102                      CAPPA = XI*(1+2*RATIO*GAMMA+RATIO**2)/(2*EMASS*
 103      +               ETA**2)
 104                      IF (CAPPA.GE.10.) THEN
 105 *     +-----------------------------------+
 106 *     I Sample from Gaussian distribution I
 107 *     +-----------------------------------+
 108                         SIGMA = XI*SQRT((1.-0.5*BET2)/CAPPA)
 109                         CALL GRNDM(RNDM,2)
 110                         F1 = -2.*LOG(RNDM(1))
 111                         DELAND = SIGMA*SQRT(F1)*COS(TWOPI*RNDM(2))
 112                      ELSE
 113                         XMEAN = -BET2-LOG(CAPPA)+GAM1
 114                         IF (CAPPA.LT.0.01) THEN
 115                            XLAMX = FLAND(XMEAN)
 116 *     +---------------------------------------------------------------+
 117 *     I Sample lambda variable from Kolbig/Schorr Landau distribution I
 118 *     +---------------------------------------------------------------+
 119 *  10                   CALL GRNDM(RNDM,1)
 120 *                       IF( RNDM(1) .GT. 0.980 ) GO TO 10
 121 *                       XLAMB = GLANDR(RNDM(1))
 122 *     +---------------------------------------------------------------+
 123 *     I Sample lambda variable from James/Hancock Landau distribution I
 124 *     +---------------------------------------------------------------+
 125    10                      CALL GLANDG(XLAMB)
 126                            IF(XLAMB.GT.XLAMX) GO TO 10
 127                         ELSE
 128 *            +---------------------------------------------------------+
 129 *            I Sample lambda variable (Landau not Vavilov) from        I
 130 *            I Rotondi&Montagna&Kolbig Vavilov distribution            I
 131 *            +---------------------------------------------------------+
 132                            CALL GRNDM(RNDM,1)
 133                            XLAMB = GVAVIV(CAPPA,BET2,RNDM(1))
 134                         ENDIF
 135 *     Calculate DE/DX - <DE/DX>
 136                         DELAND = XI*(XLAMB-XMEAN)
 137                      ENDIF
 138                   ENDIF
 139                   DE = DEMEAN + DELAND
 140                ELSE
 141                   CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
 142      +            Q(JPROB+ 10))
 143                ENDIF
 144             ELSE IF (ISTRA.LE.2) THEN
 145                IF(DNMIN.GE.DNLIM) THEN
 146                   CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
 147      +            Q(JPROB+ 10))
 148                ELSE
 149                   NMEC = NMEC+1
 150                   LMEC(NMEC)=109
 151                   DE = GSTREN(GAMMA,DCUTE,STEP)
 152 * ***   Add brem losses to ionisation
 153                   IF(ITRTYP.EQ.2) THEN
 154                      JBASE = LQ(JMA-1)+2*NEK1+IEKBIN
 155                      DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
 156                   ELSEIF(ITRTYP.EQ.5) THEN
 157                      JBASE = LQ(JMA-2)+NEK1+IEKBIN
 158                      DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
 159                   ENDIF
 160                ENDIF
 161             ENDIF
 162          ENDIF
 163       ENDIF
 164       END
 165 *CMZ :          16/02/99  19.24.38  by  Federico Carminati
 166 *CMZ :  2.03/01 28/08/98  09.33.11  by  Federico Carminati
 167 *CMZ :  2.01/00 18/06/98  17.34.13  by  Federico Carminati
 168 *CMZ :  2.00/05 28/05/98  19.04.13  by  Federico Carminati
 169 *-- Author :
 170       SUBROUTINE GTREVE
 171 C.
 172 C.    ******************************************************************
 173 C.    *                                                                *
 174 C.    *    SUBR. GTREVE                                                *
 175 C.    *                                                                *
 176 C.    *   Controls tracking of all particles belonging to the current  *
 177 C.    *    event.                                                      *
 178 C.    *                                                                *
 179 C.    *   Called by : GUTREV, called by GTRIG                          *
 180 C.    *   Authors   : R.Brun, F.Bruyant                                *
 181 C.    *                                                                *
 182 C.    ******************************************************************
 183 C.
 184 #include "geant321/pilot.h"
 185 #undef CERNLIB_GEANT321_GCBANK_INC
 186 #undef CERNLIB_GEANT321_GCLINK_INC
 187 #include "geant321/gcbank.inc"
 188 #undef CERNLIB_GEANT321_GCFLAG_INC
 189 #include "geant321/gcflag.inc"
 190 #undef CERNLIB_GEANT321_GCKINE_INC
 191 #include "geant321/gckine.inc"
 192 #undef CERNLIB_GEANT321_GCNUM_INC
 193 #include "geant321/gcnum.inc"
 194 #undef CERNLIB_GEANT321_GCSTAK_INC
 195 #include "geant321/gcstak.inc"
 196 #undef CERNLIB_GEANT321_GCTMED_INC
 197 #include "geant321/gctmed.inc"
 198 #undef CERNLIB_GEANT321_GCTRAK_INC
 199 #include "geant321/gctrak.inc"
 200 #undef CERNLIB_GEANT321_GCUNIT_INC
 201 #include "geant321/gcunit.inc"
 202 #include "sckine.inc"
 203 *KEND.
 204 *
 205       REAL UBUF(2)
 206       EQUIVALENCE (UBUF(1),WS(1))
 207       LOGICAL   BTEST
 208       DIMENSION PMOM(3),VPOS(3)
 209 C.
 210 C.    ------------------------------------------------------------------
 211       NTMSTO = 0
 212       NSTMAX = 0
 213       NALIVE = 0
 214 *         Kick start the creation of the vertex
 215       VPOS(1)=0
 216       VPOS(2)=0
 217       VPOS(3)=0
 218       PMOM(1)=0
 219       PMOM(2)=0
 220       PMOM(3)=0
 221       IPART=1
 222       CALL GSVERT(VPOS,0,0,UBUF,0,NVTX)
 223       CALL GSKINE(PMOM,IPART,NVTX,UBUF,0,NT)
 224 *
 225       MTRACK=-999
 226  10   MTROLD=MTRACK
 227       CALL RXGTRAK(MTRACK,IPART,PMOM,E,VPOS,TTOF)
 228       IF(MTROLD.LT.0) THEN
 229          MPRIMA=MTRACK
 230       ENDIF
 231       IF(MTRACK.LE.MPRIMA) THEN
 232          IF(ISWIT(4).GT.0.AND.MTRACK.GT.0) THEN
 233             IF(ISWIT(4).EQ.1.OR.MOD(MTRACK,ISWIT(4)).EQ.0) THEN
 234                WRITE(CHMAIL,10200) MTRACK
 235                CALL GMAIL(0,0)
 236             ENDIF
 237          ENDIF
 238          IF(MTROLD.GT.0) THEN
 239 C --- Output root hits tree only for each primary MTRACK
 240             CALL RXOUTH
 241          ENDIF
 242       ENDIF
 243       IF(MTRACK.LE.0) GOTO 999
 244 * Set the vertex
 245       JV=LQ(JVERTX-1)
 246       Q(JV + 1) = VPOS(1)
 247       Q(JV + 2) = VPOS(2)
 248       Q(JV + 3) = VPOS(3)
 249       Q(JV + 4) = TTOF
 250       Q(JV + 5) = 0
 251       Q(JV + 6) = 0
 252 * Set the track
 253       JK=LQ(JKINE-1)
 254       Q(JK + 1) = PMOM(1)
 255       Q(JK + 2) = PMOM(2)
 256       Q(JK + 3) = PMOM(3)
 257       Q(JK + 4) = E
 258       Q(JK + 5) = IPART
 259       Q(JK + 6) = 1
 260 * Now transport
 261 C      CALL GPVERT(0)
 262 C      CALL GPKINE(0)
 263 * Normal Gtreve code
 264       NV = NVERTX
 265       DO 40  IV = 1,NV
 266 * ***   For each vertex in turn ..
 267          JV = LQ(JVERTX-IV)
 268          NT = Q(JV+7)
 269          IF (NT.LE.0) GO TO 40
 270          TOFG   = Q(JV+4)
 271          SAFETY = 0.
 272          IF (NJTMAX.GT.0) THEN
 273             CALL GMEDIA (Q(JV+1), NUMED)
 274             IF (NUMED.EQ.0) THEN
 275                WRITE (CHMAIL, 10000) (Q(JV+I), I=1,3)
 276                CALL GMAIL (0, 0)
 277                GO TO 40
 278             ENDIF
 279             CALL GLSKLT
 280          ENDIF
 281 *  **   Loop over tracks attached to current vertex
 282          DO 20  IT = 1,NT
 283             JV   = LQ(JVERTX-IV)
 284             ITRA = Q(JV+7+IT)
 285             IF (BTEST(IQ(LQ(JKINE-ITRA)),0)) GO TO 20
 286             CALL GFKINE (ITRA, VERT, PVERT, IPART, IVERT, UBUF, NWBUF)
 287             IF (IVERT.NE.IV) THEN
 288                WRITE (CHMAIL, 10100) IV, IVERT
 289                CALL GMAIL (0, 0)
 290                GO TO 999
 291             ENDIF
 292 *   *      Store current track parameters in stack JSTAK
 293             CALL GSSTAK (2)
 294    20    CONTINUE
 295 *  **   Start tracking phase
 296    30    IF (NALIVE.NE.0) THEN
 297             NALIVE = NALIVE -1
 298 *   *      Pick-up next track in stack JSTAK, if any
 299 *   *         Initialize tracking parameters
 300             CALL GLTRAC
 301             IF (NUMED.EQ.0) GO TO 30
 302 *   *       Resume tracking
 303             CALL GUTRAK
 304             IF (IEOTRI.NE.0) GO TO 999
 305             GO TO 30
 306          ENDIF
 307 *
 308    40 CONTINUE
 309       GOTO 10
 310 *
 311 10000 FORMAT (' GTREVE : Vertex outside setup, XYZ=',3G12.4)
 312 10100 FORMAT (' GTREVE : Abnormal track/vertex connection',2I8)
 313 10200 FORMAT (' GTREVE : Transporting primary track No ',I10)
 314 *                                                             END GTREVE
 315   999 END