pythiaeRHIC/radgen/radgen_init.f

0001
0002       subroutine radgen_init(bUseLUT,bGenLUT)
0003
0004       implicit none
0005
0006        include "mc_set.inc"
0007        include "density.inc"
0008        include "cmpcom.inc"
0009        include "tailcom.inc"
0010        include "radgen.inc"
0011        include "radgenkeys.inc"
0012
0013       logical bUseLUT,bGenLUT
0014       CHARACTER*256 tname
0015       COMMON/PYNUCL/INUMOD,CHANUM,ORDER
0016       DOUBLE PRECISION INUMOD,CHANUM
0017       INTEGER ORDER
0018
0019 ! ... force block data modules to be read
0020       external radata
0021
0022 *     kill_elas_res =0, all rad corrections
0023 *     kill_elas_res =1, ???
0024 *     kill_elas_res =2, no rad corrections
0025       kill_elas_res=0
0026
0027 *     initialise the rad. correction part !
0028 *     =  2 : fint cernlib extrapolation routine
0029 *     =  1 : 2dim spline
0030 *     =  0 : make lookuptable
0031 *     = -1 : do not lookup table , calculate all events
0032       if( bUseLUT ) then
0033         if( bGenLUT ) then
0034           ixytest = 0
0035         else
0036           ixytest = 2
0037         endif
0038       else
0039         ixytest = -1
0040       endif
0041
0042       if ((INUMOD.gt.1).and.(CHANUM.gt.1)) then
0043          mcSet_TarA = INUMOD
0044          mcSet_TarZ = CHANUM
0045       endif
0046
0047       write(*,*) mcSet_TarA, mcSet_TarZ
0048
0049 *----------------------------
0050 *     ire is target 1,2,3
0051 *     add neutron target
0052       if(     mcSet_TarA .eq. 1 .and. mcSet_TarZ .eq. 0 ) then
0053         ire = 0
0054       elseif( mcSet_TarA .eq. 1 .and. mcSet_TarZ .eq. 1 ) then
0055         ire = 1
0056       elseif( mcSet_TarA .eq. 2 .and. mcSet_TarZ .eq. 1 ) then
0057         ire = 2
0058       elseif( mcSet_TarA .eq. 3 .and. mcSet_TarZ .eq. 2 ) then
0059         ire = 3
0060       elseif( mcSet_TarA .eq. 4 .and. mcSet_TarZ .eq. 2 ) then
0061         ire = 4
0062       elseif( mcSet_TarA .eq. 14 .and. mcSet_TarZ .eq. 7 ) then
0063         ire = 14
0064       elseif( mcSet_TarA .eq. 20 .and. mcSet_TarZ .eq. 10 ) then
0065         ire = 20
0066       elseif( mcSet_TarA .eq. 40 .and. mcSet_TarZ .eq. 20 ) then
0067         ire = 40
0068       elseif( mcSet_TarA .eq. 84 .and. mcSet_TarZ .eq. 36 ) then
0069         ire = 84
0070       elseif( mcSet_TarA .eq. 131 .and. mcSet_TarZ .eq. 54 ) then
0071         ire = 131
0072       elseif( mcSet_TarA .eq. 197 .and. mcSet_TarZ .eq. 79 ) then
0073         ire = 197
0074       elseif( mcSet_TarA .eq. 238 .and. mcSet_TarZ .gt. 92 ) then
0075         ire = 238
0076       else
0077         write(*,*)( 'RADGEN_INIT: invalid target selection' )
0078       endif
0079
0080 C... Using radgen with Pythia the target and beam polarisation can be
0081 C    hardcoded to zero
0082       plrun = 0.
0083       pnrun = 0.
0084       if(ire.lt.10) then
0085          tname='radgen/xytab0unp.dat'
0086       elseif (ire.lt.100) then
0087          tname='radgen/xytab00unp.dat'
0088       else
0089          tname='radgen/xytab000unp.dat'
0090       endif
0091
0092       if (ire.lt.10) then
0093           write(tname(13:13),'(i1)')ire
0094       elseif (ire.lt.100) then
0095           write(tname(13:14),'(i2)')ire
0096       else
0097           write(tname(13:15),'(i3)')ire
0098       endif
0099
0100 *----------------------------
0101 * grid of important regions in theta (7*ntk)
0102       ntk = 35
0103 *----------------------------
0104 * photonic energy grid
0105       nrr = 100
0106 *----------------------------
0107 * min. energy in the calo (resolution parameter)
0108       demin=0.10
0109
0110 *----------------------------
0111       ap=2.*amp
0112       amp2=amp**2
0113       ap2=2.*amp**2
0114       if(kill_elas_res.eq.1) amc2=4d0
0115
0116       if(ire.eq.0) then
0117         amt=.939565d0
0118         rtara=1d0
0119         rtarz=0d0
0120         fermom=0d0
0121       elseif(ire.eq.1)then
0122         amt=.938272d0
0123         rtara=1d0
0124         rtarz=1d0
0125         fermom=0d0
0126       elseif(ire.eq.2)then
0127         amt=1.87561d0
0128         rtara=2d0
0129         rtarz=1d0
0130         fermom=.07d0
0131       elseif(ire.eq.3)then
0132         amt=2.8161d0
0133         rtara=3d0
0134         rtarz=2d0
0135         fermom=.087d0
0136       elseif(ire.eq.4)then
0137         amt=3.75567d0
0138         rtara=4d0
0139         rtarz=2d0
0140         fermom=.087d0
0141       elseif(ire.eq.14)then
0142         amt=13.1447d0
0143         rtara=14d0
0144         rtarz=7d0
0145         fermom=.100d0
0146         call fordop
0147       elseif(ire.eq.20)then
0148         amt=18.7783d0
0149         rtara=20d0
0150         rtarz=10d0
0151         fermom=.100d0
0152         call fordop
0153       elseif(ire.eq.40)then
0154         amt=37.5567d0
0155         rtara=40d0
0156         rtarz=20d0
0157         fermom=.100d0
0158         call fordop
0159       elseif(ire.eq.84)then
0160         amt=78.8769d0
0161         rtara=84d0
0162         rtarz=36d0
0163         fermom=.105d0
0164         call fordop
0165       elseif(ire.eq.131)then
0166         amt=123.0132d0
0167         rtara=131d0
0168         rtarz=54d0
0169         fermom=.105d0
0170         call fordop
0171       elseif(ire.eq.197)then
0172         amt=184.9922d0
0173         rtara=197d0
0174         rtarz=79d0
0175         fermom=.105d0
0176         call fordop
0177       elseif(ire.eq.207)then
0178         amt=194.3839d0
0179         rtara=207d0
0180         rtarz=82d0
0181         fermom=.105d0
0182         call fordop
0183       elseif(ire.eq.238)then
0184         amt=223.4975d0
0185         rtara=238d0
0186         rtarz=92d0
0187         fermom=.105d0
0188         call fordop
0189       endif
0190
0191 *-----------------------------
0192       if (mcSet_XMin.ge. 1.e-05) then
0193         ntx=ntdis
0194         write(*,*)
0195      &  ('Radiative corrections for DIS kinematics (xmin=1.e-05')
0196       elseif (mcSet_XMin.ge. 1.e-09) then
0197         ntx=ntpho
0198         write(*,*)
0199      &  ('Radiative corrections for photoproduction (xmin=1.e-09)')
0200         write(*,*)
0201      &  ('Elastic and quesielastric contributions set to ZERO !')
0202         if(ixytest.gt.0 .and. mcSet_XMin.lt.1.e-07) then
0203           write(*,*)
0204      &  ('Xmin in lookup table = 1.e-07 --> extrapolation to 1.e-09')
0205         endif
0206         if(ixytest.eq.0 .and. mcSet_XMin.lt.1.e-07) then
0207           write(*,*)
0208      &  ('Xmin in lookup table = 1.e-07 --> change minimum x')
0209           stop
0210         endif
0211       else
0212         write(*,*)
0213      &  ('Minimum x value below minimum value of 10**-9 ! ')
0214         stop
0215       endif
0216
0217 * initialize lookup table
0218       if(ixytest.ge.0) then
0219
0220 * number of x bins depending on x range desired
0221 * finer grid at lower x values
0222 * number of y bins = number of x bins (equidistant in y)
0223 * (needed to avoid double loop - speed optimisation)
0224       write(*,*) ('*********************************************')
0225       write(*,*)
0226      &  ('Make sure that you did create the correct lookup table')
0227       write(6,*) 'number of x bins in RC table = ',ntx
0228       write(6,*) 'size of x bins in RC table depending on x'
0229       nty=ntx
0230       write(6,*) 'number of y bins in RC table = ',nty
0231       write(6,*) 'size of y bins in RC table = ',
0232      &  (radgen_ymax-radgen_ymin)/(nty-1)
0233       write(*,*) ('*********************************************')
0234
0235         if(ixytest.eq.0) then
0236           write(*,*)('RADGEN_INIT: Creating lookup table '//tname)
0237         elseif (ixytest.eq.2) then
0238            write(*,*)
0239      +       ('RADGEN_INIT: Loading lookup table '//tname)
0240         endif
0241         call xytabl(tname,mcSet_EneBeam,plrun,pnrun,ixytest,ire)
0242       endif
0243
0244       end