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 '(MAPA)'
  39       INCLUDE '(MULBOU)'
  40       INCLUDE '(OPPHCM)'
  41       INCLUDE '(OPPHST)'
  42       INCLUDE '(PAPROP)'
  43       INCLUDE '(STARS)'
  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  1000 CONTINUE
 152 D        IF ( SIGMCK .LT. ZERZER ) WRITE (77,*)
 153 D    &      ' ^^^CRNKVP:SIGMCK,BTNFCR',SIGMCK,BTNFCR
 154 *D        IF ( DEDXCR .LT. AZRZRZ ) WRITE (77,*)' ###CRNKVP:',
 155 *D    &  'MMAT,JTRACK,DEDXCR,NTRKCR,DTRACK(NTRKCR),CTRACK,TRNRSD,DTRCKT',
 156 *D    &   MMAT,JTRACK,DEDXCR,NTRKCR,DTRACK(NTRKCR),CTRACK,TRNRSD,DTRCKT
 157 *  |  Empty the stack if too full:
 158 D        IF ( LSTOPP .EQ. MOSTCK ) WRITE (77,*)
 159 D    &      ' ###CRNKVP:LSTOPP:',LSTOPP
 160 *  |  +----------------------------------------------------------------*
 161 *  |  |  Empty the stack if too full: compute and save the normal in
 162 *  |  |  case it is required
 163          IF ( LSTOPP .EQ. MOSTCK ) THEN
 164             LDUMMY = LBXRFL ( NEWREG, MREG, ZERZER, ZERZER, ONEONE )
 165             CALL KASOPH ( .TRUE. )
 166          END IF
 167 *  |  |
 168 *  |  +----------------------------------------------------------------*
 169          DSTPRD = - LOG ( ONEONE - FLRNDM (DSTPRD) ) / SIGMCK
 170          TRUSED = TRUSED + DSTPRD
 171 *dbgD        WRITE (77,*) ' CRNKVP:DSTPRD,TRUSED,CTRACK,TRNRSD',
 172 *dbgD    &                         DSTPRD,TRUSED,CTRACK,TRNRSD
 173          IF ( TRUSED .GE. CTRACK ) GO TO 7000
 174 *  |  +----------------------------------------------------------------*
 175 *  |  |  Sub-step loop:
 176  1500    CONTINUE
 177 D           IF ( NTRKCR .GT. NTRACK ) WRITE (77,*)
 178 D    &         ' ^^^CRNKVP:NTRKCR,NTRACK',NTRKCR,NTRACK
 179 *  |  |  +-------------------------------------------------------------*
 180 *  |  |  |  The production point is inside the current sub-step:
 181             IF ( DSTPRD .LE. TRNRSD ) THEN
 182 *  |  |  |  Update the distance from the previous interaction:
 183                TRINTR = TRINTR + DSTPRD
 184 *  |  |  |  Update the residual distance available in the current
 185 *  |  |  |  sub-step
 186                TRNRSD = TRNRSD - DSTPRD
 187                ETRKCR = ETRKCR - DSTPRD * DEDXCR
 188                LNWSUB = .FALSE.
 189 *  |  |  |
 190 *  |  |  +-------------------------------------------------------------*
 191 *  |  |  |  The production point is outside the current sub-step:
 192             ELSE
 193 *  |  |  |  Update the distance from the previous interaction:
 194                TRINTR = TRINTR + TRNRSD
 195                ETRKCR = ETRKCR - TRNRSD * DEDXCR
 196                DSTPRD = DSTPRD - TRNRSD
 197 *  |  |  |  Change sub-step:
 198                NTRKCR = NTRKCR + 1
 199 *  |  |  |  Update the residual distance available in the current
 200 *  |  |  |  sub-step
 201                TRNRSD = TTRACK (NTRKCR) * CRVCRR
 202 *  |  |  |  Update the current dE/dx:
 203                IF ( LMEQNT ) DEDXCR = DTRACK (NTRKCR) / TRNRSD
 204                LNWSUB = .TRUE.
 205             END IF
 206 *  |  |  |
 207 *  |  |  +-------------------------------------------------------------*
 208          IF ( LNWSUB ) GO TO 1500
 209 *  |  |  end sub-step loop:
 210 *  |  +----------------------------------------------------------------*
 211          PTRKCR = SQRT ( ( ETRKCR - AMPRTC ) * ( ETRKCR + AMPRTC ) )
 212 *  |  +----------------------------------------------------------------*
 213 *  |  |  Update the particle "age" (beta=p/E, dp/dE=E/p):
 214 *  |  |           Etrkcr
 215 *  |  |         /
 216 *  |  |  Dt = - | dE dx/dE E/(pc) = 1/c dx/dE [ Ptintr - Ptrkcr ]
 217 *  |  |         /
 218 *  |  |          Etintr
 219          IF ( ETINTR - ETRKCR .GT. CSNNRM * ETRKCR ) THEN
 220 *  |  |  Average dE/dx from the previous (fictitious) interaction:
 221             DEDXAV = ( ETINTR - ETRKCR ) / TRINTR
 222             ATRKCR = ATRKCR + ( PTINTR - PTRKCR ) / DEDXAV / CLIGHT
 223 *  |  |
 224 *  |  +----------------------------------------------------------------*
 225 *  |  |
 226          ELSE
 227             ATRKCR = ATRKCR + TRINTR * ETRKCR / PTRKCR / CLIGHT
 228          END IF
 229 *  |  |
 230 *  |  +----------------------------------------------------------------*
 231          PTINTR = PTRKCR
 232          ETINTR = ETRKCR
 233          TRINTR = ZERZER
 234          BETNLD = BETNCR
 235          BETNCR = PTRKCR / ETRKCR * RMXCER (MMAT)
 236 *dbgD        WRITE (77,*) ' CRNKVP:BETNCR,RMXCER,LSTOPP',
 237 *dbgD    &                         BETNCR,RMXCER(MMAT),LSTOPP
 238 *  |  No longer any chance to emit Cerenkov photons:
 239          IF ( BETNCR .LE. ONEONE ) GO TO 7000
 240          BTNFLD = BTNFCR
 241          BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
 242          FREJE  = BTNFCR / BTNFLD
 243 D        IF ( FREJE .GT. ONEPLS ) WRITE (77,*) ' ^^^CRNKVP:',
 244 D    &       'FREJE,BETNCR,BETNLD,DEDXCR,PTRKCR,ETRKCR,PTRACK,ETRACK',
 245 D    &        FREJE,BETNCR,BETNLD,DEDXCR,PTRKCR,ETRKCR,PTRACK,ETRACK
 246 *  |  Update the macroscopic sigma: beta can only decrease
 247          SIGMCK = SIGMCK * FREJE
 248          RNDREJ = FLRNDM (RNDREJ)
 249          IF ( RNDREJ .GE. FREJE ) GO TO 1000
 250 *  |--<--<--<--<--< Interaction rejected because of the decrease in beta
 251          RNDREJ = RNDREJ / FREJE
 252 *  |  Sample the emitted photon energy:
 253          EPHSMP = RNDREJ * ( EMXCER (MMAT) - EMNCER (MMAT) )
 254      &          + EMNCER (MMAT)
 255 *  |  2pi x freq.:
 256          OMGSMP = GEVOMG * EPHSMP
 257 *  |  Wavelength:
 258          WVLSMP = TWOPIP * CLIGHT / OMGSMP
 259 *  |  Compute the refraction index at the right wavelength/frequency:
 260          RFISMP = FOPTPR ( 1, WVLSMP, OMGSMP, MMAT )
 261          BETNSM = PTRKCR / ETRKCR * RFISMP
 262 *dbgD        WRITE (77,*) ' CRNKVP:BETNSM,RFISMP,EPHSMP,WVLSMP,OMGSMP',
 263 *dbgD    &                         BETNSM,RFISMP,EPHSMP,WVLSMP,OMGSMP
 264          IF ( BETNSM .LE. ONEONE ) GO TO 1000
 265 *  |--<--<--<--<--< No chance to emit Cerenkov photons at this frequency
 266          BTNFSM = ( BETNSM - ONEONE ) * ( BETNSM + ONEONE ) / BETNSM**2
 267 *  |  Compute the quantum efficiency for this emitted energy:
 268          OPSENS = FOPTSN ( WVLSMP, OMGSMP )
 269          FREJE  = BTNFSM / BTNFCR * OPSENS / OPSNMX
 270 D        IF ( FREJE .GT. ONEPLS ) WRITE (77,*)
 271 D    &      ' ^^^CRNKVP:FREJE,RFISMP,RMXCER(MMAT),OPSENS,MMAT',
 272 D    &                  FREJE,RFISMP,RMXCER(MMAT),OPSENS,MMAT
 273          RNDREJ = FLRNDM (RNDREJ)
 274 *dbgD        WRITE (77,*) ' CRNKVP:BTNFSM,BTNFCR,RNDREJ',
 275 *dbgD    &                         BTNFSM,BTNFCR,RNDREJ
 276          IF ( RNDREJ .GE. FREJE ) GO TO 1000
 277 *  |--<--<--<--<--< Interaction rejected because of the decrease in n
 278 *  |  Eventually a Cerenkov photon is going to be emitted:
 279 *  |  +----------------------------------------------------------------*
 280 *  |  |  New production sub-step:
 281          IF ( NTRKCR .NE. NTRKLD ) THEN
 282             UTRKCR = XTRACK (NTRKCR) - XTRACK (NTRKCR-1)
 283             VTRKCR = YTRACK (NTRKCR) - YTRACK (NTRKCR-1)
 284             WTRKCR = ZTRACK (NTRKCR) - ZTRACK (NTRKCR-1)
 285             TUVWCR = SQRT ( UTRKCR**2 + VTRKCR**2 + WTRKCR**2 )
 286             UTRKCR = UTRKCR / TUVWCR
 287             VTRKCR = VTRKCR / TUVWCR
 288             WTRKCR = WTRKCR / TUVWCR
 289             SINT02 = UTRKCR**2 + VTRKCR**2
 290             LSMPAN = SINT02 .LT. ANGLSQ
 291             COSTH0 = WTRKCR
 292             IF ( .NOT. LSMPAN ) THEN
 293                SINTH0 = SQRT (SINT02)
 294                COSPH0 = UTRKCR / SINTH0
 295                SINPH0 = VTRKCR / SINTH0
 296             END IF
 297             NTRKLD = NTRKCR
 298          END IF
 299 *  |  |
 300 *  |  +----------------------------------------------------------------*
 301          SRNRSD = TUVWCR * TRNRSD / CRVCRR / TTRACK (NTRKCR)
 302 D        IF ( SRNRSD .LT. ZERZER .OR. SRNRSD .GT. TUVWCR )
 303 D    &      WRITE (77,*)' ^^^CRNKVP:SRNRSD,TUVWCR,NTRKCR',
 304 D    &                              SRNRSD,TUVWCR,NTRKCR
 305          XTRKCR = XTRACK (NTRKCR) - UTRKCR * SRNRSD
 306          YTRKCR = YTRACK (NTRKCR) - VTRKCR * SRNRSD
 307          ZTRKCR = ZTRACK (NTRKCR) - WTRKCR * SRNRSD
 308 *  |  +----------------------------------------------------------------*
 309 *  |  |  Put here the check on lattice change:
 310          IF ( NULTTC .GT. 0 ) THEN
 311 *  |  |  Current (possibly curved in mag. field) non wiggled path:
 312             TRLATT = TRUSED / CRVCRR
 313             IF ( TRLATT .LE. TSLTTC (IRLTCR) ) THEN
 314                LATCCR = MULTTC (IRLTCR)
 315             ELSE
 316                DO 2000 IL = IRLTCR + 1, NULTTC
 317                   IF ( TRLATT .LE. TSLTTC (IL) ) THEN
 318                      IRLTCR = IL
 319                      LATCCR = MULTTC (IRLTCR)
 320                      GO TO 2010
 321                   END IF
 322  2000          CONTINUE
 323 D              CALL FLABRT ( 'CRNKVP', 'STOP:CRNKVP-NO-VALID-LATTICE' )
 324  2010          CONTINUE
 325             END IF
 326          END IF
 327 *  |  |
 328 *  |  +----------------------------------------------------------------*
 329 *  |  +----------------------------------------------------------------*
 330 *  |  |  Put here the check on the end point in magnetic field:
 331          IF ( LEMAGN .AND. NTRKCR .EQ. NTRACK .AND. MREG .NE. NEWREG )
 332      &      THEN
 333             NEWLSV = NEWLAT
 334             MLATSV = MLATTC
 335             NROLD  = NEWREG
 336             CALL GEOMAG ( XTRKCR, YTRKCR, ZTRKCR, UTRKCR, VTRKCR,
 337      &                    WTRKCR, NWREG , NROLD , IDSCCK )
 338             NEWLAT = NEWLSV
 339             MLATTC = MLATSV
 340             IF ( NWREG .EQ. NROLD ) GO TO 7000
 341 *  |  |-->-->--> already beyond the real crossing point:
 342          END IF
 343 *  |  |
 344 *  |  +----------------------------------------------------------------*
 345          LSTOPP = LSTOPP + 1
 346          NUMOPH = NUMOPH + 1
 347          IF ( NUMOPH .GT. 1000000000 ) THEN
 348             MUMOPH = MUMOPH + 1
 349             NUMOPH = NUMOPH + 1000000000
 350          END IF
 351          WOPTPH = WOPTPH + WTRACK
 352          POPTPH (LSTOPP) = EPHSMP
 353 *  |  +----------------------------------------------------------------*
 354 *  |  |  No easy safe way to reconstruct Dnear along a magnetic field
 355 *  |  |  step:
 356          IF ( LEMAGN ) THEN
 357             DONEAR (LSTOPP) = ZERZER
 358 *  |  |
 359 *  |  +----------------------------------------------------------------*
 360 *  |  |  Questionable: ddnear actually refers to the end track point
 361          ELSE
 362             DONEAR (LSTOPP) = MAX ( DDNEAR - CTRACK + TRUSED, ZERZER )
 363          END IF
 364 *  |  |
 365 *  |  +----------------------------------------------------------------*
 366          XOPTPH (LSTOPP) = XTRKCR
 367          YOPTPH (LSTOPP) = YTRKCR
 368          ZOPTPH (LSTOPP) = ZTRKCR
 369          NREGOP (LSTOPP) = MREG
 370          NLATOP (LSTOPP) = LATCCR
 371          COSTHE = ONEONE / BETNSM
 372          SINTHE = SQRT ( ( ONEONE - COSTHE ) * ( ONEONE + COSTHE ) )
 373          CALL SFECFE ( SINPHI, COSPHI )
 374          UPRIME = COSPHI * SINTHE
 375          VPRIME = SINPHI * SINTHE
 376          WPRIME = COSTHE
 377          IF ( LSMPAN ) THEN
 378             TXOPPH (LSTOPP) = UPRIME
 379             TYOPPH (LSTOPP) = VPRIME
 380             TZOPPH (LSTOPP) = WPRIME * COSTH0
 381          ELSE
 382             TXOPPH (LSTOPP) = UNDOXR ( UPRIME, VPRIME, WPRIME, SINPH0,
 383      &                                 COSPH0, SINTH0, COSTH0 )
 384             TYOPPH (LSTOPP) = UNDOYR ( UPRIME, VPRIME, WPRIME, SINPH0,
 385      &                                 COSPH0, SINTH0, COSTH0 )
 386             TZOPPH (LSTOPP) = UNDOZR ( UPRIME, VPRIME, WPRIME, SINPH0,
 387      &                                 COSPH0, SINTH0, COSTH0 )
 388          END IF
 389 *  |  Set-up the polarization vector: the following for linear pola-
 390 *  |  rization lying in the plane defined by the photon and the
 391 *  |  original particle direction:
 392          TXPOPP (LSTOPP) = UTRKCR / SINTHE - COSTHE / SINTHE
 393      &                   * TXOPPH (LSTOPP)
 394          TYPOPP (LSTOPP) = VTRKCR / SINTHE - COSTHE / SINTHE
 395      &                   * TYOPPH (LSTOPP)
 396          TZPOPP (LSTOPP) = WTRKCR / SINTHE - COSTHE / SINTHE
 397      &                   * TZOPPH (LSTOPP)
 398          SCADOT = TXOPPH (LSTOPP) * TXPOPP (LSTOPP)
 399      &          + TYOPPH (LSTOPP) * TYPOPP (LSTOPP)
 400      &          + TZOPPH (LSTOPP) * TZPOPP (LSTOPP)
 401          IF ( SCADOT .GT. CSNPRN ) THEN
 402             WRITE (LUNERR,'(A,1PG23.15,A)')
 403      &           ' *** Crnkvp: bad polarization',
 404      &              SCADOT, ' ***'
 405 D           WRITE (77,'(A,1PG23.15)')
 406 D    &            ' ^^^Crnkvp: bad polarization',
 407 D    &              SCADOT
 408          END IF
 409          WTOPPH (LSTOPP) = WTRACK
 410          AGOPPH (LSTOPP) = ATRKCR
 411          CMPOPP (LSTOPP) = ZERZER
 412          LOOPPH (LSTOPP) = LTRACK + 1
 413 *
 414 *
 415 *        Hook to TFluka
 416 *
 417          PXCR =  EPHSMP * TXOPPH (LSTOPP)
 418          PYCR =  EPHSMP * TYOPPH (LSTOPP)
 419          PZCR =  EPHSMP * TZOPPH (LSTOPP)
 420          POX  = TXPOPP(LSTOPP)
 421          POY  = TYPOPP(LSTOPP)
 422          POZ  = TZPOPP(LSTOPP)
 423          CALL PushCerenkovPhoton(PXCR, PYCR, PZCR, EPHSMP, XTRKCR,
 424      &        YTRKCR , ZTRKCR, ATRKCR, POX, POY, POZ, WTRACK, ITFL)
 425 *
 426 *
 427 *
 428 *  |  !!!!!! Here Stuprf should be used !!!!!!
 429          LOUOPP (LSTOPP) = ITFL
 430          DO 2100 ISPR = 1, MKBMX1
 431             SPAROK (ISPR,LSTOPP) = SPAUSR (ISPR)
 432  2100    CONTINUE
 433          DO 2200 ISPR = 1, MKBMX2
 434             ISPORK (ISPR,LSTOPP) = ISPUSR (ISPR)
 435  2200    CONTINUE
 436 *  |  Save the parent features:
 437          TPROPP (LSTOPP) = ETRACK - AMPRTC
 438          APROPP (LSTOPP) = ATRKCR
 439          LPROPP (LSTOPP) = LTRACK
 440          IPROPP (LSTOPP) = JTRACK
 441          NPROPP (LSTOPP) = 0
 442 *  |  Update the "current" particle values:
 443          ETRKCR = ETRKCR - EPHSMP
 444          DEDXCK = DEDXCK + EPHSMP
 445          PTRKCR = SQRT ( ( ETRKCR - AMPRTC ) * ( ETRKCR + AMPRTC ) )
 446          PTINTR = PTRKCR
 447          ETINTR = ETRKCR
 448          TRINTR = ZERZER
 449          BETNLD = BETNCR
 450          BETNCR = PTRKCR / ETRKCR * RMXCER (MMAT)
 451 *  |  No longer any chance to emit Cerenkov photons:
 452          IF ( BETNCR .LE. ONEONE ) GO TO 7000
 453          BTNFLD = BTNFCR
 454          BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
 455 *  |  Update the macroscopic sigma: beta can only decrease
 456          SIGMCK = SIGMCK * BTNFCR / BTNFLD
 457       GO TO 1000
 458 *  |
 459 *  +-------------------------------------------------------------------*
 460  7000 CONTINUE
 461       RETURN
 462 *=== End of subroutine Crnkvp =========================================*
 463       END
 464