TFluka/crnkvp.f

   1 *$ CREATE CRNKVP.FOR
   2 *COPY CRNKVP
   3 *
   4 *=== Crnkvp ===========================================================*
   5 *
   6       SUBROUTINE CRNKVP ( MREG, NEWREG, DDNEAR, LEMAGN, DEDXCK )
   7
   8       INCLUDE '(DBLPRC)'
   9       INCLUDE '(DIMPAR)'
  10       INCLUDE '(IOUNIT)'
  11 *
  12 *----------------------------------------------------------------------*
  13 *                                                                      *
  14 *     Copyright (C) 1997-2003      by    Alfredo Ferrari & Paola Sala  *
  15 *     All Rights Reserved.                                             *
  16 *                                                                      *
  17 *                                                                      *
  18 *     CeReNKoV Production:                                             *
  19 *                                                                      *
  20 *     Created on 21 september 1997 by    Alfredo Ferrari & Paola Sala  *
  21 *                                                   Infn - Milan       *
  22 *                                                                      *
  23 *     Last change on 30-oct-98     by    Alfredo Ferrari               *
  24 *                                                                      *
  25 *    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    *
  26 *    !!!!!! Important: Etrack, Ptrack, Atrack,  are supposed !!!!!!    *
  27 *    !!!!!! to be the INITIAL values (not yet updated at the !!!!!!    *
  28 *    !!!!!! end of the step).                                !!!!!!    *
  29 *    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    *
  30 *                                                                      *
  31 *    dN/dEdx = alpha^2 z^2 / [ r_e m_e c^2 ] [ 1 - 1 / (beta n)^2 ]    *
  32 *                                                                      *
  33 *    Note the wavelength is always assumed to be the vacuum one        *
  34 *                                                                      *
  35 *----------------------------------------------------------------------*
  36 *
  37       INCLUDE '(FHEAVY)'
  38       INCLUDE '(FLKMAT)'
  39       INCLUDE '(MULBOU)'
  40       INCLUDE '(OPPHCM)'
  41       INCLUDE '(OPPHST)'
  42       INCLUDE '(PAPROP)'
  43       INCLUDE '(SUMCOU)'
  44       INCLUDE '(TRACKR)'
  45 *
  46       PARAMETER ( CSNPRN = 100.D+00 * CSNNRM )
  47       PARAMETER ( CRNKC0 = ALPFSC * ALPFSC / RCLSEL / AMELCT )
  48       LOGICAL LEMAGN, LMEQNT, LNWSUB, LSMPAN, LDUMMY, LBXRFL
  49 *  Statement functions:
  50 *     INCLUDE '(DORTSF)'
  51       INCLUDE '(UNRTSF)'
  52 *
  53 *  No magnetic field end step Newreg check implemented for the moment:
  54 *
  55 *
  56 *      IF ( LEMAGN )
  57 *     &   CALL FLABRT ( 'CRNKVP', ' STOP:LEMAGN-NOT-YET-IMPLEMENTED' )
  58 *  No change of lattice check implemented for the moment:
  59 *      IF ( LT1TRK .NE. LT2TRK ) CALL FLABRT ( 'CRNKVP',
  60 *     &    ' STOP:LT1TRK.NE.LT2TRK-NOT-YET-IMPLEMENTED' )
  61       DEDXCK = ZERZER
  62       MMAT   = MEDIUM (MREG)
  63       IF ( JTRACK .GE. -6 ) THEN
  64          AMPRTC = AM (JTRACK)
  65       ELSE
  66          AMPRTC = AMNHEA (-JTRACK)
  67       END IF
  68 *  +-------------------------------------------------------------------*
  69 *  |  Find the maximum refractive index over the interval allowed
  70 *  |  for Cerenkov production:
  71       IF ( RMXCER (MMAT) .LT. ZERZER ) THEN
  72          RMXCER (MMAT) = RFNDMX (MMAT)
  73       END IF
  74 *  |
  75 *  +-------------------------------------------------------------------*
  76       BETNCR = PTRACK / ETRACK * RMXCER (MMAT)
  77 *  Maximum beta x n already below 1: no chance to emit any photon
  78       IF ( BETNCR .LE. ONEONE ) RETURN
  79       LMEQNT = MTRACK .EQ. NTRACK
  80       DTRCKT = ZERZER
  81 *  +-------------------------------------------------------------------*
  82 *  |  Loop on energy deposition sub-steps:
  83       DO 300 ID = 1, MTRACK
  84          DTRCKT = DTRCKT + DTRACK (ID)
  85   300 CONTINUE
  86 *  |
  87 *  +-------------------------------------------------------------------*
  88       TTRCKT = ZERZER
  89 *  +-------------------------------------------------------------------*
  90 *  |  Loop on track sub-steps (be careful with magnetic field they are
  91 *  |  not the straight distance between the end points):
  92       DO 500 IT = 1, NTRACK
  93          TTRCKT = TTRCKT + TTRACK (IT)
  94   500 CONTINUE
  95 *  |
  96 *  +-------------------------------------------------------------------*
  97       CRVCRR = CTRACK / TTRCKT
  98       IF ( JTRACK .GE. -6 ) THEN
  99          CRNKCS = CRNKC0 * ( DBLE ( ICHRGE (JTRACK) ) )**2 * OPSNMX
 100       ELSE
 101          CRNKCS = CRNKC0 * ( DBLE ( ICHEAV(-JTRACK) ) )**2 * OPSNMX
 102       END IF
 103       BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
 104 *  Macroscopic sigma for production (cm^-1)
 105       SIGMCK = CRNKCS * ( EMXCER (MMAT) - EMNCER (MMAT) ) * BTNFCR
 106 *  Average number of emitted photons:
 107       NEMPHO = NINT ( CTRACK * SIGMCK )
 108 *  Take a three sigma clearance:
 109       NEMPHO = NEMPHO + 10 + 3 * NINT ( SQRT ( CTRACK * SIGMCK ) )
 110       NEMPHO = MIN ( NEMPHO, MOSTCK )
 111 *D     IF ( LSTOPP .GT. MOSTCK - NEMPHO ) WRITE (77,*)
 112 *D    &   ' ###CRNKVP:LSTOPP:',LSTOPP,NEMPHO
 113 *  +-------------------------------------------------------------------*
 114 *  |  Empty the stack if too full: compute and save the normal in case
 115 *  |  it is required
 116       IF ( LSTOPP .GT. MOSTCK - NEMPHO ) THEN
 117          LDUMMY = LBXRFL ( NEWREG, MREG, ZERZER, ZERZER, ONEONE )
 118          CALL KASOPH ( .TRUE. )
 119       END IF
 120 *  |
 121 *  +-------------------------------------------------------------------*
 122 *  Total wiggled and curved track used up to now:
 123       TRUSED = ZERZER
 124       ETINTR = ETRACK
 125       PTINTR = PTRACK
 126       TRINTR = ZERZER
 127       NTRKLD = 0
 128 *  Current lattice number and index inside Mulbou:
 129       IRLTCR = 0
 130       LATCCR = MULTTC (IRLTCR)
 131 *  Variables in the current Ntrckr bin:
 132       NTRKCR = 1
 133       PTRKCR = PTRACK
 134       ETRKCR = ETRACK
 135       ATRKCR = ATRACK
 136       TRNRSD = TTRACK (NTRKCR) * CRVCRR
 137 *  +-------------------------------------------------------------------*
 138 *  |  Ntrack = Mtrack
 139       IF ( LMEQNT ) THEN
 140          DEDXCR = DTRACK (NTRKCR) / TRNRSD
 141 *  |
 142 *  +-------------------------------------------------------------------*
 143 *  |  Ntrack >< Mtrack
 144       ELSE
 145          DEDXCR = DTRCKT / CTRACK
 146       END IF
 147 *  |
 148 *  +-------------------------------------------------------------------*
 149 *  +-------------------------------------------------------------------*
 150 *  |  Stacking loop for Cerenkov photons:
 151       NPROD = 0
 152  1000 CONTINUE
 153 *D        IF ( SIGMCK .LT. ZERZER ) WRITE (77,*)
 154 *D    &      ' ^^^CRNKVP:SIGMCK,BTNFCR',SIGMCK,BTNFCR
 155 *D        IF ( DEDXCR .LT. AZRZRZ ) WRITE (77,*)' ###CRNKVP:',
 156 *D    &  'MMAT,JTRACK,DEDXCR,NTRKCR,DTRACK(NTRKCR),CTRACK,TRNRSD,DTRCKT',
 157 *D    &   MMAT,JTRACK,DEDXCR,NTRKCR,DTRACK(NTRKCR),CTRACK,TRNRSD,DTRCKT
 158 *  |  Empty the stack if too full:
 159 *D        IF ( LSTOPP .EQ. MOSTCK ) WRITE (77,*)
 160 *D    &      ' ###CRNKVP:LSTOPP:',LSTOPP
 161 *  |  +----------------------------------------------------------------*
 162 *  |  |  Empty the stack if too full: compute and save the normal in
 163 *  |  |  case it is required
 164          IF ( LSTOPP .EQ. MOSTCK ) THEN
 165             LDUMMY = LBXRFL ( NEWREG, MREG, ZERZER, ZERZER, ONEONE )
 166             CALL KASOPH ( .TRUE. )
 167          END IF
 168 *  |  |
 169 *  |  +----------------------------------------------------------------*
 170          DSTPRD = - LOG ( ONEONE - FLRNDM (DSTPRD) ) / SIGMCK
 171          TRUSED = TRUSED + DSTPRD
 172 *dbgD        WRITE (77,*) ' CRNKVP:DSTPRD,TRUSED,CTRACK,TRNRSD',
 173 *dbgD    &                         DSTPRD,TRUSED,CTRACK,TRNRSD
 174          IF ( TRUSED .GE. CTRACK ) GO TO 7000
 175 *  |  +----------------------------------------------------------------*
 176 *  |  |  Sub-step loop:
 177  1500    CONTINUE
 178 *D           IF ( NTRKCR .GT. NTRACK ) WRITE (77,*)
 179 *D    &         ' ^^^CRNKVP:NTRKCR,NTRACK',NTRKCR,NTRACK
 180 *  |  |  +-------------------------------------------------------------*
 181 *  |  |  |  The production point is inside the current sub-step:
 182             IF ( DSTPRD .LE. TRNRSD ) THEN
 183 *  |  |  |  Update the distance from the previous interaction:
 184                TRINTR = TRINTR + DSTPRD
 185 *  |  |  |  Update the residual distance available in the current
 186 *  |  |  |  sub-step
 187                TRNRSD = TRNRSD - DSTPRD
 188                ETRKCR = ETRKCR - DSTPRD * DEDXCR
 189                LNWSUB = .FALSE.
 190 *  |  |  |
 191 *  |  |  +-------------------------------------------------------------*
 192 *  |  |  |  The production point is outside the current sub-step:
 193             ELSE
 194 *  |  |  |  Update the distance from the previous interaction:
 195                TRINTR = TRINTR + TRNRSD
 196                ETRKCR = ETRKCR - TRNRSD * DEDXCR
 197                DSTPRD = DSTPRD - TRNRSD
 198 *  |  |  |  Change sub-step:
 199                NTRKCR = NTRKCR + 1
 200 *  |  |  |  Update the residual distance available in the current
 201 *  |  |  |  sub-step
 202                TRNRSD = TTRACK (NTRKCR) * CRVCRR
 203 *  |  |  |  Update the current dE/dx:
 204                IF ( LMEQNT ) DEDXCR = DTRACK (NTRKCR) / TRNRSD
 205                LNWSUB = .TRUE.
 206             END IF
 207 *  |  |  |
 208 *  |  |  +-------------------------------------------------------------*
 209          IF ( LNWSUB ) GO TO 1500
 210 *  |  |  end sub-step loop:
 211 *  |  +----------------------------------------------------------------*
 212          PTRKCR = SQRT ( ( ETRKCR - AMPRTC ) * ( ETRKCR + AMPRTC ) )
 213 *  |  +----------------------------------------------------------------*
 214 *  |  |  Update the particle "age" (beta=p/E, dp/dE=E/p):
 215 *  |  |           Etrkcr
 216 *  |  |         /
 217 *  |  |  Dt = - | dE dx/dE E/(pc) = 1/c dx/dE [ Ptintr - Ptrkcr ]
 218 *  |  |         /
 219 *  |  |          Etintr
 220          IF ( ETINTR - ETRKCR .GT. CSNNRM * ETRKCR ) THEN
 221 *  |  |  Average dE/dx from the previous (fictitious) interaction:
 222             DEDXAV = ( ETINTR - ETRKCR ) / TRINTR
 223             ATRKCR = ATRKCR + ( PTINTR - PTRKCR ) / DEDXAV / CLIGHT
 224 *  |  |
 225 *  |  +----------------------------------------------------------------*
 226 *  |  |
 227          ELSE
 228             ATRKCR = ATRKCR + TRINTR * ETRKCR / PTRKCR / CLIGHT
 229          END IF
 230 *  |  |
 231 *  |  +----------------------------------------------------------------*
 232          PTINTR = PTRKCR
 233          ETINTR = ETRKCR
 234          TRINTR = ZERZER
 235          BETNLD = BETNCR
 236          BETNCR = PTRKCR / ETRKCR * RMXCER (MMAT)
 237 *dbgD        WRITE (77,*) ' CRNKVP:BETNCR,RMXCER,LSTOPP',
 238 *dbgD    &                         BETNCR,RMXCER(MMAT),LSTOPP
 239 *  |  No longer any chance to emit Cerenkov photons:
 240          IF ( BETNCR .LE. ONEONE ) GO TO 7000
 241          BTNFLD = BTNFCR
 242          BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
 243          FREJE  = BTNFCR / BTNFLD
 244 *D        IF ( FREJE .GT. ONEPLS ) WRITE (77,*) ' ^^^CRNKVP:',
 245 *D    &       'FREJE,BETNCR,BETNLD,DEDXCR,PTRKCR,ETRKCR,PTRACK,ETRACK',
 246 *D    &        FREJE,BETNCR,BETNLD,DEDXCR,PTRKCR,ETRKCR,PTRACK,ETRACK
 247 *  |  Update the macroscopic sigma: beta can only decrease
 248          SIGMCK = SIGMCK * FREJE
 249          RNDREJ = FLRNDM (RNDREJ)
 250          IF ( RNDREJ .GE. FREJE ) GO TO 1000
 251 *  |--<--<--<--<--< Interaction rejected because of the decrease in beta
 252          RNDREJ = RNDREJ / FREJE
 253 *  |  Sample the emitted photon energy:
 254          EPHSMP = RNDREJ * ( EMXCER (MMAT) - EMNCER (MMAT) )
 255      &          + EMNCER (MMAT)
 256 *  |  2pi x freq.:
 257          OMGSMP = GEVOMG * EPHSMP
 258 *  |  Wavelength:
 259          WVLSMP = TWOPIP * CLIGHT / OMGSMP
 260 *  |  Compute the refraction index at the right wavelength/frequency:
 261          RFISMP = FOPTPR ( 1, WVLSMP, OMGSMP, MMAT )
 262          BETNSM = PTRKCR / ETRKCR * RFISMP
 263 *dbgD        WRITE (77,*) ' CRNKVP:BETNSM,RFISMP,EPHSMP,WVLSMP,OMGSMP',
 264 *dbgD    &                         BETNSM,RFISMP,EPHSMP,WVLSMP,OMGSMP
 265          IF ( BETNSM .LE. ONEONE ) GO TO 1000
 266 *  |--<--<--<--<--< No chance to emit Cerenkov photons at this frequency
 267          BTNFSM = ( BETNSM - ONEONE ) * ( BETNSM + ONEONE ) / BETNSM**2
 268 *  |  Compute the quantum efficiency for this emitted energy:
 269          OPSENS = FOPTSN ( WVLSMP, OMGSMP )
 270          FREJE  = BTNFSM / BTNFCR * OPSENS / OPSNMX
 271 *D        IF ( FREJE .GT. ONEPLS ) WRITE (77,*)
 272 *D    &      ' ^^^CRNKVP:FREJE,RFISMP,RMXCER(MMAT),OPSENS,MMAT',
 273 *D    &                  FREJE,RFISMP,RMXCER(MMAT),OPSENS,MMAT
 274          RNDREJ = FLRNDM (RNDREJ)
 275 *dbgD        WRITE (77,*) ' CRNKVP:BTNFSM,BTNFCR,RNDREJ',
 276 *dbgD    &                         BTNFSM,BTNFCR,RNDREJ
 277          IF ( RNDREJ .GE. FREJE ) GO TO 1000
 278 *  |--<--<--<--<--< Interaction rejected because of the decrease in n
 279 *  |  Eventually a Cerenkov photon is going to be emitted:
 280 *  |  +----------------------------------------------------------------*
 281 *  |  |  New production sub-step:
 282          IF ( NTRKCR .NE. NTRKLD ) THEN
 283             UTRKCR = XTRACK (NTRKCR) - XTRACK (NTRKCR-1)
 284             VTRKCR = YTRACK (NTRKCR) - YTRACK (NTRKCR-1)
 285             WTRKCR = ZTRACK (NTRKCR) - ZTRACK (NTRKCR-1)
 286             TUVWCR = SQRT ( UTRKCR**2 + VTRKCR**2 + WTRKCR**2 )
 287             UTRKCR = UTRKCR / TUVWCR
 288             VTRKCR = VTRKCR / TUVWCR
 289             WTRKCR = WTRKCR / TUVWCR
 290             SINT02 = UTRKCR**2 + VTRKCR**2
 291             LSMPAN = SINT02 .LT. ANGLSQ
 292             COSTH0 = WTRKCR
 293             IF ( .NOT. LSMPAN ) THEN
 294                SINTH0 = SQRT (SINT02)
 295                COSPH0 = UTRKCR / SINTH0
 296                SINPH0 = VTRKCR / SINTH0
 297             END IF
 298             NTRKLD = NTRKCR
 299          END IF
 300 *  |  |
 301 *  |  +----------------------------------------------------------------*
 302          SRNRSD = TUVWCR * TRNRSD / CRVCRR / TTRACK (NTRKCR)
 303 *D        IF ( SRNRSD .LT. ZERZER .OR. SRNRSD .GT. TUVWCR )
 304 *D    &      WRITE (77,*)' ^^^CRNKVP:SRNRSD,TUVWCR,NTRKCR',
 305 *D    &                              SRNRSD,TUVWCR,NTRKCR
 306          XTRKCR = XTRACK (NTRKCR) - UTRKCR * SRNRSD
 307          YTRKCR = YTRACK (NTRKCR) - VTRKCR * SRNRSD
 308          ZTRKCR = ZTRACK (NTRKCR) - WTRKCR * SRNRSD
 309 *  |  +----------------------------------------------------------------*
 310 *  |  |  Put here the check on lattice change:
 311          IF ( NULTTC .GT. 0 ) THEN
 312 *  |  |  Current (possibly curved in mag. field) non wiggled path:
 313             TRLATT = TRUSED / CRVCRR
 314             IF ( TRLATT .LE. TSLTTC (IRLTCR) ) THEN
 315                LATCCR = MULTTC (IRLTCR)
 316             ELSE
 317                DO 2000 IL = IRLTCR + 1, NULTTC
 318                   IF ( TRLATT .LE. TSLTTC (IL) ) THEN
 319                      IRLTCR = IL
 320                      LATCCR = MULTTC (IRLTCR)
 321                      GO TO 2010
 322                   END IF
 323  2000          CONTINUE
 324 *D              CALL FLABRT ( 'CRNKVP', 'STOP:CRNKVP-NO-VALID-LATTICE' )
 325  2010          CONTINUE
 326             END IF
 327          END IF
 328 *  |  |
 329 *  |  +----------------------------------------------------------------*
 330 *  |  +----------------------------------------------------------------*
 331 *  |  |  Put here the check on the end point in magnetic field:
 332          IF ( LEMAGN .AND. NTRKCR .EQ. NTRACK .AND. MREG .NE. NEWREG )
 333      &      THEN
 334             NEWLSV = NEWLAT
 335             MLATSV = MLATTC
 336             NROLD  = NEWREG
 337             CALL GEOMAG ( XTRKCR, YTRKCR, ZTRKCR, UTRKCR, VTRKCR,
 338      &                    WTRKCR, NWREG , NROLD , IDSCCK )
 339             NEWLAT = NEWLSV
 340             MLATTC = MLATSV
 341             IF ( NWREG .EQ. NROLD ) GO TO 7000
 342 *  |  |-->-->--> already beyond the real crossing point:
 343          END IF
 344 *  |  |
 345 *  |  +----------------------------------------------------------------*
 346          LSTOPP = LSTOPP + 1
 347          NUMOPH = NUMOPH + 1
 348          IF ( NUMOPH .GT. 1000000000 ) THEN
 349             MUMOPH = MUMOPH + 1
 350             NUMOPH = NUMOPH + 1000000000
 351          END IF
 352          WOPTPH = WOPTPH + WTRACK
 353          POPTPH (LSTOPP) = EPHSMP
 354 *  |  +----------------------------------------------------------------*
 355 *  |  |  No easy safe way to reconstruct Dnear along a magnetic field
 356 *  |  |  step:
 357          IF ( LEMAGN ) THEN
 358             DONEAR (LSTOPP) = ZERZER
 359 *  |  |
 360 *  |  +----------------------------------------------------------------*
 361 *  |  |  Questionable: ddnear actually refers to the end track point
 362          ELSE
 363             DONEAR (LSTOPP) = MAX ( DDNEAR - CTRACK + TRUSED, ZERZER )
 364          END IF
 365 *  |  |
 366 *  |  +----------------------------------------------------------------*
 367          XOPTPH (LSTOPP) = XTRKCR
 368          YOPTPH (LSTOPP) = YTRKCR
 369          ZOPTPH (LSTOPP) = ZTRKCR
 370          NREGOP (LSTOPP) = MREG
 371          NLATOP (LSTOPP) = LATCCR
 372          COSTHE = ONEONE / BETNSM
 373          SINTHE = SQRT ( ( ONEONE - COSTHE ) * ( ONEONE + COSTHE ) )
 374          CALL SFECFE ( SINPHI, COSPHI )
 375          UPRIME = COSPHI * SINTHE
 376          VPRIME = SINPHI * SINTHE
 377          WPRIME = COSTHE
 378          IF ( LSMPAN ) THEN
 379             TXOPPH (LSTOPP) = UPRIME
 380             TYOPPH (LSTOPP) = VPRIME
 381             TZOPPH (LSTOPP) = WPRIME * COSTH0
 382          ELSE
 383             TXOPPH (LSTOPP) = UNDOXR ( UPRIME, VPRIME, WPRIME, SINPH0,
 384      &                                 COSPH0, SINTH0, COSTH0 )
 385             TYOPPH (LSTOPP) = UNDOYR ( UPRIME, VPRIME, WPRIME, SINPH0,
 386      &                                 COSPH0, SINTH0, COSTH0 )
 387             TZOPPH (LSTOPP) = UNDOZR ( UPRIME, VPRIME, WPRIME, SINPH0,
 388      &                                 COSPH0, SINTH0, COSTH0 )
 389          END IF
 390 *  |  Set-up the polarization vector: the following for linear pola-
 391 *  |  rization lying in the plane defined by the photon and the
 392 *  |  original particle direction:
 393          TXPOPP (LSTOPP) = UTRKCR / SINTHE - COSTHE / SINTHE
 394      &                   * TXOPPH (LSTOPP)
 395          TYPOPP (LSTOPP) = VTRKCR / SINTHE - COSTHE / SINTHE
 396      &                   * TYOPPH (LSTOPP)
 397          TZPOPP (LSTOPP) = WTRKCR / SINTHE - COSTHE / SINTHE
 398      &                   * TZOPPH (LSTOPP)
 399          SCADOT = TXOPPH (LSTOPP) * TXPOPP (LSTOPP)
 400      &          + TYOPPH (LSTOPP) * TYPOPP (LSTOPP)
 401      &          + TZOPPH (LSTOPP) * TZPOPP (LSTOPP)
 402          IF ( SCADOT .GT. CSNPRN ) THEN
 403             WRITE (LUNERR,'(A,1PG23.15,A)')
 404      &           ' *** Crnkvp: bad polarization',
 405      &              SCADOT, ' ***'
 406 *D           WRITE (77,'(A,1PG23.15)')
 407 *D    &            ' ^^^Crnkvp: bad polarization',
 408 *D    &              SCADOT
 409          END IF
 410          WTOPPH (LSTOPP) = WTRACK
 411          AGOPPH (LSTOPP) = ATRKCR
 412          CMPOPP (LSTOPP) = ZERZER
 413          LOOPPH (LSTOPP) = LTRACK + 1
 414 *
 415 *
 416 *        Hook to TFluka
 417 *
 418          PXCR =  EPHSMP * TXOPPH (LSTOPP)
 419          PYCR =  EPHSMP * TYOPPH (LSTOPP)
 420          PZCR =  EPHSMP * TZOPPH (LSTOPP)
 421          POX  = TXPOPP(LSTOPP)
 422          POY  = TYPOPP(LSTOPP)
 423          POZ  = TZPOPP(LSTOPP)
 424          CALL pshckp(PXCR, PYCR, PZCR, EPHSMP, XTRKCR,
 425      &        YTRKCR , ZTRKCR, ATRKCR, POX, POY, POZ, WTRACK, ITFL)
 426          NPROD = NPROD + 1
 427          CALL ustckv(NPROD, MREG, XTRKCR, YTRKCR, ZTRKCR)
 428 *
 429 *
 430 *
 431 *  |  !!!!!! Here Stuprf should be used !!!!!!
 432          LOUOPP (LSTOPP) = ITFL
 433          DO 2100 ISPR = 1, MKBMX1
 434             SPAROK (ISPR,LSTOPP) = SPAUSR (ISPR)
 435  2100    CONTINUE
 436          DO 2200 ISPR = 1, MKBMX2
 437             ISPORK (ISPR,LSTOPP) = ISPUSR (ISPR)
 438  2200    CONTINUE
 439 *  |  Save the parent features:
 440          TPROPP (LSTOPP) = ETRACK - AMPRTC
 441          APROPP (LSTOPP) = ATRKCR
 442          LPROPP (LSTOPP) = LTRACK
 443          IPROPP (LSTOPP) = JTRACK
 444          NPROPP (LSTOPP) = 0
 445 *  |  Update the "current" particle values:
 446          ETRKCR = ETRKCR - EPHSMP
 447          DEDXCK = DEDXCK + EPHSMP
 448          PTRKCR = SQRT ( ( ETRKCR - AMPRTC ) * ( ETRKCR + AMPRTC ) )
 449          PTINTR = PTRKCR
 450          ETINTR = ETRKCR
 451          TRINTR = ZERZER
 452          BETNLD = BETNCR
 453          BETNCR = PTRKCR / ETRKCR * RMXCER (MMAT)
 454 *  |  No longer any chance to emit Cerenkov photons:
 455          IF ( BETNCR .LE. ONEONE ) GO TO 7000
 456          BTNFLD = BTNFCR
 457          BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
 458 *  |  Update the macroscopic sigma: beta can only decrease
 459          SIGMCK = SIGMCK * BTNFCR / BTNFLD
 460       GO TO 1000
 461 *  |
 462 *  +-------------------------------------------------------------------*
 463  7000 CONTINUE
 464       RETURN
 465 *=== End of subroutine Crnkvp =========================================*
 466       END
 467