* * $Id$ * * $Log$ * Revision 1.1.1.1 1995/10/24 10:19:58 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.44 by S.Giani *-- Author : *$ CREATE RBKEKV.FOR *COPY RBKEKV * * *=== rbkekv ===========================================================* * * SUBROUTINE RBKEKV ( IT, EXSOP, TO, AMSS, TKIN, TSTRCK, & PSTRCK, ARECL, TKRECL, COD, SID ) #include "geant321/dblprc.inc" #include "geant321/dimpar.inc" #include "geant321/iounit.inc" * *----------------------------------------------------------------------* * version by Alfredo Ferrari * INFN - Milan * last change 19 april 93 by Alfredo Ferrari * INFN - Milan * * To be called from the high energy production * * this is a subroutine of fluka to sample intranuclear cascade * particles: it is based on the old Rakeka from J. Ranft * * input variables: * it = type of the secondary requested; 1=proton, 2=neutron * bbtar = atomic weight of the medium * * output variables: * tkin = kinetic energy of the secondary in GeV before applying * the nuclear well (and eventually the Coulomb barreer) * tstrck= kinetic energy of the secondary in GeV * pstrck= momentum of the secondary in GeV/c * sid,cod = sine and cosine of the angle between * the directions of the primary * and secondary particles * ********************************************************************* * PARAMETER ( PI = PIPIPI ) PARAMETER ( PIO2 = PI / 2.D+00 ) #include "geant321/nucdat.inc" LOGICAL LINIT, LZEROI DIMENSION EXSOP (2), AMMOLD (2) DIMENSION PFINIT (50,2), TKINIT (50,2), TSINIT (50,2), & TRINIT (50,2), NINI (2) REAL RNDM(4) SAVE AMMOLD, PFINIT, TKINIT, TSINIT, TRINIT, NINI DATA AMMOLD / 0.9382796D+00, 0.9395731D+00 / * LINIT = .FALSE. IF ( NINI (IT) .GT. 0 ) THEN PKFRMI = PFINIT (NINI(IT),IT) TKIN = TKINIT (NINI(IT),IT) TSTRCK = TSINIT (NINI(IT),IT) TKRECL = TRINIT (NINI(IT),IT) NINI (IT) = NINI (IT) - 1 EXSOP (IT) = 0.D+00 GO TO 3500 ELSE GO TO 200 END IF ENTRY RBKINI ( IT, LZEROI, EXSOP, TKIN, TSTRCK, & PSTRCK, ARECL, TKRECL ) LINIT = .TRUE. IF ( LZEROI ) THEN NINI (1) = 0 NINI (2) = 0 RETURN END IF 200 CONTINUE * In this version the low energy component has been suppressed, since * they are now produced by the evaporation model, A. Ferrari. * The parameters needed for the sampling have been already set in * Corrin: EXSOP (IT) = 0.D+00 * Sample the Fermi momentum ESLPFF = ESLOPE (IT) EXUPFF = EXUPNU (IT) EKUPFF = EKUPNU (IT) EXDWFF = EXMNNU (IT) EKDWFF = EKMNNU (IT) ERCLFF = ERCLAV (IT) IRECNT = 0 * +-------------------------------------------------------------------* * | 100 CONTINUE * | Sample the Fermi momentum CALL GRNDM(RNDM,4) PKFRMI = PFRMMX (IT) * MAX ( RNDM (1), RNDM (2), & RNDM (3) ) PKFRSQ = PKFRMI * PKFRMI TKIN = - ESLPFF * LOG ( EXDWFF - RNDM (4) * ( & EXDWFF - EXUPFF ) ) TKRECL = 0.5D+00 * PKFRSQ / ( AMUC12 * ARECL ) * ( 1.D+00 - & 0.25D+00 * PKFRSQ / ( ARECL**2 * AMUCSQ ) ) TSTRCK = SQRT ( AMNUSQ (IT) + PKFRSQ ) + TKIN - AMNUCL (IT) & - V0WELL (IT) - TKRECL - EBNDNG (IT) * | +----------------------------------------------------------------* * | | IF ( TSTRCK .LE. VEFFNU (IT) - V0WELL (IT) ) THEN EXSOP (IT) = EXSOP (IT) + TKIN IRECNT = IRECNT + 1 * | | +-------------------------------------------------------------* * | | | IF ( IRECNT .GT. 10 ) THEN * | | | +----------------------------------------------------------* * | | | | IF ( TKRECL - ERCLFF .GT. EKUPFF - EKDWFF .AND. & IRECNT .LT. 20 ) THEN ERCLFF = ERCLFF + TKRECL - ERCLFF EKUPFF = EKUPFF + TKRECL - ERCLFF EKDWFF = EKDWFF + TKRECL - ERCLFF AHELP = EXP ( - ( TKRECL - ERCLFF ) / ESLPFF ) EXUPFF = EXUPFF * AHELP EXDWFF = EXDWFF * AHELP * | | | | * | | | +----------------------------------------------------------* * | | | | ELSE IF ( IRECNT .GT. 15 ) THEN TKRECL = MAX ( TKRECL, ERCLFF ) ERCLFF = ERCLFF + TKRECL EKUPFF = EKUPFF + TKRECL - ERCLFF EKDWFF = EKDWFF + TKRECL - ERCLFF AHELP = EXP ( - TKRECL / ESLPFF ) EXUPFF = EXUPFF * AHELP EXDWFF = EXDWFF * AHELP END IF * | | | | * | | | +----------------------------------------------------------* END IF * | | | * | | +-------------------------------------------------------------* GO TO 100 * +-<|--<--<--<--<--< go to resampling END IF * | | * | +----------------------------------------------------------------* * +-------------------------------------------------------------------* * | IF ( LINIT ) THEN NINI (IT) = NINI (IT) + 1 PFINIT (NINI(IT),IT) = PKFRMI TKINIT (NINI(IT),IT) = TKIN TSINIT (NINI(IT),IT) = TSTRCK TRINIT (NINI(IT),IT) = TKRECL RETURN END IF * | * +-------------------------------------------------------------------* 3500 CONTINUE * | Masses have been updated PSTRCK = SQRT ( TSTRCK * ( TSTRCK + 2.D+00 * AMNUCL (IT) ) ) * ********************* Sample the angle ******************************** * Polar angle selection ADE=0.090D0*(1.D0+0.081D0*ATO1O3)/TKIN DEX=EXP(- PIO2 * PIO2 / ADE) AN1=(1.D0-DEX)*ADE/2.D0 AN2=DEX*PIO2 AN=AN1+AN2 AN1=AN1/AN CALL GRNDM(RNDM,1) IF(RNDM(1).GT.AN1) GO TO 3 2 CONTINUE CALL GRNDM(RNDM,1) DE=SQRT(-ADE*LOG(1.D0-RNDM(1)*(1.D0-DEX))) IF(DE.GT.PIO2) GO TO 2 C WRITE(LUNOUT,*)IT,TO,AMSS,SQAMSS,T,P,DE COD = COS (DE) SID = SIN (DE) RETURN 3 CONTINUE CALL GRNDM(RNDM,1) COD = - RNDM(1) SID = SQRT ( (1.D0 + COD ) * ( 1.D0 - COD ) ) C WRITE(LUNOUT,*)IT,TO,AMSS,SQAMSS,T,P,DE RETURN ENTRY RBKMIN (IT) NINI(IT) = NINI(IT)-1 RETURN END