GEANT321/fluka/nucrel.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:04  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.44  by  S.Giani
  11 *-- Author :
  12 *$ CREATE NUCREL.FOR
  13 *COPY NUCREL
  14 *
  15 *=== nucrel ===========================================================*
  16 *
  17 C
  18 C   HJM 17/11/88
  19 C
  20 C                MODIFIED ELASTIC SCATTERING ROUTINE
  21 C                (COMP. KMU-HEP INTERNAL NOTE 88-01)
  22 C
  23 C------------------------------------------------------------------
  24 C
  25       SUBROUTINE NUCREL(IT,PLAB,EKIN,CX,CY,CZ,ANUC,WEE)
  26
  27 #include "geant321/dblprc.inc"
  28 #include "geant321/dimpar.inc"
  29 #include "geant321/iounit.inc"
  30 C********************************************************************
  31 C     VERSION JULY 81 BY             J. RANFT
  32 C                                    LEIPZIG
  33 C     LAST CHANGE A. FERRARI 26-9-89: RECOIL ENERGY SCORING ADDED
  34 C
  35 C     THIS IS A SUBROUTINE OF FLUKA82 TO SAMPLE ELASTIC INTERACTIONS
  36 C
  37 C     INPUT VARIABLES:
  38 C     IT       = TYPE OF THE PRIMARY
  39 C     PLAB     = PRIMARY PARTICLE LAB MOMENTUM (GEV/C)
  40 C     CX,CY,CZ = PRIMARY PARTICLE DIRECTION COSINES
  41 C     ANUC     = ATOMIC WEIGHT OF THE TARGET NUCLEUS
  42 C
  43 C     OUTPUT VARIABLES:
  44 C     CCX,CCY,CCZ = DIRECTION COSINES OF SCATTERED PARTICLE
  45 C
  46 C********************************************************************
  47 C
  48       PARAMETER (AMUC12 = 0.9314943228D0)
  49 #include "geant321/aadat.inc"
  50 #include "geant321/finuc.inc"
  51 *------ Common part has been added by A. Ferrari ----------------------*
  52 #include "geant321/part2.inc"
  53 C----------------------------------------------------------------
  54       REAL RNDM(1)
  55       DIMENSION ITT(NALLWP),SHP(6),BP(6),XNN(6,6),BA(6,6),XNA(6),AM(6)
  56       SAVE ITT,SHP,BP,XNN,BA,XNA,AM
  57       DATA ITT/1,2,0,0,0,0,0,1,2,0,0,5,3,4,5,6,1,2,5,1,1,1,3,5,6,
  58      &         3,0,0,0,1,2,2,2,1,2,1,2,1,2/
  59 C     ITT(IT) - INTERNAL PARTICLE NUMBER
  60       DATA AM/ 9.01D0, 12.01D0, 26.98D0, 63.55D0, 118.69D0, 207.19D0/
  61 C     AM - ATOMIC WEIGHTS OF SEVERAL TARGET MATERIALS
  62       DATA XNA/3.5D0, 3.4D0, 4.5D0, 6.7D0, 8.2D0, 9.5D0/
  63 C     XNA(AM) - EFFECTIVE NUCLEON NUMBER
  64       DATA BP/12.D0, 12.D0, 4*10.D0/
  65 C     BP(IIT) - INCOHERENT SLOPE (GEV/C)**-2
  66       DATA SHP/2*40.D0, 2*25.D0, 2*20.D0/
  67 C     SHP(IIT) - TOTAL HADRON - NUCLEON CROSS SECTION  (MB)
  68 C(((((((((((((((      BA(IIT,A)    XNN(IIT,A)       (((((((((((((((((
  69 C     BA - COHERENT SLOPE   (GEV/C)**-2
  70 C     XNN - NORMALIZATION FOR PROJECTILE -NUCLEUS SCATTERING  (MB)
  71       DATA XNN/
  72      ( 256.6D0,  240.5D0,  181.0D0,  167.1D0,  144.0D0,  139.0D0,
  73      ( 335.3D0,  355.0D0,  241.0D0,  229.7D0,  181.5D0,  196.0D0,
  74      ( 755.1D0, 1344.0D0,  506.0D0,  480.8D0,  431.0D0,  403.0D0,
  75      (1656.0D0, 1886.0D0, 1104.0D0, 1144.0D0,  914.0D0,  939.0D0,
  76      (3005.0D0, 2559.0D0, 2112.0D0, 2068.0D0, 1723.0D0, 1791.0D0,
  77      (5123.0D0, 5444.0D0, 4071.0D0, 3391.0D0, 3368.0D0, 3168.0D0/
  78       DATA BA/
  79      ( 72.0D0,  69.3D0,  65.2D0,  63.4D0,  60.0D0,  65.6D0,
  80      ( 72.3D0,  75.9D0,  65.5D0,  63.1D0,  59.5D0,  68.1D0,
  81      (119.4D0, 128.7D0, 107.3D0, 105.0D0, 105.6D0,  98.8D0,
  82      (201.0D0, 212.0D0, 188.9D0, 183.0D0, 179.0D0, 178.0D0,
  83      (311.0D0, 315.0D0, 286.0D0, 277.0D0, 270.0D0, 266.0D0,
  84      (454.0D0, 448.0D0, 435.0D0, 397.0D0, 416.0D0, 403.0D0/
  85 C
  86 C********************************************************************
  87 C     PARAMETRIZATION
  88 C     DSIGMA/DOMEGA=C*EXP(-DS*THETA**2)+E*EXP(-FS*THETA**2)
  89 C                   COHERENT PART             INCOHERENT PART
  90 C********************************************************************
  91 C
  92 C***
  93       NP=1
  94       TV=0.D0
  95 C***
  96       IF(IT.EQ.30) GO TO 301
  97       IF(ANUC.LT.0.99D0) THEN
  98          CXR(1)=CX
  99          CYR(1)=CY
 100          CZR(1)=CZ
 101          TKI(1)=EKIN
 102          PLR(1)=PLAB
 103          KPART(1)=IT
 104          RETURN
 105       ENDIF
 106 C
 107 C--------------------------------------------------------------------
 108 C         HJM 10/88       ELASTIC SCATTERING INTO 2-BODY FINALSTATE
 109 C                         FOR NUCLEON-PROTON INITIAL STATE
 110 C         AF  9/91        ALSO PBAR,NBAR-PROTON
 111 C                         FINAL STATE PARTICLES IN /FINUC/
 112       IF( (ANUC.LT.1.5D0) .AND. (IT.LE.2 .OR. IT.EQ.8 .OR. IT.EQ.9)
 113      &    .AND. (EKIN.LT.3.5D0) ) THEN
 114          CALL NUPREL(IT,EKIN,PLAB,CX,CY,CZ)
 115          WEI(1)=WEE
 116          WEI(2)=WEE
 117 C
 118          RETURN
 119       ENDIF
 120 C-------------------------------------------------------------------
 121 C
 122 C********************************************************************
 123 C     1=P,N, 2=AP,AN, 3=PI+, 4=PI-, 5=K+,K0, 6=K-,K0 BAR
 124 C********************************************************************
 125 C
 126       IIT=ITT(IT)
 127       IF(IIT.EQ.0)RETURN
 128       AP=PLAB
 129       AP2=AP**2
 130       A=ANUC
 131       IF(IIT.GE.2) GO TO 101
 132       IF(AP.GT.20.D0) GO TO 101
 133 C
 134 C********************************************************************
 135 C     FOR PROTONS BELOW 20 GEV/C
 136 C     DSIGMA/DOMEGA=12.5*A**1.63*EXP(-14.5*A**0.66*T)+
 137 C                   17.5*A**0.33*EXP(-10*T)    FOR A<62
 138 C
 139 C     DSIGMA/DOMEGA=50*A**1.33*EXP(-60*A**0.33*T)+
 140 C     20*A**0.4*EXP(-10*T)  FOR A>62
 141 C********************************************************************
 142 C
 143       A3=A**0.3333333333333333D0
 144       ATAR=A
 145       IF (ATAR .GT. 62.D0) GO TO 1
 146       C=12.5D0*A**1.63D0
 147       DS=14.5D0*A3*A3*AP2
 148       E=17.5D0*A3
 149       GO TO 2
 150   1   CONTINUE
 151       C=50.D0*A3*A
 152       DS=60.D0*A3*AP2
 153       E=20.D0*A**0.4D0
 154   2   CONTINUE
 155       FS=10.D0*AP2
 156       GO TO 3
 157 C
 158 C********************************************************************
 159 C     FOR PROTONS OVER 20 GEV/C AND OTHER PARTICLES
 160 C********************************************************************
 161 C
 162   101 CONTINUE
 163       DO 200 I=1,6
 164       IF(A.LE.AM(I)) GO TO 201
 165   200 CONTINUE
 166       K=6
 167       GO TO 202
 168   201 K=I
 169   202 KK=K+1
 170       IF(KK.GT.6) KK=K-1
 171       XNEL=XNN(IIT,K)+(A-AM(K))*(XNN(IIT,KK)-XNN(IIT,K))/(AM(KK)-AM(K))
 172       C=XNEL**2
 173       DS=LOG(BA(IIT,KK)/BA(IIT,K))/LOG(AM(KK)/AM(K))
 174       DS=BA(IIT,K)*(A/AM(K))**DS
 175       DS=DS*AP2
 176       E=XNA(K)+(A-AM(K))*(XNA(KK)-XNA(K))/(AM(KK)-AM(K))
 177       E=E*SHP(IIT)**2
 178       FS=BP(IIT)*AP2
 179       GO TO 3
 180 C
 181 C********************************************************************
 182 C     FOR HEAVY IONS:
 183 C
 184 C       FOR APROJ>25 OR (APROJ>10 AND ATARG>100)
 185 C DSIG/DOM = 216 * PTOT**1.86 * (AP**0.7 + AT**0.7)**2.2 *
 186 C          EXP(-16.1 * (AP**1.2 + AT**0.9)**0.8 * T) +
 187 C          0.3 * PTOT**1.86 * (AP**0.7 + AT**0.7)**2.2 * EXP(-23 * T)
 188 C
 189 C       FOR LIGHTER NUCLIDES:
 190 C DSIG/DOM = 78 * PTOT**1.78 * (AP + AT**0.9)**2.1 *
 191 C          EXP(-16.1 * (AP**1.2 + AT**0.9)**0.8 * T) +
 192 C          0.5 * PTOT**1.78 * (AP + AT**0.9)**2.1 * EXP(-30 * T)
 193 C
 194 C       THE PARAMETRIZATION IS OBTAINED BY FITTING TO
 195 C       DATA GENERATED BY THE SOFT-SPHERES MODEL.
 196 C       (REF. CHAUVIN & AL. PHYS.REV.C. 28 1970 (1983))
 197 C       THE PARAMETRIZATION IS QUITE CRUDE AND FAILS AT LARGE
 198 C       ANGLES.
 199 C       THERE IS NO EXPERIMENTAL VERIFICATION FOR THE MODEL AT
 200 C       ENERGIES ABOVE 100 MEV/A
 201 C                                   MIKA HUHTINEN 25.8.88
 202 C********************************************************************
 203 C
 204   301 CONTINUE
 205       ITARA=NINT(ANUC)
 206 C     ALWAYS LET THE LIGHTER PARTICLE BE THE PROJECTILE (THEN THE
 207 C     PARAMETRIZATION WILL GIVE BETTER RESULTS)
 208       IF (ITARA.GE.IPROA) THEN
 209         PROA=IPROA
 210         TARA=ITARA
 211       ELSE
 212         PROA=ITARA
 213         TARA=IPROA
 214       ENDIF
 215 C     TOTAL LAB. MOMENTUM
 216       PLTOT=PLAB*PROA
 217 C     SQUARE OF CMS MOMENTUM
 218       AP2=((PLTOT*TARA)/(TARA+PROA))**2
 219
 220       DS=16.1D0*(PROA**1.2D0 + TARA**0.9D0)**0.8D0*AP2
 221
 222       IF (PROA.GT.25.D0) GO TO 302
 223       IF ((PROA.GT.10.D0).AND.(TARA.GT.100.D0)) GO TO 302
 224
 225       CFF=PLTOT**1.78D0*(PROA + TARA**0.9D0)**2.1D0
 226       C=78.D0*CFF
 227       E=0.5D0*CFF
 228       FS=30.D0*AP2
 229       GO TO 3
 230
 231   302 CONTINUE
 232       CFF=PLTOT**1.86D0*(PROA**0.7D0 + TARA**0.7D0)**2.2D0
 233       C=216.D0*CFF
 234       E=0.3D0*CFF
 235       FS=23.D0*AP2
 236 C
 237 C********************************************************************
 238 C     SAMPLE THE ANGLE AND ROTATE TO GET NEW DIRECTION COSINES
 239 C********************************************************************
 240 C
 241     3 CONTINUE
 242       EDS=1.D-9
 243       IF(DS.LT.10.5D0) EDS=EXP(-2.D0*DS)
 244       EDS2=EDS**2
 245       DD=C*(1.D0-EDS2)/DS
 246       EFS=1.D-9
 247       IF(FS.LT.10.5D0) EFS=EXP(-2.D0*FS)
 248       EFS2=EFS**2
 249       FF=E*(1.D0-EFS2)/FS
 250       DF=2.D0*DS
 251       CALL GRNDM(RNDM,1)
 252       VV=RNDM(1)*(DD+FF)
 253       IF (VV.GT.DD) DF=2.D0*FS
 254   31  CONTINUE
 255       CALL GRNDM(RNDM,1)
 256       V2=LOG(RNDM(1))/DF
 257       IF (ABS(V2).GT.2.D0) GO TO 31
 258       V1=V2
 259       COD=V1+1.D0
 260       SID=SQRT(-V1*(2.D0+V1))
 261 *  ******* Computing the recoil energy!!!! A. Ferrari, 26-9-89 ******* *
 262       KP    = IPTOKP (IT)
 263       APP2  = AAM(KP)**2
 264       EPROJ = EKIN + AAM(KP)
 265 *  +-------------------------------------------------------------------*
 266 *  |  Hydrogen: generate the recoil proton
 267       IF ( ANUC .LT. 1.2D+00 ) THEN
 268          ATT  = AAM (1)
 269          ATT2 = ATT**2
 270          ETTT = EPROJ + ATT
 271          PTTL = PLAB
 272          PTTT = 0.D+00
 273          HELP = ATT2 - APP2 * ( 1.D+00 - COD ) * ( 1.D+00 + COD )
 274 *  | Now: it can happen for amproj > amprot on hydrogen.
 275          IF ( HELP .LT. 0.D+00 ) GO TO 3
 276          PLAB = PLAB * ( ( APP2  + ATT*EPROJ ) * COD + (EPROJ + ATT)
 277      &        * SQRT (HELP) ) / ( (EPROJ + ATT)**2 - (PLAB*COD)**2 )
 278          NP = 2
 279          TV = 0.D+00
 280          WEI   (1) = WEE
 281          WEI   (2) = WEE
 282          KPART (1) = IT
 283          KPART (2) = 1
 284          TKI   (1) = SQRT ( PLAB**2 + APP2 ) - AAM (KP)
 285          TKI   (2) = EKIN - TKI (1)
 286          EKIN      = TKI (1)
 287          PLR   (1) = PLAB
 288          PLR   (2) = SQRT ( TKI (2) * ( TKI (2) + 2.D+00 * AAM (1) ) )
 289          CALL SFECFE ( SFE, CFE )
 290          CALL TTRANS ( CX, CY, CZ, COD, SID, SFE, CFE, CCX, CCY, CCZ )
 291          CXR (1) = CCX
 292          CYR (1) = CCY
 293          CZR (1) = CCZ
 294          SINTHE  = PLAB / PLR (2) * SID
 295          COSTHE  = SQRT ( ( 1.D+00 - SINTHE ) * ( 1.D+00 + SINTHE ) )
 296          CFE     = - CFE
 297          SFE     = - SFE
 298          CALL TTRANS ( CX, CY, CZ, COSTHE, SINTHE, SFE, CFE, CXR (2),
 299      &                 CYR (2), CZR (2) )
 300          RETURN
 301 *  |
 302 *  +-------------------------------------------------------------------*
 303 *  |  "Heavy" nuclei
 304 *  |  the following line to patch el. scatt. of an, ap, lambdas, sigmas
 305 *  |  etc on hydrogen ( Now it should be useless )
 306       ELSE
 307          ATT  = MAX ( ANUC * AMUC12, AAM (KP) )
 308          ATT2 = ATT**2
 309          HELP  = ATT2 - APP2 * ( 1.D+00 - COD ) * ( 1.D+00 + COD )
 310          PLAB  = PLAB*( (APP2  + ATT*EPROJ)*COD + (EPROJ + ATT)
 311      &         * SQRT (HELP) ) / ( (EPROJ + ATT)**2 - (PLAB*COD)**2 )
 312          TV = EPROJ - SQRT ( PLAB**2 + APP2 )
 313          EKIN = EKIN - TV
 314 *  | ***** The following lines only for debugging *****
 315          IF ( PLAB .LE. 0.D+00 ) THEN
 316             WRITE (LUNOUT,*)' *** Nucrel: plab',PLAB,' ***'
 317             WRITE (LUNERR,*)' *** Nucrel: plab',PLAB,' ***'
 318          END IF
 319 *  | ***** End debugging lines *****
 320       END IF
 321 *  |              Now the recoil energy is put into Tv
 322 *  +-------------------------------------------------------------------*
 323       CALL SFECFE(SFE,CFE)
 324       CALL TTRANS(CX,CY,CZ,COD,SID,SFE,CFE,CCX,CCY,CCZ)
 325 C---------------------------------------------------------
 326 C   MODIFICATION HJM 28/10/88
 327 C   NOTE:  HEAVY ION OPTION  NOT CONSISTENTLY TREATED NOW]]]
 328 C
 329       WEI(1)=WEE
 330       CXR(1)=CCX
 331       CYR(1)=CCY
 332       CZR(1)=CCZ
 333       TKI(1)=EKIN
 334       PLR(1)=PLAB
 335       KPART(1)=IT
 336       RETURN
 337       END