+++ /dev/null
-*$ CREATE CRNKVP.FOR
-*COPY CRNKVP
-*
-*=== Crnkvp ===========================================================*
-*
- SUBROUTINE CRNKVP ( MREG, NEWREG, DDNEAR, LEMAGN, DEDXCK )
-
- INCLUDE '(DBLPRC)'
- INCLUDE '(DIMPAR)'
- INCLUDE '(IOUNIT)'
-*
-*----------------------------------------------------------------------*
-* *
-* Copyright (C) 1997-2003 by Alfredo Ferrari & Paola Sala *
-* All Rights Reserved. *
-* *
-* *
-* CeReNKoV Production: *
-* *
-* Created on 21 september 1997 by Alfredo Ferrari & Paola Sala *
-* Infn - Milan *
-* *
-* Last change on 30-oct-98 by Alfredo Ferrari *
-* *
-* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
-* !!!!!! Important: Etrack, Ptrack, Atrack, are supposed !!!!!! *
-* !!!!!! to be the INITIAL values (not yet updated at the !!!!!! *
-* !!!!!! end of the step). !!!!!! *
-* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
-* *
-* dN/dEdx = alpha^2 z^2 / [ r_e m_e c^2 ] [ 1 - 1 / (beta n)^2 ] *
-* *
-* Note the wavelength is always assumed to be the vacuum one *
-* *
-*----------------------------------------------------------------------*
-*
- INCLUDE '(FHEAVY)'
- INCLUDE '(FLKMAT)'
- INCLUDE '(MULBOU)'
- INCLUDE '(OPPHCM)'
- INCLUDE '(OPPHST)'
- INCLUDE '(PAPROP)'
- INCLUDE '(SUMCOU)'
- INCLUDE '(TRACKR)'
-*
- PARAMETER ( CSNPRN = 100.D+00 * CSNNRM )
- PARAMETER ( CRNKC0 = ALPFSC * ALPFSC / RCLSEL / AMELCT )
- LOGICAL LEMAGN, LMEQNT, LNWSUB, LSMPAN, LDUMMY, LBXRFL
-* Statement functions:
-* INCLUDE '(DORTSF)'
- INCLUDE '(UNRTSF)'
-*
-* No magnetic field end step Newreg check implemented for the moment:
-*
-*
-* IF ( LEMAGN )
-* & CALL FLABRT ( 'CRNKVP', ' STOP:LEMAGN-NOT-YET-IMPLEMENTED' )
-* No change of lattice check implemented for the moment:
-* IF ( LT1TRK .NE. LT2TRK ) CALL FLABRT ( 'CRNKVP',
-* & ' STOP:LT1TRK.NE.LT2TRK-NOT-YET-IMPLEMENTED' )
- DEDXCK = ZERZER
- MMAT = MEDIUM (MREG)
- IF ( JTRACK .GE. -6 ) THEN
- AMPRTC = AM (JTRACK)
- ELSE
- AMPRTC = AMNHEA (-JTRACK)
- END IF
-* +-------------------------------------------------------------------*
-* | Find the maximum refractive index over the interval allowed
-* | for Cerenkov production:
- IF ( RMXCER (MMAT) .LT. ZERZER ) THEN
- RMXCER (MMAT) = RFNDMX (MMAT)
- END IF
-* |
-* +-------------------------------------------------------------------*
- BETNCR = PTRACK / ETRACK * RMXCER (MMAT)
-* Maximum beta x n already below 1: no chance to emit any photon
- IF ( BETNCR .LE. ONEONE ) RETURN
- LMEQNT = MTRACK .EQ. NTRACK
- DTRCKT = ZERZER
-* +-------------------------------------------------------------------*
-* | Loop on energy deposition sub-steps:
- DO 300 ID = 1, MTRACK
- DTRCKT = DTRCKT + DTRACK (ID)
- 300 CONTINUE
-* |
-* +-------------------------------------------------------------------*
- TTRCKT = ZERZER
-* +-------------------------------------------------------------------*
-* | Loop on track sub-steps (be careful with magnetic field they are
-* | not the straight distance between the end points):
- DO 500 IT = 1, NTRACK
- TTRCKT = TTRCKT + TTRACK (IT)
- 500 CONTINUE
-* |
-* +-------------------------------------------------------------------*
- CRVCRR = CTRACK / TTRCKT
- IF ( JTRACK .GE. -6 ) THEN
- CRNKCS = CRNKC0 * ( DBLE ( ICHRGE (JTRACK) ) )**2 * OPSNMX
- ELSE
- CRNKCS = CRNKC0 * ( DBLE ( ICHEAV(-JTRACK) ) )**2 * OPSNMX
- END IF
- BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
-* Macroscopic sigma for production (cm^-1)
- SIGMCK = CRNKCS * ( EMXCER (MMAT) - EMNCER (MMAT) ) * BTNFCR
-* Average number of emitted photons:
- NEMPHO = NINT ( CTRACK * SIGMCK )
-* Take a three sigma clearance:
- NEMPHO = NEMPHO + 10 + 3 * NINT ( SQRT ( CTRACK * SIGMCK ) )
- NEMPHO = MIN ( NEMPHO, MOSTCK )
-*D IF ( LSTOPP .GT. MOSTCK - NEMPHO ) WRITE (77,*)
-*D & ' ###CRNKVP:LSTOPP:',LSTOPP,NEMPHO
-* +-------------------------------------------------------------------*
-* | Empty the stack if too full: compute and save the normal in case
-* | it is required
- IF ( LSTOPP .GT. MOSTCK - NEMPHO ) THEN
- LDUMMY = LBXRFL ( NEWREG, MREG, ZERZER, ZERZER, ONEONE )
- CALL KASOPH ( .TRUE. )
- END IF
-* |
-* +-------------------------------------------------------------------*
-* Total wiggled and curved track used up to now:
- TRUSED = ZERZER
- ETINTR = ETRACK
- PTINTR = PTRACK
- TRINTR = ZERZER
- NTRKLD = 0
-* Current lattice number and index inside Mulbou:
- IRLTCR = 0
- LATCCR = MULTTC (IRLTCR)
-* Variables in the current Ntrckr bin:
- NTRKCR = 1
- PTRKCR = PTRACK
- ETRKCR = ETRACK
- ATRKCR = ATRACK
- TRNRSD = TTRACK (NTRKCR) * CRVCRR
-* +-------------------------------------------------------------------*
-* | Ntrack = Mtrack
- IF ( LMEQNT ) THEN
- DEDXCR = DTRACK (NTRKCR) / TRNRSD
-* |
-* +-------------------------------------------------------------------*
-* | Ntrack >< Mtrack
- ELSE
- DEDXCR = DTRCKT / CTRACK
- END IF
-* |
-* +-------------------------------------------------------------------*
-* +-------------------------------------------------------------------*
-* | Stacking loop for Cerenkov photons:
- NPROD = 0
- 1000 CONTINUE
-*D IF ( SIGMCK .LT. ZERZER ) WRITE (77,*)
-*D & ' ^^^CRNKVP:SIGMCK,BTNFCR',SIGMCK,BTNFCR
-*D IF ( DEDXCR .LT. AZRZRZ ) WRITE (77,*)' ###CRNKVP:',
-*D & 'MMAT,JTRACK,DEDXCR,NTRKCR,DTRACK(NTRKCR),CTRACK,TRNRSD,DTRCKT',
-*D & MMAT,JTRACK,DEDXCR,NTRKCR,DTRACK(NTRKCR),CTRACK,TRNRSD,DTRCKT
-* | Empty the stack if too full:
-*D IF ( LSTOPP .EQ. MOSTCK ) WRITE (77,*)
-*D & ' ###CRNKVP:LSTOPP:',LSTOPP
-* | +----------------------------------------------------------------*
-* | | Empty the stack if too full: compute and save the normal in
-* | | case it is required
- IF ( LSTOPP .EQ. MOSTCK ) THEN
- LDUMMY = LBXRFL ( NEWREG, MREG, ZERZER, ZERZER, ONEONE )
- CALL KASOPH ( .TRUE. )
- END IF
-* | |
-* | +----------------------------------------------------------------*
- DSTPRD = - LOG ( ONEONE - FLRNDM (DSTPRD) ) / SIGMCK
- TRUSED = TRUSED + DSTPRD
-*dbgD WRITE (77,*) ' CRNKVP:DSTPRD,TRUSED,CTRACK,TRNRSD',
-*dbgD & DSTPRD,TRUSED,CTRACK,TRNRSD
- IF ( TRUSED .GE. CTRACK ) GO TO 7000
-* | +----------------------------------------------------------------*
-* | | Sub-step loop:
- 1500 CONTINUE
-*D IF ( NTRKCR .GT. NTRACK ) WRITE (77,*)
-*D & ' ^^^CRNKVP:NTRKCR,NTRACK',NTRKCR,NTRACK
-* | | +-------------------------------------------------------------*
-* | | | The production point is inside the current sub-step:
- IF ( DSTPRD .LE. TRNRSD ) THEN
-* | | | Update the distance from the previous interaction:
- TRINTR = TRINTR + DSTPRD
-* | | | Update the residual distance available in the current
-* | | | sub-step
- TRNRSD = TRNRSD - DSTPRD
- ETRKCR = ETRKCR - DSTPRD * DEDXCR
- LNWSUB = .FALSE.
-* | | |
-* | | +-------------------------------------------------------------*
-* | | | The production point is outside the current sub-step:
- ELSE
-* | | | Update the distance from the previous interaction:
- TRINTR = TRINTR + TRNRSD
- ETRKCR = ETRKCR - TRNRSD * DEDXCR
- DSTPRD = DSTPRD - TRNRSD
-* | | | Change sub-step:
- NTRKCR = NTRKCR + 1
-* | | | Update the residual distance available in the current
-* | | | sub-step
- TRNRSD = TTRACK (NTRKCR) * CRVCRR
-* | | | Update the current dE/dx:
- IF ( LMEQNT ) DEDXCR = DTRACK (NTRKCR) / TRNRSD
- LNWSUB = .TRUE.
- END IF
-* | | |
-* | | +-------------------------------------------------------------*
- IF ( LNWSUB ) GO TO 1500
-* | | end sub-step loop:
-* | +----------------------------------------------------------------*
- PTRKCR = SQRT ( ( ETRKCR - AMPRTC ) * ( ETRKCR + AMPRTC ) )
-* | +----------------------------------------------------------------*
-* | | Update the particle "age" (beta=p/E, dp/dE=E/p):
-* | | Etrkcr
-* | | /
-* | | Dt = - | dE dx/dE E/(pc) = 1/c dx/dE [ Ptintr - Ptrkcr ]
-* | | /
-* | | Etintr
- IF ( ETINTR - ETRKCR .GT. CSNNRM * ETRKCR ) THEN
-* | | Average dE/dx from the previous (fictitious) interaction:
- DEDXAV = ( ETINTR - ETRKCR ) / TRINTR
- ATRKCR = ATRKCR + ( PTINTR - PTRKCR ) / DEDXAV / CLIGHT
-* | |
-* | +----------------------------------------------------------------*
-* | |
- ELSE
- ATRKCR = ATRKCR + TRINTR * ETRKCR / PTRKCR / CLIGHT
- END IF
-* | |
-* | +----------------------------------------------------------------*
- PTINTR = PTRKCR
- ETINTR = ETRKCR
- TRINTR = ZERZER
- BETNLD = BETNCR
- BETNCR = PTRKCR / ETRKCR * RMXCER (MMAT)
-*dbgD WRITE (77,*) ' CRNKVP:BETNCR,RMXCER,LSTOPP',
-*dbgD & BETNCR,RMXCER(MMAT),LSTOPP
-* | No longer any chance to emit Cerenkov photons:
- IF ( BETNCR .LE. ONEONE ) GO TO 7000
- BTNFLD = BTNFCR
- BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
- FREJE = BTNFCR / BTNFLD
-*D IF ( FREJE .GT. ONEPLS ) WRITE (77,*) ' ^^^CRNKVP:',
-*D & 'FREJE,BETNCR,BETNLD,DEDXCR,PTRKCR,ETRKCR,PTRACK,ETRACK',
-*D & FREJE,BETNCR,BETNLD,DEDXCR,PTRKCR,ETRKCR,PTRACK,ETRACK
-* | Update the macroscopic sigma: beta can only decrease
- SIGMCK = SIGMCK * FREJE
- RNDREJ = FLRNDM (RNDREJ)
- IF ( RNDREJ .GE. FREJE ) GO TO 1000
-* |--<--<--<--<--< Interaction rejected because of the decrease in beta
- RNDREJ = RNDREJ / FREJE
-* | Sample the emitted photon energy:
- EPHSMP = RNDREJ * ( EMXCER (MMAT) - EMNCER (MMAT) )
- & + EMNCER (MMAT)
-* | 2pi x freq.:
- OMGSMP = GEVOMG * EPHSMP
-* | Wavelength:
- WVLSMP = TWOPIP * CLIGHT / OMGSMP
-* | Compute the refraction index at the right wavelength/frequency:
- RFISMP = FOPTPR ( 1, WVLSMP, OMGSMP, MMAT )
- BETNSM = PTRKCR / ETRKCR * RFISMP
-*dbgD WRITE (77,*) ' CRNKVP:BETNSM,RFISMP,EPHSMP,WVLSMP,OMGSMP',
-*dbgD & BETNSM,RFISMP,EPHSMP,WVLSMP,OMGSMP
- IF ( BETNSM .LE. ONEONE ) GO TO 1000
-* |--<--<--<--<--< No chance to emit Cerenkov photons at this frequency
- BTNFSM = ( BETNSM - ONEONE ) * ( BETNSM + ONEONE ) / BETNSM**2
-* | Compute the quantum efficiency for this emitted energy:
- OPSENS = FOPTSN ( WVLSMP, OMGSMP )
- FREJE = BTNFSM / BTNFCR * OPSENS / OPSNMX
-*D IF ( FREJE .GT. ONEPLS ) WRITE (77,*)
-*D & ' ^^^CRNKVP:FREJE,RFISMP,RMXCER(MMAT),OPSENS,MMAT',
-*D & FREJE,RFISMP,RMXCER(MMAT),OPSENS,MMAT
- RNDREJ = FLRNDM (RNDREJ)
-*dbgD WRITE (77,*) ' CRNKVP:BTNFSM,BTNFCR,RNDREJ',
-*dbgD & BTNFSM,BTNFCR,RNDREJ
- IF ( RNDREJ .GE. FREJE ) GO TO 1000
-* |--<--<--<--<--< Interaction rejected because of the decrease in n
-* | Eventually a Cerenkov photon is going to be emitted:
-* | +----------------------------------------------------------------*
-* | | New production sub-step:
- IF ( NTRKCR .NE. NTRKLD ) THEN
- UTRKCR = XTRACK (NTRKCR) - XTRACK (NTRKCR-1)
- VTRKCR = YTRACK (NTRKCR) - YTRACK (NTRKCR-1)
- WTRKCR = ZTRACK (NTRKCR) - ZTRACK (NTRKCR-1)
- TUVWCR = SQRT ( UTRKCR**2 + VTRKCR**2 + WTRKCR**2 )
- UTRKCR = UTRKCR / TUVWCR
- VTRKCR = VTRKCR / TUVWCR
- WTRKCR = WTRKCR / TUVWCR
- SINT02 = UTRKCR**2 + VTRKCR**2
- LSMPAN = SINT02 .LT. ANGLSQ
- COSTH0 = WTRKCR
- IF ( .NOT. LSMPAN ) THEN
- SINTH0 = SQRT (SINT02)
- COSPH0 = UTRKCR / SINTH0
- SINPH0 = VTRKCR / SINTH0
- END IF
- NTRKLD = NTRKCR
- END IF
-* | |
-* | +----------------------------------------------------------------*
- SRNRSD = TUVWCR * TRNRSD / CRVCRR / TTRACK (NTRKCR)
-*D IF ( SRNRSD .LT. ZERZER .OR. SRNRSD .GT. TUVWCR )
-*D & WRITE (77,*)' ^^^CRNKVP:SRNRSD,TUVWCR,NTRKCR',
-*D & SRNRSD,TUVWCR,NTRKCR
- XTRKCR = XTRACK (NTRKCR) - UTRKCR * SRNRSD
- YTRKCR = YTRACK (NTRKCR) - VTRKCR * SRNRSD
- ZTRKCR = ZTRACK (NTRKCR) - WTRKCR * SRNRSD
-* | +----------------------------------------------------------------*
-* | | Put here the check on lattice change:
- IF ( NULTTC .GT. 0 ) THEN
-* | | Current (possibly curved in mag. field) non wiggled path:
- TRLATT = TRUSED / CRVCRR
- IF ( TRLATT .LE. TSLTTC (IRLTCR) ) THEN
- LATCCR = MULTTC (IRLTCR)
- ELSE
- DO 2000 IL = IRLTCR + 1, NULTTC
- IF ( TRLATT .LE. TSLTTC (IL) ) THEN
- IRLTCR = IL
- LATCCR = MULTTC (IRLTCR)
- GO TO 2010
- END IF
- 2000 CONTINUE
-*D CALL FLABRT ( 'CRNKVP', 'STOP:CRNKVP-NO-VALID-LATTICE' )
- 2010 CONTINUE
- END IF
- END IF
-* | |
-* | +----------------------------------------------------------------*
-* | +----------------------------------------------------------------*
-* | | Put here the check on the end point in magnetic field:
- IF ( LEMAGN .AND. NTRKCR .EQ. NTRACK .AND. MREG .NE. NEWREG )
- & THEN
- NEWLSV = NEWLAT
- MLATSV = MLATTC
- NROLD = NEWREG
- CALL GEOMAG ( XTRKCR, YTRKCR, ZTRKCR, UTRKCR, VTRKCR,
- & WTRKCR, NWREG , NROLD , IDSCCK )
- NEWLAT = NEWLSV
- MLATTC = MLATSV
- IF ( NWREG .EQ. NROLD ) GO TO 7000
-* | |-->-->--> already beyond the real crossing point:
- END IF
-* | |
-* | +----------------------------------------------------------------*
- LSTOPP = LSTOPP + 1
- NUMOPH = NUMOPH + 1
- IF ( NUMOPH .GT. 1000000000 ) THEN
- MUMOPH = MUMOPH + 1
- NUMOPH = NUMOPH + 1000000000
- END IF
- WOPTPH = WOPTPH + WTRACK
- POPTPH (LSTOPP) = EPHSMP
-* | +----------------------------------------------------------------*
-* | | No easy safe way to reconstruct Dnear along a magnetic field
-* | | step:
- IF ( LEMAGN ) THEN
- DONEAR (LSTOPP) = ZERZER
-* | |
-* | +----------------------------------------------------------------*
-* | | Questionable: ddnear actually refers to the end track point
- ELSE
- DONEAR (LSTOPP) = MAX ( DDNEAR - CTRACK + TRUSED, ZERZER )
- END IF
-* | |
-* | +----------------------------------------------------------------*
- XOPTPH (LSTOPP) = XTRKCR
- YOPTPH (LSTOPP) = YTRKCR
- ZOPTPH (LSTOPP) = ZTRKCR
- NREGOP (LSTOPP) = MREG
- NLATOP (LSTOPP) = LATCCR
- COSTHE = ONEONE / BETNSM
- SINTHE = SQRT ( ( ONEONE - COSTHE ) * ( ONEONE + COSTHE ) )
- CALL SFECFE ( SINPHI, COSPHI )
- UPRIME = COSPHI * SINTHE
- VPRIME = SINPHI * SINTHE
- WPRIME = COSTHE
- IF ( LSMPAN ) THEN
- TXOPPH (LSTOPP) = UPRIME
- TYOPPH (LSTOPP) = VPRIME
- TZOPPH (LSTOPP) = WPRIME * COSTH0
- ELSE
- TXOPPH (LSTOPP) = UNDOXR ( UPRIME, VPRIME, WPRIME, SINPH0,
- & COSPH0, SINTH0, COSTH0 )
- TYOPPH (LSTOPP) = UNDOYR ( UPRIME, VPRIME, WPRIME, SINPH0,
- & COSPH0, SINTH0, COSTH0 )
- TZOPPH (LSTOPP) = UNDOZR ( UPRIME, VPRIME, WPRIME, SINPH0,
- & COSPH0, SINTH0, COSTH0 )
- END IF
-* | Set-up the polarization vector: the following for linear pola-
-* | rization lying in the plane defined by the photon and the
-* | original particle direction:
- TXPOPP (LSTOPP) = UTRKCR / SINTHE - COSTHE / SINTHE
- & * TXOPPH (LSTOPP)
- TYPOPP (LSTOPP) = VTRKCR / SINTHE - COSTHE / SINTHE
- & * TYOPPH (LSTOPP)
- TZPOPP (LSTOPP) = WTRKCR / SINTHE - COSTHE / SINTHE
- & * TZOPPH (LSTOPP)
- SCADOT = TXOPPH (LSTOPP) * TXPOPP (LSTOPP)
- & + TYOPPH (LSTOPP) * TYPOPP (LSTOPP)
- & + TZOPPH (LSTOPP) * TZPOPP (LSTOPP)
- IF ( SCADOT .GT. CSNPRN ) THEN
- WRITE (LUNERR,'(A,1PG23.15,A)')
- & ' *** Crnkvp: bad polarization',
- & SCADOT, ' ***'
-*D WRITE (77,'(A,1PG23.15)')
-*D & ' ^^^Crnkvp: bad polarization',
-*D & SCADOT
- END IF
- WTOPPH (LSTOPP) = WTRACK
- AGOPPH (LSTOPP) = ATRKCR
- CMPOPP (LSTOPP) = ZERZER
- LOOPPH (LSTOPP) = LTRACK + 1
-*
-*
-* Hook to TFluka
-*
- PXCR = EPHSMP * TXOPPH (LSTOPP)
- PYCR = EPHSMP * TYOPPH (LSTOPP)
- PZCR = EPHSMP * TZOPPH (LSTOPP)
- POX = TXPOPP(LSTOPP)
- POY = TYPOPP(LSTOPP)
- POZ = TZPOPP(LSTOPP)
- CALL pshckp(PXCR, PYCR, PZCR, EPHSMP, XTRKCR,
- & YTRKCR , ZTRKCR, ATRKCR, POX, POY, POZ, WTRACK, ITFL)
- NPROD = NPROD + 1
- CALL ustckv(NPROD, MREG, XTRKCR, YTRKCR, ZTRKCR)
-*
-*
-*
-* | !!!!!! Here Stuprf should be used !!!!!!
- LOUOPP (LSTOPP) = ITFL
- DO 2100 ISPR = 1, MKBMX1
- SPAROK (ISPR,LSTOPP) = SPAUSR (ISPR)
- 2100 CONTINUE
- DO 2200 ISPR = 1, MKBMX2
- ISPORK (ISPR,LSTOPP) = ISPUSR (ISPR)
- 2200 CONTINUE
-* | Save the parent features:
- TPROPP (LSTOPP) = ETRACK - AMPRTC
- APROPP (LSTOPP) = ATRKCR
- LPROPP (LSTOPP) = LTRACK
- IPROPP (LSTOPP) = JTRACK
- NPROPP (LSTOPP) = 0
-* | Update the "current" particle values:
- ETRKCR = ETRKCR - EPHSMP
- DEDXCK = DEDXCK + EPHSMP
- PTRKCR = SQRT ( ( ETRKCR - AMPRTC ) * ( ETRKCR + AMPRTC ) )
- PTINTR = PTRKCR
- ETINTR = ETRKCR
- TRINTR = ZERZER
- BETNLD = BETNCR
- BETNCR = PTRKCR / ETRKCR * RMXCER (MMAT)
-* | No longer any chance to emit Cerenkov photons:
- IF ( BETNCR .LE. ONEONE ) GO TO 7000
- BTNFLD = BTNFCR
- BTNFCR = ( BETNCR - ONEONE ) * ( BETNCR + ONEONE ) / BETNCR**2
-* | Update the macroscopic sigma: beta can only decrease
- SIGMCK = SIGMCK * BTNFCR / BTNFLD
- GO TO 1000
-* |
-* +-------------------------------------------------------------------*
- 7000 CONTINUE
- RETURN
-*=== End of subroutine Crnkvp =========================================*
- END
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff