Minor changes to comply with the new version of AliAltroRawStream (C.Cheshkov)
[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 '(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) = LLOUSE
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