Add Cerenkov photons to TVirtualMCStack.
[u/mrichter/AliRoot.git] / 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