bin/HARDCOL/bfkl.f

0001 C
0002 C     Returns the interpolated value of the BFKL
0003 C     amplitude for hard color singlet exchange.
0004 C     It needs a file in which the exact values
0005 C     are stored in the chosen form.
0006 C
0007 C     Modifications for PYTHIA:
0008 C        -gives result in GeV-2 instead of mb
0009 C        -not function of controll anymore
0010 C        -checks range of (y,q)
0011
0012
0013 C
0014 C    initialisation: reads tabularized values of cross section
0015 C
0016
0017 C    Modified by Sheila Amaral, 16 Oct 2009
0018 C
0019 C    Implemented the file path using environment variable CMSSW_BASE
0020 C    The .dat files should be in data/ directory
0021
0022       SUBROUTINE read_hcs_file
0023       IMPLICIT NONE
0024
0025 C   Pythia variables
0026       INTEGER MSTP,MSTI
0027       DOUBLE PRECISION PARP,PARI
0028       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0029
0030 C   Local variables
0031       REAL*8 hcstab(12,25,3)
0032       COMMON / hcstab / hcstab
0033
0034       CHARACTER hcs_name*120   !!
0035       COMMON/ hcsfile / hcs_name
0036
0037       REAL*8 x1,x2
0038
0039       INTEGER iq,iy
0040
0041 C  Choose value of s_0
0042 C  The *.dat files must reside in the directory where the program is run.
0043 C
0044       character cmsdir*82
0045       call getenv('CMSSW_BASE',cmsdir)
0046
0047 C
0048       IF(MSTP(198).EQ.1) THEN
0049          hcs_name = cmsdir(1:index(cmsdir,' ')-1)
0050      $            //'/src/GeneratorInterface/GenExtensions/data/'
0051      $            //'dJYs05rc1.dat'
0052       ELSEIF(MSTP(198).EQ.2) THEN
0053          hcs_name = cmsdir(1:index(cmsdir,' ')-1)
0054      $            //'/src/GeneratorInterface/GenExtensions/data/'
0055      $            //'dJYs1rc1.dat'
0056       ELSEIF(MSTP(198).EQ.3) THEN
0057          hcs_name = cmsdir(1:index(cmsdir,' ')-1)
0058      $            //'/src/GeneratorInterface/GenExtensions/data/'
0059      $            //'dJYs2rc1.dat'
0060       ELSE
0061          WRITE(*,*)
0062          WRITE(*,*) 'ERROR: MSTP(198) must be 1,2 or 3!'
0063          WRITE(*,*) 'STOPPING EXECUTION...'
0064          WRITE(*,*)
0065          STOP
0066       ENDIF
0067
0068       WRITE(*,*) 'read_hcs_file: chose ',hcs_name
0069
0070       OPEN( UNIT = 5, FILE = hcs_name, STATUS = 'OLD')
0071 C      REWIND(5)
0072       READ(5,*)
0073
0074       DO 10, iq = 1, 12
0075         DO 20, iy = 1,25
0076           READ(5,*) hcstab(iq,iy,2),hcstab(iq,iy,1),x1,x2,
0077      >              hcstab(iq,iy,3)
0078  20     CONTINUE
0079  10   CONTINUE
0080
0081       CLOSE(5)
0082       RETURN
0083
0084       END
0085
0086
0087 C
0088 C    returns the  value based on direct quadratic interpolation
0089 C    table of hcs XS values (multiplied by some factors) for
0090 C    0<Y<11.5, 5<Q<144
0091 C
0092       FUNCTION HCS_XS(y,q)
0093       IMPLICIT NONE
0094
0095       REAL*8  HCS_XS
0096       REAL*8  y,q, controll
0097
0098       REAL*8 Pi, CONV
0099
0100       PARAMETER ( Pi = 3.141592654D0 )
0101       PARAMETER ( CONV = 0.3894D6 )
0102
0103       REAL*8 hcstab(12,25,3)
0104       COMMON / hcstab / hcstab
0105
0106       INTEGER NQ,NY
0107       PARAMETER ( NY = 25 )
0108       PARAMETER ( NQ = 12 )
0109
0110       REAL*8 YMAX, YSTEP, QRAT, QMIN
0111       PARAMETER ( YMAX  = 12.D0 )
0112       PARAMETER ( YSTEP = 5.D-1 )
0113       PARAMETER ( QMIN =  5.D0  )
0114       PARAMETER ( QRAT =  1.4D0 )
0115
0116       INTEGER iq,iy
0117
0118       REAL*8 q1,q2,q3,y1,y2,y3,z1,zq2,zq3,zy2,zy3,zqy
0119       REAL*8 aq,ay,aqy,bq,by,c, t2
0120
0121 C      WRITE(*,*)
0122 C      WRITE(*,*) 'y=',y
0123 C      WRITE(*,*) 'q=',q
0124
0125       IF(y.GT.11.5D0.OR.q.LT.5D0.OR.q.GT.144D0) THEN
0126         HCS_XS=0D0
0127 C         WRITE(*,*)
0128 C         WRITE(*,*) 'WARNING! Outside of kinematical region where'
0129 C         WRITE(*,*) 'calculation of d(sigma-hat)/d(t-hat) is valid!'
0130 C         WRITE(*,*) 'y=',y
0131 C         WRITE(*,*) 'q=',q
0132 C         WRITE(*,*)
0133         GOTO 30
0134       ENDIF
0135
0136       iy = IDINT(y/YSTEP)+1
0137       iq = IDINT(DLOG(q/QMIN) / DLOG(QRAT)) + 1
0138
0139       q1 = hcstab(iq,iy,1)
0140       q2 = hcstab(iq+1,iy,1)
0141       q3 = hcstab(iq+2,iy,1)
0142
0143       y1 = hcstab(iq,iy,2)
0144       y2 = hcstab(iq,iy+1,2)
0145       y3 = hcstab(iq,iy+2,2)
0146
0147       z1 =  hcstab(iq,iy,3)
0148       zq2 = hcstab(iq+1,iy,3)
0149       zq3 = hcstab(iq+2,iy,3)
0150       zy2 = hcstab(iq,iy+1,3)
0151       zy3 = hcstab(iq,iy+2,3)
0152       zqy = hcstab(iq+1,iy+1,3)
0153
0154
0155 C  010619: aq and ay changed according to Leszek
0156       aq = ((q3-q1)*(zq2-z1)-(q2-q1)*(zq3-z1))/
0157      >     ((q2-q1)*(q3-q1)*(q2-q3))
0158
0159       ay =  ((y3-y1)*(zy2-z1)-(y2-y1)*(zy3-z1))/
0160      >     ((y2-y1)*(y3-y1)*(y2-y3))
0161
0162       aqy = (z1 + zqy - zq2 - zy2)/(q2*y2+q1*y1-q1*y2-q2*y1)
0163
0164       bq = (zq2 - z1)/(q2-q1) - aq*(q1+q2) - aqy*y1
0165       by = (zy2 - z1)/(y2-y1) - ay*(y1+y2) - aqy*q1
0166
0167       c = zqy - aq*q2*q2 - bq*q2 - ay*y2*y2 - by*y2 - aqy*y2*q2
0168
0169       t2 = Pi*32.D0/9.D0 / q**4
0170       t2 = t2* (12.D0*Pi / 25.D0 / DLOG(q**2/0.04D0))**2
0171
0172       controll = aq*q**2+ay*y**2+aqy*q*y+bq*q+by*y+c
0173
0174       HCS_XS = t2*(aq*q**2+ay*y**2+aqy*q*y+bq*q+by*y+c)
0175
0176 C      WRITE(6,*), q1, q2, q3
0177 C      WRITE(6,*), y1, y2, y3
0178 C      WRITE(6,*), z1,zq2,zy2,zqy
0179 C      WRITE(6,*), controll
0180
0181  30   RETURN
0182       END
0183
0184
0185