pythiaeRHIC/pythia/pypdfu.f

0001
0002 C*********************************************************************
0003
0004 C...PYPDFU
0005 C...Gives electron, muon, tau, photon, pi+, neutron, proton and hyperon
0006 C...parton distributions according to a few different parametrizations.
0007 C...Note that what is coded is x times the probability distribution,
0008 C...i.e. xq(x,Q2) etc.
0009
0010       SUBROUTINE PYPDFU(KF,X,Q2,XPQ)
0011
0012 C...Double precision and integer declarations.
0013       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0014       IMPLICIT INTEGER(I-N)
0015       INTEGER PYK,PYCHGE,PYCOMP
0016 C...Commonblocks.
0017       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0018       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0019       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0020       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0021       COMMON/PYINT1/MINT(400),VINT(400)
0022       COMMON/PYINT8/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
0023      &XPDIR(-6:6)
0024       COMMON/PYINT9/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
0025       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0026      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0027      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0028       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/,/PYINT8/,
0029      &/PYINT9/,/PYINTM/
0030 C...Local arrays.
0031       DIMENSION XPQ(-25:25),XPEL(-25:25),XPGA(-6:6),VXPGA(-6:6),
0032      &XPPI(-6:6),XPPR(-6:6),XPVAL(-6:6),PPAR(6,2)
0033       SAVE PPAR
0034
0035 C...Interface to PDFLIB.
0036       COMMON/W50513/XMIN,XMAX,Q2MIN,Q2MAX
0037       SAVE /W50513/
0038       DOUBLE PRECISION XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU,
0039      &VALUE(20),XMIN,XMAX,Q2MIN,Q2MAX
0040       CHARACTER*20 PARM(20)
0041       DATA VALUE/20*0D0/,PARM/20*' '/
0042
0043 C...Added by liang 1/6/12
0044 C...Switches for nuclear PDF correction
0045       COMMON/PYNUCL/INUMOD,CHANUM,ORDER
0046       SAVE /PYNUCL/
0047       DOUBLE PRECISION INUMOD,CHANUM
0048       INTEGER ORDER
0049       CHARACTER*20 chpset
0050       CHARACTER*20 Nprm
0051
0052
0053 C...Data related to Schuler-Sjostrand photon distributions.
0054       DATA ALAMGA/0.2D0/, PMCGA/1.3D0/, PMBGA/4.6D0/
0055
0056 C...Valence PDF momentum integral parametrizations PER PARTON!
0057       DATA (PPAR(1,IPAR),IPAR=1,2) /0.385D0,1.60D0/
0058       DATA (PPAR(2,IPAR),IPAR=1,2) /0.480D0,1.56D0/
0059       PAVG(IFL,Q2)=PPAR(IFL,1)/(1D0+PPAR(IFL,2)*
0060      &LOG(LOG(MAX(Q2,1D0)/0.04D0)))
0061
0062 C...Reset parton distributions.
0063       MINT(92)=0
0064       DO 100 KFL=-25,25
0065         XPQ(KFL)=0D0
0066   100 CONTINUE
0067       DO 110 KFL=-6,6
0068         XPVAL(KFL)=0D0
0069   110 CONTINUE
0070
0071 C...Check x and particle species.
0072       IF(X.LE.0D0.OR.X.GE.1D0) THEN
0073         WRITE(MSTU(11),5000) X
0074         GOTO 9999
0075       ENDIF
0076       KFA=IABS(KF)
0077       IF(KFA.NE.11.AND.KFA.NE.13.AND.KFA.NE.15.AND.KFA.NE.22.AND.
0078      &KFA.NE.211.AND.KFA.NE.2112.AND.KFA.NE.2212.AND.KFA.NE.3122.AND.
0079      &KFA.NE.3112.AND.KFA.NE.3212.AND.KFA.NE.3222.AND.KFA.NE.3312.AND.
0080      &KFA.NE.3322.AND.KFA.NE.3334.AND.KFA.NE.111.AND.KFA.NE.321.AND.
0081      &KFA.NE.310.AND.KFA.NE.130) THEN
0082         WRITE(MSTU(11),5100) KF
0083         GOTO 9999
0084       ENDIF
0085
0086 C...Electron (or muon or tau) parton distribution call.
0087       IF(KFA.EQ.11.OR.KFA.EQ.13.OR.KFA.EQ.15) THEN
0088         CALL PYPDEL(KFA,X,Q2,XPEL)
0089         DO 120 KFL=-25,25
0090           XPQ(KFL)=XPEL(KFL)
0091   120   CONTINUE
0092
0093 C...Photon parton distribution call (VDM+anomalous).
0094       ELSEIF(KFA.EQ.22.AND.MINT(109).LE.1) THEN
0095         IF(MSTP(56).EQ.1.AND.MSTP(55).EQ.1) THEN
0096           CALL PYPDGA(X,Q2,XPGA)
0097           DO 130 KFL=-6,6
0098             XPQ(KFL)=XPGA(KFL)
0099   130     CONTINUE
0100           XPVU=4D0*(XPQ(2)-XPQ(1))/3D0
0101           XPVAL(1)=XPVU/4D0
0102           XPVAL(2)=XPVU
0103           XPVAL(3)=MIN(XPQ(3),XPVU/4D0)
0104           XPVAL(4)=MIN(XPQ(4),XPVU)
0105           XPVAL(5)=MIN(XPQ(5),XPVU/4D0)
0106           XPVAL(-1)=XPVAL(1)
0107           XPVAL(-2)=XPVAL(2)
0108           XPVAL(-3)=XPVAL(3)
0109           XPVAL(-4)=XPVAL(4)
0110           XPVAL(-5)=XPVAL(5)
0111         ELSEIF(MSTP(56).EQ.1.AND.MSTP(55).GE.5.AND.MSTP(55).LE.8) THEN
0112           Q2MX=Q2
0113           P2MX=0.36D0
0114           IF(MSTP(55).GE.7) P2MX=4.0D0
0115           IF(MSTP(57).EQ.0) Q2MX=P2MX
0116           P2=0D0
0117           IF(VINT(120).LT.0D0) P2=VINT(120)**2
0118           CALL PYGGAM(MSTP(55)-4,X,Q2MX,P2,MSTP(60),F2GAM,XPGA)
0119           DO 140 KFL=-6,6
0120             XPQ(KFL)=XPGA(KFL)
0121             XPVAL(KFL)=VXPDGM(KFL)
0122   140     CONTINUE
0123           VINT(231)=P2MX
0124         ELSEIF(MSTP(56).EQ.1.AND.MSTP(55).GE.9.AND.MSTP(55).LE.12) THEN
0125           Q2MX=Q2
0126           P2MX=0.36D0
0127           IF(MSTP(55).GE.11) P2MX=4.0D0
0128           IF(MSTP(57).EQ.0) Q2MX=P2MX
0129           P2=0D0
0130           IF(VINT(120).LT.0D0) P2=VINT(120)**2
0131           CALL PYGGAM(MSTP(55)-8,X,Q2MX,P2,MSTP(60),F2GAM,XPGA)
0132           DO 150 KFL=-6,6
0133             XPQ(KFL)=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
0134             XPVAL(KFL)=VXPVMD(KFL)+VXPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
0135   150     CONTINUE
0136           VINT(231)=P2MX
0137         ELSEIF(MSTP(56).EQ.2) THEN
0138 C...Call PDFLIB parton distributions.
0139           PARM(1)='NPTYPE'
0140           VALUE(1)=3
0141           PARM(2)='NGROUP'
0142           VALUE(2)=MSTP(55)/1000
0143           PARM(3)='NSET'
0144           VALUE(3)=MOD(MSTP(55),1000)
0145           IF(MINT(93).NE.3000000+MSTP(55)) THEN
0146             CALL PDFSET(PARM,VALUE)
0147             MINT(93)=3000000+MSTP(55)
0148           ENDIF
0149           XX=X
0150           QQ2=MAX(0D0,Q2MIN,Q2)
0151           IF(MSTP(57).EQ.0) QQ2=Q2MIN
0152           P2=0D0
0153           IF(VINT(120).LT.0D0) P2=VINT(120)**2
0154           IP2=MSTP(60)
0155           IF(MSTP(55).EQ.5004) THEN
0156             IF(5D0*P2.LT.QQ2.AND.
0157      &      QQ2.GT.0.6D0.AND.QQ2.LT.5D4.AND.
0158      &      P2.GE.0D0.AND.P2.LT.10D0.AND.
0159      &      XX.GT.1D-4.AND.XX.LT.1D0) THEN
0160               CALL STRUCTP(XX,QQ2,P2,IP2,UPV,DNV,USEA,DSEA,STR,CHM,
0161      &        BOT,TOP,GLU)
0162             ELSE
0163               UPV=0D0
0164               DNV=0D0
0165               USEA=0D0
0166               DSEA=0D0
0167               STR=0D0
0168               CHM=0D0
0169               BOT=0D0
0170               TOP=0D0
0171               GLU=0D0
0172             ENDIF
0173           ELSE
0174             IF(P2.LT.QQ2) THEN
0175               CALL STRUCTP(XX,QQ2,P2,IP2,UPV,DNV,USEA,DSEA,STR,CHM,
0176      &        BOT,TOP,GLU)
0177             ELSE
0178               UPV=0D0
0179               DNV=0D0
0180               USEA=0D0
0181               DSEA=0D0
0182               STR=0D0
0183               CHM=0D0
0184               BOT=0D0
0185               TOP=0D0
0186               GLU=0D0
0187             ENDIF
0188           ENDIF
0189           VINT(231)=Q2MIN
0190           XPQ(0)=GLU
0191           XPQ(1)=DNV
0192           XPQ(-1)=DNV
0193           XPQ(2)=UPV
0194           XPQ(-2)=UPV
0195           XPQ(3)=STR
0196           XPQ(-3)=STR
0197           XPQ(4)=CHM
0198           XPQ(-4)=CHM
0199           XPQ(5)=BOT
0200           XPQ(-5)=BOT
0201           XPQ(6)=TOP
0202           XPQ(-6)=TOP
0203           XPVU=4D0*(XPQ(2)-XPQ(1))/3D0
0204           XPVAL(1)=XPVU/4D0
0205           XPVAL(2)=XPVU
0206           XPVAL(3)=MIN(XPQ(3),XPVU/4D0)
0207           XPVAL(4)=MIN(XPQ(4),XPVU)
0208           XPVAL(5)=MIN(XPQ(5),XPVU/4D0)
0209           XPVAL(-1)=XPVAL(1)
0210           XPVAL(-2)=XPVAL(2)
0211           XPVAL(-3)=XPVAL(3)
0212           XPVAL(-4)=XPVAL(4)
0213           XPVAL(-5)=XPVAL(5)
0214         ELSE
0215           WRITE(MSTU(11),5200) KF,MSTP(56),MSTP(55)
0216         ENDIF
0217
0218 C...Pion/gammaVDM parton distribution call.
0219       ELSEIF(KFA.EQ.211.OR.KFA.EQ.111.OR.KFA.EQ.321.OR.KFA.EQ.130.OR.
0220      &KFA.EQ.310.OR.(KFA.EQ.22.AND.MINT(109).EQ.2)) THEN
0221         IF(KFA.EQ.22.AND.MSTP(56).EQ.1.AND.MSTP(55).GE.5.AND.
0222      &  MSTP(55).LE.12) THEN
0223           ISET=1+MOD(MSTP(55)-1,4)
0224           Q2MX=Q2
0225           P2MX=0.36D0
0226           IF(ISET.GE.3) P2MX=4.0D0
0227           IF(MSTP(57).EQ.0) Q2MX=P2MX
0228           P2=0D0
0229           IF(VINT(120).LT.0D0) P2=VINT(120)**2
0230           CALL PYGGAM(ISET,X,Q2MX,P2,MSTP(60),F2GAM,XPGA)
0231           DO 160 KFL=-6,6
0232             XPQ(KFL)=XPVMD(KFL)
0233             XPVAL(KFL)=VXPVMD(KFL)
0234   160     CONTINUE
0235           VINT(231)=P2MX
0236         ELSEIF(MSTP(54).EQ.1.AND.MSTP(53).GE.1.AND.MSTP(53).LE.3) THEN
0237           CALL PYPDPI(X,Q2,XPPI)
0238           DO 170 KFL=-6,6
0239             XPQ(KFL)=XPPI(KFL)
0240   170     CONTINUE
0241           XPVAL(2)=XPQ(2)-XPQ(-2)
0242           XPVAL(-1)=XPQ(-1)-XPQ(1)
0243         ELSEIF(MSTP(54).EQ.2) THEN
0244 C...Call PDFLIB parton distributions.
0245           PARM(1)='NPTYPE'
0246           VALUE(1)=2
0247           PARM(2)='NGROUP'
0248           VALUE(2)=MSTP(53)/1000
0249           PARM(3)='NSET'
0250           VALUE(3)=MOD(MSTP(53),1000)
0251           IF(MINT(93).NE.2000000+MSTP(53)) THEN
0252             CALL PDFSET(PARM,VALUE)
0253             MINT(93)=2000000+MSTP(53)
0254           ENDIF
0255           XX=X
0256           QQ=SQRT(MAX(0D0,Q2MIN,Q2))
0257           IF(MSTP(57).EQ.0) QQ=SQRT(Q2MIN)
0258           CALL STRUCTM(XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU)
0259           VINT(231)=Q2MIN
0260           XPQ(0)=GLU
0261           XPQ(1)=DSEA
0262           XPQ(-1)=UPV+DSEA
0263           XPQ(2)=UPV+USEA
0264           XPQ(-2)=USEA
0265           XPQ(3)=STR
0266           XPQ(-3)=STR
0267           XPQ(4)=CHM
0268           XPQ(-4)=CHM
0269           XPQ(5)=BOT
0270           XPQ(-5)=BOT
0271           XPQ(6)=TOP
0272           XPQ(-6)=TOP
0273           XPVAL(2)=UPV
0274           XPVAL(-1)=UPV
0275         ELSE
0276           WRITE(MSTU(11),5200) KF,MSTP(54),MSTP(53)
0277         ENDIF
0278
0279 C...Anomalous photon parton distribution call.
0280       ELSEIF(KFA.EQ.22.AND.MINT(109).EQ.3) THEN
0281         Q2MX=Q2
0282         P2MX=PARP(15)**2
0283         IF(MSTP(56).EQ.1.AND.MSTP(55).LE.8) THEN
0284           IF(MSTP(55).EQ.5.OR.MSTP(55).EQ.6) P2MX=0.36D0
0285           IF(MSTP(55).EQ.7.OR.MSTP(55).EQ.8) P2MX=4.0D0
0286           IF(MSTP(57).EQ.0) Q2MX=P2MX
0287           P2=0D0
0288           IF(VINT(120).LT.0D0) P2=VINT(120)**2
0289           CALL PYGGAM(MSTP(55)-4,X,Q2MX,P2,MSTP(60),F2GM,XPGA)
0290           DO 180 KFL=-6,6
0291             XPQ(KFL)=XPANL(KFL)+XPANH(KFL)
0292             XPVAL(KFL)=VXPANL(KFL)+VXPANH(KFL)
0293   180     CONTINUE
0294           VINT(231)=P2MX
0295         ELSEIF(MSTP(56).EQ.1) THEN
0296           IF(MSTP(55).EQ.9.OR.MSTP(55).EQ.10) P2MX=0.36D0
0297           IF(MSTP(55).EQ.11.OR.MSTP(55).EQ.12) P2MX=4.0D0
0298           IF(MSTP(57).EQ.0) Q2MX=P2MX
0299           P2=0D0
0300           IF(VINT(120).LT.0D0) P2=VINT(120)**2
0301           CALL PYGGAM(MSTP(55)-8,X,Q2MX,P2,MSTP(60),F2GM,XPGA)
0302           DO 190 KFL=-6,6
0303             XPQ(KFL)=MAX(0D0,XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL))
0304             XPVAL(KFL)=MAX(0D0,VXPANL(KFL)+XPBEH(KFL)+XPDIR(KFL))
0305   190     CONTINUE
0306           VINT(231)=P2MX
0307         ELSEIF(MSTP(56).EQ.2) THEN
0308           IF(MSTP(57).EQ.0) Q2MX=P2MX
0309           CALL PYGANO(0,X,Q2MX,P2MX,ALAMGA,XPGA,VXPGA)
0310           DO 200 KFL=-6,6
0311             XPQ(KFL)=XPGA(KFL)
0312             XPVAL(KFL)=VXPGA(KFL)
0313   200     CONTINUE
0314           VINT(231)=P2MX
0315         ELSEIF(MSTP(55).GE.1.AND.MSTP(55).LE.5) THEN
0316           IF(MSTP(57).EQ.0) Q2MX=P2MX
0317           CALL PYGVMD(0,MSTP(55),X,Q2MX,P2MX,PARP(1),XPGA,VXPGA)
0318           DO 210 KFL=-6,6
0319             XPQ(KFL)=XPGA(KFL)
0320             XPVAL(KFL)=VXPGA(KFL)
0321   210     CONTINUE
0322           VINT(231)=P2MX
0323         ELSE
0324   220     RKF=11D0*PYR(0)
0325           KFR=1
0326           IF(RKF.GT.1D0) KFR=2
0327           IF(RKF.GT.5D0) KFR=3
0328           IF(RKF.GT.6D0) KFR=4
0329           IF(RKF.GT.10D0) KFR=5
0330           IF(KFR.EQ.4.AND.Q2.LT.PMCGA**2) GOTO 220
0331           IF(KFR.EQ.5.AND.Q2.LT.PMBGA**2) GOTO 220
0332           IF(MSTP(57).EQ.0) Q2MX=P2MX
0333           CALL PYGVMD(0,KFR,X,Q2MX,P2MX,PARP(1),XPGA,VXPGA)
0334           DO 230 KFL=-6,6
0335             XPQ(KFL)=XPGA(KFL)
0336             XPVAL(KFL)=VXPGA(KFL)
0337   230     CONTINUE
0338           VINT(231)=P2MX
0339         ENDIF
0340
0341 C...Proton parton distribution call.
0342       ELSE
0343         IF(MSTP(52).EQ.1.AND.MSTP(51).GE.1.AND.MSTP(51).LE.20) THEN
0344           CALL PYPDPR(X,Q2,XPPR)
0345           DO 240 KFL=-6,6
0346             XPQ(KFL)=XPPR(KFL)
0347   240     CONTINUE
0348           XPVAL(1)=XPQ(1)-XPQ(-1)
0349           XPVAL(2)=XPQ(2)-XPQ(-2)
0350         ELSEIF(MSTP(52).EQ.2) THEN
0351 C...Call PDFLIB parton distributions.
0352           PARM(1)='NPTYPE'
0353           VALUE(1)=1
0354           PARM(2)='NGROUP'
0355           VALUE(2)=MSTP(51)/1000
0356           PARM(3)='NSET'
0357           VALUE(3)=MOD(MSTP(51),1000)
0358           IF(MINT(93).NE.1000000+MSTP(51)) THEN
0359 C...Uncommented by liang to link LHAPDF 1/6/12
0360             CALL PDFSET(PARM,VALUE)
0361             MINT(93)=1000000+MSTP(51)
0362           ENDIF
0363           XX=X
0364           QQ=SQRT(MAX(0D0,Q2MIN,Q2))
0365           IF(MSTP(57).EQ.0) QQ=SQRT(Q2MIN)
0366 C...Modified by liang to add nuclear correction 1/6/12*************
0367           IF(INUMOD.GT.1D0) THEN
0368             IF(ORDER/100.EQ.1) THEN
0369 C               print*,'Now run LO nuclear pdf'
0370                ipset=ORDER-100
0371                write(chpset,'(I2)') ipset
0372                Nprm='EPS09LO,'//chpset
0373                CALL SETLHAPARM(Nprm)
0374             ELSE IF(ORDER/100.EQ.2) THEN
0375 C               print*,'Now run NLO nuclear pdf'
0376                ipset=ORDER-200
0377                write(chpset,'(I2)') ipset
0378                Nprm='EPS09NLO,'//chpset
0379                CALL SETLHAPARM(Nprm)
0380             ENDIF
0381             CALL STRUCTA(XX,QQ,INUMOD,UPV,DNV,USEA,DSEA,STR,CHM,BOT,
0382      &      TOP,GLU)
0383 C...Isospin transformation for nucleus A with Z protons************
0384             A=INUMOD
0385             Z=CHANUM
0386             UPVp=UPV
0387             UPVn=DNV
0388             USEAp=USEA
0389             USEAn=DSEA
0390             DNVp=DNV
0391             DNVn=UPV
0392             DSEAp=DSEA
0393             DSEAn=USEA
0394             UPV=Z/A*UPVp+(A-Z)/A*UPVn
0395             USEA=Z/A*USEAp+(A-Z)/A*USEAn
0396             DNV=Z/A*DNVp+(A-Z)/A*DNVn
0397             DSEA=Z/A*DSEAp+(A-Z)/A*DSEAn
0398           ELSE
0399             CALL STRUCTM(XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU)
0400           ENDIF
0401 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0402           VINT(231)=Q2MIN
0403           XPQ(0)=GLU
0404           XPQ(1)=DNV+DSEA
0405           XPQ(-1)=DSEA
0406           XPQ(2)=UPV+USEA
0407           XPQ(-2)=USEA
0408           XPQ(3)=STR
0409           XPQ(-3)=STR
0410           XPQ(4)=CHM
0411           XPQ(-4)=CHM
0412           XPQ(5)=BOT
0413           XPQ(-5)=BOT
0414           XPQ(6)=TOP
0415           XPQ(-6)=TOP
0416           XPVAL(1)=DNV
0417           XPVAL(2)=UPV
0418         ELSE
0419           WRITE(MSTU(11),5200) KF,MSTP(52),MSTP(51)
0420         ENDIF
0421       ENDIF
0422
0423 C...Isospin average for pi0/gammaVDM.
0424       IF(KFA.EQ.111.OR.(KFA.EQ.22.AND.MINT(109).EQ.2)) THEN
0425         IF(KFA.EQ.22.AND.MSTP(55).GE.5.AND.MSTP(55).LE.12) THEN
0426           XPV=XPQ(2)-XPQ(1)
0427           XPQ(2)=XPQ(1)
0428           XPQ(-2)=XPQ(-1)
0429         ELSE
0430           XPS=0.5D0*(XPQ(1)+XPQ(-2))
0431           XPV=0.5D0*(XPQ(2)+XPQ(-1))-XPS
0432           XPQ(2)=XPS
0433           XPQ(-1)=XPS
0434         ENDIF
0435         XPVL=0.5D0*(XPVAL(1)+XPVAL(2)+XPVAL(-1)+XPVAL(-2))+
0436      &  XPVAL(3)+XPVAL(4)+XPVAL(5)
0437         DO 250 KFL=-6,6
0438           XPVAL(KFL)=0D0
0439   250   CONTINUE
0440         IF(KFA.EQ.22.AND.MINT(105).LE.223) THEN
0441           XPQ(1)=XPQ(1)+0.2D0*XPV
0442           XPQ(2)=XPQ(2)+0.8D0*XPV
0443           XPVAL(1)=0.2D0*XPVL
0444           XPVAL(2)=0.8D0*XPVL
0445         ELSEIF(KFA.EQ.22.AND.MINT(105).EQ.333) THEN
0446           XPQ(3)=XPQ(3)+XPV
0447           XPVAL(3)=XPVL
0448         ELSEIF(KFA.EQ.22.AND.MINT(105).EQ.443) THEN
0449           XPQ(4)=XPQ(4)+XPV
0450           XPVAL(4)=XPVL
0451           IF(MSTP(55).GE.9) THEN
0452             DO 260 KFL=-6,6
0453               XPQ(KFL)=0D0
0454   260       CONTINUE
0455           ENDIF
0456         ELSE
0457           XPQ(1)=XPQ(1)+0.5D0*XPV
0458           XPQ(2)=XPQ(2)+0.5D0*XPV
0459           XPVAL(1)=0.5D0*XPVL
0460           XPVAL(2)=0.5D0*XPVL
0461         ENDIF
0462         DO 270 KFL=1,6
0463           XPQ(-KFL)=XPQ(KFL)
0464           XPVAL(-KFL)=XPVAL(KFL)
0465   270   CONTINUE
0466
0467 C...Rescale for gammaVDM by effective gamma -> rho coupling.
0468 C+++Do not rescale?
0469         IF(KFA.EQ.22.AND.MINT(109).EQ.2.AND..NOT.(MSTP(56).EQ.1
0470      &  .AND.MSTP(55).GE.5.AND.MSTP(55).LE.12)) THEN
0471           DO 280 KFL=-6,6
0472             XPQ(KFL)=VINT(281)*XPQ(KFL)
0473             XPVAL(KFL)=VINT(281)*XPVAL(KFL)
0474   280     CONTINUE
0475           VINT(232)=VINT(281)*XPV
0476         ENDIF
0477
0478 C...Simple recipes for kaons.
0479       ELSEIF(KFA.EQ.321) THEN
0480         XPQ(-3)=XPQ(-3)+XPQ(-1)-XPQ(1)
0481         XPQ(-1)=XPQ(1)
0482         XPVAL(-3)=XPVAL(-1)
0483         XPVAL(-1)=0D0
0484       ELSEIF(KFA.EQ.130.OR.KFA.EQ.310) THEN
0485         XPS=0.5D0*(XPQ(1)+XPQ(-2))
0486         XPV=0.5D0*(XPQ(2)+XPQ(-1))-XPS
0487         XPQ(2)=XPS
0488         XPQ(-1)=XPS
0489         XPQ(1)=XPQ(1)+0.5D0*XPV
0490         XPQ(-1)=XPQ(-1)+0.5D0*XPV
0491         XPQ(3)=XPQ(3)+0.5D0*XPV
0492         XPQ(-3)=XPQ(-3)+0.5D0*XPV
0493         XPV=0.5D0*(XPVAL(2)+XPVAL(-1))
0494         XPVAL(2)=0D0
0495         XPVAL(-1)=0D0
0496         XPVAL(1)=0.5D0*XPV
0497         XPVAL(-1)=0.5D0*XPV
0498         XPVAL(3)=0.5D0*XPV
0499         XPVAL(-3)=0.5D0*XPV
0500
0501 C...Isospin conjugation for neutron.
0502       ELSEIF(KFA.EQ.2112) THEN
0503         XPSV=XPQ(1)
0504         XPQ(1)=XPQ(2)
0505         XPQ(2)=XPSV
0506         XPSV=XPQ(-1)
0507         XPQ(-1)=XPQ(-2)
0508         XPQ(-2)=XPSV
0509         XPSV=XPVAL(1)
0510         XPVAL(1)=XPVAL(2)
0511         XPVAL(2)=XPSV
0512
0513 C...Simple recipes for hyperon (average valence parton distribution).
0514       ELSEIF(KFA.EQ.3122.OR.KFA.EQ.3112.OR.KFA.EQ.3212.OR.KFA.EQ.3222
0515      &  .OR.KFA.EQ.3312.OR.KFA.EQ.3322.OR.KFA.EQ.3334) THEN
0516         XPV=(XPQ(1)+XPQ(2)-XPQ(-1)-XPQ(-2))/3D0
0517         XPS=0.5D0*(XPQ(-1)+XPQ(-2))
0518         XPQ(1)=XPS
0519         XPQ(2)=XPS
0520         XPQ(-1)=XPS
0521         XPQ(-2)=XPS
0522         XPQ(KFA/1000)=XPQ(KFA/1000)+XPV
0523         XPQ(MOD(KFA/100,10))=XPQ(MOD(KFA/100,10))+XPV
0524         XPQ(MOD(KFA/10,10))=XPQ(MOD(KFA/10,10))+XPV
0525         XPV=(XPVAL(1)+XPVAL(2))/3D0
0526         XPVAL(1)=0D0
0527         XPVAL(2)=0D0
0528         XPVAL(KFA/1000)=XPVAL(KFA/1000)+XPV
0529         XPVAL(MOD(KFA/100,10))=XPVAL(MOD(KFA/100,10))+XPV
0530         XPVAL(MOD(KFA/10,10))=XPVAL(MOD(KFA/10,10))+XPV
0531       ENDIF
0532
0533 C...Charge conjugation for antiparticle.
0534       IF(KF.LT.0) THEN
0535         DO 290 KFL=1,25
0536           IF(KFL.EQ.21.OR.KFL.EQ.22.OR.KFL.EQ.23.OR.KFL.EQ.25) GOTO 290
0537           XPSV=XPQ(KFL)
0538           XPQ(KFL)=XPQ(-KFL)
0539           XPQ(-KFL)=XPSV
0540   290   CONTINUE
0541         DO 300 KFL=1,6
0542           XPSV=XPVAL(KFL)
0543           XPVAL(KFL)=XPVAL(-KFL)
0544           XPVAL(-KFL)=XPSV
0545   300  CONTINUE
0546       ENDIF
0547
0548 C...MULTIPLE INTERACTIONS - PDF RESHAPING.
0549 C...Set side.
0550       JS=MINT(30)
0551 C...Only reshape PDFs for the non-first interactions;
0552 C...But need valence/sea separation already from first interaction.
0553       IF ((JS.EQ.1.OR.JS.EQ.2).AND.MINT(35).GE.2) THEN
0554         KFVSEL=KFIVAL(JS,1)
0555 C...If valence quark kicked out of pi0 or gamma then that decides
0556 C...whether we should consider state as d dbar, u ubar, s sbar, etc.
0557         IF(KFVSEL.NE.0.AND.(KFA.EQ.111.OR.KFA.EQ.22)) THEN
0558           XPVL=0D0
0559           DO 310 KFL=1,6
0560             XPVL=XPVL+XPVAL(KFL)
0561             XPQ(KFL)=MAX(0D0,XPQ(KFL)-XPVAL(KFL))
0562             XPVAL(KFL)=0D0
0563   310     CONTINUE
0564           XPQ(IABS(KFVSEL))=XPQ(IABS(KFVSEL))+XPVL
0565           XPVAL(IABS(KFVSEL))=XPVL
0566           DO 320 KFL=1,6
0567             XPQ(-KFL)=XPQ(KFL)
0568             XPVAL(-KFL)=XPVAL(KFL)
0569   320     CONTINUE
0570
0571 C...If valence quark kicked out of K0S or K0S then that decides whether
0572 C...we should consider state as d sbar or s dbar.
0573         ELSEIF(KFVSEL.NE.0.AND.(KFA.EQ.130.OR.KFA.EQ.310)) THEN
0574           KFS=1
0575           IF(KFVSEL.EQ.-1.OR.KFVSEL.EQ.3) KFS=-1
0576           XPQ(KFS)=XPQ(KFS)+XPVAL(-KFS)
0577           XPVAL(KFS)=XPVAL(KFS)+XPVAL(-KFS)
0578           XPQ(-KFS)=MAX(0D0,XPQ(-KFS)-XPVAL(-KFS))
0579           XPVAL(-KFS)=0D0
0580           KFS=-3*KFS
0581           XPQ(KFS)=XPQ(KFS)+XPVAL(-KFS)
0582           XPVAL(KFS)=XPVAL(KFS)+XPVAL(-KFS)
0583           XPQ(-KFS)=MAX(0D0,XPQ(-KFS)-XPVAL(-KFS))
0584           XPVAL(-KFS)=0D0
0585         ENDIF
0586
0587 C...XPQ distributions are nominal for a (signed) beam particle
0588 C...of KF type, with 1-Sum(x_prev) rescaled to 1.
0589         CMPFAC=1D0
0590         NRESC=0
0591  345    NRESC=NRESC+1
0592         PVCTOT(JS,-1)=0D0
0593         PVCTOT(JS, 0)=0D0
0594         PVCTOT(JS, 1)=0D0
0595         DO 350 IFL=-6,6
0596           IF(IFL.EQ.0) GOTO 350
0597
0598 C...Count up number of original IFL valence quarks.
0599           IVORG=0
0600           IF(KFIVAL(JS,1).EQ.IFL) IVORG=IVORG+1
0601           IF(KFIVAL(JS,2).EQ.IFL) IVORG=IVORG+1
0602           IF(KFIVAL(JS,3).EQ.IFL) IVORG=IVORG+1
0603 C...For pi0/gamma/K0S/K0L without valence flavour decided yet, here
0604 C...bookkeep as if d dbar (for total momentum sum in valence sector).
0605           IF(KFIVAL(JS,1).EQ.0.AND.IABS(IFL).EQ.1) IVORG=1
0606 C...Count down number of remaining IFL valence quarks. Skip current
0607 C...interaction initiator.
0608           IVREM=IVORG
0609           DO 330 I1=1,NMI(JS)
0610             IF (I1.EQ.MINT(36)) GOTO 330
0611             IF (K(IMI(JS,I1,1),2).EQ.IFL.AND.IMI(JS,I1,2).EQ.0)
0612      &           IVREM=IVREM-1
0613   330     CONTINUE
0614
0615 C...Separate out original VALENCE and SEA content.
0616           VAL=XPVAL(IFL)
0617           SEA=MAX(0D0,XPQ(IFL)-VAL)
0618           XPSVC(IFL,0)=VAL
0619           XPSVC(IFL,-1)=SEA
0620
0621 C...Rescale valence content if changed.
0622           IF (IVORG.NE.0.AND.IVREM.NE.IVORG) XPSVC(IFL,0)=
0623      &    (VAL*IVREM)/IVORG
0624
0625 C...Momentum integrals of original and removed valence quarks.
0626           IF(IVORG.NE.0) THEN
0627 C...For p/n/pbar/nbar beams can split into d_val and u_val.
0628 C...Isospin conjugation for neutrons
0629             IF(KFA.EQ.2212.OR.KFA.EQ.2112) THEN
0630               IAFLP=IABS(IFL)
0631               IF (KFA.EQ.2112) IAFLP=3-IAFLP
0632               VPAVG=PAVG(IAFLP,Q2)
0633 C...For other baryons average d_val and u_val, like for PDFs.
0634             ELSEIF(KFA.GT.1000) THEN
0635               VPAVG=(PAVG(1,Q2)+2D0*PAVG(2,Q2))/3D0
0636 C...For mesons and photon average d_val and u_val and scale by 3/2.
0637 C...Very crude, especially for photon.
0638             ELSE
0639               VPAVG=0.5D0*(PAVG(1,Q2)+2D0*PAVG(2,Q2))
0640             ENDIF
0641             PVCTOT(JS,-1)=PVCTOT(JS,-1)+IVORG*VPAVG
0642             PVCTOT(JS, 0)=PVCTOT(JS, 0)+(IVORG-IVREM)*VPAVG
0643           ENDIF
0644
0645 C...Now add companions (at X with partner having been at Z=XASSOC).
0646 C...NOTE: due to the assumed simple x scaling, the partner was at what
0647 C...corresponds to a higher Z than XASSOC, if there were intermediate
0648 C...scatterings. Nothing done about that for the moment.
0649           DO 340 IVC=1,NVC(JS,IFL)
0650 C...Skip companions that have been kicked out
0651             IF (XASSOC(JS,IFL,IVC).LE.0D0) THEN
0652               XPSVC(IFL,IVC)=0D0
0653               GOTO 340
0654             ELSE
0655 C...Momentum fraction of the partner quark.
0656 C...Use rescaled YS = XS/(1-Sum_rest) where X and XS are not in "rest".
0657               XS=XASSOC(JS,IFL,IVC)
0658               XREM=VINT(142+JS)
0659               YS=XS/(XREM+XS)
0660 C...Momentum fraction of the companion quark.
0661 C...Rescale from X = x/XREM to Y = x/(1-Sum_rest) -> factor (1-YS).
0662               Y=X*(1D0-YS)
0663               XPSVC(IFL,IVC)=PYFCMP(Y/CMPFAC,YS/CMPFAC,MSTP(87))
0664 C...Add to momentum sum, with rescaling compensation factor.
0665               XCFAC=(XREM+XS)/XREM*CMPFAC
0666               PVCTOT(JS,1)=PVCTOT(JS,1)+XCFAC*PYPCMP(YS/CMPFAC,MSTP(87))
0667             ENDIF
0668   340     CONTINUE
0669   350   CONTINUE
0670
0671 C...Wait until all flavours treated, then rescale seas and gluon.
0672         XPSVC(0,-1)=XPQ(0)
0673         XPSVC(0,0)=0D0
0674         RSFAC=1D0+(PVCTOT(JS,0)-PVCTOT(JS,1))/(1D0-PVCTOT(JS,-1))
0675         IF (RSFAC.LE.0D0) THEN
0676 C...First calculate factor needed to exactly restore pz cons.
0677           IF (NRESC.EQ.1) CMPFAC =
0678      &         (1D0-(PVCTOT(JS,-1)-PVCTOT(JS,0)))/PVCTOT(JS,1)
0679 C...Add a bit of headroom
0680           CMPFAC=0.99*CMPFAC
0681 C...Try a few times if more headroom is needed, then print error message.
0682           IF (NRESC.LE.10) GOTO 345
0683           CALL PYERRM(15,
0684      &         '(PYPDFU:) Negative reshaping factor persists!')
0685           WRITE(MSTU(11),5300) (PVCTOT(JS,ITMP),ITMP=-1,1), RSFAC
0686           RSFAC=0D0
0687         ENDIF
0688         DO 370 IFL=-6,6
0689           XPSVC(IFL,-1)=RSFAC*XPSVC(IFL,-1)
0690 C...Also store resulting distributions in XPQ
0691           XPQ(IFL)=0D0
0692           DO 360 ISVC=-1,NVC(JS,IFL)
0693             XPQ(IFL)=XPQ(IFL)+XPSVC(IFL,ISVC)
0694   360     CONTINUE
0695   370   CONTINUE
0696 C...Save companion reweighting factor for PYPTIS.
0697         VINT(140)=CMPFAC
0698       ENDIF
0699
0700
0701 C...Allow gluon also in position 21.
0702       XPQ(21)=XPQ(0)
0703
0704 C...Check positivity and reset above maximum allowed flavour.
0705       DO 380 KFL=-25,25
0706         XPQ(KFL)=MAX(0D0,XPQ(KFL))
0707         IF(IABS(KFL).GT.MSTP(58).AND.IABS(KFL).LE.8) XPQ(KFL)=0D0
0708   380 CONTINUE
0709
0710 C...Formats for error printouts.
0711  5000 FORMAT(' Error: x value outside physical range; x =',1P,D12.3)
0712  5100 FORMAT(' Error: illegal particle code for parton distribution;',
0713      &' KF =',I5)
0714  5200 FORMAT(' Error: unknown parton distribution; KF, library, set =',
0715      &3I5)
0716  5300 FORMAT(' Original valence momentum fraction : ',F6.3/
0717      &       ' Removed valence momentum fraction  : ',F6.3/
0718      &       ' Added companion momentum fraction  : ',F6.3/
0719      &       ' Resulting rescale factor           : ',F6.3)
0720
0721 C...Reset side pointer and return
0722  9999 MINT(30)=0
0723
0724       RETURN
0725       END