4 *=== Crnkvp ===========================================================*
6 SUBROUTINE CRNKVP ( MREG, NEWREG, DDNEAR, LEMAGN, DEDXCK )
12 *----------------------------------------------------------------------*
14 * Copyright (C) 1997-2003 by Alfredo Ferrari & Paola Sala *
15 * All Rights Reserved. *
18 * CeReNKoV Production: *
20 * Created on 21 september 1997 by Alfredo Ferrari & Paola Sala *
23 * Last change on 30-oct-98 by Alfredo Ferrari *
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 * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
31 * dN/dEdx = alpha^2 z^2 / [ r_e m_e c^2 ] [ 1 - 1 / (beta n)^2 ] *
33 * Note the wavelength is always assumed to be the vacuum one *
35 *----------------------------------------------------------------------*
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:
53 * No magnetic field end step Newreg check implemented for the moment:
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' )
63 IF ( JTRACK .GE. -6 ) THEN
66 AMPRTC = AMNHEA (-JTRACK)
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)
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
81 * +-------------------------------------------------------------------*
82 * | Loop on energy deposition sub-steps:
84 DTRCKT = DTRCKT + DTRACK (ID)
87 * +-------------------------------------------------------------------*
89 * +-------------------------------------------------------------------*
90 * | Loop on track sub-steps (be careful with magnetic field they are
91 * | not the straight distance between the end points):
93 TTRCKT = TTRCKT + TTRACK (IT)
96 * +-------------------------------------------------------------------*
97 CRVCRR = CTRACK / TTRCKT
98 IF ( JTRACK .GE. -6 ) THEN
99 CRNKCS = CRNKC0 * ( DBLE ( ICHRGE (JTRACK) ) )**2 * OPSNMX
101 CRNKCS = CRNKC0 * ( DBLE ( ICHEAV(-JTRACK) ) )**2 * OPSNMX
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
116 IF ( LSTOPP .GT. MOSTCK - NEMPHO ) THEN
117 LDUMMY = LBXRFL ( NEWREG, MREG, ZERZER, ZERZER, ONEONE )
118 CALL KASOPH ( .TRUE. )
121 * +-------------------------------------------------------------------*
122 * Total wiggled and curved track used up to now:
128 * Current lattice number and index inside Mulbou:
130 LATCCR = MULTTC (IRLTCR)
131 * Variables in the current Ntrckr bin:
136 TRNRSD = TTRACK (NTRKCR) * CRVCRR
137 * +-------------------------------------------------------------------*
140 DEDXCR = DTRACK (NTRKCR) / TRNRSD
142 * +-------------------------------------------------------------------*
145 DEDXCR = DTRCKT / CTRACK
148 * +-------------------------------------------------------------------*
149 * +-------------------------------------------------------------------*
150 * | Stacking loop for Cerenkov photons:
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. )
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 * | +----------------------------------------------------------------*
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
187 TRNRSD = TRNRSD - DSTPRD
188 ETRKCR = ETRKCR - DSTPRD * DEDXCR
191 * | | +-------------------------------------------------------------*
192 * | | | The production point is outside the current sub-step:
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:
200 * | | | Update the residual distance available in the current
202 TRNRSD = TTRACK (NTRKCR) * CRVCRR
203 * | | | Update the current dE/dx:
204 IF ( LMEQNT ) DEDXCR = DTRACK (NTRKCR) / TRNRSD
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):
217 * | | Dt = - | dE dx/dE E/(pc) = 1/c dx/dE [ Ptintr - Ptrkcr ]
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
225 * | +----------------------------------------------------------------*
228 ATRKCR = ATRKCR + TRINTR * ETRKCR / PTRKCR / CLIGHT
231 * | +----------------------------------------------------------------*
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
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) )
257 OMGSMP = GEVOMG * EPHSMP
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
293 IF ( .NOT. LSMPAN ) THEN
294 SINTH0 = SQRT (SINT02)
295 COSPH0 = UTRKCR / SINTH0
296 SINPH0 = VTRKCR / SINTH0
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)
317 DO 2000 IL = IRLTCR + 1, NULTTC
318 IF ( TRLATT .LE. TSLTTC (IL) ) THEN
320 LATCCR = MULTTC (IRLTCR)
324 *D CALL FLABRT ( 'CRNKVP', 'STOP:CRNKVP-NO-VALID-LATTICE' )
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 )
337 CALL GEOMAG ( XTRKCR, YTRKCR, ZTRKCR, UTRKCR, VTRKCR,
338 & WTRKCR, NWREG , NROLD , IDSCCK )
341 IF ( NWREG .EQ. NROLD ) GO TO 7000
342 * | |-->-->--> already beyond the real crossing point:
345 * | +----------------------------------------------------------------*
348 IF ( NUMOPH .GT. 1000000000 ) THEN
350 NUMOPH = NUMOPH + 1000000000
352 WOPTPH = WOPTPH + WTRACK
353 POPTPH (LSTOPP) = EPHSMP
354 * | +----------------------------------------------------------------*
355 * | | No easy safe way to reconstruct Dnear along a magnetic field
358 DONEAR (LSTOPP) = ZERZER
360 * | +----------------------------------------------------------------*
361 * | | Questionable: ddnear actually refers to the end track point
363 DONEAR (LSTOPP) = MAX ( DDNEAR - CTRACK + TRUSED, ZERZER )
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
379 TXOPPH (LSTOPP) = UPRIME
380 TYOPPH (LSTOPP) = VPRIME
381 TZOPPH (LSTOPP) = WPRIME * COSTH0
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 )
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
395 TYPOPP (LSTOPP) = VTRKCR / SINTHE - COSTHE / SINTHE
397 TZPOPP (LSTOPP) = WTRKCR / SINTHE - COSTHE / SINTHE
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',
406 *D WRITE (77,'(A,1PG23.15)')
407 *D & ' ^^^Crnkvp: bad polarization',
410 WTOPPH (LSTOPP) = WTRACK
411 AGOPPH (LSTOPP) = ATRKCR
412 CMPOPP (LSTOPP) = ZERZER
413 LOOPPH (LSTOPP) = LTRACK + 1
418 PXCR = EPHSMP * TXOPPH (LSTOPP)
419 PYCR = EPHSMP * TYOPPH (LSTOPP)
420 PZCR = EPHSMP * TZOPPH (LSTOPP)
424 CALL pshckp(PXCR, PYCR, PZCR, EPHSMP, XTRKCR,
425 & YTRKCR , ZTRKCR, ATRKCR, POX, POY, POZ, WTRACK, ITFL)
427 CALL ustckv(NPROD, MREG, XTRKCR, YTRKCR, ZTRKCR)
431 * | !!!!!! Here Stuprf should be used !!!!!!
432 LOUOPP (LSTOPP) = ITFL
433 DO 2100 ISPR = 1, MKBMX1
434 SPAROK (ISPR,LSTOPP) = SPAUSR (ISPR)
436 DO 2200 ISPR = 1, MKBMX2
437 ISPORK (ISPR,LSTOPP) = ISPUSR (ISPR)
439 * | Save the parent features:
440 TPROPP (LSTOPP) = ETRACK - AMPRTC
441 APROPP (LSTOPP) = ATRKCR
442 LPROPP (LSTOPP) = LTRACK
443 IPROPP (LSTOPP) = JTRACK
445 * | Update the "current" particle values:
446 ETRKCR = ETRKCR - EPHSMP
447 DEDXCK = DEDXCK + EPHSMP
448 PTRKCR = SQRT ( ( ETRKCR - AMPRTC ) * ( ETRKCR + AMPRTC ) )
453 BETNCR = PTRKCR / ETRKCR * RMXCER (MMAT)
454 * | No longer any chance to emit Cerenkov photons:
455 IF ( BETNCR .LE. ONEONE ) GO TO 7000
457 BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
458 * | Update the macroscopic sigma: beta can only decrease
459 SIGMCK = SIGMCK * BTNFCR / BTNFLD
462 * +-------------------------------------------------------------------*
465 *=== End of subroutine Crnkvp =========================================*