pythiaeRHIC/erhic/pythia_radgen_extras.f

0001 C A collection of routines needed by RADGEN
0002
0003 !---------------------------------------------------------------
0004 ! Calculate the F2 structure function
0005
0006       subroutine mkf2 (DQ2,DX,A,Z,DF2,DF1)
0007
0008       implicit none
0009       include "py6strf.inc"
0010       include "mc_set.inc"
0011
0012         double precision dQ2, dx, dF2, dF1, DF2NF2P
0013         DOUBLE PRECISION F2ALLM, DNP(5), DR, XPQ
0014         DOUBLE PRECISION gamma2, dnu, w2, pmass2
0015         integer A, Z
0016         DIMENSION XPQ(-25:25)
0017
0018 C parameters for ratio F2(n)/F2(p)
0019 C measured at NMC, (NMC Amaudruz et al. CERN-PPE/91-167))
0020         data DNP   / 0.976D0,    -1.34D0,      1.319D0,
0021      &              -2.133D0,     1.533D0/
0022
0023         pmass2=massp**2
0024         w2=pmass2+(dq2*(1/dx-1))
0025         dnu=(w2-pmass2+dq2)/(2.*massp)
0026         gamma2=dq2/(dnu**2)
0027
0028         if ((A.eq.1).and.(Z.eq.1)) then
0029 *:        ALLM:
0030          IF(genSet_FStruct(1:4).EQ.'ALLM') THEN
0031             Call MKR(DQ2,DX,DR)
0032             DF2=F2ALLM(dx,dq2)
0033             DF1=(1.D0+gamma2)/(2.D0*dx)/(1.D0+DR)*DF2
0034          ELSEIF(genSet_FStruct(1:4).EQ.'F2PY') THEN
0035             call F2PYTH(dx,dq2,df1,df2,z)
0036          ELSEIF(genSet_FStruct(1:4).EQ.'PDF') THEN
0037             call PYPDFU(2212,dX,dQ2,XPQ)
0038             DF2=1D0/9D0*(XPQ(1)+XPQ(-1)+XPQ(3)+XPQ(-3))+
0039      &          4D0/9D0*(XPQ(2)+XPQ(-2))
0040             DF1=(1.D0+gamma2)/(2.D0*dx)*DF2
0041          ELSE
0042 *:           error:
0043             write(*,*)('invalid parametrisation choice in mkf2')
0044          ENDIF
0045       ELSEIF(A.EQ.2.and.z.eq.1)THEN
0046 *:        ALLM:
0047          IF(genSet_FStruct(1:4).EQ.'ALLM') THEN
0048             Call MKR(DQ2,DX,DR)
0049             DF2=F2ALLM(dx,dq2)
0050             DF2NF2P=DNP(1)+dx*(DNP(2)+dx*(DNP(3)+dx*(DNP(4)+dx*DNP(5))))
0051             DF2=DF2*0.5*(df2nf2p+1.)
0052             DF1=(1.D0+gamma2)/(2.D0*dx)/(1.D0+DR)*DF2
0053          ELSE
0054 *:           error
0055             write(*,*)('MKF2: invalid parametrisation choice FStruct')
0056          ENDIF
0057 *
0058 ! ... neutron = 2*(deuterium_per_nucleon) - proton:
0059       ELSEIF(A.EQ.1.and.z.eq.0)THEN
0060 *:        ALLM:
0061          IF(genSet_FStruct(1:4).EQ.'ALLM') THEN
0062             DF2=F2ALLM(dx,dq2)
0063             Call MKR(DQ2,DX,DR)
0064             DF2NF2P=DNP(1)+dx*(DNP(2)+dx*(DNP(3)+dx*(DNP(4)+dx*DNP(5))))
0065             DF2=DF2*df2nf2p
0066             DF1=(1.D0+gamma2)/(2.D0*dx)/(1.D0+DR)*DF2
0067          ELSEIF(genSet_FStruct(1:4).EQ.'F2PY') THEN
0068             call F2PYTH(dx,dq2,df1,df2,z)
0069          ELSEIF(genSet_FStruct(1:4).EQ.'PDF') THEN
0070             call PYPDFU(2112,dX,dQ2,XPQ)
0071             DF2=1D0/9D0*(XPQ(1)+XPQ(-1)+XPQ(3)+XPQ(-3))+
0072      &          4D0/9D0*(XPQ(2)+XPQ(-2))
0073             DF1=(1.D0+gamma2)/(2.D0*dx)*DF2
0074          ELSE
0075 *:           error
0076             write(*,*)('MKF2: invalid parametrisation choice FStruct')
0077          ENDIF
0078       ELSEIF((A.gt.1).and.(Z.gt.1)) then
0079          IF(genSet_FStruct(1:4).EQ.'PDF') THEN
0080             call PYPDFU(2212,dX,dQ2,XPQ)
0081             DF2=1D0/9D0*(XPQ(1)+XPQ(-1)+XPQ(3)+XPQ(-3))+
0082      &          4D0/9D0*(XPQ(2)+XPQ(-2))
0083             DF1=(1.D0+gamma2)/(2.D0*dx)*DF2
0084          ELSE
0085 *:           error:
0086             write(*,*)('invalid parametrisation choice in mkf2')
0087          ENDIF
0088       ELSE
0089 *:        error:
0090          write(*,*)('MKF2: invalid target type'), A, Z
0091       ENDIF
0092
0093        return
0094        end
0095
0096 !---------------------------------------------------------------
0097 ! Calculate R = sigma_L/sigma_T
0098
0099       SUBROUTINE MKR(DQ2,DX,DR)
0100       IMPLICIT NONE
0101       include "mc_set.inc"
0102       include "py6strf.inc"
0103
0104       DOUBLE PRECISION DQ2, DX
0105       DOUBLE PRECISION DR,DELTAR
0106
0107       IF ( genSet_R .EQ. '1990' ) THEN
0108 *        Whitlow et al.,  Phys.Lett.B 250(1990),193
0109          CALL R1990(DQ2,DX,DR,DELTAR)
0110          py6R=DR
0111       ELSE IF ( genSet_R .EQ. '1998' ) THEN
0112 *        E143, hep-ex/9808028
0113          CALL R1998(DQ2,DX,DR,DELTAR)
0114          py6R=DR
0115       ELSE IF ( genSet_R .eq. '0' ) THEN
0116 * pure transverse (sigma_L=0)
0117          DR = 0.d0
0118          py6R=0.d0
0119       ELSE
0120          write(*,*)( 'MKR: invalid choice for R parametrization' )
0121       ENDIF
0122
0123       RETURN
0124       END
0125
0126 C------------------------------------------------------------------
0127
0128       SUBROUTINE R1990(DQ2,DX,DR,DELTAR)
0129
0130       IMPLICIT NONE
0131
0132       DOUBLE PRECISION DQ2, DX
0133       DOUBLE PRECISION DR, DELTAR
0134
0135       REAL R
0136       REAL QQ35, XX
0137       REAL FAC, RLOG, Q2THR
0138       REAL R_A, R_B, R_C
0139 C
0140 C Data-Definition of R-Calculation, see
0141 C            L.W.WHITLOW, SLAC-REPORT-357,
0142 C            PH.D. THESIS, STANFORD UNIVERSITY,
0143 C            MARCH 1990.
0144       REAL AR1990(3), BR1990(3), CR1990(3)
0145       DATA AR1990  / .06723, .46714, 1.89794 /
0146       DATA BR1990  / .06347, .57468, -.35342 /
0147       DATA CR1990  / .05992, .50885, 2.10807 /
0148
0149       DELTAR = 0.
0150
0151       XX=real(DX)
0152       IF (DQ2.LT.0.35) THEN
0153         QQ35=0.35
0154       ELSE
0155         QQ35=real(DQ2)
0156       ENDIF
0157 C
0158 C *** If Q2 < 0.35 then variable "R" is calculated at the fixed Q2 of 0.35
0159 C
0160       FAC   = 1+12.*(QQ35/(1.+QQ35))*(.125**2/(XX**2+.125**2))
0161       RLOG  = FAC/LOG(QQ35/.04)
0162       Q2THR = 5.*(1.-XX)**5
0163
0164       R_A   = AR1990(1)*RLOG +
0165      &        AR1990(2)/SQRT(SQRT(QQ35**4+AR1990(3)**4))
0166       R_B   = BR1990(1)*RLOG +
0167      &        BR1990(2)/QQ35 + BR1990(3)/(QQ35**2+.3**2)
0168       R_C   = CR1990(1)*RLOG +
0169      &        CR1990(2)/SQRT((QQ35-Q2THR)**2+CR1990(3)**2)
0170       R     = (R_A+R_B+R_C)/3.
0171
0172       IF (DQ2.GE.0.35) THEN
0173         DR=dble(R)
0174       ELSE
0175         DR=dble(R)*DQ2/0.35
0176       ENDIF
0177
0178 c      print*,'R:',R
0179
0180       END
0181
0182 C-----------------------------------------------------------------------
0183       SUBROUTINE R1998(DQ2,DX,DR,DELTAR)
0184
0185 C new fit to R  hep-ex/9808028 E143 Collab.
0186 C it is based on the old 3 paramter forms
0187 C 0.005<x<0.86, 0.5<Q2<130 GeV2
0188 C E143 x-section measurement 0.03<x<0.4
0189 C with  overall norm uncertainty 2.5%
0190 C
0191 C U. Stoesslein, October 1998
0192 C
0193
0194       IMPLICIT NONE
0195
0196       DOUBLE PRECISION DQ2,DX,DR,DELTAR
0197       DOUBLE PRECISION Q2,Q2MAX,FAC,RLOG,Q2THR
0198       DOUBLE PRECISION R_A_NEW,R_A,R_B_NEW,R_B,R_C
0199
0200       DOUBLE PRECISION A(6),B(6),C(6)
0201
0202       DATA A/ .0485, 0.5470, 2.0621, -.3804,  0.5090, -.0285/
0203       DATA B/ .0481, 0.6114, -.3509, -.4611,  0.7172, -.0317/
0204       DATA C/ .0577, 0.4644, 1.8288,12.3708,-43.1043,41.7415/
0205
0206       DOUBLE PRECISION DR1998
0207       EXTERNAL DR1998
0208
0209 * use R(0.35) if q2 is below 0.35
0210       Q2=DQ2
0211       Q2MAX=0.35
0212       IF(Q2.LT.Q2MAX) Q2=Q2MAX
0213
0214       FAC   = 1+12.*(Q2/(1.+Q2))*(.125**2/(DX**2+.125**2))
0215       RLOG  = FAC/LOG(Q2/.04)
0216       Q2thr = C(4)*DX+C(5)*DX*DX+C(6)*DX*DX*DX
0217
0218 * new additional terms
0219       R_A_NEW = (1.+A(4)*DX+A(5)*DX*DX)*DX**(A(6))
0220       R_A   = A(1)*RLOG + A(2)/SQRT(SQRT(Q2**4+A(3)**4))*R_A_NEW
0221       R_B_NEW = (1.+B(4)*DX+B(5)*DX*DX)*DX**(B(6))
0222       R_B   = B(1)*RLOG + (B(2)/Q2 + B(3)/(Q2**2+0.3**2))*R_B_NEW
0223       R_C   = C(1)*RLOG + C(2)/SQRT((Q2-Q2thr)**2+C(3)**2)
0224       DR    = (R_A+R_B+R_C)/3.
0225
0226 * straight line fit extrapolation to R(Q2=0)=0
0227       if (dq2.lt.q2max) DR = DR*DQ2/Q2MAX
0228
0229 * I assume the fit uncertainty only for measured Q2 range
0230       if (Q2 .GT. 0.5) then
0231          DELTAR = DR1998(DX,Q2)
0232       else
0233          DELTAR=DR
0234       endif
0235
0236       RETURN
0237       END
0238
0239 C--------------------------------------------------------------------
0240       DOUBLE PRECISION FUNCTION DR1998(DX,DQ2)
0241
0242 * Parameterize uncertainty in R1998
0243 * associated to the fitting procedure only
0244
0245       IMPLICIT NONE
0246       DOUBLE PRECISION DX,DQ2
0247
0248       DR1998 = 0.0078-0.013*DX+(0.070-0.39*DX+0.7*DX*DX)/(1.7+DQ2)
0249
0250       RETURN
0251       END
0252
0253 !---------------------------------------------------------------
0254 ! Calculate the asymmetries A1 and A2. Routine is empty because Pythia
0255 ! is unpolarised, but radgen expects it
0256
0257        subroutine mkasym (dQ2, dX, A, Z, dA1, dA2)
0258
0259        implicit none
0260
0261        double precision dQ2, dx, dA1, dA2
0262        integer A, Z
0263
0264         da1 = 0.D0
0265         da2 = 0.D0
0266
0267        return
0268        end
0269
0270 !---------------------------------------------------------------
0271 ! Calculate the dilution factor
0272
0273        double precision function fdilut (dQ2, dx, A)
0274
0275        implicit none
0276
0277        DOUBLE PRECISION DQ2, DX, DF
0278        DOUBLE PRECISION DNP, DFN, DFP
0279        DOUBLE PRECISION DZ, DF2NF2P
0280        INTEGER A
0281        dimension dnp(7)
0282 *
0283 C ... fit to NMC  F2n/F2p data (86/87+89 T1,T14)
0284        data dnp/
0285      +   0.67225D+00,
0286      +   0.16254D+01,
0287      +  -0.15436D+00,
0288      +   0.48301D-01,
0289      +        0.41979D+00,
0290      +   0.47331D-01,
0291      +  -0.17816D+00/
0292
0293 ! Definitions of 'intrinsic' dilutions for neutron and proton (GNOME confused)
0294
0295        dfn=1.
0296        dfp=1.
0297
0298 ! Only 3He has a dilution, and we determine it as F2n/(2*F2p + F2n)
0299
0300        if (A.ne.3) then
0301          df = 1.
0302        else
0303          dZ = 1./2.*DLOG(1.+DEXP(2.0-1000.*dX))
0304          df2nf2p = dnp(1)*(1.0-dx)**dnp(2)+dnp(3)*dx**dnp(4)
0305      1            +(dnp(5)+dnp(6)*dz+dnp(7)*dz**2)
0306          df = dfn*(1./((2./df2nf2p)+1))
0307        endif
0308
0309        FDILUT = DF
0310        return
0311
0312        END
0313
0314 !---------------------------------------------------------------
0315 ! Function to calculate F2 from ALLM parametrisation
0316
0317       double precision FUNCTION f2allm(x,q2)
0318
0319       implicit double precision (a-h,o-z)
0320
0321       REAL M02,M12,LAM2,M22
0322       COMMON/ALLM/SP,AP,BP,SR,AR,BR,S,XP,XR,F2P,F2R
0323 C  POMERON
0324       PARAMETER (
0325      , S11   =   0.28067, S12   =   0.22291, S13   =   2.1979,
0326      , A11   =  -0.0808 , A12   =  -0.44812, A13   =   1.1709,
0327      , B11   =   0.60243**2, B12   =   1.3754**2, B13   =   1.8439,
0328      , M12   =  49.457 )
0329
0330 C  REGGEON
0331       PARAMETER (
0332      , S21   =   0.80107, S22   =   0.97307, S23   =   3.4942,
0333      , A21   =   0.58400, A22   =   0.37888, A23   =   2.6063,
0334      , B21   =   0.10711**2, B22   =   1.9386**2, B23   =   0.49338,
0335      , M22   =   0.15052 )
0336 C
0337       PARAMETER ( M02=0.31985, LAM2=0.065270, Q02=0.46017 +LAM2 )
0338       PARAMETER ( ALFA=112.2, XMP2=0.8802)
0339 C
0340       W2=q2*(1./x -1.)+xmp2
0341       W=dsqrt(w2)
0342 C
0343       IF(Q2.EQ.0.) THEN
0344        S=0.
0345        Z=1.
0346 C
0347 C   POMERON
0348 C
0349        XP=1./(1.+(W2-XMP2)/(Q2+M12))
0350        AP=A11
0351        BP=B11
0352        SP=S11
0353        F2P=SP*XP**AP
0354 C
0355 C   REGGEON
0356 C
0357        XR=1./(1.+(W2-XMP2)/(Q2+M22))
0358        AR=A21
0359        BR=B21
0360        SR=S21
0361        F2R=SR*XR**AR
0362 C
0363       ELSE
0364        S=DLOG(DLOG((Q2+Q02)/LAM2)/DLOG(Q02/LAM2))
0365        Z=1.-X
0366 C
0367 C   POMERON
0368 C
0369        XP=1./(1.+(W2-XMP2)/(Q2+M12))
0370        AP=A11+(A11-A12)*(1./(1.+S**A13)-1.)
0371        BP=B11+B12*S**B13
0372        SP=S11+(S11-S12)*(1./(1.+S**S13)-1.)
0373        F2P=SP*XP**AP*Z**BP
0374 C
0375 C   REGGEON
0376 C
0377        XR=1./(1.+(W2-XMP2)/(Q2+M22))
0378        AR=A21+A22*S**A23
0379        BR=B21+B22*S**B23
0380        SR=S21+S22*S**S23
0381        F2R=SR*XR**AR*Z**BR
0382
0383 C
0384       ENDIF
0385
0386 c      CIN=ALFA/(Q2+M02)*(1.+4.*XMP2*Q2/(Q2+W2-XMP2)**2)/Z
0387 c      SIGal=CIN*(F2P+F2R)
0388 c      f2allm=sigal/alfa*(q2**2*(1.-x))/(q2+4.*xmp2*x**2)
0389       f2allm = q2/(q2+m02)*(F2P+F2R)
0390
0391       RETURN
0392       END
0393 C--------------------------------------------------------------------
0394
0395 C.....Total inclusive Pythia cross section according to the model in
0396 C.....C.Friberg, T.Sjoestrand, J. High Energy Phys. JHEP 0009, 010 (2000)
0397 C.....(neglecting GVMD cross sections)
0398 C.....P.Liebing, 08/16/2003
0399
0400       SUBROUTINE F2PYTH(x,q2,f1,f2,z)
0401
0402       IMPLICIT NONE
0403       COMMON/PYINT1/MINT(400),VINT(400)
0404       INTEGER MINT
0405       DOUBLE PRECISION VINT
0406       SAVE/PYINT1/
0407       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0408       INTEGER MSTP,MSTI
0409       DOUBLE PRECISION PARP,PARI
0410       SAVE/PYPARS/
0411       DOUBLE PRECISION f2dis,f1dis,sdis,svmd1,svmd2,sigvm,x,q2,dipol
0412       DOUBLE PRECISION mrho2,alpha,eps,eta,pmass,pmass2,gevmb,rvmd
0413       DOUBLE PRECISION w2,gamma2,nu,convf2,convf1,conv,sigh,pi
0414       DOUBLE PRECISION f1,f2,df2allm,df2allml,df1allm,df1allml,r,rl
0415       DOUBLE PRECISION w2l,gamma2l,nul,xl,sw2
0416       INTEGER z
0417 C...Local arrays.
0418       DOUBLE PRECISION XPQ,XPAR,YPAR
0419       INTEGER i
0420       DIMENSION XPQ(-25:25),XPAR(4),YPAR(4)
0421
0422 C...X and Y parameters of sigmatot = X * s**epsilon + Y * s**(-eta).
0423       DATA XPAR/2*13.63D0,10.01D0,0.970D0/
0424       DATA YPAR/2*31.79D0,-1.51D0,-0.146D0/
0425       DATA EPS/0.0808D0/,ETA/-0.4525D0/
0426       DATA MRHO2/0.591822D0/,ALPHA/7.297352533D-3/,PMASS/0.93827D0/,
0427      &     GEVMB/0.389379292D0/,pi/3.14159265358979324D0/
0428
0429       EXTERNAL PYPDFU
0430       f2dis=0D0
0431       f1dis=0D0
0432       sdis=0D0
0433       svmd1=0D0
0434       svmd2=0D0
0435       sigvm=0D0
0436       rvmd=0D0
0437       dipol=0D0
0438       w2=0D0
0439       gamma2=0D0
0440       nu=0D0
0441       convf2=0D0
0442       convf1=0D0
0443       conv=0D0
0444       sigh=0D0
0445       f1=0D0
0446       f2=0D0
0447       df1allm=0D0
0448       df2allm=0D0
0449       df1allml=0D0
0450       df2allml=0D0
0451       r=0D0
0452       rl=0D0
0453       sw2=0D0
0454 c....Kinematic factors and constants
0455       pmass2=pmass**2
0456       nul=-1D0
0457       gamma2l=0D0
0458       w2l=-1D0
0459       xl=0D0
0460
0461       if ((x.gt.0D0).and.(x.le.1D0)) then
0462        w2=pmass2+(q2*(1./x-1.))
0463       else
0464        f1=0D0
0465        f2=0D0
0466        return
0467       endif
0468
0469       if (w2.lt.4D0) then
0470        w2l=w2
0471        w2=4D0
0472        xl=x
0473        nul=(w2l-pmass2+q2)/(2D0*pmass)
0474       endif
0475
0476       nu=(w2-pmass2+q2)/(2D0*pmass)
0477
0478       if (nu.gt.0D0) then
0479        gamma2=q2/(nu**2)
0480        if (nul.gt.0D0) then
0481         gamma2l=q2/(nul**2)
0482        endif
0483       else
0484        f1=0D0
0485        f2=0D0
0486        return
0487       endif
0488
0489 c....For W2<4, we don't trust the PYTHIA F2, so we calculate
0490 C....F2-ALLM(W2,Q2). The real kinematics have and "l" at the end,
0491 C....the kinematics without the "l" are the ones we get by setting W2 to 4
0492 c....Output: f2allm=F2-ALLM(W2=4,Q2),f2allml=F2-ALLM(W2=w2l,Q2)
0493       if (w2l.gt.0D0) then
0494        sw2=((w2l-pmass2)/(4D0-pmass2))**10
0495       endif
0496 c....This factor is needed to convert the Pythia virtual gamma cross
0497 c....section for VMD to the same level as F2
0498 c....The kinematic factors making the (ep) cross section out of F2 are
0499 c....provided by RADGEN
0500       conv=q2*(1D0-x)/(4D0*pi**2*alpha)/gevmb
0501 c....Pythia PDF call, sum PDFs to F2
0502       if (z.eq.1) then
0503       call PYPDFU(2212,X,Q2,XPQ)
0504       elseif (z.eq.0) then
0505       call PYPDFU(2112,X,Q2,XPQ)
0506       endif
0507        f2dis=1D0/9D0*(XPQ(1)+XPQ(-1)+XPQ(3)+XPQ(-3))+
0508      &       4D0/9D0*(XPQ(2)+XPQ(-2))
0509 c....Suppression factor for DIS
0510       if (MSTP(19).eq.0) then
0511        sdis=1.
0512       else
0513        sdis=q2/(q2+mrho2)
0514        if (MSTP(19).gt.1) then
0515         sdis=sdis**2
0516        endif
0517       endif
0518 C....Sum of Hadronic (Vector Meson) cross sections * Photon couplings
0519 C....const.
0520       sigh=0.
0521       do 10 i=1,4
0522        sigh=sigh+alpha/PARP(160+i)*(XPAR(i)*w2**eps+YPAR(i)*w2**eta)
0523    10 continue
0524 C....W2/Q2 suppression of VMD and (1+epsilon R_VMD)
0525       svmd1=(w2/(w2+q2))**MSTP(20)
0526       if (MSTP(20).eq.0) then
0527        dipol=2.575D0
0528       else
0529        dipol=2D0
0530       endif
0531       if (MSTP(17).eq.6) then
0532        rvmd=PARP(165)*(q2/mrho2)**PARP(166)
0533       else
0534 c    ...Attention: This is only good for MSTP(17)=4, i.e., the Pythia
0535 c    ...default
0536        rvmd=PARP(165)*(4.*mrho2*q2)/(mrho2+q2)**2
0537       endif
0538 C  .... Dipole factor for VMD
0539       svmd2=(mrho2/(mrho2+q2))**dipol
0540 C.....virtual photon xsec for VMD
0541       sigvm=svmd1*svmd2*sigh
0542       convf2=(1D0+rvmd)/(1D0+gamma2)
0543       convf1=1D0/(2D0*x)
0544 c.....Total "F2"
0545       f2=sdis*f2dis+conv*convf2*sigvm
0546       f1dis=(1.D0+gamma2)/(2.D0*x)*f2dis
0547       f1=sdis*f1dis+conv*convf1*sigvm
0548       if (w2l.gt.0D0) then
0549 C.....Here we scale F2-ALLM(W2=w2l,Q2) by the factor
0550 C.....F2-PYTH(W2=4,Q2)/F2-ALLM(W2=4,Q2) (normalize ALLM to PYTHIA model)
0551         f2=sw2*f2
0552         f1=sw2*f1
0553       endif
0554       RETURN
0555       END