]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TFluka/crnkvp.f
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Add Cerenkov photons to TVirtualMCStack.
[u/mrichter/AliRoot.git] / TFluka / crnkvp.f
... / ...
CommitLineData
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 )
111D IF ( LSTOPP .GT. MOSTCK - NEMPHO ) WRITE (77,*)
112D & ' ###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
152D IF ( SIGMCK .LT. ZERZER ) WRITE (77,*)
153D & ' ^^^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:
158D IF ( LSTOPP .EQ. MOSTCK ) WRITE (77,*)
159D & ' ###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
177D IF ( NTRKCR .GT. NTRACK ) WRITE (77,*)
178D & ' ^^^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
243D IF ( FREJE .GT. ONEPLS ) WRITE (77,*) ' ^^^CRNKVP:',
244D & 'FREJE,BETNCR,BETNLD,DEDXCR,PTRKCR,ETRKCR,PTRACK,ETRACK',
245D & 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
270D IF ( FREJE .GT. ONEPLS ) WRITE (77,*)
271D & ' ^^^CRNKVP:FREJE,RFISMP,RMXCER(MMAT),OPSENS,MMAT',
272D & 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)
302D IF ( SRNRSD .LT. ZERZER .OR. SRNRSD .GT. TUVWCR )
303D & WRITE (77,*)' ^^^CRNKVP:SRNRSD,TUVWCR,NTRKCR',
304D & 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
323D 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, ' ***'
405D WRITE (77,'(A,1PG23.15)')
406D & ' ^^^Crnkvp: bad polarization',
407D & 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