Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:54

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