Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:30

0001 C=======================================================================
0002 C          SSSSSS   IIIIIII  BBBBB   YY      YY   L        L
0003 C         S            I     B    B    YY  YY     L        L
0004 C          SSSSS       I     BBBBB       YY       L        L
0005 C               S      I     B    B      YY       L        L
0006 C         SSSSSS    IIIIIII  BBBBB       YY       LLLLLLL  LLLLLLL
0007 C=======================================================================
0008 C  Code for SIBYLL:  hadronic interaction Monte Carlo event generator
0009 C=======================================================================
0010 C
0011 C   Version 2.1     (28-Sep-2001)
0012 C
0013 C       By   Ralph Engel
0014 C            R.S. Fletcher
0015 C            T.K. Gaisser
0016 C            Paolo Lipari
0017 C            Todor Stanev
0018 C
0019 C-----------------------------------------------------------------------
0020 C***  Please  have people who want this code contact one of the authors.
0021 C***  Please report any problems.       ****
0022 C
0023 C      For a correct copy contact:
0024 C                REngel@bartol.udel.edu
0025 C                Gaisser@bartol.udel.edu
0026 C                Stanev@bartol.udel.edu
0027 C                Lipari@roma1.infn.it
0028 C
0029 c Sept 15, 2008:  all COMMONS aligned for double precision   by D. Heck
0030 C-----------------------------------------------------------------------
0031 
0032 
0033 
0034       SUBROUTINE SIBYLL (K_beam, IATARG, Ecm)
0035 C-----------------------------------------------------------------------
0036 C...Main routine for the production of hadronic events,
0037 C.  generates an inelastic hadronic interaction of
0038 C.  a `projectile particle' of code K_beam with a
0039 C.  target nucleus of mass number A = IATARG (integer)
0040 C.  IATARG = 0 is an "air" nucleus  (superposition of oxygen and nitrogen)
0041 C.  with c.m. energy for the hadron-nucleon system Ecm (GeV)
0042 C.
0043 C.  Allowed values of K_beam: 7,8,9,10,11,12,13,14,-13,-14
0044 C.                            pi+-,K+-,KL,KS,p,n,pbar,nbar
0045 C.
0046 C.  The output is contained in COMMON /S_PLIST/ that contains:
0047 C.
0048 C.     NP           number of final particles
0049 C.     P(1:NP, 1:5) 4-momenta + masses of the final particles
0050 C.     LLIST (1:NP) codes of final particles.
0051 C.  the reaction is studied in the c.m. of  hadron-nucleon system
0052 C.
0053 C.  The COMMON block /S_CHIST/ contains information about the
0054 C.  the structure of the  generated event:
0055 C.    NW   = number of wounded nucleons
0056 C.    NJET = total number of hard interactions
0057 C.    NSOF = total number of soft interactions
0058 C.    NNSOF (1:NW) = number of soft pomeron cuts in each interaction
0059 C.    NNJET (1:NW) = number of minijets produced in each interaction
0060 C.    XJ1 (1:Index) = x1  for each string
0061 C.    XJ2 (1:Index) = x2   "   "     "
0062 C.    PTJET (1:Index) = pT   "   "     "
0063 C.    NNPJET (1:Index) = total number of particles in each string
0064 C.    NNPSTR (1:2*NW) = number of particles in each `beam string'
0065 C.    JDIF(1:NW) = diffraction code
0066 C----------------------------------------------------------------------
0067       SAVE
0068 
0069       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
0070       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
0071       COMMON /S_DEBUG/ Ncall, Ndebug
0072       PARAMETER (NW_max = 20)
0073       PARAMETER (NS_max = 20, NH_max = 50)
0074       PARAMETER (NJ_max = (NS_max+NH_max)*NW_max)
0075       COMMON /S_CHIST/ X1J(NJ_max),X2J(NJ_max),
0076      &    X1JSUM(NW_max),X2JSUM(NW_max),PTJET(NJ_max),PHIJET(NJ_max),
0077      &    NNPJET(NJ_max),NNPSTR(2*NW_max),NNSOF(NW_max),NNJET(NW_max),
0078      &    JDIF(NW_max),NW,NJET,NSOF
0079       COMMON /S_CCSTR/ X1(2*NW_max),X2(2*NW_max),
0080      &    PXB(2*NW_max),PYB(2*NW_max),PXT(2*NW_max),PYT(2*NW_max),
0081      &    IFLB(2*NW_max),IFLT(2*NW_max)
0082       COMMON /S_CLDIF/ LDIFF
0083       COMMON /S_CQDIS/ PPT0 (33),ptflag
0084       COMMON /S_CUTOFF/ STR_mass_val, STR_mass_sea
0085 
0086       DIMENSION LL(6:14)
0087       DATA LL /7*2,2*1/
0088       DATA FOX /0.257/
0089 
0090       if(Ndebug.gt.1)
0091      &  print *,' SIBYLL: called with (K_beam,IATARG,Ecm):',
0092      &  K_beam,IATARG,Ecm
0093 
0094       kb = K_beam
0095       SQS = Ecm
0096       S = SQS*SQS
0097 
0098       Ncall = Ncall+1
0099 
0100  100  CONTINUE
0101 
0102       NP = 0
0103       NJET = 0
0104       NSOF = 0
0105       IATARGET = IATARG
0106 
0107 C...Generate an 'air' interaction by choosing Nitrogen or Oxygen
0108 
0109       IF (IATARGET .EQ. 0) THEN
0110           R = S_RNDM(NP)
0111           IATARGET = 14
0112           IF (R .LT. FOX)  IATARGET = 16
0113       ENDIF
0114       L = LL(IABS(KB))
0115 
0116 C...Generate number NW wounded nucleons, and diffraction code.
0117 
0118 1000  CALL SIB_START_EV (Ecm, L, IATARGET, NW, JDIF)
0119 
0120 C...limits on simulation of pure diffraction dissociation
0121       IF((LDIFF.NE.0).and.(NW.EQ.1)) THEN
0122          IF((LDIFF.EQ.-1) .AND. (JDIF(1).NE.0) ) GOTO 1000
0123          IF((LDIFF.EQ. 1) .AND. ((JDIF(1).NE.0).AND.(JDIF(1).NE.3)))
0124      +     GOTO 1000
0125          IF((LDIFF.EQ. 5) .AND. (JDIF(1).EQ.2)) GOTO 1000
0126          IF((LDIFF.GE. 2) .AND. (LDIFF.LE.4)) THEN
0127            JDIF(1) = LDIFF-1
0128          ENDIF
0129       ENDIF
0130 
0131 C...Diffractive/non-diffractive interactions
0132 
0133       IF((NW.EQ.1).and.(JDIF(1).NE.0)) THEN
0134         CALL SIB_DIFF (KB, JDIF(1), Ecm, 1, IREJ)
0135       ELSE
0136         CALL SIB_NDIFF (KB, IATARGET, Ecm, 1, IREJ)
0137       ENDIF
0138 
0139       IF (IREJ.NE.0) THEN
0140         if(Ndebug.gt.2) print *,
0141      &    'SIBYLL: rejection (Ecm,Ncall,Nw,JDIF):',Ecm,Ncall,NW,JDIF(1)
0142         GOTO 100
0143       ENDIF
0144 
0145 C...Check energy-momentum conservation
0146 
0147       CALL PFsum(1,NP,Esum,PXsum,PYsum,PZsum,NF)
0148       IF (ABS(Esum/(0.5*Ecm*FLOAT(NW+1)) - 1.) .GT. 1.E-03)  THEN
0149 ctp         WRITE(*,*) ' SIBYLL: energy not conserved (L,call): ',L,Ncall
0150 ctp         WRITE(*,*) ' sqs_inp = ', Ecm, ' sqs_out = ', Esum
0151 ctp         CALL SIB_LIST(6)
0152 ctp         WRITE(*,*) ' SIBYLL: event rejected'
0153          goto 100
0154       ENDIF
0155 
0156 C...list final state particles
0157       if(Ndebug.gt.10) call sib_list(6)
0158 
0159       RETURN
0160       END
0161 
0162 
0163       SUBROUTINE SIB_NDIFF (K_beam, IATARGET, Ecm, Irec, IREJ)
0164 C----------------------------------------------------------------------
0165 C...Non-diffractive or multiple non-diff./diff. interactions
0166 C.    Irec  flag to avoid recursive calls of SIB_DIFF and SIB_NDIFF
0167 C----------------------------------------------------------------------
0168       SAVE
0169 
0170       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
0171       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
0172       COMMON /S_DEBUG/ Ncall, Ndebug
0173       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
0174 
0175       PARAMETER (NW_max = 20)
0176       PARAMETER (NS_max = 20, NH_max = 50)
0177       PARAMETER (NJ_max = (NS_max+NH_max)*NW_max)
0178       COMMON /S_CHIST/ X1J(NJ_max),X2J(NJ_max),
0179      &    X1JSUM(NW_max),X2JSUM(NW_max),PTJET(NJ_max),PHIJET(NJ_max),
0180      &    NNPJET(NJ_max),NNPSTR(2*NW_max),NNSOF(NW_max),NNJET(NW_max),
0181      &    JDIF(NW_max),NW,NJET,NSOF
0182       COMMON /S_CCSTR/ X1(2*NW_max),X2(2*NW_max),
0183      &    PXB(2*NW_max),PYB(2*NW_max),PXT(2*NW_max),PYT(2*NW_max),
0184      &    IFLB(2*NW_max),IFLT(2*NW_max)
0185 
0186       COMMON /S_CLDIF/ LDIFF
0187       COMMON /S_CQDIS/ PPT0 (33),ptflag
0188       COMMON /S_CUTOFF/ STR_mass_val, STR_mass_sea
0189 
0190       DIMENSION X2JET(NW_max),BET(2*NW_max),GAM(2*NW_max),EE(2*NW_max)
0191 
0192       DIMENSION QMAS(33)
0193       DIMENSION LL(6:14)
0194       DATA QMAS
0195      &  /2*0.35,0.6,7*0.,2*1.1,1.25,7*0.,1.25,1.1,1.25,7*0,2*1.25,1.5/
0196       DATA LL /7*2,2*1/
0197 
0198       if(Ndebug.gt.1)
0199      &  print *,' SIB_NDIFF: called with (K_beam,IATARGET,Ecm,Irec):',
0200      &  K_beam,IATARGET,Ecm,Irec
0201 
0202       IREJ = 1
0203 
0204       NP_0    = NP
0205       SQS_0   = SQS
0206 
0207       SQS   = Ecm
0208       S     = SQS*SQS
0209 
0210 *        print *,' current NP,QSQ (-3) ',NP,SQS,NP_0
0211 
0212 C...`soft increase of pT'
0213 
0214 C Setting ptflag = 0 will result in
0215 C underestimating the P_t at high energies.
0216       if (ptflag.gt.0.0) then
0217             ptu=.3+.08*log10(sqs/30.)
0218             pts=.45+.08*log10(sqs/30.)
0219             ptqq=.6+.08*log10(sqs/30.)
0220             PPT0 (1) = PTU
0221             PPT0 (2) = PTU
0222             PPT0 (3) = PTS
0223             PPT0 (10) = PTQQ
0224             DO J=11,33
0225                 PPT0(J) = PTQQ
0226             ENDDO
0227       endif
0228 
0229 C...energy-dependent transverse momentum cutoff
0230 
0231       PTmin = PAR(10)+PAR(11)*EXP(PAR(12)*SQRT(LOG(SQS)))
0232       XMIN = 4.*PTmin**2/S
0233       ZMIN = LOG(XMIN)
0234 
0235 2000  CONTINUE
0236 
0237 C...sample multiple interaction configuration
0238 *        print *,' current NP,QSQ (-2a) ',NP,SQS,NP_0
0239 
0240       L = LL(IABS(K_beam))
0241       DO I=1,NW
0242         if(JDIF(I).eq.0) then
0243           CALL CUT_PRO(L, SQS, PTmin, NNSOF(I), NNJET(I))
0244         else
0245           NNSOF(I) = 1
0246           NNJET(I) = 0
0247         endif
0248       ENDDO
0249 
0250 *        print *,' current NP,QSQ (-2b) ',NP,SQS,NP_0
0251 
0252 C...sample x values
0253 
0254       ITRY = 0
0255 3000  CONTINUE
0256       ITRY = ITRY+1
0257       IF(ITRY.GT.5) GOTO 2000
0258       NP = NP_0
0259       NJET = 0
0260       NSOF = 0
0261       Nall = 0
0262       X1JET = 0.
0263       DO JW=1,NW
0264 C...hard sea-sea interactions
0265          X2JET(JW) = 0.
0266          X1JSUM(JW) = 0.
0267          X2JSUM(JW) = 0.
0268          DO JJ=1,NNJET(JW)
0269            Nall = Nall+1
0270            NJET = NJET+1
0271 *          print *,' Ncall,JW,NW,Njet,NNJET(JW),Nall',
0272 *    &       Ncall,JW,NW,Njet,NNJET(JW),Nall
0273            CALL SAMPLE_hard (L,X1J(Nall),X2J(Nall),PTJET(Nall))
0274 *        print *,' current NP,QSQ (-2c) ',NP,SQS,NP_0
0275            X1JET = X1JET + X1J(Nall)
0276            X2JET(JW) = X2JET(JW)+X2J(Nall)
0277            if(Ndebug.gt.2)
0278      &       print *,' SIB_NDIFF: hard JJ,JW,X1JET,X2JET(JW):',
0279      &       JJ,JW,X1JET,X2JET(JW)
0280            IF ((X2JET(JW).GT.0.9).OR.(X1JET.GT.0.9)) then
0281              if(Ndebug.gt.2) print *,
0282      &         ' SIB_NDIFF: not enough phase space (Ncall,Njet):',
0283      &         Ncall,Njet
0284              GOTO 3000
0285            ENDIF
0286            X1JSUM(JW) = X1JSUM(JW)+X1J(Nall)
0287            X2JSUM(JW) = X2JSUM(JW)+X2J(Nall)
0288          ENDDO
0289 C...soft sea-sea interactions
0290          NSOF_JW = 0
0291          DO JJ=1,NNSOF(JW)-1
0292 ctp060203           CALL SAMPLE_soft (L,STR_mass_sea,X1S,X2S,PTSOF)
0293            CALL SAMPLE_soft (STR_mass_sea,X1S,X2S,PTSOF)
0294 *        print *,' current NP,QSQ (-2d) ',NP,SQS,NP_0
0295            IF ((X2JET(JW)+X2S.LT.0.9).AND.(X1JET+X1S.LT.0.9)) THEN
0296              NSOF = NSOF+1
0297              Nall = Nall+1
0298 *            print *,' Ncall,JW,NW,Nsof,NNSOF(JW),Nall',
0299 *    &         Ncall,JW,NW,Nsof,NNSOF(JW),Nall
0300              NSOF_JW = NSOF_JW+1
0301              X1J(Nall) = X1S
0302              X2J(Nall) = X2S
0303              PTjet(Nall) = PTsof
0304              X1JSUM(JW) = X1JSUM(JW)+X1S
0305              X2JSUM(JW) = X2JSUM(JW)+X2S
0306              X1JET = X1JET + X1S
0307              X2JET(JW) = X2JET(JW)+X2S
0308            ENDIF
0309            if(Ndebug.gt.2)
0310      &       print *,' SIB_NDIFF: soft JJ,JW,X1JET,X2JET(JW):',
0311      &       JJ,JW,X1JET,X2JET(JW)
0312          ENDDO
0313          NNSOF(JW) = NSOF_JW+1
0314 ctp060203 3500    CONTINUE
0315       ENDDO
0316 
0317 *        print *,' current NP,QSQ (-1) ',NP,SQS,NP_0
0318 
0319 C...Prepare 2*NW valence/sea color strings.
0320 
0321       CALL BEAM_SPLIT (K_beam, NW, X1, IFLB, X1JET, LXBAD, STR_mass_val)
0322       IF (LXBAD .EQ. 1) then
0323         if(Ndebug.gt.2) print *,' BEAM_SPLIT: rejection (Ncall):',Ncall
0324         NP    = NP_0
0325         SQS   = SQS_0
0326         S     = SQS*SQS
0327         return
0328       ENDIF
0329 *        print *,' current NP,QSQ (-1a) ',NP,SQS,NP_0
0330       DO J=1,NW
0331          J1=2*(J-1)+1
0332          J2=J1+1
0333 *        print *,' J,J1,J2,NW ',J,J1,J2,NW
0334          KT=13
0335          IF (IATARGET .GT. 1)  KT = 13+INT(2.*S_RNDM(0))
0336          CALL HSPLI (KT,IFLT(J2),IFLT(J1))
0337 *        XMINA = 2.*STR_mass_val/(SQS*(1.-X2JET(J)))
0338          XMINA = 1./(SQS*(1.-X2JET(J)))**2
0339 C        XMINA = 2.*0.20/(SQS*(1.-X2JET(J)))  ! change RSF. 5-92
0340          CHI=CHIDIS (KT,IFLT(J2),IFLT(J1))
0341          XVAL=1.-X2JET(J)
0342          IF (XVAL.LT.XMINA) GOTO 3000
0343          X2(J2) = MAX(CHI*XVAL,XMINA)
0344          X2(J2) = MIN(X2(J2),XVAL-XMINA)
0345          X2(J1) = XVAL-X2(J2)
0346       ENDDO
0347 
0348 C...Generates primordial pT for the partons
0349 *        print *,' current NP,QSQ (-1b) ',NP,SQS,NP_0
0350 
0351       DO J=1,NW
0352          J1 = 2*(J-1)+1
0353          J2 = J1+1
0354          CALL PTDIS (10,PXT(J1),PYT(J1))
0355          if (j.eq.1) then
0356             CALL PTDIS (10,PXB(J2),PYB(J2))
0357          else
0358             CALL PTDIS (IFLB(J2),PXB(J2),PYB(J2))
0359          endif
0360          PXB(J1) = -PXB(J2)
0361          PYB(J1) = -PYB(J2)
0362          PXT(J2) = -PXT(J1)
0363          PYT(J2) = -PYT(J1)
0364       ENDDO
0365 
0366 *        print *,' current NP,QSQ (-1c) ',NP,SQS,NP_0
0367 C...Check consistency of kinematics
0368 
0369       DO J=1,2*NW
0370          EE(J) = SQS*SQRT(X1(J)*X2(J))
0371          XM1 = SQRT(PXB(J)**2+PYB(J)**2+QMAS(IABS(IFLB(J)))**2)
0372          XM2 = SQRT(PXT(J)**2+PYT(J)**2+QMAS(IABS(IFLT(J)))**2)
0373 *        print *,' current NP,QSQ (-1d) ',NP,SQS,NP_0
0374 *        print *,' J,IFLB(J),IFLT(J),NW ',J,IFLB(J),IFLT(J),NW
0375          IF (EE(J) .LT. XM1+XM2+0.3)  GOTO 2000
0376       ENDDO
0377 
0378 C...Fragmentation of soft/hard sea color strings
0379 
0380 *     print *,' current NP,SQS (0)',NP,SQS,NP_0
0381 
0382       DO I=1,Nall
0383         NOLD=NP
0384         CALL JET_FRAG (I)
0385         NNPJET (I) = NP-NOLD
0386 *       print *,' current NP,SQS (1)',NP,SQS,NP_0
0387       ENDDO
0388 
0389 C...Fragment the 2*NW valence/sea color strings
0390 
0391       DO JW=1,NW
0392         if((Irec.eq.1).and.(JDIF(JW).ne.0)) then
0393           J1 = 2*JW-1
0394           J2 = J1+1
0395           X1D = X1(J1)+X1(J2)
0396           X2D = X2(J1)+X2(J2)
0397           EE (J1) = SQS*SQRT(X1D*X2D)
0398           BET(J1) = (X1D-X2D)/(X1D+X2D)
0399           GAM(J1) = (X1D+X2D)/(2.*SQRT(X1D*X2D))
0400           if(JW.eq.1) then
0401             KD = K_beam
0402           else
0403             KD = 9
0404           endif
0405           Nold = NP
0406           call SIB_DIFF(KD, JDIF(JW), EE(J1), 0, IREJ)
0407           if(IREJ.ne.0) print *,' SIB_NDIFF: SIB_DIFF rejection:',Ncall
0408           DO K=NOLD+1,NP
0409             PZ = P(K,3)
0410             P(K,3) = GAM(J1)*(PZ+BET(J1)*P(K,4))
0411             P(K,4) = GAM(J1)*(P(K,4)+BET(J1)*PZ)
0412           ENDDO
0413           NNPSTR(J1) = NP-Nold
0414           NNPSTR(J2) = 0
0415         else
0416           DO J=2*JW-1,2*JW
0417             EE (J) = SQS*SQRT(X1(J)*X2(J))
0418             BET(J) = (X1(J)-X2(J))/(X1(J)+X2(J))
0419             GAM(J) = (X1(J)+X2(J))/(2.*SQRT(X1(J)*X2(J)))
0420             NOLD=NP
0421             CALL STRING_FRAG
0422      &        (EE(J),IFLB(J),IFLT(J),PXB(J),PYB(J),PXT(J),PYT(J),IFBAD)
0423             IF (IFBAD .EQ. 1) then
0424               if(Ndebug.gt.2)
0425      &          print *,' STRING_FRAG: rejection (Ncall):',Ncall
0426               GOTO 2000
0427             ENDIF
0428             DO K=NOLD+1,NP
0429               PZ = P(K,3)
0430               P(K,3) = GAM(J)*(PZ+BET(J)*P(K,4))
0431               P(K,4) = GAM(J)*(P(K,4)+BET(J)*PZ)
0432             ENDDO
0433             NNPSTR(J) = NP-NOLD
0434           ENDDO
0435         endif
0436       ENDDO
0437 
0438       IREJ = 0
0439       SQS   = SQS_0
0440       S     = SQS*SQS
0441 
0442       if(Ndebug.gt.2)
0443      &  print *,'SIB_NDIFF: generated interactions (Ns,Nh):',
0444      &  NSOF+NW,NJET
0445 
0446       RETURN
0447       END
0448 
0449 
0450       SUBROUTINE SIBNUC (IAB, IAT, SQS)
0451 C-----------------------------------------------------------------------
0452 C.  Routine that generates the interaction of a nucleus of
0453 C.  mass number IAB with a  target nucleus  of mass IAT
0454 C.  (IAT=0 : air).
0455 C.  SQS (GeV) is the  center of mass energy of each
0456 C.  nucleon - nucleon cross section
0457 C-----------------------------------------------------------------------
0458       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
0459       COMMON /S_PLNUC/ PA(5,40000), LLA(40000), NPA
0460       COMMON /S_MASS1/ AM(49), AM2(49)
0461       COMMON /CKFRAG/ KODFRAG
0462       PARAMETER (IAMAX=56)
0463       COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
0464      +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
0465      +         ,JJAEL(IAMAX), JJBEL(IAMAX)
0466       COMMON /FRAGMENTS/ PPP(3,60)
0467       DIMENSION SIGDIF(3)
0468       DIMENSION IAF(60)
0469       SAVE
0470       DATA RPOX /0.3624/
0471 
0472 C...Target mass
0473       IF (IAT .EQ. 0) THEN
0474          IATARGET = 14 + 2*INT((1.+RPOX)*S_RNDM(0))
0475       ELSE
0476          IATARGET = IAT
0477       ENDIF
0478 
0479 C...Single nucleon (proton) case
0480 
0481       IF (IAB .EQ. 1)  THEN
0482          NPA = 0
0483          CALL SIBYLL (13,IATARGET, SQS)
0484          CALL DECSIB
0485          DO J=1,NP
0486             LA = IABS(LLIST(J))
0487             IF (LA .LT. 10000)  THEN
0488                NPA = NPA + 1
0489                LLA(NPA) = LLIST(J)
0490                DO K=1,5
0491                   PA(K,NPA) = P(J,K)
0492                ENDDO
0493             ENDIF
0494          ENDDO
0495          RETURN
0496       ENDIF
0497 
0498 
0499 C...Nuclei
0500 
0501       CALL SIB_SIGMA_HP(1,SQS,SIGT,SIGEL,SIG0,SIGDIF,SLOPE,RHO)
0502       CALL INT_NUC (IATARGET, IAB, SIG0, SIGEL)
0503 
0504 C...fragment spectator nucleons
0505       NBT = NB + NBEL
0506       IF (KODFRAG .EQ. 1)  THEN
0507           CALL FRAGM1(IAB,NBT, NF, IAF)
0508       ELSE IF(KODFRAG .EQ. 2)  THEN
0509           CALL FRAGM2(IAB,NBT, NF, IAF)
0510       ELSE
0511           CALL FRAGM (IATARGET, IAB, NBT,B, NF, IAF)
0512       ENDIF
0513 
0514 C...Spectator fragments
0515       NPA = 0
0516       DO J=1,NF
0517          NPA = NPA+1
0518          if(NPA.gt.40000) then
0519            write(6,'(1x,a,2i8)')
0520      &       'SIBNUC: no space left in S_PLNUC (NPA,NF)',NPA,NF
0521            NPA = NPA-1
0522            return
0523          endif
0524          LLA(NPA) = 1000+IAF(J)
0525          PA(1,NPA) = 0.
0526          PA(2,NPA) = 0.
0527          PA(3,NPA) = SQS/2.
0528          PA(4,NPA) = SQS/2.
0529          PA(5,NPA) = FLOAT(IAF(J))*0.5*(AM(13)+AM(14))
0530       ENDDO
0531 
0532 C...Elastically scattered fragments
0533       DO J=1,NBEL
0534          NPA = NPA+1
0535          if(NPA.gt.40000) then
0536            write(6,'(1x,a,2i8)')
0537      &       'SIBNUC: no space left in S_PLNUC (NPA,NBEL)',NPA,NBEL
0538            NPA = NPA-1
0539            return
0540          endif
0541          LLA(NPA) = 1001
0542          PA(1,NPA) = 0.
0543          PA(2,NPA) = 0.
0544          PA(3,NPA) = SQS/2.
0545          PA(4,NPA) = SQS/2.
0546          PA(5,NPA) = 0.5*(AM(13)+AM(14))
0547       ENDDO
0548 
0549 C...Superimpose NB  nucleon interactions
0550       DO JJ=1,NB
0551           CALL SIBYLL (13,IATARGET, SQS)
0552           CALL DECSIB
0553           DO J=1,NP
0554              LA = IABS(LLIST(J))
0555              IF (LA .LT. 10000)   THEN
0556                 NPA = NPA + 1
0557                 if(NPA.gt.40000) then
0558                   write(6,'(1x,a,2i8)')
0559      &              'SIBNUC: no space left in S_PLNUC (NPA,NP)',NPA,NP
0560                   NPA = NPA-1
0561                   return
0562                 endif
0563                 LLA(NPA) = LLIST(J)
0564                 DO K=1,5
0565                     PA(K,NPA) = P(J,K)
0566                 ENDDO
0567             ENDIF
0568          ENDDO
0569       ENDDO
0570 
0571       RETURN
0572       END
0573 
0574 
0575 
0576       FUNCTION CHIDIS (KPARTin, IFL1, IFL2)
0577 C...Generate CHI (fraction of energy of a hadron carried by
0578 C.                the valence quark, or diquark, as specified by IFL1)
0579 C.  INPUT KPART = code of particle
0580 C.        IFL1, IFL2 = codes of partons (3, 3bar of color)
0581 C.........................................................
0582       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
0583       COMMON /S_DEBUG/ Ncall, Ndebug
0584       COMMON /S_CPSPL/ CCHIK(3,6:14)
0585       COMMON /S_CUTOFF/ STR_mass_val, STR_mass_sea
0586       SAVE
0587 
0588       kpart=IABS(kpartin)
0589       IFQ=IABS(IFL1)
0590       IF (IFQ.GT.10) IFQ=IABS(IFL2)
0591       CUT=2.*STR_mass_val/SQS
0592 100   CHIDIS=S_RNDM(0)**2
0593       if (chidis.lt.cut) goto 100
0594       if (chidis.gt.(1.-cut)) goto 100
0595       IF((CHIDIS**2/(CHIDIS**2+CUT**2))**0.5
0596      +   *(1.-CHIDIS)**CCHIK(IFQ,KPART).LT.S_RNDM(0)) GOTO 100
0597       CHIDIS = MAX(0.5*CUT,CHIDIS)
0598       CHIDIS = MIN(1.-CUT,CHIDIS)
0599       IF (IABS(IFL1).GT.10)  CHIDIS=1.-CHIDIS
0600       RETURN
0601       END
0602 
0603 
0604 
0605       SUBROUTINE HSPLI (KF, KP1,KP2)
0606 C...This subroutine splits one hadron of code KF
0607 C.  into 2 partons of code KP1 and KP2
0608 C.  KP1 refers to a color triplet [q or (qq)bar]
0609 C.  KP2 to a a color anti-triplet [qbar or (qq)]
0610 C.  allowed inputs:
0611 C.  KF = 6:14 pi0,pi+-,k+-,k0L,k0s, p,n
0612 C.     = -13,-14  pbar,nbar
0613 C-------------------------------------------------
0614       SAVE
0615       L = IABS(KF)-5
0616 C...Test for good input
0617       IF ( (L .LE. 0) .OR. (L.GT. 9) ) THEN
0618          WRITE(6,*)
0619      &      'HSPLI : Routine entered with illegal particle code ',KF
0620       ENDIF
0621       GOTO (50,100,200,300,400,500,500,600,700), L
0622 
0623 50    R = S_RNDM(0)              ! pi0
0624       IF (R.LE.0.)  THEN
0625          KP1 = 1
0626          KP2 = -1
0627       ELSE
0628         KP1 = 2
0629         KP2 = -2
0630       ENDIF
0631       RETURN
0632 100   KP1 = 1                  ! pi+
0633       KP2 = -2
0634       RETURN
0635 200   KP1 = 2                  ! pi-
0636       KP2 = -1
0637       RETURN
0638 300   KP1 = 1                  ! k+
0639       KP2 = -3
0640       RETURN
0641 400   KP1 = 3                  ! k-
0642       KP2 = -1
0643       RETURN
0644 500   KP1 = 2                  ! k0l, k0s
0645       KP2 = -3
0646       IF (S_RNDM(0).GT. 0.5)  THEN
0647         KP1 = 3
0648         KP2 = -2
0649       ENDIF
0650       RETURN
0651 600   R = 6.*S_RNDM(0)            ! p/pbar
0652       IF (R .LT.3.)       THEN
0653         KP1 = 1
0654         KP2 = 12
0655       ELSEIF (R .LT. 4.)  THEN
0656         KP1 = 1
0657         KP2 = 21
0658       ELSE
0659         KP1 = 2
0660         KP2 = 11
0661       ENDIF
0662       IF (KF .LT. 0)      THEN
0663         KPP = KP1
0664         KP1 = -KP2
0665         KP2 = -KPP
0666       ENDIF
0667       RETURN
0668 700   R = 6.*S_RNDM(0)                  ! n/nbar
0669       IF (R .LT.3.)       THEN
0670          KP1 = 2
0671          KP2 = 12
0672       ELSEIF (R .LT. 4.)  THEN
0673         KP1 = 2
0674         KP2 = 21
0675       ELSE
0676         KP1 = 1
0677         KP2 = 22
0678       ENDIF
0679       IF (KF .LT. 0)      THEN
0680         KPP = KP1
0681         KP1 = -KP2
0682         KP2 = -KPP
0683       ENDIF
0684       RETURN
0685       END
0686 
0687 
0688 
0689       SUBROUTINE JET_FRAG (Index)
0690 C-----------------------------------------------------------------------
0691 C.   Fragmentation of a jet-jet system
0692 C.   Input : kinematical variables of a jet-jet system,
0693 C.           taken from /S_CHIST/
0694 C-----------------------------------------------------------------------
0695 
0696       REAL*8 DX1J, DX2J, DBETJ
0697       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
0698       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
0699       COMMON /S_DEBUG/ Ncall, Ndebug
0700       PARAMETER (NW_max = 20)
0701       PARAMETER (NS_max = 20, NH_max = 50)
0702       PARAMETER (NJ_max = (NS_max+NH_max)*NW_max)
0703       COMMON /S_CHIST/ X1J(NJ_max),X2J(NJ_max),
0704      &    X1JSUM(NW_max),X2JSUM(NW_max),PTJET(NJ_max),PHIJET(NJ_max),
0705      &    NNPJET(NJ_max),NNPSTR(2*NW_max),NNSOF(NW_max),NNJET(NW_max),
0706      &    JDIF(NW_max),NW,NJET,NSOF
0707       SAVE
0708       DATA PGG /1./
0709 
0710       if(Ndebug.gt.2) then
0711         print *,' JET_FRAG: called for entry (I,NP):',Index,NP
0712         print *,' JET_FRAG: (X1J,X2J,PTjet):',X1J(Index),X2J(Index),
0713      &    PTjet(Index)
0714       endif
0715 
0716       E0 = SQRT(S*X1J(Index)*X2J(Index))
0717       TH = ASIN(MIN(0.999999,2.*PTJET(Index)/E0))
0718       FI = 6.283185*S_RNDM(0)
0719       NOLD = NP
0720       IF ( (E0.LT.8.) .OR. (S_RNDM(0).GT.PGG)) THEN
0721          IS = -1 + 2.*INT(1.9999*S_RNDM(0))
0722  100     IFL1 = IS*(INT((2.+0.3)*S_RNDM(0))+1)
0723          XM = 2.*QMASS(IFL1)+0.3
0724          if(E0.LE.XM) GOTO 100
0725          CALL STRING_FRAG (E0,IFL1,-IFL1,0.,0.,0.,0.,IFBAD)
0726          if(IFBAD.ne.0) print *,
0727      &     ' JET_FRAG: rejection in STRING_FRAG (IFL,E0):',IFL1,E0
0728       ELSE
0729          CALL GG_FRAG(E0)
0730       ENDIF
0731       DX1J = X1J(Index)
0732       DX2J = X2J(Index)
0733       DBETJ = (DX1J-DX2J)/(DX1J+DX2J)
0734       CALL SIROBO (NOLD+1,NP,TH,FI,0.D0,0.D0,DBETJ)
0735 
0736       if(Ndebug.gt.2) print *,' JET_FRAG: particles produced:',NP-NOLD
0737 
0738       RETURN
0739       END
0740 
0741 
0742 
0743       SUBROUTINE STRING_FRAG(E0,IFL1,IFL2,PX1,PY1,PX2,PY2,IFBAD)
0744 C-----------------------------------------------------------------------
0745 C.  This routine fragments a string of energy E0
0746 C.  the ends of the strings  have flavors IFL1 and IFL2
0747 C.  the particles produced are in the  jet-jet frame
0748 C.  with IFL1 going in the +z direction
0749 C.     E0 = total energy in jet-jet system
0750 C.  This version consider also a primordial pT attached
0751 C.  to the ends of the string PX1,PY1,  PX2,PY2
0752 C.  OUTPUT:  IFBAD =1  kinematically impossible decay
0753 c
0754 c      Modified Nov. 91.  RSF and TSS to fragment symmetrically
0755 c      ie forward and backward are fragmented as leading.
0756 c      Change- Dec. 92  RSF.  call to ptdis moved- to use flavor
0757 c      of NEW quark in fragmentation.
0758 C-----------------------------------------------------------------------
0759 
0760       COMMON /S_DEBUG/ Ncall, Ndebug
0761       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
0762       COMMON /S_MASS1/ AM(49), AM2(49)
0763       DIMENSION WW(2,2), PTOT(4), PX(3),PY(3),IFL(3)
0764       DIMENSION LPOINT(3000), PMQ(3)
0765       LOGICAL LRANK
0766       SAVE
0767       DATA LRANK/.true./
0768 
0769       if(Ndebug.gt.2) then
0770         print *,
0771      &    ' STRING_FRAG: called with (E0,IFL1,IFL2,PX1,PY1,PX2,PY2)',
0772      &    E0,IFL1,IFL2,PX1,PY1,PX2,PY2
0773         print *,' STRING_FRAG: NP before fragmentation:',NP
0774       endif
0775 
0776 C...initialise
0777       NTRY = 0
0778       IFBAD = 0
0779 200      NTRY = NTRY + 1
0780       IF (NTRY .GT. 50)  THEN
0781          IFBAD = 1
0782          RETURN
0783       ENDIF
0784       I = NP
0785       DO K=1,2
0786          WW(K,1) = 1.
0787          WW(K,2) = 0.
0788       ENDDO
0789       PX(1) = PX1
0790       PY(1) = PY1
0791       PX(2) = PX2
0792       PY(2) = PY2
0793       PX(3) = 0.
0794       PY(3) = 0.
0795       PTOT (1) = PX1+PX2
0796       PTOT (2) = PY1+PY2
0797       PTOT (3) = 0.
0798       PTOT (4) = E0
0799       IFL(1) = IFL1
0800       IFL(2) = IFL2
0801       PMQ(1) = QMASS(IFL(1))
0802       PMQ(2) = QMASS(IFL(2))
0803 
0804       IBLEAD = 0
0805 C
0806 C      SET FLAG FOR GENERATION OF LEADING PARTICLES.
0807 C      "AND" IS FOR PPBAR ( DIQUARK AT BOTH ENDS)
0808 C      "OR" IS FOR PP, PPI, ( DIQUARK AT ONE END.)
0809 C
0810       IF (IABS(IFL1) .GT. 10 .AND. IABS(IFL2) .GT. 10)  THEN
0811          IBLEAD = 2
0812          I = I+1
0813          JT = 1.5+S_RNDM(0)
0814          GOTO 350
0815       ENDIF
0816       IF (IABS(IFL1) .GT. 10 .OR. IABS(IFL2) .GT. 10)  THEN
0817          IBLEAD = 1
0818          I = I+1
0819          JT = 2
0820          IF (IABS(IFL2) .GT. 10) JT = 1
0821          GOTO 350
0822       ENDIF
0823 
0824 C...produce new particle: side, pT
0825 300   continue
0826       I=I+1
0827       if(i.gt.8000) then
0828         write(6,'(1x,a,i8)')
0829      &    'STRING_FRAG: no space left in S_PLIST:',I
0830         stop
0831       endif
0832       IF (IBLEAD .GT. 0)  THEN
0833            JT = 3 - JT
0834            GO TO 350
0835        ENDIF
0836 c
0837 ctp060203 349     continue
0838          JT=1.5+S_RNDM(0)
0839  350      JR=3-JT
0840       LPOINT(I) = JT
0841 
0842 C...particle ID and pt.
0843 ctp060203 999        continue
0844       CALL SIB_IFLAV (IFL(JT), 0, IFL(3), LLIST(I))
0845 ctp060302 991    continue
0846       PMQ(3) = QMASS(IFL(3))
0847       P(I,5) = AM(IABS(LLIST(I)))
0848       CALL PTDIS (IFL(3), PX(3),PY(3))
0849 C...fill transverse momentum
0850       P(I,1) = PX(JT) + PX(3)
0851       P(I,2) = PY(JT) + PY(3)
0852       XMT2 = P(I,5)**2+P(I,1)**2+P(I,2)**2
0853 
0854 
0855 C...test end of fragmentation
0856 
0857       WREM2 = PTOT(4)**2-PTOT(1)**2-PTOT(2)**2-PTOT(3)**2
0858       IF (WREM2 .LT. 0.1)  GOTO 200
0859 *     WMIN = PMQ(1)+PMQ(2)+2.*PMQ(3)+ 0.6 + (2.*S_RNDM(0)-1.)*0.2
0860       WMIN = PMQ(1)+PMQ(2)+2.*PMQ(3)+ 1.1 + (2.*S_RNDM(0)-1.)*0.2
0861 c      WMIN = PMQ(jr)+sqrt(xmt2)+pmq(3)+ 1.1 +(2.*S_RNDM(0)-1.)*0.2
0862 c      IF (WREM2 .LT. WMIN**2) goto 400
0863       IF (WREM2 .LT. WMIN**2)    Then!   goto 400
0864          if (abs(ifl(3)).ne.3) GOTO 400
0865           goto 200
0866       endif
0867 
0868 C...Choose z
0869       IF (IBLEAD .GT. 0.and.abs(ifl(jt)).gt.10)  THEN
0870 C        Special frag. for leading Baryon only
0871          Z = ZBLEAD (IABS(LLIST(I)))
0872          IBLEAD = IBLEAD - 1
0873       ELSE
0874          Z = ZDIS (IFL(3),ifl(jt),XMT2)
0875       ENDIF
0876 
0877       WW(JT,2) = Z*WW(JT,1)
0878       WW(JR,2) = XMT2/(WW(JT,2)*E0**2)
0879 
0880       P(I,3) = WW(1,2)*0.5*E0 - WW(2,2)*0.5*E0
0881       P(I,4) = WW(1,2)*0.5*E0 + WW(2,2)*0.5*E0
0882 
0883       DO J=1,4
0884          PTOT (J) = PTOT(J) - P(I,J)
0885       ENDDO
0886       DO K=1,2
0887          WW(K,1) = WW(K,1) - WW(K,2)
0888       ENDDO
0889 
0890 C...Reset pT and flavor at ends of the string
0891       PX(JT) = -PX(3)
0892       PY(JT) = -PY(3)
0893       IFL(JT) =-IFL(3)
0894       PMQ(JT) = PMQ(3)
0895       GOTO 300
0896 
0897 C...Final two hadrons
0898 400      IF (IFL(JR)*IFL(3) .GT. 100)  GOTO 200
0899       CALL SIB_IFLAV (IFL(JR), -IFL(3), IFLA, LLIST(I+1))
0900       P(I+1,5) = AM(IABS(LLIST(I+1)))
0901       P(I,1)   = PX(JT)+PX(3)
0902       P(I,2)   = PY(JT)+PY(3)
0903       I1 = I+1
0904       P(I+1,1) = PX(JR)-PX(3)
0905       P(I+1,2) = PY(JR)-PY(3)
0906       XM1 = P(I,5)**2+P(I,1)**2+P(I,2)**2
0907       XM2 = P(I1,5)**2+P(I1,1)**2+P(I1,2)**2
0908       IF (SQRT(XM1)+SQRT(XM2) .GT. SQRT(WREM2)) GOTO 200
0909       WREM = SQRT(WREM2)
0910       EA1 = (WREM2+XM1-XM2)/(2.*WREM)
0911       PA2 = (EA1**2-XM1)
0912       if (pa2.gt.0)  then
0913             PA = SQRT(PA2)
0914       else
0915             goto 200
0916       endif
0917       BA = PTOT(3)/PTOT(4)
0918       GA = PTOT(4)/WREM
0919       S = FLOAT(3-2*JT)
0920       P(I,3) = GA*(BA*EA1+S*PA)
0921       P(I,4) = GA*(EA1+BA*S*PA)
0922       P(I+1,3) = PTOT(3)-P(I,3)
0923       P(I+1,4) = PTOT(4)-P(I,4)
0924       NA= NP+1
0925       NP=I+1
0926 
0927 C...reorder  particles along chain (in rank)
0928       IF (LRANK)  THEN
0929       N1 = NA-1
0930       N2 = 0
0931       DO J=NA,NP
0932          IF(LPOINT(J) .EQ. 2)  THEN
0933             N2=N2+1
0934             LLIST (NP+N2) = LLIST(J)
0935             DO K=1,5
0936                P(NP+N2,K)=P(J,K)
0937             ENDDO
0938          ELSE
0939             N1= N1+1
0940             IF (N1.LT.J)   THEN
0941                LLIST(N1) = LLIST(J)
0942                DO K=1,5
0943                   P(N1,K) = P(J,K)
0944                ENDDO
0945             ENDIF
0946          ENDIF
0947       ENDDO
0948       JJ=N1
0949       DO J=NP+N2,NP+1,-1
0950          JJ= JJ+1
0951          LLIST(JJ) = LLIST(J)
0952          DO K=1,5
0953              P(JJ,K) = P(J,K)
0954          ENDDO
0955       ENDDO
0956       ENDIF
0957 
0958       if(Ndebug.gt.2)
0959      &  print *,' STRING_FRAG: NP after fragmentation:',NP
0960 
0961       RETURN
0962       END
0963 
0964 
0965 
0966       FUNCTION ZDIS (IFL1,ifl2, XMT2)
0967 C...z distribution
0968       COMMON /S_CZDIS/ FAin, FB0in
0969       COMMON /S_CZDISs/ FAs1, fAs2
0970       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
0971       COMMON /S_DEBUG/ Ncall, Ndebug
0972       SAVE
0973 
0974       fa=fain
0975       fb0=fb0in
0976 CDH   correction  may 10-1996
0977       if (iabs(kb).ge.13) then   ! baryons only
0978           if (abs(ifl2).eq.3)  fa=fain+fas2
0979           if (abs(ifl1).eq.3)  fa=fain+fas1
0980       endif
0981       FB = FB0*XMT2
0982       IF(FA.GT.0.01.AND.ABS(FA-1.)/FB.LE.0.01) ZMAX=FB/(1.+FB)+
0983      +  (1.-FA)*FB**2/(1.+FB)**3
0984       IF(FA.GT.0.01.AND.ABS(FA-1.)/FB.GT.0.01) ZMAX=0.5*(1.+FB-
0985      +  SQRT((1.-FB)**2+4.*FA*FB))/(1.-FA)
0986       IF(ZMAX.LT.0.1)  ZDIV=2.75*ZMAX
0987       IF(ZMAX.GT.0.85)
0988      +     ZDIV=ZMAX-0.6/FB**2+(FA/FB)*ALOG((0.01+FA)/FB)
0989 C...Choice if z, preweighted for peaks at low or high z
0990 100   Z=S_RNDM(0)
0991       IDIV=1
0992       FPRE=1.
0993       IF (ZMAX.LT.0.1)  THEN
0994          IF(1..LT.S_RNDM(0)*(1.-ALOG(ZDIV)))  IDIV=2
0995          IF (IDIV.EQ.1)  Z=ZDIV*Z
0996          IF (IDIV.EQ.2)  Z=ZDIV**Z
0997          IF (IDIV.EQ.2)  FPRE=ZDIV/Z
0998       ELSEIF (ZMAX.GT.0.85)  THEN
0999          IF(1..LT.S_RNDM(0)*(FB*(1.-ZDIV)+1.)) IDIV=2
1000          IF (IDIV.EQ.1)  Z=ZDIV+ALOG(Z)/FB
1001          IF (IDIV.EQ.1)  FPRE=EXP(FB*(Z-ZDIV))
1002          IF (IDIV.EQ.2)  Z=ZDIV+Z*(1.-ZDIV)
1003       ENDIF
1004 C...weighting according to the correct formula
1005       IF (Z.LE.FB/(50.+FB).OR.Z.GE.1.)  GOTO 100
1006       FVAL=(ZMAX/Z)*EXP(FB*(1./ZMAX-1./Z))
1007       IF(FA.GT.0.01)  FVAL=((1.-Z)/(1.-ZMAX))**FA*FVAL
1008       IF(FVAL.LT.S_RNDM(0)*FPRE)  GOTO 100
1009       ZDIS=Z
1010       RETURN
1011       END
1012 
1013 
1014 
1015       FUNCTION ZBLEAD (LB)
1016 C...fragmentation function for leading baryon
1017 C.  simple form:  f(z) = a + x**b
1018 C   INPUT : LB = particle code.
1019 C..................................................
1020       COMMON /S_CZLEAD/ CLEAD, FLEAD
1021 c      COMMON /S_SZLEAD/ CLEADs, FLEADs
1022       COMMON /S_CHP/ ICHP(49), ISTR(49), IBAR(49)
1023       SAVE
1024 
1025             IC = ICHP(Lb)*ISIGN(1,Lb)
1026 
1027       if (lb.ge.34.and.lb.le.39)  then  ! Lambda's and Sigma's
1028   665               ZBLEAD = S_RNDM(0)
1029                 if (zblead.le..01) goto 665
1030 c          zblead=zdisn(1) ! blead**2   ! soft
1031       else if (ic.eq.0)     then
1032           zblead=zdisn(1)   ! blead**2   !soft
1033       else if (ic.eq.1)  then  ! fast protons only
1034             if (abs(lb).eq.13) then
1035               IF (S_RNDM(0) .LT. CLEAD)  THEN
1036   666               ZBLEAD = S_RNDM(0)
1037                 if (zblead.le..01) goto 666
1038               ELSE
1039                   zblead=1.-zdisn(1)  ! zblead**2   !hard
1040               ENDIF
1041             continue
1042            else
1043                zblead=zdisn(1)  ! zblead**2   !hard
1044            endif
1045       else if (ic.eq.2)  then  ! fast delta++
1046           zblead=1.- zdisn(1)  ! (zblead)**.3333
1047       else
1048                zblead=S_RNDM(0) ! zdisn(1)     !hard
1049       endif
1050        RETURN
1051       END
1052 
1053 
1054 
1055       FUNCTION ZDISN (n)
1056 C...Generate (1-x)**n
1057       SAVE
1058 666   rmin=1.1
1059       do i=1,n+1
1060 c  use dummy argument to prevent compiler optimization
1061          R1=S_RNDM(i)
1062          IF (R1.LE.RMIN) RMIN=R1
1063       ENDDO
1064       ZDISn=RMIN
1065       if (zdisn.le..01) goto 666
1066       if (zdisn.ge..99) goto 666
1067       RETURN
1068       END
1069 
1070 
1071 
1072       SUBROUTINE GG_FRAG (E0)
1073 C...This routine fragments a  gluon-gluon system
1074 C.  of mass E0 (GeV)
1075 C.  the particles produced are in the  jet-jet frame
1076 C.  oriented along the z axis
1077 C...........................................................
1078       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
1079       COMMON /S_MASS1/ AM(49), AM2(49)
1080       DIMENSION WW(2,2),PTOT(4),PX(3),PY(3),IFL(3),PMQ(3)
1081       SAVE
1082 
1083 C...Generate the 'forward' leading particle.
1084 100   I = NP+1
1085       I0 = -1 + 2.*INT(1.9999*S_RNDM(0))
1086       CALL SIB_IFLAV(I0,0,IFL1, LDUM)
1087       CALL SIB_IFLAV(IFL1,0,IFL2, LLIST(I))
1088       CALL PTDIS(IFL1,PX1,PY1)
1089       CALL PTDIS(IFL2,PX2,PY2)
1090       P(I,1) = PX1+PX2
1091       P(I,2) = PY1+PY2
1092       P(I,5) = AM(IABS(LLIST(I)))
1093       XM1 = P(I,5)**2+P(I,1)**2+P(I,2)**2
1094       Z1 = ZDIS (IFL1,1,0.25*XM1)
1095       Z2 = ZDIS (IFL2,1,0.25*XM1)
1096       T1  = 4.*XM1/(E0*E0*(Z1+Z2))
1097       P(I,4) = 0.25*E0*(Z1+Z2 + T1)
1098       P(I,3) = 0.25*E0*(Z1+Z2 - T1)
1099 
1100 C...Generate the 'backward' leading particle.
1101       I = I+1
1102       CALL SIB_IFLAV(-I0,0,IFL3, LDUM)
1103       CALL SIB_IFLAV(IFL3,0,IFL4, LLIST(I))
1104       CALL PTDIS(IFL3,PX3,PY3)
1105       CALL PTDIS(IFL4,PX4,PY4)
1106       P(I,1) = PX3+PX4
1107       P(I,2) = PY3+PY4
1108       P(I,5) = AM(IABS(LLIST(I)))
1109       XM2 = P(I,5)**2+P(I,1)**2+P(I,2)**2
1110       Z3 = ZDIS (IFL3,1,0.25*XM2)
1111       Z4 = ZDIS (IFL4,1,0.25*XM2)
1112       T2  = 4.*XM2/(E0*E0*(Z3+Z4))
1113       P(I,4) = 0.25*E0*( Z3+Z4 + T2)
1114       P(I,3) = 0.25*E0*(-Z3-Z4 + T2)
1115 
1116 C...Fragment the two remaning strings
1117       N0 = 0
1118       DO KS=1,2
1119 
1120       NTRY = 0
1121 200      NTRY = NTRY+1
1122       I = NP+2+N0
1123       IF (NTRY .GT. 30)  GOTO 100
1124 
1125       IF (KS .EQ. 1)  THEN
1126          WW(1,1) = 0.5 * (1 - Z1 - 0.5*T2)
1127          WW(2,1) = 0.5 * (1 - Z3 - 0.5*T1)
1128          PX(1) = -PX1
1129          PY(1) = -PY1
1130          PX(2) = -PX3
1131          PY(2) = -PY3
1132          IFL(1) = -IFL1
1133          IFL(2) = -IFL3
1134       ELSE
1135          WW(1,1) = 0.5 * (1 - Z2 - 0.5*T2)
1136          WW(2,1) = 0.5 * (1 - Z4 - 0.5*T1)
1137          PX(1) = -PX2
1138          PY(1) = -PY2
1139          PX(2) = -PX4
1140          PY(2) = -PY4
1141          IFL(1) = -IFL2
1142          IFL(2) = -IFL4
1143       ENDIF
1144       PX(3) = 0.
1145       PY(3) = 0.
1146       PTOT (1) = PX(1)+PX(2)
1147       PTOT (2) = PY(1)+PY(2)
1148       PTOT (3) = 0.5*E0*(WW(1,1)-WW(2,1))
1149       PTOT (4) = 0.5*E0*(WW(1,1)+WW(2,1))
1150 
1151       PMQ(1) = QMASS(IFL(1))
1152       PMQ(2) = QMASS(IFL(2))
1153 
1154 C...produce new particle: side, pT
1155 300      I=I+1
1156       if(i.gt.8000) then
1157         write(6,'(1x,a,i8)')
1158      &    'GG_FRAG: no space left in S_PLIST:',I
1159         stop
1160       endif
1161       JT=1.5+S_RNDM(0)
1162       JR=3-JT
1163 c      CALL PTDIS (IFL(JT), PX(3),PY(3))
1164 
1165 C...particle ID
1166       CALL SIB_IFLAV (IFL(JT), 0, IFL(3), LLIST(I))
1167       PMQ(3) = QMASS(IFL(3))
1168       P(I,5) = AM(IABS(LLIST(I)))
1169 
1170       CALL PTDIS (IFL(3), PX(3),PY(3))
1171 
1172 C...test end of fragmentation
1173       WREM2 = PTOT(4)**2-PTOT(1)**2-PTOT(2)**2-PTOT(3)**2
1174       IF (WREM2 .LT. 0.1)  GOTO 200
1175       WMIN = PMQ(1)+PMQ(2)+2.*PMQ(3)+1.1 + (2.*S_RNDM(0)-1.)*0.2
1176       IF (WREM2 .LT. WMIN**2)  GOTO 400
1177 
1178 C...fill transverse momentum
1179       P(I,1) = PX(JT) + PX(3)
1180       P(I,2) = PY(JT) + PY(3)
1181 
1182 C...Choose z
1183       XMT2 = P(I,5)**2+P(I,1)**2+P(I,2)**2
1184       Z = ZDIS (ifl(3),IFL(JT), XMT2)
1185 
1186       WW(JT,2) = Z*WW(JT,1)
1187       WW(JR,2) = XMT2/(WW(JT,2)*E0**2)
1188 
1189       P(I,3) = WW(1,2)*0.5*E0 - WW(2,2)*0.5*E0
1190       P(I,4) = WW(1,2)*0.5*E0 + WW(2,2)*0.5*E0
1191 
1192       DO J=1,4
1193          PTOT (J) = PTOT(J) - P(I,J)
1194       ENDDO
1195       DO K=1,2
1196          WW(K,1) = WW(K,1) - WW(K,2)
1197       ENDDO
1198 
1199 C...Reset pT and flavor at ends of the string
1200       PX(JT) = -PX(3)
1201       PY(JT) = -PY(3)
1202       IFL(JT) =-IFL(3)
1203       PMQ(JT) = PMQ(3)
1204       GOTO 300
1205 
1206 C...Final two hadrons
1207 400   IF (IFL(JR)*IFL(3) .GT. 100)  GOTO 200
1208       CALL SIB_IFLAV (IFL(JR), -IFL(3), IFLA, LLIST(I+1))
1209       P(I+1,5) = AM(IABS(LLIST(I+1)))
1210       P(I,1)   = PX(JT)+PX(3)
1211       P(I,2)   = PY(JT)+PY(3)
1212       I1 = I+1
1213       P(I1,1) = PX(JR)-PX(3)
1214       P(I1,2) = PY(JR)-PY(3)
1215       XM1 = P(I,5)**2+P(I,1)**2+P(I,2)**2
1216       XM2 = P(I1,5)**2+P(I1,1)**2+P(I1,2)**2
1217       IF (SQRT(XM1)+SQRT(XM2) .GT. SQRT(WREM2)) GOTO 200
1218       if (ptot(4).le.0) goto 200
1219       WREM = SQRT(WREM2)
1220       EA1 = (WREM2+XM1-XM2)/(2.*WREM)
1221       PA2 = (EA1**2-XM1)
1222       if (pa2.ge.0.0) then
1223         PA = SQRT(pa2)
1224       else
1225        goto 200
1226       endif
1227       BA = PTOT(3)/PTOT(4)
1228       GA = PTOT(4)/WREM
1229       S = FLOAT(3-2*JT)
1230       P(I,3) = GA*(BA*EA1+S*PA)
1231       P(I,4) = GA*(EA1+BA*S*PA)
1232       P(I+1,3) = PTOT(3)-P(I,3)
1233       P(I+1,4) = PTOT(4)-P(I,4)
1234       N0 = I-NP-1
1235       ENDDO                  ! loop on two `remaining strings'
1236       NP = I+1
1237       RETURN
1238       END
1239 
1240 
1241 
1242       FUNCTION QMASS(IFL)
1243 C-----------------------------------------------------------------------
1244 C...Return quark or diquark constituent masses
1245 C-----------------------------------------------------------------------
1246       DIMENSION QMAS(3)
1247       SAVE
1248       DATA QMAS /0.325,0.325,0.5/
1249 
1250       IFLA = IABS(IFL)
1251       IF (IFLA .LE. 3)       THEN
1252          QMASS = QMAS(IFLA)
1253       ELSE
1254          QMA = QMAS(IFLA/10)
1255          QMB = QMAS(MOD(IFLA,10))
1256          QMASS = QMA+QMB
1257       ENDIF
1258       RETURN
1259       END
1260 
1261 
1262 
1263       SUBROUTINE SIB_IFLAV (IFL1,IFL2_A, IFL2, KF)
1264 C-----------------------------------------------------------------------
1265 C.  This subroutine receives as input IFL1 the flavor code
1266 C.  of a quark (antiquark) and  generates the antiquark (quark)
1267 C.  of flavor code IFL2 that combine with the original parton
1268 C.  to compose an hadron of code KF. ONLY 3 FLAVORS
1269 C.  If (IFL2_A.NE.0) returns an hadron KF composed of IFL1 and IFL2_A
1270 C-----------------------------------------------------------------------
1271 
1272       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
1273 
1274       DIMENSION KFLA(3,3,2), CDIAG(12), KDIAG(6)
1275       DIMENSION KBAR(30), CFR(12), KFR(80)
1276       SAVE
1277       DATA KFLA /0,8,10,7,0,22,9,21,0,0,26,29,25,0,31,28,30,0/
1278       DATA CDIAG /0.5,0.25,0.5,0.25,1.,0.5,0.5,0.,0.5,0.,1.,1./
1279       DATA KDIAG /6,23,24,27,32,33/
1280       DATA KBAR /13,14,34,35,36,37,38,9*0,39,3*0,40,41,42,43,44,
1281      +             45,46,47,48,49/
1282       DATA CFR /0.75,0.,0.5,0.,0.,1.,0.1667,0.3333,0.0833,0.6667,
1283      +            0.1667,0.3333/
1284       DATA KFR/0,16,17,19,100,104,109,115,0,26,27,29,122,126,131,137
1285      +  ,0,40,42,47,144,158,178,205,0,1,3,6,10,15,21,28,0,0,56,57,240,
1286      +  246,256,271,0,0,1,3,6,10,15,21,60,61,64,70,292,307,328,356,
1287      +  0,1,3,6,10,15,21,28,16*0/
1288 
1289 
1290       IFLA = IABS(IFL1)
1291       IFL2A = IFL2_A
1292       IF (IFL2A .NE. 0)  THEN
1293          IFL2A = MOD(IFL2A,100)
1294          IFL2 = IFL2A
1295          IFLB = IABS(IFL2A)
1296          MB = 0
1297          IF (IFLB .GT. 10)   MB=1
1298          IF (IFLA .GT. 10)   MB=2
1299       ELSE
1300           MB = 2
1301          IF (IFLA .LT. 10)   THEN
1302              MB = 1
1303              IF ((1.+PAR(1))*S_RNDM(0).LT. 1.)  MB=0
1304          ENDIF
1305       ENDIF
1306 
1307       IF (MB .EQ. 0)  THEN
1308          IF (IFL2A.EQ.0)
1309      +        IFL2=ISIGN(1+INT((2.+PAR(2))*S_RNDM(0)),-IFL1)
1310          IFLD = MAX(IFL1,IFL2)
1311          IFLE = MIN(IFL1,IFL2)
1312          GOTO 100
1313       ENDIF
1314 
1315 C...Decide if the diquark must be split
1316       IF (MB .EQ. 2 .AND. IFLA .GT. 100)   THEN
1317          IFLA = MOD(IFLA,100)
1318            GOTO 200
1319       ENDIF
1320       IF (MB .EQ. 2 .AND. IFLA .EQ. 0)   THEN
1321           IF (S_RNDM(0) .LT. PAR(8))  THEN
1322              MB = 0
1323              IFLG = MOD(IFL1,10)
1324              IFLH =(IFL1-IFLG)/10
1325              IF (S_RNDM(0) .GT. 0.5)  THEN
1326                 IFLDUM = IFLG
1327                 IFLG = IFLH
1328                 IFLH = IFLDUM
1329              ENDIF
1330              IFL11=IFLG
1331              IFL22=ISIGN(1+INT((2.+PAR(2))*S_RNDM(0)),-IFL1)
1332              IFLD = MAX(IFL11,IFL22)
1333              IFLE = MIN(IFL11,IFL22)
1334              IFL2 = -IFLH*10+IFL22
1335              IF (S_RNDM(0) .GT. 0.5)  IFL2 = IFL22*10-IFLH
1336              IFL2 = IFL2+ISIGN(100,IFL2)
1337           ENDIF
1338       ENDIF
1339 
1340 C...Form a meson: consider spin and flavor mixing for the diagonal states
1341 100      IF (MB .EQ. 0)  THEN
1342          IF1 = IABS(IFLD)
1343          IF2 = IABS(IFLE)
1344          IFLC = MAX(IF1,IF2)
1345          KSP = INT(PAR(5)+S_RNDM(0))
1346          KSP = MIN(KSP,1)
1347          IF (IFLC.EQ.3)  KSP = INT(PAR(6)+S_RNDM(0))
1348          IF (IF1 .NE. IF2)   THEN
1349             KF = KFLA(IF1,IF2,KSP+1)
1350          ELSE
1351             R = S_RNDM(0)
1352             JF=1+INT(R+CDIAG(6*KSP+2*IF1-1))+
1353      +             INT(R+CDIAG(6*KSP+2*IF1))
1354             JF = MIN(JF,3)
1355             KF=KDIAG(JF+3*KSP)
1356          ENDIF
1357          RETURN
1358       ENDIF
1359 
1360 C...Form a baryon
1361 200      IF (IFL2A .NE. 0)   THEN
1362           IF (MB .EQ. 1)  THEN
1363              IFLD = IFLA
1364              IFLE = IFLB/10
1365              IFLF = MOD(IFLB,10)
1366           ELSE
1367              IFLD = IFLB
1368              IFLE = IFLA/10
1369              IFLF = MOD(IFLA,10)
1370           ENDIF
1371           LFR = 3+2*((2*(IFLE-IFLF))/(1+IABS(IFLE-IFLF)))
1372           IF(IFLD.NE.IFLE.AND.IFLD.NE.IFLF)  LFR=LFR+1
1373       ELSE
1374 110          CONTINUE
1375           IF(MB.EQ.1)   THEN            ! generate diquark
1376              IFLD = IFLA
1377 120             IFLE = 1+INT((2.+PAR(2)*PAR(3))*S_RNDM(0))
1378              IFLF = 1+INT((2.+PAR(2)*PAR(3))*S_RNDM(0))
1379              IF(IFLE.GE.IFLF.AND.PAR(4).LT.S_RNDM(0))    GOTO 120
1380              IF(IFLE.LT.IFLF.AND.PAR(4)*S_RNDM(0).GT.1.) GOTO 120
1381              IFL2=ISIGN(10*IFLE+IFLF,IFL1)
1382           ELSE                  ! generate quark
1383              IFL2=ISIGN(1+INT((2.+PAR(2))*S_RNDM(0)),IFL1)
1384              IFLD=IABS(IFL2)
1385              IFLE=IFLA/10
1386              IFLF=MOD(IFLA,10)
1387           ENDIF
1388 C...SU(6) factors for baryon formation
1389              LFR=3+2*((2*(IFLE-IFLF))/(1+IABS(IFLE-IFLF)))
1390           IF(IFLD.NE.IFLE.AND.IFLD.NE.IFLF)  LFR=LFR+1
1391           WT = CFR(2*LFR-1)+PAR(7)*CFR(2*LFR)
1392           IF(IFLE.LT.IFLF)   WT=WT/3.
1393           IF (WT.LT.S_RNDM(0)) GOTO 110
1394       ENDIF
1395 
1396 C...Form Baryon
1397       IFLG=MAX(IFLD,IFLE,IFLF)
1398       IFLI=MIN(IFLD,IFLE,IFLF)
1399       IFLH=IFLD+IFLE+IFLF-IFLG-IFLI
1400       KSP=2+2*INT(1.-CFR(2*LFR-1)+(CFR(2*LFR-1)+PAR(7)*
1401      1       CFR(2*LFR))*S_RNDM(0))
1402 
1403 C...Distinguish Lambda- and Sigma- like particles
1404       IF (KSP.EQ.2.AND.IFLG.GT.IFLH.AND.IFLH.GT.IFLI)  THEN
1405       IF(IFLE.GT.IFLF.AND.IFLD.NE.IFLG) KSP=2+INT(0.75+S_RNDM(0))
1406        IF(IFLE.LT.IFLF.AND.IFLD.EQ.IFLG) KSP=3
1407        IF(IFLE.LT.IFLF.AND.IFLD.NE.IFLG) KSP=2+INT(0.25+S_RNDM(0))
1408       ENDIF
1409       KF=KFR(16*KSP-16+IFLG)+KFR(16*KSP-8+IFLH)+IFLI
1410       KF=ISIGN(KBAR(KF-40),IFL1)
1411 
1412       RETURN
1413       END
1414 
1415 
1416 
1417       SUBROUTINE PTDIS (IFL,PX,PY)
1418 C...Generate pT
1419       COMMON /S_CQDIS/ PPT0(33),ptflag
1420       SAVE
1421 
1422       PT = PPT0(IABS(IFL))*SQRT(-ALOG(MAX(1E-10,S_RNDM(0))))
1423       PHI= 6.2831853*S_RNDM(0)
1424       PX=PT*COS(PHI)
1425       PY=PT*SIN(PHI)
1426       RETURN
1427       END
1428 
1429 
1430 
1431       SUBROUTINE SIROBO( NBEG, NEND, THE, PHI, DBEX, DBEY, DBEZ)
1432 C **********************************************************************
1433 C   THIS IS A SLIGHTLY ALTERED VERSION OF "LUROBO" [JETSET63.PYTHIA]   *
1434 C SET TO WORK IN THE SIBYL ENVIROMENT. THE TRANSFORMATION IS PERFORMED *
1435 C ON PARTICLES NUMBER FROM NBEG TO NEND. COMMON BLOCKS CHANGED.        *
1436 C                                      TSS,   Oct '87                  *
1437 C  modification  use directly BETA in double precision in input (PL)   *
1438 C **********************************************************************
1439       COMMON /S_PLIST/ PLIST(8000,5), LLIST(8000), NP
1440       DIMENSION ROT(3,3),PV(3)
1441       DOUBLE PRECISION DP(4),DBEX,DBEY,DBEZ,DGA,DBEP,DGABEP
1442       SAVE
1443 
1444       IF(THE**2+PHI**2 .LE. 1E-20) GO TO 131
1445 C...ROTATE (TYPICALLY FROM Z AXIS TO DIRECTION THETA,PHI)
1446        ROT(1,1)=COS(THE)*COS(PHI)
1447        ROT(1,2)=-SIN(PHI)
1448        ROT(1,3)=SIN(THE)*COS(PHI)
1449        ROT(2,1)=COS(THE)*SIN(PHI)
1450        ROT(2,2)=COS(PHI)
1451        ROT(2,3)=SIN(THE)*SIN(PHI)
1452        ROT(3,1)=-SIN(THE)
1453        ROT(3,2)=0.
1454        ROT(3,3)=COS(THE)
1455        DO 120 I=NBEG,NEND
1456        DO 100 J=1,3
1457  100   PV(J)=PLIST(I,J)
1458        DO 110 J=1,3
1459  110   PLIST(I,J)=ROT(J,1)*PV(1)+ROT(J,2)*PV(2)+ROT(J,3)*PV(3)
1460  120   CONTINUE
1461  131    IF(DBEX**2+DBEY**2+DBEZ**2 .LE. 1D-20) GO TO 151
1462 C...LORENTZ BOOST (TYPICALLY FROM REST TO MOMENTUM/ENERGY=BETA)
1463        DGA=1D0/DSQRT(1D0-DBEX**2-DBEY**2-DBEZ**2)
1464        DO 140 I=NBEG, NEND
1465        DO 130 J=1,4
1466  130   DP(J)=PLIST(I,J)
1467        DBEP=DBEX*DP(1)+DBEY*DP(2)+DBEZ*DP(3)
1468        DGABEP=DGA*(DGA*DBEP/(1D0+DGA)+DP(4))
1469        PLIST(I,1)=DP(1)+DGABEP*DBEX
1470        PLIST(I,2)=DP(2)+DGABEP*DBEY
1471        PLIST(I,3)=DP(3)+DGABEP*DBEZ
1472        PLIST(I,4)=DGA*(DP(4)+DBEP)
1473  140   CONTINUE
1474  151   RETURN
1475       END
1476 
1477 
1478       SUBROUTINE BEAM_SPLIT (L, NW, XX, IFL, XJET, LXBAD, STR_mass)
1479 C...This subroutine split a hadron of code L
1480 C.  into 2*NW partons, each of energy XX(j) and
1481 C.  flavor IFL.  The minimum fractional energy of
1482 C.  each parton is X_min = 2*STR_mass/sqrt(s)
1483 C.
1484 C.  Variable qmas changed to STR_mass to agree with name in SIBYLL
1485 C.      and added to calling sequenceto insure symmetry.
1486 C.     Also a factor of (1-xjet) is added to the def. of xmin for nw=1
1487 C.                               RSF  Apr-2-92
1488 C---------------------------------------------------------------------
1489 
1490       PARAMETER (NW_max = 20)
1491       DIMENSION XX(2*NW_max), IFL(2*NW_max)
1492 
1493       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
1494       COMMON /S_DEBUG/ Ncall, Ndebug
1495       SAVE
1496 
1497       DATA AC /-0.2761856692/             ! log(2) - gamma(Eulero)
1498       DATA GAMMA /2./
1499       DATA NBAD / 0 /
1500 
1501 c-------
1502 c  New code to handle low energy p nuc problem.
1503 c------
1504       LXBAD = 0
1505       XMIN = 2.*STR_mass/SQS
1506       IF (1.-XJET .LT. FLOAT(2*NW)*XMIN)  THEN
1507          NBAD = NBAD + 1
1508          LXBAD = 1
1509 ctp         IF (NBAD .LE. 20) THEN
1510 ctp           WRITE (6, *) 'BEAM_SPLIT: kinematically forbidden situation'
1511 ctp           WRITE (6, 5)  NBAD, SQS, XJET, NW
1512 ctp         ENDIF
1513 ctp 5       FORMAT(1X,'NBAD = ',I3,3X,'sqs = ',E10.3,
1514 ctp     &            3X, 'x_jet = ', F9.3, 3X, ' NW = ',I2)
1515 ctp         IF (NBAD .eq. 20) THEN
1516 ctp           WRITE (6, *)
1517 ctp     &     ' BEAM_SPLIT : Last warning about bad splittings '
1518 ctp           WRITE (6, *) ' The energy threshold is probably too low.'
1519 ctp         ENDIF
1520          RETURN
1521       ENDIF
1522 
1523       IF (NW .EQ. 1)  THEN
1524          XVAL = 1.-XJET
1525          GOTO 200
1526       ENDIF
1527 
1528 C...Choose total energy of sea partons
1529       N = 2*(NW-1)
1530       Z1 = LOG(FLOAT(N))
1531       Z2 = LOG(0.5*SQS*(1.-XJET)/STR_mass-2.)
1532 100   R=S_RNDM(0)
1533       Z=(Z1+AC)*(1.+R*(((Z2+AC)/(Z1+AC))**N-1.))**(1./FLOAT(N))-AC
1534       XSEA = XMIN*EXP(Z)
1535       IF ( (1.-XSEA)**GAMMA .LT. S_RNDM(0)) GOTO 100
1536 C...Split the energy  of sea partons among the different partons
1537       XREM = XSEA - FLOAT(N)*XMIN
1538       DO J=3,N+1
1539 c  use dummy argument to prevent compiler from optimizing
1540          XA = XREM*S_RNDM(j)
1541          XREM = XREM - XA
1542 *        print *,' BEAM_SPLIT: XX index ',J
1543          XX(J) = XMIN + XA
1544       ENDDO
1545 *     print *,' BEAM_SPLIT: XX index ',N+2
1546       XX(N+2) = XMIN + XREM
1547       XVAL = 1.-XSEA-XJET
1548 C...Flavor of sea partons
1549       DO J=1,N/2
1550          J1 =  3 + (J-1)*2
1551 *        print *,' BEAM_SPLIT: flavour indices ',J1,J1+1
1552 c  use dummy argument to prevent compiler from optimizing
1553          IFL(J1) = INT(1.+1.99*S_RNDM(j))
1554          IFL(J1+1) = -IFL(J1)
1555       ENDDO
1556 C...Prepare the valence partons
1557 200   CALL HSPLI (L,IFL(1),IFL(2))
1558       CHI = CHIDIS(L,IFL(1),IFL(2))
1559       XX(1) = MAX(CHI*XVAL,XMIN)
1560       XX(1) = MIN(XX(1),XVAL-XMIN)
1561 C      FOR MESONS, SPLIT ENERGY SYMETRICALLY.
1562 C????? SPLIT K'S WITH ENERGY TO S QUARK?
1563 C
1564       if (abs(l).le.12.and.S_RNDM(0).le.0.5) xx(1)=XVAL-XX(1)
1565       XX(2) = XVAL-XX(1)
1566 
1567       RETURN
1568       END
1569 
1570 C--------------------------------------------------------------------------
1571 C    CODE OF ANALYSIS (not needed to generate events)
1572 C--------------------------------------------------------------------------
1573 
1574       SUBROUTINE PFsum(N1,N2,ETOT,PXT,PYT,PZT,NF)
1575 C...Return the energy,px,py,pz and the number of stable
1576 C.  particles in the list between N1 and N2
1577       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
1578       SAVE
1579 
1580       NF=0
1581       ETOT=0.
1582       PXT=0.
1583       PYT=0.
1584       PZT=0.
1585       DO J=N1,N2
1586          L = LLIST(J)
1587          IF (IABS(L) .LT. 10000)  THEN
1588            NF = NF+1
1589            ETOT = ETOT + ABS( P(J,4) )
1590            PXT = PXT + P(J,1)
1591            PYT = PYT + P(J,2)
1592            PZT = PZT + P(J,3)
1593          ENDIF
1594       ENDDO
1595       RETURN
1596       END
1597 
1598 
1599       SUBROUTINE QNUM (JQ,JS,JB,JBA, NC, NF)
1600 C...Return the quantum numbers of one event
1601 C.  JQ = charge, JB = baryon number, JS = strangeness
1602 C.  JBA = (number of baryons+antibaryons)
1603 C.  NC  = number of charged particles
1604 C.  NF  = number of final particles
1605 C..................................................
1606       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
1607       COMMON /S_CHP/ ICHP(49), ISTR(49), IBAR(49)
1608       SAVE
1609 
1610       JQ = 0
1611       JB = 0
1612       JS = 0
1613       JBA= 0
1614       NC = 0
1615       NF = 0
1616       DO J=1,NP
1617           L = LLIST(J)
1618           LL = IABS(L)
1619           IF (LL .LT. 10000)  THEN
1620               IF(ICHP(LL) .NE. 0) NC = NC + 1
1621               NF = NF + 1
1622               JQ = JQ + ICHP(LL)*ISIGN(1,L)
1623               JB = JB + IBAR(LL)*ISIGN(1,L)
1624               JBA= JBA+ IBAR(LL)
1625               JS = JS + ISTR(LL)*ISIGN(1,L)
1626           ENDIF
1627       ENDDO
1628       RETURN
1629       END
1630 
1631 
1632       SUBROUTINE SIB_LIST(LUN)
1633 C-----------------------------------------------------------------------
1634 C...This routine prints the event record for the
1635 C.  current event on unit LUN
1636 C-----------------------------------------------------------------------
1637 
1638       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
1639       COMMON /S_PLIST1/ LLIST1(8000)
1640       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
1641       COMMON /S_DEBUG/ Ncall, Ndebug
1642       PARAMETER (NW_max = 20)
1643       PARAMETER (NS_max = 20, NH_max = 50)
1644       PARAMETER (NJ_max = (NS_max+NH_max)*NW_max)
1645       COMMON /S_CHIST/ X1J(NJ_max),X2J(NJ_max),
1646      &    X1JSUM(NW_max),X2JSUM(NW_max),PTJET(NJ_max),PHIJET(NJ_max),
1647      &    NNPJET(NJ_max),NNPSTR(2*NW_max),NNSOF(NW_max),NNJET(NW_max),
1648      &    JDIF(NW_max),NW,NJET,NSOF
1649       COMMON /S_CCSTR/ X1(2*NW_max),X2(2*NW_max),
1650      &    PXB(2*NW_max),PYB(2*NW_max),PXT(2*NW_max),PYT(2*NW_max),
1651      &    IFLB(2*NW_max),IFLT(2*NW_max)
1652       CHARACTER*6 NAMP
1653       COMMON /S_CNAM/ NAMP (0:49)
1654       COMMON /S_CHP/ ICHP(49), ISTR(49), IBAR(49)
1655 
1656       CHARACTER CODE*18
1657       CHARACTER*18 NAMDIF(0:3)
1658       SAVE
1659       DATA NAMDIF /'Non-diff. event   ',
1660      &  'Beam diffraction  ','Target diffraction','Double diffraction'/
1661 
1662       WRITE (LUN,*)
1663       WRITE (LUN, *) ' Event record '
1664       if(NW.eq.1) WRITE (LUN,*) '  ',NAMDIF(JDIF(1))
1665       WRITE (LUN,*) '  N_w/N_s/N_j = ', NW, NSOF, NJET
1666       WRITE (LUN,100)
1667 
1668 C...Print particle list
1669       ichar = 0
1670       ibary = 0
1671       DO J=1,NP
1672         L = MOD(LLIST(J),10000)
1673         CODE = '                  '
1674         CODE(1:6) = NAMP(IABS(L))
1675         IF (L .LT. 0) CODE(7:9) = 'bar'
1676         IF(IABS(LLIST(J)) .GT. 10000)   CODE(10:10) = '*'
1677         WRITE (LUN,120) J, CODE, LLIST1(J), (P(J,K),K=1,4)
1678         if(abs(LLIST(J)).LT.10000) then
1679           ichar = ichar+sign(1,l)*ICHP(iabs(l))
1680           ibary = ibary+sign(1,l)*IBAR(iabs(l))
1681         endif
1682       ENDDO
1683       CALL PFsum(1,NP,Esum,PXsum,PYsum,PZsum,NF)
1684       WRITE(LUN,140) PXsum,PYsum,PZsum,Esum
1685 100      FORMAT(3X,'N  Particle',12X,'Ori',6x,'PX',9x,'PY',9x,'PZ'
1686      +         ,9x,'E', /, 3X,70('-'))
1687 120      FORMAT(1X,I4,1X,A18,1X,I4,2X,2(F9.3,2X),2(E9.3,2X))
1688 140      FORMAT(1X,'Tot = ',24X,2(F9.3,2X),G9.3,2X,E9.3)
1689       write(LUN,'(1x,a,i3,3x,a,i3)') 'Total charge:',ichar,
1690      &  'baryon number:',ibary
1691 
1692       RETURN
1693       END
1694 
1695 
1696 
1697       SUBROUTINE KCODE (J,CODE,NC)
1698 C...Produce the code for parton J
1699 C.  Input K, Output CODE, NC=number of characters
1700 C..................................................
1701       CHARACTER*5 CODE
1702       CHARACTER*1 NAMQ(3)
1703       SAVE
1704       DATA NAMQ /'U','D','S'/
1705 
1706       CODE = '     '
1707       IF(J.EQ.0)  THEN
1708          CODE(1:3) = 'GLU'
1709          NC = 3
1710          RETURN
1711       ENDIF
1712       JA = IABS(J)
1713       J1 = MOD(JA,10)
1714       J2 = (JA-J1)/10
1715       IF(JA .GT. 10) THEN
1716          CODE(1:1) = NAMQ(J2)
1717          CODE(2:2) = NAMQ(J1)
1718          NC = 2
1719       ELSE
1720          CODE(1:1) = NAMQ(J1)
1721          NC = 1
1722       ENDIF
1723       IF (J .LT. 0)  THEN
1724          CODE(NC+1:NC+3) = 'bar'
1725          NC = NC+3
1726       ENDIF
1727       RETURN
1728       END
1729 
1730 
1731 
1732 C----------------------------------------------------------------------------
1733 C  Code for sampling
1734 C-----------------------------------------------------------------------------
1735 
1736 ctp060302      SUBROUTINE SAMPLE_soft (L, STR_mass_min, X1,X2,PT)
1737       SUBROUTINE SAMPLE_soft (STR_mass_min, X1,X2,PT)
1738 C-----------------------------------------------------------------------
1739 C...Routine for the sampling the kinematical variables
1740 C.  that characterize a soft cut pomeron (x1,x2, pT)
1741 C.  from the differential cross section:
1742 C.     d3sigma/(dx1 dx2 dpT)
1743 C.  INPUT:  L=1 incident proton, L=2  incident pi
1744 C.          (soft strings identical for pi and p interactions)
1745 C.  OUTPUT:  X1, X2, PT (GeV)
1746 C-----------------------------------------------------------------------
1747 
1748       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
1749       COMMON /S_DEBUG/ Ncall, Ndebug
1750       COMMON /S_CQDIS/ PPT0(33),ptflag
1751       SAVE
1752 
1753       ZSOF = 2.*LOG(STR_mass_min/SQS)
1754  100  Z1=-ZSOF*S_RNDM(0)+ZSOF
1755       Z2=-ZSOF*S_RNDM(0)+ZSOF
1756       IF(Z1+Z2.LE.ZSOF) GOTO 100
1757       X1=EXP(Z1)
1758       X2=EXP(Z2)
1759       STR_mass2 = sqrt(X1*X2*S)/2.
1760  150  PT = PPT0(10)*SQRT(-ALOG(MAX(1E-10,S_RNDM(0))))
1761       IF(PT.GT.PTmin) GOTO 150
1762       IF(PT.GE.STR_mass2) GOTO 150
1763 
1764       RETURN
1765       END
1766 
1767 
1768 
1769       SUBROUTINE SAMPLE_hard (L, X1,X2,PT)
1770 C-----------------------------------------------------------------------
1771 C...Routine for the sampling the kinematical variables
1772 C.  that determine a  jet-jet  system (x1,x2, pT)
1773 C.  from the differential cross section:
1774 C.     d3sigma/(dx1 dx2 dpT)
1775 C.  This version assumes the `single parton approximation'
1776 C.  INPUT:  L=1 incident proton, L=2  incident pi
1777 C.  OUTPUT:  X1, X2, PT (GeV)
1778 C-----------------------------------------------------------------------
1779 
1780       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
1781       COMMON /S_DEBUG/ Ncall, Ndebug
1782       SAVE
1783 
1784 100   Z1=ZSAMPLE (ZMIN,L)
1785       Z2=ZSAMPLE (ZMIN,1)
1786       SIG=1.-XMIN*EXP(-Z1-Z2)
1787       IF (SIG .LT. S_RNDM(0))  GOTO 100
1788       X1=EXP(Z1)
1789       X2=EXP(Z2)
1790       Q2=PTmin**2/(1.-S_RNDM(0)*SIG)
1791       PT=SQRT(Q2*(1.-Q2/(S*X1*X2)))
1792 
1793       RETURN
1794       END
1795 
1796 
1797 
1798       FUNCTION ZSAMPLE (ZMIN,L)
1799 C...This function returns as output a value z=log(x)
1800 C.  distributed as f(x) = g(x) + 4/9 *(q(x) + qbar(x))
1801 C.  from a minimum value ZMIN to 0,
1802 C.  for a proton (L=1) or a pi (L=2)
1803 C.  needs to be initialised with: CALL ZSAMPLE_INI
1804 C.....................................................
1805       COMMON /S_CZGEN/ XA,XB,XMAX,ZA,ZB,ZMAX,DX,DZ,      APART(2),
1806      +   FFA(2),FFB(2),
1807      +   DFX(2),DFZ(2),XX(200,2),ZZ(200,2),FFX(200,2),FFZ(200,2),
1808      +   NX,NZ
1809       SAVE
1810 
1811       F = PART_INT(ZMIN,L)*S_RNDM(0)
1812       IF (F .GE. FFA(L))  THEN
1813          ZSAMPLE = ZA - (F-FFA(L))/APART(L)
1814       ELSE IF (F .GE. FFB(L))  THEN
1815          JF = (F-FFB(L))/DFZ(L) + 1
1816          F0 = FFB(L) + DFZ(L)*FLOAT(JF-1)
1817          T = (F-F0)/DFZ(L)
1818          ZSAMPLE = ZZ(JF,L)*(1.-T)+ZZ(JF+1,L)*T
1819       ELSE
1820          JF = F/DFX(L)+1
1821          F0 = DFX(L)*FLOAT(JF-1)
1822          T = (F-F0)/DFX(L)
1823          X = XX(JF,L)*(1.-T)+XX(JF+1,L)*T
1824          ZSAMPLE = LOG(X)
1825       ENDIF
1826 
1827       RETURN
1828       END
1829 
1830 
1831 
1832       FUNCTION PART_INT (ZMIN,L)
1833 C...This function returns as output the integral of
1834 C.  the parton structure function:
1835 C.     f(x) = g(x) + 4/9 *(q(x) + qbar(x))
1836 C.  from xmin = exp(zmin) to 1
1837 C.  for a proton (L=1) or a pi (L=2)
1838 C.  needs to be initialised with: CALL ZSAMPLE_INI
1839 C.....................................................
1840       COMMON /S_CZGEN/ XA,XB,XMAX,ZA,ZB,ZMAX,DX,DZ,      APART(2),
1841      +   FFA(2),FFB(2),
1842      +   DFX(2),DFZ(2),XX(200,2),ZZ(200,2),FFX(200,2),FFZ(200,2),
1843      +   NX,NZ
1844       SAVE
1845 
1846       IF (ZMIN .LT. ZA)  THEN
1847          PART_INT = FFA(L) + APART(L)*(ZA-ZMIN)
1848       ELSE IF (ZMIN .LT. ZB) THEN
1849          JZ = (ZB-ZMIN)/DZ+1
1850          JZ = min(JZ,199)
1851          Z0 = ZB-DZ*FLOAT(JZ-1)
1852          T = (Z0-ZMIN)/DZ
1853          PART_INT = FFZ(JZ,L)*(1.-T) + FFZ(JZ+1,L)*T
1854       ELSE
1855          X = EXP(ZMIN)
1856          JX = (XMAX-X)/DX+1
1857          JX = min(JX,199)
1858          X0 = XMAX-DX*FLOAT(JX-1)
1859          T = (X0-X)/DX
1860          PART_INT = FFX(JX,L)*(1.-T) + FFX(JX+1,L)*T
1861       ENDIF
1862       RETURN
1863       END
1864 
1865 
1866 
1867       SUBROUTINE ZSAMPLE_INI
1868 C...This subroutine initialise the generation of
1869 C.  z = log(x)  for the generation  of z according
1870 C.  to the structure functions
1871 C..................................................
1872       COMMON /S_CZGEN/ XA,XB,XMAX,ZA,ZB,ZMAX,DX,DZ,      APART(2),
1873      +   FFA(2),FFB(2),
1874      +   DFX(2),DFZ(2),XX(200,2),ZZ(200,2),FFX(200,2),FFZ(200,2),
1875      +   NX,NZ
1876       SAVE
1877 
1878       XA = 1.E-04
1879       XB = 1.E-01
1880       XMAX = 0.80
1881       ZA = LOG(XA)
1882       ZB = LOG(XB)
1883       ZMAX = LOG(XMAX)
1884       NX = 200
1885       NZ = 200
1886       DX = (XMAX-XB)/FLOAT(NX-1)
1887       DZ = (ZB-ZA)/FLOAT(NZ-1)
1888 
1889       DO L=1,2
1890 C         very small x:  f(x) = A/x
1891          APART(L) = PARTON(0.,L)
1892 
1893 C         large x: interpolation in x
1894          FFX(1,L) = 0.
1895          DO J=2,NX
1896             X = XMAX - DX*(FLOAT(J)-0.5)
1897              G = PARTON(X,L)/X
1898             FFX(J,L) = FFX(J-1,L)+G*DX
1899          ENDDO
1900          CALL INVERT_ARRAY (FFX(1,L),XMAX,-DX,NX,XX(1,L),FMIN,
1901      +                        DFX(L))
1902 
1903 C         small x: interpolation in log(x)
1904          FFZ(1,L) = FFX(NX,L)
1905          DO J=2,NZ
1906             Z = ZB - DZ*(FLOAT(J)-0.5)
1907             X = EXP(Z)
1908             G = PARTON(X,L)
1909             FFZ(J,L) = FFZ(J-1,L)+G*DZ
1910          ENDDO
1911          CALL INVERT_ARRAY (FFZ(1,L),ZB,-DZ,NZ,ZZ(1,L),FMIN,DFZ(L))
1912          FFA(L) = FFZ(NZ,L)
1913          FFB(L) = FFX(NX,L)
1914       ENDDO
1915       RETURN
1916       END
1917 
1918 
1919 
1920       FUNCTION PARTON(X,L)
1921 C...This function returns the structure function
1922 C.   f(x) = x * [ g(x) + 4/9 *(q(x) + qbar(x)) ]
1923 C.  for a proton.
1924 C................................................
1925 
1926       parameter (beta=1.925978)
1927       SAVE
1928 
1929       IF (L .EQ. 2)  GOTO 1000
1930 
1931 C...Eichten et al.  (set 1)
1932 ctp060203 100      uv = 1.78 * x**0.5 * (1.-x**1.51)**3.5
1933       uv = 1.78 * x**0.5 * (1.-x**1.51)**3.5
1934       dv = 0.67 * x**0.4 * (1.-x**1.51)**4.5
1935       us = 0.182 * (1.-x)**8.54
1936       ss = 0.081 * (1.-x)**8.54
1937       qq0 = uv + dv + 4.*us + 2.*ss
1938       glu0 = (2.62 + 9.17*x)* (1.-x)**5.90
1939       parton = glu0 + 4./9.*qq0
1940       return
1941 
1942 1000      continue
1943 
1944 C...Owens set 1   from STRF from Wisc. Pheno. group. for q2=q2_min
1945       AV=.4
1946       BV=.7
1947 c      BETA=GGAMMA(AV)*GGAMMA(BV+1.)/GGAMMA(AV+BV+1.)  =1.925978
1948       uv=X**(AV)*(1.-X)**BV/BETA
1949       dv=uv
1950 
1951       A=.9
1952       BET=5.
1953       us=(A*(1.-X)**BET)/6.
1954 
1955       A=.888
1956       BET=3.11
1957       GA1=6.0
1958       glu0=A*(1.-X)**BET*(1.+GA1*X)
1959 c   Bug Fix thanks to Sue Kashahara- correct factor in front of
1960 c   sea quarks for Owens S.F.  5-94
1961       qq0 = uv + dv + 6.*us
1962       parton = (glu0 + 4./9.*qq0)
1963       return
1964 
1965       end
1966 
1967 
1968 
1969       BLOCK DATA PARAM_INI
1970 C-----------------------------------------------------------------------
1971 C....This block data contains default values
1972 C.   of the parameters used in fragmentation
1973 C-----------------------------------------------------------------------
1974 
1975       COMMON /S_DEBUG/ Ncall, Ndebug
1976       COMMON /S_CZDIS/ FA, FB0
1977       COMMON /S_CZDISs/ FAs1, fAs2
1978       COMMON /S_CZLEAD/ CLEAD, FLEAD
1979       COMMON /S_CPSPL/ CCHIK(3,6:14)
1980       COMMON /S_CQDIS/ PPT0 (33),ptflag
1981       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
1982       COMMON /S_CUTOFF/ STR_mass_val, STR_mass_sea
1983       COMMON /CKFRAG/ KODFRAG
1984       SAVE
1985 
1986 C...mass cutoff for soft strings
1987       data STR_mass_val /.35/
1988       data STR_mass_sea /1./
1989 C...Longitudinal Fragmentation function
1990       DATA FA /0.5/, FB0 /0.8/
1991 C...Longitudinal Fragmentation function for leading baryons
1992        DATA CLEAD  /0.6/, FLEAD  /0.6/
1993 c      strange fragmentation
1994       data FAs1 /3./, fAs2 /3./
1995 c      data FAs1 /0./, fAs2 /0./
1996 C...pT of sea partons
1997       DATA PTFLAG /1./
1998       DATA PPT0 /0.30,0.30,0.450,30*0.60/
1999 C...Splitting parameters
2000       DATA CCHIK /21*2.,6*3./
2001 C...Parameters of flavor formation
2002       DATA PAR /0.04,0.25,0.25,0.14,0.3,0.3,0.15,0.,
2003      &          7.0, 11*0. /
2004 C...Fragmentation of nuclei
2005       DATA KODFRAG /0/
2006 C...Debug label and event counter
2007       DATA Ndebug /0/
2008       DATA Ncall /0/
2009 
2010       END
2011 
2012 
2013 
2014       SUBROUTINE PARAM_PRINT(LUN)
2015 
2016       COMMON /S_CZDIS/ FA, FB0
2017       COMMON /S_CZLEAD/ CLEAD, FLEAD
2018       COMMON /S_CPSPL/ CCHIK(3,6:14)
2019       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
2020       COMMON /S_DEBUG/ Ncall, Ndebug
2021       COMMON /S_CQDIS/ PPT0 (33),ptflag
2022       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
2023       SAVE
2024 
2025       WRITE (LUN, 25)
2026 25      FORMAT( //,1x,40('-'), /
2027      +   ' SIBYLL MONTE CARLO PROGRAM. Version 2.1',/,1x,40('-'),/
2028      +   ' List of parameters: ' )
2029 
2030       WRITE (LUN, 31) FA, FB0
2031 31      FORMAT (' Parameters of longitudinal fragmentation: ', /,
2032      +          '  f(z) = (1-z)**a * exp(-b * mt**2/z) ', /,
2033      +          '  a = ', f9.3, 3x, ' b = ', f9.3, ' GeV**-2' )
2034       WRITE (LUN, 32) CLEAD, 1./FLEAD-1.
2035 32      FORMAT (' Parameters of leading fragmentation: ', /,
2036      +   '  f(z) = c + (1-z)**a ', /,
2037      +   '  c = ',f9.3,3x,' a = ',f9.3)
2038 
2039       WRITE (LUN, 35) PPT0(1), PPT0(3), PPT0(11),ppt0(10)
2040 35      FORMAT (' <pT> of sea partons ', /,
2041      +   2x,'<pT>(u/d) ',F8.3,2x,'<pT>(s) ',f8.3,2x,'<pT>(qq) ',f8.3,
2042      +     2x,'<pT>(val) ',f8.3)
2043 
2044       WRITE (LUN, 120) (PAR(K),K=1,12)
2045 120      FORMAT (1x, 'Parameters of flavor formation: ',/,
2046      +   3x,'PAR(1) = Prob(qq)/Prob(q) =              ',F10.2,/,
2047      +   3x,'PAR(2) = Prob(s)/Prob(u)  =              ',F10.2,/,
2048      +   3x,'PAR(3) = Prob(us)/Prob(ud) =             ',F10.2,/,
2049      +   3x,'PAR(4) = Prob(ud_0)/Prob(ud_1) =         ',F10.2,/,
2050      +   3x,'PAR(5) = Prob(Vector)/Prob(Scalar) =     ',F10.2,/,
2051      +   3x,'PAR(6) = Prob(K*)/Prob(K) =              ',F10.2,/,
2052      +   3x,'PAR(7) = Prob(spin 3/2)/Prob(spin=1/2) = ',F10.2,/,
2053      +   3x,'PAR(8) = Prob(B-M-Bbar)/Prob(B-Bbar) =   ',F10.2,/,
2054      +   3x,'PAR(9) = Phase space suppression of MI = ',F10.2,/,
2055      +   3x,'PAR(10)= Low-energy limit for pt cutoff= ',F10.2,/,
2056      +   3x,'PAR(11)= Pt cutoff factor for exp      = ',F10.2,/,
2057      +   3x,'PAR(12)= Pt cutoff factor in exp       = ',F10.2)
2058 
2059       WRITE (LUN, 40)
2060       WRITE (LUN, 41) CCHIK (1,13), CCHIK(2,13)
2061 40      FORMAT(' Parameters of hadron splitting ' )
2062 41      FORMAT('   p -> [(ud) u] splitting: alpha = ', F10.3, /,
2063      +         '   p -> [(uu) d] splitting: alpha = ', F10.3 )
2064 
2065       RETURN
2066       END
2067 
2068 
2069 
2070 C-----------------------------------------------------------------------
2071 C  Code for diffraction
2072 C-----------------------------------------------------------------------
2073 
2074 
2075       SUBROUTINE SIB_DIFF (L0, JDIF1, Ecm, Irec, IREJ)
2076 C-----------------------------------------------------------------------
2077 C...diffraction dissociation
2078 C.  INPUT L0 = index of "beam particle"
2079 C.             the target is assumed to be a proton.
2080 C.    JDIF1 = 1  "beam diffraction"
2081 C.          = 2  "target diffraction"
2082 C.          = 3  "double diffraction"
2083 C     Irec  flag to avoid recursive calls of SIB_DIFF and SIB_NDIFF
2084 C-----------------------------------------------------------------------
2085 
2086       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
2087       COMMON /S_RUN/ SQS, S, PTmin, XMIN, ZMIN, kb ,kt
2088       COMMON /S_DEBUG/ Ncall, Ndebug
2089       COMMON /S_MASS1/ AM(49), AM2(49)
2090       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
2091       DIMENSION XM2MIN(3), ALXMIN(3)
2092       DIMENSION P0(5)
2093       DIMENSION KK(6:14)
2094       SAVE
2095 
2096       DATA PI /3.1415926/
2097       DATA KK /3*2,4*3,2*1/
2098       DATA XM2MIN /1.5, 0.2, 0.6/                  ! M_x**2(min) GeV**2
2099       DATA ALXMIN /0.405465,-1.6094379,-0.5108256/ ! log[M_x**2(min)]
2100       DATA SLOP0 /6.5/                 ! b (slope_ for Mx**2 > 5 GeV**2
2101       DATA ASLOP /31.10362/            ! fit to the slope parameter.
2102       DATA BSLOP /-15.29012/
2103 
2104       if(Ndebug.gt.1)
2105      &  print *,' SIB_DIFF: called with (L0,JDIF1,Ecm):',
2106      &  L0,JDIF1,Ecm
2107 
2108       IREJ = 1
2109       LA = IABS(L0)
2110       XM2MAX = PAR(13)*Ecm*Ecm
2111 
2112 C...Double diffraction
2113       IF (JDIF1 .EQ. 3)   THEN
2114          K = KK(LA)
2115          AL = LOG(XM2MAX/XM2MIN(K))
2116          ALX = ALXMIN(K) + AL*S_RNDM(0)
2117          XMB2 = EXP(ALX)
2118          XMB = SQRT (XMB2)
2119          AL = LOG(XM2MAX/XM2MIN(1))
2120          ALX = ALXMIN(1) + AL*S_RNDM(0)
2121          XMT2 = EXP(ALX)
2122          XMT = SQRT (XMT2)
2123          X1 = 1.+(XMB2-XMT2)/(Ecm*Ecm)
2124          X2 = 2.-X1
2125          SLOPE = MAX(SLOP0, ASLOP+BSLOP*ALX)
2126 50       T = -LOG(S_RNDM(0))/SLOPE
2127          PT = SQRT(T)
2128          PZ1 = 0.25*Ecm*Ecm*X1*X1-XMB2-PT*PT
2129          PZ2 = 0.25*Ecm*Ecm*X2*X2-XMT2-PT*PT
2130          IF (PZ1.LT.0. .OR. PZ2.LT.0.)   GOTO 50
2131          PHI = PI*S_RNDM(0)
2132          P0(5) = XMB
2133          P0(4) = 0.5*Ecm*X1
2134          P0(1) = PT*COS(PHI)
2135          P0(2) = PT*SIN(PHI)
2136          P0(3) = SQRT(PZ1)
2137          CALL DIFDEC (L0, Irec, P0)
2138          P0(5) = XMT
2139          P0(4) = 0.5*Ecm*X2
2140          P0(1) = -P0(1)
2141          P0(2) = -P0(2)
2142          P0(3) = -SQRT(PZ2)
2143          CALL DIFDEC (13, Irec, P0)
2144          IREJ = 0
2145          RETURN
2146       ENDIF
2147 
2148 C...Single diffraction
2149       IF (JDIF1.EQ. 1)  THEN
2150          K = KK(LA)
2151          EM  = AM(13)
2152          EM2 = AM2(13)
2153          L = 13
2154          ZD = -1.
2155       ELSE
2156          K = 1
2157          EM  = AM(LA)
2158          EM2 = AM2(LA)
2159          L = L0
2160          ZD = +1.
2161       ENDIF
2162 C...Generate the mass of the diffracted system Mx (1/Mx**2 distribution)
2163       AL = LOG(XM2MAX/XM2MIN(K))
2164       ALX = ALXMIN(K) + AL*S_RNDM(0)
2165       XM2 = EXP(ALX)
2166       XM = SQRT (XM2)
2167       XMB = XM
2168       XMT = XM
2169 C...Generate the Kinematics of the pseudoelastic hadron
2170       X = 1.-(XM2-EM2)/(Ecm*Ecm)
2171       NP = NP+1
2172       P(NP,4) = 0.5*Ecm*X
2173       SLOPE = MAX(SLOP0, ASLOP+BSLOP*ALX)
2174 60    T = -LOG(MAX(1.E-10,S_RNDM(0)))/SLOPE
2175       PT = SQRT(T*X)
2176       PZ2 = P(NP,4)**2-EM2 - PT*PT
2177       IF (PZ2 .LT.0.)   GOTO 60
2178       PHI = PI*S_RNDM(0)
2179       P(NP,3) = SQRT(PZ2)*ZD
2180       P(NP,1) = PT*COS(PHI)
2181       P(NP,2) = PT*SIN(PHI)
2182       P(NP,5) = EM
2183       LLIST(NP) = L
2184 C...Generating the hadronic system recoling against the produced particle
2185       P0(5) = SQRT(XM2)
2186       P0(4) = 0.5*Ecm*(2.-X)
2187       DO J=1,3
2188          P0(J) = -P(NP,J)
2189       ENDDO
2190       CALL DIFDEC (L0, Irec, P0)
2191       IREJ = 0
2192 
2193       RETURN
2194       END
2195 
2196 
2197 
2198       SUBROUTINE DIFDEC (L0, Irec, P0)
2199 C-----------------------------------------------------------------------
2200 C..."decay" of an excited state with the quantum numbers
2201 C.   of particle L0 and the 5-momentum P0
2202 C.   - low energy: phase space decay (fire ball model)
2203 C.   - intermediate energy: one-string decay (longitudinal phase space)
2204 C.   - high energy: pomeron-hadron scattering (multi-string model)
2205 C-----------------------------------------------------------------------
2206 
2207       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
2208       COMMON /S_MASS1/ AM(49), AM2(49)
2209       COMMON /S_CHP/ ICHP(49), ISTR(49), IBAR(49)
2210       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
2211       DIMENSION P0(5), LL(10), PD(10,5), BE(3), LCON(6:14)
2212       SAVE
2213       DATA EMIN /0.7/
2214       DATA EMIN2 /10./
2215       DATA LCON /7,6,6,11,11,9,9,14,13/
2216       DATA PCHEX /0.33/            ! probability of charge exchange
2217 
2218       LA = IABS(L0)
2219       DELTAE = P0(5) - AM(LA)
2220 
2221 C...pomeron-hadron scattering (pi0 is used instead of pomeron)
2222       IF ((IPAR(10).gt.0).and.(Irec.gt.0).and.(DELTAE.gt.EMIN2))  THEN
2223          N1 = NP+1
2224 
2225  50      CONTINUE
2226            CALL SIB_NDIFF(LA, 6, P0(5), 0, IREJ)
2227          IF(IREJ.NE.0) THEN
2228            NP = N1-1
2229            GOTO 50
2230          ENDIF
2231 
2232          DO J=1,3
2233             BE(J)=P0(J)/P0(4)
2234          ENDDO
2235          GA=P0(4)/P0(5)
2236          if(P0(3).gt.0.) then
2237            do i=N1,NP
2238              P(I,3) = -P(I,3)
2239            enddo
2240          endif
2241          DO I=N1,NP
2242             BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
2243             DO J=1,3
2244                P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
2245             ENDDO
2246             P(I,4)=GA*(P(I,4)+BEP)
2247          ENDDO
2248 
2249 C..."string-like" decay
2250       ELSE IF (DELTAE .GT. EMIN)  THEN
2251            N1 = NP+1
2252          CALL HSPLI(L0,IFL1,IFL2)
2253          IF (P0(3) .GT. 0.)  THEN
2254             IFLA = IFL2
2255             IFL2 = IFL1
2256             IFL1 = IFLA
2257          ENDIF
2258 10         CALL STRING_FRAG (P0(5),IFL1,IFL2,0.,0.,0.,0.,IFBAD)
2259          IF (IFBAD .EQ. 1)  GOTO 10
2260          DO J=1,3
2261             BE(J)=P0(J)/P0(4)
2262          ENDDO
2263          GA=P0(4)/P0(5)
2264          DO I=N1,NP
2265             BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
2266             DO J=1,3
2267                P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
2268             ENDDO
2269             P(I,4)=GA*(P(I,4)+BEP)
2270          ENDDO
2271 
2272 C...Phase space decay of the excited state
2273       ELSE
2274         AV = 2.*SQRT(DELTAE)
2275 ctp060203 100     NPI = AV*(1.+0.5*GASDEV(0))
2276 100     NPI = AV*(1.+0.5*GASDEV(LA))
2277         IF(NPI.LE.0.OR.NPI.GT.9.OR.AM(LA)+NPI*AM(7)+0.02
2278      .            .GT.P0(5))  GOTO 100
2279         IF (S_RNDM(0).LT.PCHEX)  THEN
2280             LL(NPI+1) = LCON(LA)*ISIGN(1,L0)
2281             IF( (L0 .EQ. 6) .OR. (L0 .EQ. 11) )
2282      .             LL(NPI+1) = LL(NPI+1)+INT(1.99999*S_RNDM(0))
2283         ELSE
2284             LL(NPI+1) = L0
2285         ENDIF
2286         JQQ = ICHP(LA)*ISIGN(1,L0)-
2287      .            ICHP(IABS(LL(NPI+1)))*ISIGN(1,LL(NPI+1))
2288 120     JQTOT = 0.
2289         DO K=1,NPI-1
2290            LL(K) = 6+INT(S_RNDM(0)*2.99999)
2291            JQTOT = JQTOT + ICHP(LL(K))
2292         ENDDO
2293         JQR = JQQ-JQTOT
2294         IF (JQR.LT.-1.OR.JQR.GT.1)  GOTO 120
2295         LL(NPI) = 6+JQR
2296         IF (LL(NPI) .EQ. 5)  LL(NPI)=8
2297         CALL DECPAR (0,P0,NPI+1,LL, PD)
2298         DO J=1,NPI+1
2299            NP = NP+1
2300            LLIST(NP) = LL(J)
2301            DO K=1,5
2302               P(NP,K) = PD(J,K)
2303            ENDDO
2304         ENDDO
2305       ENDIF
2306 
2307       RETURN
2308       END
2309 
2310 
2311       SUBROUTINE CUT_PRO (L, SQS, PTmin, NSOFR, NJETR)
2312 C-----------------------------------------------------------------------
2313 C...Generate a number of soft/hard (jet-)pairs for a 'projectile'
2314 C.  (K=1:p),(K=2:pi) interacting with a nucleon at sqrt(s)=SQS(GeV)
2315 C-----------------------------------------------------------------------
2316 
2317       COMMON /S_DEBUG/ Ncall, Ndebug
2318       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
2319       PARAMETER (NS_max = 20, NH_max = 50)
2320       COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
2321      &    SSIGN(61,3), ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
2322       COMMON /S_CUTOFF/ STR_mass_val, STR_mass_sea
2323       SAVE
2324 
2325       K = L
2326       if(K.eq.3) K = 2
2327 
2328       AL = LOG10 (SQS)
2329       IF (AL .LT. ASQSMIN)  THEN
2330           WRITE(*,*)  ' CUT_PRO:  low sqrt(s) ', SQS
2331           NSOFR = 1
2332           NJETR = 0
2333           RETURN
2334       ENDIF
2335       IF (AL .GT. ASQSMAX)  THEN
2336           WRITE(*,*)  ' CUT_PRO:  sqrt(s) out of bounds ', SQS
2337           NJETR = 0
2338           RETURN
2339       ENDIF
2340 
2341       J1 = (AL - ASQSMIN)/DASQS + 1
2342       J1 = MIN(J1,60)
2343       J1 = MAX(J1,1)
2344       J2 = J1+1
2345       T = (AL-ASQSMIN)/DASQS - FLOAT(J1-1)
2346 
2347       R = 0.9999*S_RNDM(0)
2348       DO I=0,NS_max
2349         DO J=0,NH_max
2350           IF (R.LT.(1.-T)*PJETC(I,J,J1,K)+T*PJETC(I,J,J2,K)) GOTO 100
2351         ENDDO
2352       ENDDO
2353 100   CONTINUE
2354 
2355 C...phase space limitation
2356 
2357  120  CONTINUE
2358       XM = FLOAT(2*I)*STR_mass_sea + FLOAT(2*J)*PTmin
2359       PACC = EXP(PAR(9)*(2.-XM)/SQS)
2360       IF(S_RNDM(0).GT.PACC) THEN
2361         IF(I+J.GT.1) THEN
2362           IF(I.GT.0) THEN
2363             I = I-1
2364             GOTO 120
2365           ELSE IF(J.GT.0) THEN
2366             J = J-1
2367             GOTO 120
2368           ENDIF
2369         ENDIF
2370       ENDIF
2371 
2372       NSOFR = I
2373       NJETR = J
2374 
2375       if(Ndebug.gt.2)
2376      &  print *,' CUT_PRO: (L,SQS,PTmin,Ns,Nh)',K,SQS,PTmin,I,J
2377 
2378       RETURN
2379       END
2380 
2381 
2382 C===========================================================================
2383 C  Code for initialization
2384 C===========================================================================
2385 
2386 
2387       SUBROUTINE SIBYLL_INI
2388 C-----------------------------------------------------------------------
2389 C...Initialization routine for SYBILL
2390 C.
2391 C.  the routine fills the COMMON block /CCSIG/ that contains
2392 C.  important information for the generation of events
2393 C.
2394 C     PARAMETER (NS_max = 20, NH_max = 50)
2395 C     COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
2396 C    &    SSIGN(61,3), ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
2397 C.
2398 C.  NSQS = number of energy points  (61 is current version)
2399 C.  ASQSMIN = log_10 [sqrt(s) GeV]   minimum value
2400 C.  ASQSMIN = log_10 [sqrt(s) GeV]   maximum value
2401 C.  DASQS   = step  in log_10[sqrt(s)]
2402 C.            DASQS = (ASQSMAX - ASQSMIN)/(NSQS-1)
2403 C.
2404 C.  SSIG(J,1) inelastic cross section for pp interaction
2405 C.            at energy: sqrt(s)(GeV) = 10**[ASQSMIN+DASQS*(J-1)]
2406 C.  SSIG(J,2)  inelastic cross section for pi-p interaction
2407 C.  SSIGN(J,1) inelastic cross section for p-Air interaction
2408 C.  SSIGN(J,2) inelastic cross section for pi-Air interaction
2409 C.
2410 C.  PJETC(n_s,n_j,J,1) Cumulative  probability distribution
2411 C.                 for the production of n_s soft interactions and
2412 C.                 n_j (n_j=0:30) jet pairs at sqrt(s) labeled
2413 C.                 by J, for p-p interaction
2414 C.  PJETC(n_s,n_j,J,2) Same as above for pi-p interaction
2415 C.  ALINT(J,1)   proton-air  interaction length (g cm-2)
2416 C.  ALINT(J,2)   pi-air  interaction length (g cm-2)
2417 C-----------------------------------------------------------------------
2418       SAVE
2419 cdh
2420       WRITE(*,100)
2421  100  FORMAT(' ','====================================================',
2422      *     /,' ','|                                                  |',
2423      *     /,' ','|                 S I B Y L L  2.1                 |',
2424      *     /,' ','|                                                  |',
2425      *     /,' ','|         HADRONIC INTERACTION MONTE CARLO         |',
2426      *     /,' ','|                        BY                        |',
2427      *     /,' ','|                   Ralph ENGEL                    |',
2428      *     /,' ','|           R.S. FLETCHER, T.K. GAISSER            |',
2429      *     /,' ','|               P. LIPARI, T. STANEV               |',
2430      *     /,' ','|                                                  |',
2431      *     /,' ','| Publication to be cited when using this program: |',
2432      *     /,' ','| R. Engel et al., Proc. 26th ICRC, 1 (1999) 415   |',
2433      *     /,' ','|                                                  |',
2434      *     /,' ','| last modified:  28. Sept. 2001 by R. Engel       |',
2435      *     /,' ','====================================================',
2436      *     /)
2437 cdh
2438 
2439 *     WRITE(*,*) ' Initialization of SIBYLL 2.1 event generator '
2440       CALL PAR_INI
2441       CALL JET_INI
2442       CALL ZSAMPLE_INI
2443       CALL BLOCK_INI
2444       CALL NUC_GEOM_INI
2445       CALL SIG_AIR_INI
2446 
2447       RETURN
2448       END
2449 
2450 
2451 
2452       SUBROUTINE PAR_INI
2453 C------------------------------------------------------------
2454       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
2455       SAVE
2456 
2457 C...Model switches
2458 
2459 C...amplitude/cross section fit parameter set
2460       IPAR(1) = 1
2461       IPAR(2) = 0
2462 
2463 C...recursive diffraction (default=on)
2464       IPAR(10) = 1
2465 
2466 C...Model parameters
2467 
2468 C...energy dependence of PTmin
2469       PAR(10) = 1.
2470       PAR(11) = 0.065
2471       IF(IPAR(2).EQ.0) THEN
2472         PAR(12) = 0.9
2473       ELSE
2474         PAR(12) = 1.12
2475       ENDIF
2476 
2477 C...max mass in diffraction dissociation (Md_max**2/s)
2478       PAR(13) = 0.2
2479 
2480 
2481       RETURN
2482       END
2483 
2484 
2485 
2486       SUBROUTINE JET_INI
2487 C-----------------------------------------------------------------------
2488 C...Compute table of cross sections, and table of probability
2489 C.  for the production of multiple soft and hard interactions
2490 C.
2491 C.  The output of this routine  is the COMMON block /S_CCSIG/
2492 C.  that contains  the cross sections h-p, h-Air, and the
2493 C.  cumulative probability of NS soft and NH hard interactions
2494 C-----------------------------------------------------------------------
2495 
2496       PARAMETER (NS_max = 20, NH_max = 50)
2497       COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
2498      &    SSIGN(61,3), ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
2499       COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
2500      &    SSIG_DD(61,3),SSIG_B(61,3),SSIG_RHO(61,3)
2501 
2502       DIMENSION Pjet(0:NS_max,0:NH_max)
2503       DIMENSION SIG_df(3),SIGDIF(3),SIGDIF_pi(3),
2504      &          PS_tab(61),PH_tab(61),PT_tab(61)
2505       SAVE
2506 
2507 
2508 C...spacing in energy for table of cross sections.
2509 
2510       NSQS = 61
2511       ASQSMIN = 1.
2512       ASQSMAX = 7.
2513       DASQS = (ASQSMAX-ASQSMIN)/FLOAT(NSQS-1)
2514 
2515 C...initialization of proton and pion tables
2516 
2517       DO KK=1,2
2518 
2519 ctp         WRITE(6,'(2(/,1X,A,A))')
2520 ctp     &     'Table: J, sqs,  PT_cut,  SIG_tot,  SIG_inel,  B_el,  ',
2521 ctp     &     'rho,  <n_s>,  <n_h>',
2522 ctp     &     '-----------------------------------------------------',
2523 ctp     &     '-------------------'
2524 
2525          JINT = KK
2526          DO J=1, NSQS
2527            ASQS = ASQSMIN + DASQS*FLOAT(J-1)
2528            SQS = 10.**ASQS
2529 
2530            CALL SIB_SIG (JINT, SQS, PTmin,
2531      &                   SIG_tot, SIG_inel, SIG_df, B_el, Pjet)
2532 
2533 C...low-energy interpolation with data-parametrizations
2534            call SIB_HADCSL(JINT,SQS,
2535      &                     SIGTOT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2536            if(SQS.le.100.) then
2537              SIG_TOT  = SIGTOT
2538              SIG_inel = SIGINEL
2539              B_EL     = SLOPE
2540            else if(SQS.le.1000.) then
2541              Xi = log(SQS/100.)/2.30258509299405
2542              SIG_TOT  = Xi*SIG_TOT+(1.-Xi)*SIGTOT
2543              SIG_inel = Xi*SIG_inel+(1.-Xi)*SIGINEL
2544              B_EL     = Xi*B_EL+(1.-Xi)*SLOPE
2545            endif
2546 
2547            SSIG_TOT(J,KK) = SIG_TOT
2548            SSIG(J,KK)     = SIG_inel
2549            SSIG_SD1(J,KK) = SIGDIF(1)
2550            SSIG_SD2(J,KK) = SIGDIF(2)
2551            SSIG_DD(J,KK)  = SIG_df(3)
2552            SSIG_B(J,KK)   = B_EL
2553            SSIG_RHO(J,KK) = RHO
2554 
2555            PSUM = 0.
2556            PH = 0.
2557            PS = 0.
2558            DO NS=0,NS_max
2559              DO NJ=0,NH_max
2560 
2561                PS = PS+FLOAT(NS)*Pjet(NS,NJ)
2562                PH = PH+FLOAT(NJ)*Pjet(NS,NJ)
2563 
2564                PSUM = PSUM+Pjet(NS,NJ)
2565                PJETC(NS,NJ,J,KK) = PSUM
2566 
2567              ENDDO
2568            ENDDO
2569            PS_tab(J) = PS
2570            PH_tab(J) = PH
2571            PT_tab(J) = PTmin
2572 
2573 ctp           WRITE(6,'(3X,I2,1P,E12.3,0P,4F8.2,3F8.3)')
2574 ctp     &       JINT,SQS,PTmin,SIG_tot,SIG_inel,B_el,RHO,PS,PH
2575 
2576          ENDDO
2577       ENDDO
2578 
2579 C...initialization of kaon tables
2580 
2581       JINT = 3
2582 
2583 ctp      WRITE(6,'(2(/,1X,A,A))')
2584 ctp     &  'Table: J, sqs,  PT_cut,  SIG_tot,  SIG_inel,  B_el,  ',
2585 ctp     &  'rho,  <n_s>,  <n_h>',
2586 ctp     &  '-----------------------------------------------------',
2587 ctp     &  '-------------------'
2588       DO J=1, NSQS
2589         ASQS = ASQSMIN + DASQS*FLOAT(J-1)
2590         SQS = 10.**ASQS
2591 C...use pion cross section rescaled for high-energy extrapolation
2592         SIG_tot   = SSIG_TOT(J,2)
2593         SIG_inel  = SSIG(J,2)
2594         SIG_df(1) = SSIG_SD1(J,2)
2595         SIG_df(2) = SSIG_SD2(J,2)
2596         SIG_df(3) = SSIG_DD(J,2)
2597         B_el = SSIG_B(J,2)
2598         PTmin = PT_tab(J)
2599         PS = PS_tab(J)
2600         PH = PH_tab(J)
2601 
2602 C...low-energy interpolation with data-parametrizations
2603         call SIB_HADCSL(2,SQS,
2604      &                  SIGTOT_pi,SIGEL_pi,SIGINEL,SIGDIF_pi,SLOPE,RHO)
2605         call SIB_HADCSL(3,SQS,
2606      &                  SIGTOT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2607         SIG_el    = (SIGEL/SIGEL_pi)*(SIG_TOT-SIG_inel)
2608         SIG_TOT   = (SIGTOT/SIGTOT_pi)*SIG_TOT
2609         SIG_inel  = SIG_TOT-SIG_el
2610         SIG_df(3) = (SIGDIF(3)/SIGDIF_pi(3))*SIG_df(3)
2611         if(SQS.le.100.) then
2612           SIG_TOT  = SIGTOT
2613           SIG_inel = SIGINEL
2614           B_EL     = SLOPE
2615         else if(SQS.le.1000.) then
2616           Xi = log(SQS/100.)/2.30258509299405
2617           SIG_TOT  = Xi*SIG_TOT+(1.-Xi)*SIGTOT
2618           SIG_inel = Xi*SIG_inel+(1.-Xi)*SIGINEL
2619           B_EL     = Xi*B_EL+(1.-Xi)*SLOPE
2620         endif
2621 
2622         SSIG_TOT(J,3) = SIG_TOT
2623         SSIG(J,3)     = SIG_inel
2624         SSIG_SD1(J,3) = SIGDIF(1)
2625         SSIG_SD2(J,3) = SIGDIF(2)
2626         SSIG_DD(J,3)  = SIG_df(3)
2627         SSIG_B(J,3)   = B_EL
2628         SSIG_RHO(J,3) = RHO
2629 
2630 ctp        WRITE(6,'(3X,I2,1P,E12.3,0P,4F8.2,3F8.3)')
2631 ctp     &    JINT,SQS,PTmin,SIG_tot,SIG_inel,B_el,RHO,PS,PH
2632 
2633       ENDDO
2634 
2635 
2636       RETURN
2637       END
2638 
2639 
2640       SUBROUTINE INI_WRITE (LUN)
2641 C-----------------------------------------------------------------------
2642 C   This subroutine prints on unit LUN
2643 C   a table of the cross sections  used in the program
2644 C   and of the average number of hard interactions, and the average
2645 C   number of wounded nucleons in a hadron-air interaction
2646 C-----------------------------------------------------------------------
2647 
2648       PARAMETER (NS_max = 20, NH_max = 50)
2649       COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
2650      &    SSIGN(61,3), ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
2651       DIMENSION PJ(2),PS(2),PW(2)
2652       SAVE
2653 
2654       DATA ATARG /14.514/
2655 
2656 *      CALL PARAM_PRINT(LUN)
2657       WRITE (LUN, 10)
2658       WRITE (LUN, 15)
2659       WRITE (LUN, 16)
2660       WRITE (LUN, 18)
2661 10    FORMAT(//,' Table of cross sections, and average number',
2662      &         ' of minijets and wounded nucleons ')
2663 15    FORMAT('        [sqrt(s) in GeV, cross sections in mbarn]. ')
2664 16    FORMAT(' sqrt(s) sig(pp) sig(pA) <n_s> <n_j> <n_w>',
2665      &    ' sig(pip) sig(piA) <n_s> <n_j> <n_w>')
2666 18    FORMAT(1X,77('-') )
2667       DO J=1,61,1
2668          SQS = 10.**(ASQSMIN + DASQS*FLOAT(J-1))
2669 
2670          DO K=1,2
2671 
2672            PW(K) = ATARG*SSIG(J,K)/SSIGN(J,K)
2673 
2674            PJ(K) = 0.
2675            PS(K) = 0.
2676            DO NS=0,NS_max
2677              DO NJ=0,NH_max
2678                IF(NJ.GT.0) THEN
2679                  PROB = PJETC(NS,NJ,J,K) - PJETC(NS,NJ-1,J,K)
2680                ELSE IF(NS.GT.0) THEN
2681                  PROB = PJETC(NS,NJ,J,K) - PJETC(NS-1,NH_max,J,K)
2682                ELSE
2683                  PROB = 0.
2684                ENDIF
2685                PJ(K) = PJ(K)+FLOAT(NJ)*PROB
2686                PS(K) = PS(K)+FLOAT(NS)*PROB
2687              ENDDO
2688            ENDDO
2689 
2690          ENDDO
2691 
2692          WRITE(LUN,20) SQS,SSIG(J,1),SSIGN(J,1),PS(1),PJ(1),PW(1)
2693      &                      ,SSIG(J,2),SSIGN(J,2),PS(2),PJ(2),PW(2)
2694 
2695       ENDDO
2696       WRITE (LUN, 18)
2697 20    FORMAT (1X,E8.2, 2(2F7.1,1X,3F6.2,1X))
2698 
2699       RETURN
2700       END
2701 
2702 C*************************************************************************
2703 C=========================================================================
2704 C. UTILITIES ROUTINES
2705 C=========================================================================
2706 C***********************************************************************
2707 
2708 C=======================================================================
2709 C. Code for the wounded nucleon distribution
2710 C=======================================================================
2711 
2712 
2713       SUBROUTINE SIB_START_EV (SQS, L, IA, NW, JDIF)
2714 C-----------------------------------------------------------------------
2715 C...Beginning of a SIBYLL interaction
2716 C.
2717 C.  INPUT : SQS = c.m.s. energy (GeV)
2718 C.          L = 1:proton, 2:charged pion
2719 C.          IA = mass of target nucleon
2720 C.
2721 C.  OUTPUT: NW    = number of wounded nucleons
2722 C.          JDIF(JW)  = diffraction code    !!!! changed to field !!!!
2723 C.                  (0 : non-diffractive interaction)
2724 C.                  (1 : forward diffraction)
2725 C.                  (2 : backward diffraction)
2726 C.                  (3 : double diffraction)
2727 C.
2728 C-----------------------------------------------------------------------
2729 
2730       PARAMETER (NW_max = 20)
2731       DIMENSION JDIF(NW_max)
2732       COMMON /S_CNCM0/ B, BMAX, NTRY, NA
2733 
2734       DIMENSION SIGDIF(3)
2735       SAVE
2736 
2737 C...sample number of wounded nucleons
2738       CALL SIB_SIGMA_HP(L,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2739       IF (IA .GT. 1)  THEN
2740          CALL INT_H_NUC (IA, SIGT, SLOPE, RHO)
2741       ELSE
2742          NA = 1
2743       ENDIF
2744       NW = NA
2745 
2746 C...new treatment of diffraction
2747       PF = SIGDIF(1)/SIGINEL
2748       PB = SIGDIF(2)/SIGINEL
2749       PD = SIGDIF(3)/SIGINEL
2750       P0 = 1.-PF-PB-PD
2751       P1 = P0 + PF
2752       P2 = P1 + PB
2753       DO K=1, NW
2754         R = S_RNDM(0)
2755         IF (R .LT. P0)  THEN
2756           JDIF(K) = 0
2757         ELSE IF (R .LT. P1)  THEN
2758           JDIF(K) = 1
2759         ELSE IF (R .LT. P2)  THEN
2760           JDIF(K) = 2
2761         ELSE
2762           JDIF(K) = 3
2763         ENDIF
2764       ENDDO
2765 
2766       RETURN
2767       END
2768 
2769 
2770 
2771       SUBROUTINE INT_H_NUC (IA, SIGT, SLOPE, RHO)
2772 C...Compute with a montecarlo method the "multiple interaction structure"
2773 C.  of an hadron-nucleus collision.
2774 C.
2775 C.
2776 C.  INPUT : IA               = mass of target nucleus
2777 C.          SIGT (mbarn)     = total hp cross section
2778 C.          SLOPE (GeV**-2)  = slope of hp elastic scattering
2779 C.          RHO              = real/imaginary part of forward elastic
2780 C.                             scattering amplitude
2781 C.
2782 C.  OUTPUT : in COMMON block /CNCMS0/
2783 C.           B = impact parameter (fm)
2784 C.           BMAX = maximum impact parameter for generation
2785 C.           NTRY = number of "trials" before one interaction
2786 C.           NA = number of wounded nucleons in A
2787 C. Author : P.Lipari  (may 1993)
2788 C---------------------------------------------------------------------------
2789       PARAMETER (IAMAX=56)
2790       COMMON /S_CNCM0/ B, BMAX, NTRY, NA
2791       DIMENSION XA(IAMAX), YA(IAMAX)
2792       SAVE
2793       DATA PI /3.1415926/
2794       DATA CMBARN /0.389385/
2795 
2796       CC = SIGT/(4.*PI*SLOPE*CMBARN)
2797       DEN = 2.*SLOPE*CMBARN*0.1
2798       BMAX = 10.                             ! fm
2799       NTRY = 0
2800       CALL NUC_CONF (IA, XA, YA)
2801 1000  B = BMAX*SQRT(S_RNDM(0))
2802       PHI = 2.*PI*S_RNDM(0)
2803       BX = B*COS(PHI)
2804       BY = B*SIN(PHI)
2805       NTRY = NTRY+1
2806       NA = 0
2807       DO JA=1,IA
2808          S = (XA(JA)-BX)**2 + (YA(JA)-BY)**2
2809          F = EXP(-S/DEN)
2810          PEL = CC*CC*(1.+RHO*RHO)*F*F
2811          PINEL  = 2.*CC*F-PEL
2812          R = S_RNDM(0)
2813          IF (R .LT. PINEL)  THEN
2814             NA = NA + 1
2815          ENDIF
2816       ENDDO
2817       IF (NA .EQ. 0)  GOTO 1000
2818       RETURN
2819       END
2820 
2821 
2822 C==========================================================================
2823 C. Cross sections
2824 C==========================================================================
2825 
2826 
2827       SUBROUTINE SIG_AIR_INI
2828 C-----------------------------------------------------------------------
2829 C...Initialize the cross section and interaction lengths on air
2830 C.  (this version initializes p-air, pi-air, and K-air cross sections)
2831 C-----------------------------------------------------------------------
2832 
2833       PARAMETER (NS_max = 20, NH_max = 50)
2834       COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
2835      &    SSIGN(61,3), ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
2836       COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
2837      &    SSIG_DD(61,3),SSIG_B(61,3),SSIG_RHO(61,3)
2838 
2839       DIMENSION SIGDIF(3)
2840       SAVE
2841 
2842       DATA AVOG /6.0221367E-04/
2843 
2844       ATARGET = 14.514
2845 
2846 C...particle loop (p, pi, K)
2847       DO K=1,3
2848         DO J=1,NSQS
2849 
2850            ASQS = ASQSMIN + DASQS*FLOAT(J-1)
2851            SQS = 10.**ASQS
2852 
2853            CALL SIB_SIGMA_HP(K,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2854            CALL SIG_H_AIR(SIGT, SLOPE, RHO, SSIGT, SSIGEL, SSIGQE)
2855 
2856 C  particle production cross section
2857            SSIGN(J,K) = SSIGT-SSIGQE
2858            ALINT(J,K) = 1./(AVOG*SSIGn(j,K)/ATARGET)
2859 
2860         ENDDO
2861       ENDDO
2862 
2863       RETURN
2864       END
2865 
2866 
2867       SUBROUTINE SIB_SIGMA_HP(L,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2868 C-----------------------------------------------------------------------
2869 C     Hadron-proton cross sections, taken from interpolation table
2870 C     calculated by SIBYLL_INI
2871 C
2872 C     input:       L     1      proton-proton
2873 C                        2      pi-proton
2874 C                        3      K-proton
2875 C                  SQS   sqrt(s)
2876 C
2877 C     output:      SIGT       total cross section (mb)
2878 C                  SIGEL      elastic cross section (mb)
2879 C                  SIGINEL    inelastic cross section (mb)
2880 C                  SIGDIF     diffraction dissociation CS (mb)
2881 C                  SLOPE      elastic slope parameter (GeV^-2)
2882 C                  RHO        real/imaginary part of forward amplitude
2883 C-----------------------------------------------------------------------
2884 
2885       DIMENSION SIGDIF(3)
2886 
2887       PARAMETER (NS_max = 20, NH_max = 50)
2888       COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
2889      &    SSIGN(61,3), ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
2890       COMMON /S_CCSIG2/ SSIG_TOT(61,3),SSIG_SD1(61,3),SSIG_SD2(61,3),
2891      &    SSIG_DD(61,3),SSIG_B(61,3),SSIG_RHO(61,3)
2892       SAVE
2893 
2894       IF(NSQS.LE.0) THEN
2895         WRITE(6,'(//,1X,A)')
2896      &    'SIB_SIGMA_HP: interpolation table not initialized.'
2897         STOP
2898       ENDIF
2899 
2900       AL = LOG10(SQS)
2901       J1 = (AL - 1.)*10. + 1
2902       if((j1.lt.1).or.(j1.ge.NSQS)) then
2903 c        write (6,'(1x,a,i3,1p,e12.3)')
2904 c     &    'SIB_SIGMA_HP: energy out of range ',L,sqs
2905         J1 = min(J1,NSQS-1)
2906         J1 = max(J1,1)
2907       endif
2908       T = (AL-1.)*10. - FLOAT(J1-1)
2909       SIGT    = SSIG_TOT(J1,L)*(1.-T) + SSIG_TOT(J1+1,L)*T
2910       SIGINEL = SSIG(J1,L)*(1.-T) + SSIG(J1+1,L)*T
2911       SIGEL   = SIGT-SIGINEL
2912       SIGDIF(1) = SSIG_SD1(J1,L)*(1.-T) + SSIG_SD1(J1+1,L)*T
2913       SIGDIF(2) = SSIG_SD2(J1,L)*(1.-T) + SSIG_SD2(J1+1,L)*T
2914       SIGDIF(3) = SSIG_DD(J1,L)*(1.-T) + SSIG_DD(J1+1,L)*T
2915       SLOPE   = SSIG_B(J1,L) *(1.-T) + SSIG_B(J1+1,L)*T
2916       RHO     = SSIG_RHO(J1,L) *(1.-T) + SSIG_RHO(J1+1,L)*T
2917 
2918       RETURN
2919       END
2920 
2921 
2922       SUBROUTINE SIB_SIGMA_HAIR (L,SQS,SIGprod)
2923 C-----------------------------------------------------------------------
2924 C     Hadron-air cross sections, taken from interpolation table
2925 C     calculated by SIBYLL_INI
2926 C
2927 C     input:       L     1      proton-air
2928 C                        2      pi-air
2929 C                        3      K-air
2930 C                  SQS   sqrt(s)
2931 C
2932 C     output:      SIGprod    particle production cross section (mb)
2933 C-----------------------------------------------------------------------
2934 
2935       PARAMETER (NS_max = 20, NH_max = 50)
2936       COMMON /S_CCSIG/ SSIG(61,3), PJETC(0:NS_max,0:NH_max,61,2),
2937      &    SSIGN(61,3), ALINT(61,3), ASQSMIN, ASQSMAX, DASQS, NSQS
2938       SAVE
2939 
2940       IF(NSQS.LE.0) THEN
2941         WRITE(6,'(//,1X,A)')
2942      &    'SIB_SIGMA_HAIR: interpolation table not initialized.'
2943         STOP
2944       ENDIF
2945 
2946       AL = LOG10(SQS)
2947       J1 = (AL - 1.)*10. + 1
2948       if((j1.lt.1).or.(j1.ge.NSQS)) then
2949         write (6,'(1x,a,i3,1p,e12.3)')
2950      &    'SIB_SIGMA_HAIR: energy out of range ',L,sqs
2951         J1 = min(J1,NSQS-1)
2952         J1 = max(J1,1)
2953       endif
2954       T = (AL-1.)*10. - FLOAT(J1-1)
2955       SIGprod = SSIGN(J1,L)*(1.-T) + SSIGN(J1+1,L)*T
2956 
2957       RETURN
2958       END
2959 
2960 
2961       SUBROUTINE SIB_HADCSL(L,ECM,SIGTOT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
2962 C-----------------------------------------------------------------------
2963 C     low-energy cross section parametrizations (target always proton)
2964 C
2965 C     input:   L           beam particle: (1 - proton,
2966 C                                          2 - pion,
2967 C                                          3 - kaon)
2968 C                          target is always proton
2969 C              ECM         c.m. energy (GeV)
2970 C
2971 C     output:  SIGTOT      total cross section (mb)
2972 C              SIGEL       elastic cross section (mb)
2973 C              SIGDIF      diffractive cross section (sd-1,sd-2,dd, mb)
2974 C              SLOPE       forward elastic slope (GeV**-2)
2975 C              RHO         real/imaginary part of elastic amplitude
2976 C-----------------------------------------------------------------------
2977       DIMENSION SIGDIF(3)
2978 
2979       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
2980       SAVE
2981 
2982 C  proton-proton cross section as reference
2983       CALL SIB_HADCS1(1,ECM,SIGTOT,SIGEL,SIGINEL,SLOPE,RHO)
2984 
2985 C  parametrization for diffraction
2986       Xi_min = 1.5/(ECM*ECM)
2987       Xi_max = PAR(13)
2988       SIGeff = SIGEL
2989       call SIB_HADCS2(ECM,Xi_min,Xi_max,SIGeff,SIGDIF)
2990 
2991       if(L.eq.1) return
2992 
2993 C  regge motivated rescaling of diffraction dissociation
2994       sigtot_pp = SIGTOT
2995       sigel_pp  = SIGEL
2996       slope_pp  = SLOPE
2997       CALL SIB_HADCS1(L,ECM,SIGTOT,SIGEL,SIGINEL,SLOPE,RHO)
2998       SIGDIF(1) = slope_pp/SLOPE*SIGTOT/sigtot_pp*SIGDIF(1)
2999       SIGDIF(2) = slope_pp/SLOPE*SIGEL/sigel_pp*SIGDIF(2)
3000       SIGDIF(3) = SIGTOT/sigtot_pp*SIGDIF(3)
3001 
3002       RETURN
3003       END
3004 
3005 
3006       SUBROUTINE SIB_HADCS1(L,ECM,SIGTOT,SIGEL,SIGINEL,SLOPE,RHO)
3007 C-----------------------------------------------------------------------
3008 C     low-energy cross section parametrizations
3009 C
3010 C     input:   L           beam particle: (1 - proton,
3011 C                                          2 - pion,
3012 C                                          3 - kaon)
3013 C                          target is always proton
3014 C              ECM         c.m. energy (GeV)
3015 C
3016 C     output:  SIGTOT      total cross section (mb)
3017 C              SIGEL       elastic cross section (mb)
3018 C              SIGDIF      diffractive cross section (sd-1,sd-2,dd, mb)
3019 C              SLOPE       forward elastic slope (GeV**-2)
3020 C              RHO         real/imaginary part of elastic amplitude
3021 C
3022 C     comments:
3023 C     - low-energy data interpolation uses PDG fits from 1992
3024 C     - slopes from ???, new fit to pp data
3025 C     - high-energy extrapolation by Donnachie-Landshoff like fit made
3026 C       by PDG 1996
3027 C     - analytic extension of amplitude to calculate rho
3028 C-----------------------------------------------------------------------
3029 
3030       DIMENSION TPDG92(7,2,6),TPDG96(9,6),BURQ83(3,6),XMA(6)
3031       SAVE
3032 
3033       DATA TPDG92  /
3034      &  3.D0, 2100.D0, 48.D0, 0.D0, 1.D0, 0.522D0, -4.51D0,
3035      &  3.D0, 2100.D0, 11.9D0, 26.9D0, -1.21D0, 0.169D0, -1.85D0,
3036      &  5.D0, 2100.D0, 38.4D0, 77.6D0, -0.64D0, 0.26D0, -1.2D0,
3037      &  5.D0, 2100.D0, 10.2D0, 52.7D0, -1.16D0, 0.125D0, -1.28D0,
3038      &  4.D0, 340.D0,  16.4D0, 19.3D0, -0.42D0, 0.19D0, 0.D0,
3039      &  4.D0, 340.D0,  0.D0, 11.4D0, -0.4D0, 0.079D0, 0.D0,
3040      &  2.5D0, 370.D0, 33.D0, 14.D0, -1.36D0, 0.456D0, -4.03D0,
3041      &  2.5D0, 370.D0, 1.76D0, 11.2D0, -0.64D0, 0.043D0, 0.D0,
3042      &  2.D0, 310.D0,  18.1D0, 0.D0, 1.D0, 0.26D0, -1.D0,
3043      &  2.D0, 310.D0,  5.D0, 8.1D0, -1.8D0, 0.16D0, -1.3D0,
3044      &  3.D0, 310.D0,  32.1D0, 0.D0, 1.D0, 0.66D0, -5.6D0,
3045      &  3.D0, 310.D0,  7.3D0, 0.D0, 1.D0, 0.29D0, -2.4D0  /
3046 
3047       DATA TPDG96  /
3048      &  50.D0, 22.D0,0.079D0,0.25D0,0.D0,
3049      &         77.15D0,-21.05D0,0.46D0,0.9D0,
3050      &  50.D0, 22.D0,0.079D0,0.25D0,0.D0,
3051      &         77.15D0,21.05D0,0.46D0,0.9D0,
3052      &  10.D0, 13.70,0.079D0,0.25D0,0.D0,
3053      &         31.85D0,-4.05D0,0.45D0,0.9D0,
3054      &  10.D0, 13.70,0.079D0,0.25D0,0.D0,
3055      &         31.85D0,4.05D0,0.45D0,0.9D0,
3056      &  10.D0, 12.20,0.079D0,0.25D0,0.D0,
3057      &         17.35D0,-9.05D0,0.50D0,0.9D0,
3058      &  10.D0, 12.20,0.079D0,0.25D0,0.D0,
3059      &         17.35D0,9.05D0,0.50D0,0.9D0  /
3060 
3061       DATA BURQ83 /
3062      &  8.557D0,  0.00D0, 0.574D0,
3063      &  11.13D0,  7.23D0, 0.30D0,
3064      &  9.11D0,  -0.73D0, 0.28D0,
3065      &  9.11D0,   0.65D0, 0.28D0,
3066      &  8.55D0,  -5.98D0, 0.28D0,
3067      &  8.55D0,   1.60D0, 0.28D0  /
3068 
3069       DATA XMA / 2*0.93956563, 2*0.13956995, 2*0.493677 /
3070       DATA GEV2MB /0.389365/
3071       DATA PI /3.14159265358979/
3072 
3073 C  find index
3074       IF(L.eq.1) THEN
3075         K = 1                            ! p p
3076       ELSE IF(L.eq.2) THEN
3077         K = 3                            ! pi+ p
3078 *       K = 4                            ! pi- p
3079       ELSE IF(L.eq.3) THEN
3080         K = 5                            ! K+ p
3081 *       K = 6                            ! K- p
3082       ELSE
3083         GOTO 100
3084       ENDIF
3085 
3086 C  calculate lab momentum
3087       SS = ECM**2
3088       E1 = (SS-XMA(1)**2-XMA(K)**2)/(2.*XMA(1))
3089       PL = SQRT((E1-XMA(K))*(E1+XMA(K)))
3090       PLL = LOG(PL)
3091 
3092 C  check against lower limit
3093       IF(ECM.LE.XMA(1)+XMA(K)) GOTO 200
3094 
3095       XP  = TPDG96(2,K)*SS**TPDG96(3,K)
3096       YP  = TPDG96(6,K)/SS**TPDG96(8,K)
3097       YM  = TPDG96(7,K)/SS**TPDG96(8,K)
3098 
3099       PHR = TAN(PI/2.*(1.-TPDG96(8,K)))
3100       PHP = TAN(PI/2.*(1.+TPDG96(3,K)))
3101       RHO = (-YP/PHR + YM*PHR - XP/PHP)/(YP+YM+XP)
3102 
3103       SLOPE = BURQ83(1,K)+BURQ83(2,K)/SQRT(PL)+BURQ83(3,K)*PLL
3104 
3105 C  select energy range and interpolation method
3106       IF(PL.LT.TPDG96(1,K)) THEN
3107         SIGTOT = TPDG92(3,1,K)+TPDG92(4,1,K)*PL**TPDG92(5,1,K)
3108      &          + TPDG92(6,1,K)*PLL**2+TPDG92(7,1,K)*PLL
3109         SIGEL  = TPDG92(3,2,K)+TPDG92(4,2,K)*PL**TPDG92(5,2,K)
3110      &          + TPDG92(6,2,K)*PLL**2+TPDG92(7,2,K)*PLL
3111       ELSE IF(PL.LT.TPDG92(2,1,K)) THEN
3112         SIGTO1 = TPDG92(3,1,K)+TPDG92(4,1,K)*PL**TPDG92(5,1,K)
3113      &          + TPDG92(6,1,K)*PLL**2+TPDG92(7,1,K)*PLL
3114         SIGEL1 = TPDG92(3,2,K)+TPDG92(4,2,K)*PL**TPDG92(5,2,K)
3115      &          + TPDG92(6,2,K)*PLL**2+TPDG92(7,2,K)*PLL
3116         SIGTO2 = YP+YM+XP
3117         SIGEL2 = SIGTO2**2/(16.*PI*SLOPE*GEV2MB)*(1.+RHO**2)
3118         X2 = LOG(PL/TPDG96(1,K))/LOG(TPDG92(2,1,K)/TPDG96(1,K))
3119         X1 = 1. - X2
3120         SIGTOT = SIGTO2*X2 + SIGTO1*X1
3121         SIGEL  = SIGEL2*X2 + SIGEL1*X1
3122       ELSE
3123         SIGTOT = YP+YM+XP
3124         SIGEL  = SIGTOT**2/(16.*PI*SLOPE*GEV2MB)*(1.+RHO**2)
3125       ENDIF
3126       SIGINEL = SIGTOT-SIGEL
3127 
3128       RETURN
3129 
3130  100  CONTINUE
3131         WRITE(6,'(1X,2A,2I7)') 'SIB_HADCSL: ',
3132      &    'invalid beam particle: ',L
3133         RETURN
3134 
3135  200  CONTINUE
3136         WRITE(6,'(1X,2A,1P,E12.4)') 'SIB_HADCSL: ',
3137      &    'energy too small (Ecm): ',ECM
3138 
3139       RETURN
3140       END
3141 
3142 
3143       SUBROUTINE SIB_HADCS2(SQS,Xi_min,Xi_max,SIGeff,SIGDIF)
3144 C-----------------------------------------------------------------------
3145 C   cross section for diffraction dissociation
3146 C
3147 C   - single diffraction dissociation:
3148 C     Goulianos' parametrization (Ref: PL B358 (1995) 379)
3149 C   - double diffration dissociation: simple scaling model using
3150 C     single diff. cross section
3151 C
3152 C     in addition rescaling for different particles is applied using
3153 C     internal rescaling tables (not implemented yet)
3154 C
3155 C     input:     SQS         c.m. energy (GeV)
3156 C                Xi_min      min. diff mass (squared) = Xi_min*SQS**2
3157 C                Xi_max      max. diff mass (squared) = Xi_max*SQS**2
3158 C                SIGeff      effective cross section for DD scaling
3159 C
3160 C     output:    sig_sd1     cross section for diss. of particle 1 (mb)
3161 C                sig_sd2     cross section for diss. of particle 2 (mb)
3162 C                sig_dd      cross section for diss. of both particles
3163 C-----------------------------------------------------------------------
3164 
3165       DIMENSION SIGDIF(3)
3166       DOUBLE PRECISION Xpos1(96),Xwgh1(96),Xpos2(96),Xwgh2(96)
3167       DOUBLE PRECISION xil,xiu,tl,tu
3168       SAVE
3169 
3170 C  model parameters
3171       DATA delta    / 0.104 /
3172       DATA alphap   / 0.25 /
3173       DATA beta0    / 6.56 /
3174       DATA gpom0    / 1.21 /
3175       DATA xm_p     / 0.938 /
3176       DATA x_rad2   / 0.71 /
3177 
3178 C  integration precision
3179       DATA Ngau1    / 32 /
3180       DATA Ngau2    / 32 /
3181 
3182       DATA PI /3.14159265358979/
3183       DATA GEV2MB /0.389365/
3184 
3185 
3186       SIGDIF(1) = 0.
3187       SIGDIF(2) = 0.
3188       SIGDIF(3) = 0.
3189 
3190       XIL = LOG(Xi_min)
3191       XIU = LOG(Xi_max)
3192 
3193       if(XIL.ge.XIU) return
3194 
3195       SS = SQS*SQS
3196       xm4_p2 = 4.*xm_p**2
3197       fac = beta0**2/(16.*PI)
3198 
3199       t1 = -5.
3200       t2 = 0.
3201       tl = x_rad2/3./(1.-t1/x_rad2)**3
3202       tu = x_rad2/3./(1.-t2/x_rad2)**3
3203 
3204 C  flux renormalization and cross section for pp/ppbar case
3205 
3206       Xnorm  = 0.
3207 
3208       xil = log(1.5/SS)
3209       xiu = log(0.1)
3210 
3211       IF(xiu.LE.xil) goto 1000
3212 
3213       CALL SIB_GAUSET(xil,xiu,Ngau1,xpos1,xwgh1)
3214       CALL SIB_GAUSET(tl,tu,Ngau2,xpos2,xwgh2)
3215 
3216       do i1=1,Ngau1
3217 
3218         xi = exp(xpos1(i1))
3219         w_xi = Xwgh1(i1)
3220 
3221         do i2=1,Ngau2
3222 
3223           tt = x_rad2-x_rad2*(x_rad2/(3.*xpos2(i2)))**(1./3.)
3224 
3225           alpha_t =  1.+delta+alphap*tt
3226           f2_t = ((xm4_p2-2.8*tt)/(xm4_p2-tt))**2
3227 
3228           Xnorm = Xnorm
3229      &      + f2_t*xi**(2.-2.*alpha_t)*Xwgh2(i2)*w_xi
3230 
3231         enddo
3232       enddo
3233 
3234       Xnorm = Xnorm*fac
3235 
3236  1000 continue
3237 
3238       XIL = LOG(Xi_min)
3239       XIU = LOG(Xi_max)
3240 
3241       T1 = -5.
3242       T2 = 0.
3243 
3244       TL = x_rad2/3./(1.-t1/x_rad2)**3
3245       TU = x_rad2/3./(1.-t2/x_rad2)**3
3246 
3247 C  single diffraction diss. cross section
3248 
3249       CSdiff = 0.
3250 
3251       CALL SIB_GAUSET(XIL,XIU,NGAU1,XPOS1,XWGH1)
3252       CALL SIB_GAUSET(TL,TU,NGAU2,XPOS2,XWGH2)
3253 
3254       do i1=1,Ngau1
3255 
3256         xi = exp(xpos1(i1))
3257         w_xi = Xwgh1(i1)*beta0*gpom0*(xi*ss)**delta
3258 
3259         do i2=1,Ngau2
3260 
3261           tt = x_rad2-x_rad2*(x_rad2/(3.*xpos2(i2)))**(1./3.)
3262 
3263           alpha_t =  1.+delta+alphap*tt
3264           f2_t = ((xm4_p2-2.8*tt)/(xm4_p2-tt))**2
3265 
3266           CSdiff = CSdiff
3267      &      + f2_t*xi**(2.-2.*alpha_t)*Xwgh2(i2)*w_xi
3268 
3269         enddo
3270       enddo
3271 
3272       CSdiff = CSdiff*fac*GEV2MB/MAX(1.,Xnorm)
3273 
3274 *     write(6,'(1x,1p,4e14.3)')
3275 *    &  sqrt(SS),Xnorm,2.*CSdiff*MAX(1.,Xnorm),2.*CSdiff
3276 
3277       SIGDIF(1) = CSdiff
3278       SIGDIF(2) = CSdiff
3279 
3280 C  double diff. dissociation from simple probability consideration
3281 *     Pdiff = 0.5-sqrt(0.25-CSdiff/SIGeff)
3282       Pdiff = CSdiff/SIGeff
3283       SIGDIF(3) = Pdiff*Pdiff*SIGeff
3284 
3285       RETURN
3286       END
3287 
3288 
3289       SUBROUTINE DECSIB
3290 C-----------------------------------------------------------------------
3291 C...Decay all unstable particle in Sibyll
3292 C.  decayed particle have the code increased by 10000
3293 C
3294 C   changed to allow for multiple calls to DECSIB in one event
3295 C-----------------------------------------------------------------------
3296       COMMON /S_CSYDEC/ CBR(102), KDEC(612), LBARP(49), IDB(49)
3297       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
3298       COMMON /S_PLIST1/ LLIST1(8000)
3299       COMMON /S_MASS1/ AM(49), AM2(49)
3300       DIMENSION P0(5), LL(10), PD(10,5)
3301       SAVE
3302 
3303       NN = 1
3304       DO J=1,NP
3305          LLIST1(J) = 0
3306       ENDDO
3307       DO WHILE (NN .LE. NP)
3308          L= LLIST(NN)
3309          LA = IABS(L)
3310          if(LA.lt.50) then
3311            IF (IDB(LA) .GT. 0)  THEN
3312               DO K=1,5
3313                 P0(K) = P(NN,K)
3314               ENDDO
3315               CALL DECPAR (L,P0,ND,LL,PD)
3316               LLIST(NN) = LLIST(NN)+ISIGN(10000,LLIST(NN))
3317               DO J=1,ND
3318                 NP = NP+1
3319                 if(NP.gt.8000) then
3320                   write(6,'(1x,a,2i8)')
3321      &              'DECSIB: no space left in S_PLIST (NP,ND):',NP,ND
3322                   NP = NP-1
3323                   return
3324                 endif
3325                 DO K=1,5
3326                   P(NP,K) = PD(J,K)
3327                 ENDDO
3328                 LLIST(NP)=LL(J)
3329                 LLIST1(NP)=NN
3330               ENDDO
3331            ENDIF
3332          endif
3333          NN = NN+1
3334       ENDDO
3335 
3336       RETURN
3337       END
3338 
3339 
3340 
3341       SUBROUTINE DECPAR (LA,P0,ND,LL,P)
3342 C-----------------------------------------------------------------------
3343 C...This subroutine generates the decay of a particle
3344 C.  with ID = LA, and 5-momentum P0(1:5)
3345 C.  into ND particles of 5-momenta P(j,1:5) (j=1:ND)
3346 C.
3347 C.  If the initial particle code is LA=0
3348 C.  then ND and LL(1:ND) are considered as  input and
3349 C.  the routine generates a phase space decay into ND
3350 C.  particles of codes LL(1:nd)
3351 C.
3352 C.  june 1992
3353 C.  This version  contains the decay of polarized muons
3354 C.  The muon codes are  L =  4 : mu+ R
3355 C.                          -4 : mu+ L
3356 C.                           5 : mu- L
3357 C.                          -5 : mu- R
3358 C-----------------------------------------------------------------------
3359       COMMON /S_CSYDEC/ CBR(102), KDEC(612), LBARP(49), IDB(49)
3360       COMMON /S_MASS1/ AM(49), AM2(49)
3361       DIMENSION P0(5), LL(10), P(10,5)
3362       DIMENSION PV(10,5), RORD(10), UE(3),BE(3), FACN(3:10)
3363       SAVE
3364       DATA FACN /2.,5.,15.,60.,250.,1500.,12000.,120000./
3365       DATA PI /3.1415926/
3366 
3367 C...c.m.s. Momentum in two particle decays
3368 cdh  corrected 25.4.02
3369       PAWT(A,B,C) = SQRT((A**2-(B+C)**2+1.e-5)*(A**2-(B-C)**2))/(2.*A)
3370 
3371 C...Phase space decay into the particles in the list
3372       IF (LA .EQ. 0)  THEN
3373           MAT = 0
3374           MBST = 0
3375           PS = 0.
3376           DO J=1,ND
3377 CDH          following statements corrected by D.H. dec 20.,1995
3378              P (J,5) = AM(IABS(LL(J)))
3379              PV(J,5) = AM(IABS(LL(J)))
3380              PS = PS+P(J,5)
3381           ENDDO
3382           DO J=1,4
3383              PV(1,J) = P0(J)
3384           ENDDO
3385           PV(1,5) = P0(5)
3386           GOTO 140
3387       ENDIF
3388 
3389 C...Choose decay channel
3390       L = IABS(LA)
3391       ND=0
3392       IDC = IDB(L)-1
3393       IF (IDC+1 .LE.0)  RETURN
3394       RBR = S_RNDM(0)
3395 110   IDC=IDC+1
3396       IF(RBR.GT.CBR(IDC))  GOTO 110
3397 
3398       KD =6*(IDC-1)+1
3399       ND = KDEC(KD)
3400       MAT= KDEC(KD+1)
3401       MBST=0
3402       IF (MAT .GT.0 .AND. P0(4) .GT. 20*P0(5)) MBST=1
3403       IF (MAT .GT.0 .AND. MBST .EQ. 0)
3404      +        BETA = SQRT(P0(1)**2+P0(2)**2+P0(3)**2)/P0(4)
3405       PS = 0.
3406       DO J=1,ND
3407          LL(J) = KDEC(KD+1+J)
3408          P(J,5)  = AM(LL(J))
3409          PV(J,5) = AM(LL(J))
3410          PS = PS + P(J,5)
3411       ENDDO
3412       DO J=1,4
3413          PV(1,J) = 0.
3414          IF (MBST .EQ. 0)  PV(1,J) = P0(J)
3415       ENDDO
3416       IF (MBST .EQ. 1)  PV(1,4) = P0(5)
3417       PV(1,5) = P0(5)
3418 
3419 140   IF (ND .EQ. 2) GOTO 280
3420 
3421       IF (ND .EQ. 1)  THEN
3422          DO J=1,4
3423             P(1,J) = P0(J)
3424          ENDDO
3425          RETURN
3426       ENDIF
3427 
3428 C...Calculate maximum weight for ND-particle decay
3429       WWTMAX = 1./FACN(ND)
3430       PMAX=PV(1,5)-PS+P(ND,5)
3431       PMIN=0.
3432       DO IL=ND-1,1,-1
3433          PMAX = PMAX+P(IL,5)
3434          PMIN = PMIN+P(IL+1,5)
3435          WWTMAX = WWTMAX*PAWT(PMAX,PMIN,P(IL,5))
3436       ENDDO
3437 
3438 C...generation of the masses, compute weight, if rejected try again
3439 240   RORD(1) = 1.
3440       DO 260 IL1=2,ND-1
3441       RSAV = S_RNDM(il1)
3442       DO 250 IL2=IL1-1,1,-1
3443       IF(RSAV.LE.RORD(IL2))   GOTO 260
3444 250     RORD(IL2+1)=RORD(IL2)
3445 260     RORD(IL2+1)=RSAV
3446       RORD(ND) = 0.
3447       WT = 1.
3448       DO 270 IL=ND-1,1,-1
3449       PV(IL,5)=PV(IL+1,5)+P(IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
3450 270   WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(IL,5))
3451       IF (WT.LT.S_RNDM(0)*WWTMAX)   GOTO 240
3452 
3453 C...Perform two particle decays in respective cm frame
3454 280   DO 300 IL=1,ND-1
3455       PA=PAWT(PV(IL,5),PV(IL+1,5),P(IL,5))
3456       UE(3)=2.*S_RNDM(il)-1.
3457       PHI=2.*PI*S_RNDM(il)
3458       UT = SQRT(1.-UE(3)**2)
3459       UE(1) = UT*COS(PHI)
3460       UE(2) = UT*SIN(PHI)
3461       DO 290 J=1,3
3462       P(IL,J)=PA*UE(J)
3463 290   PV(IL+1,J)=-PA*UE(J)
3464       P(IL,4)=SQRT(PA**2+P(IL,5)**2)
3465 300   PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
3466 
3467 C...Lorentz transform decay products to lab frame
3468       DO 310 J=1,4
3469 310   P(ND,J)=PV(ND,J)
3470       DO 340 IL=ND-1,1,-1
3471       DO 320 J=1,3
3472 320   BE(J)=PV(IL,J)/PV(IL,4)
3473       GA=PV(IL,4)/PV(IL,5)
3474       DO 340 I=IL,ND
3475       BEP = BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
3476       DO 330 J=1,3
3477 330   P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
3478 340   P(I,4)=GA*(P(I,4)+BEP)
3479 
3480 C...Weak decays
3481       IF (MAT .EQ. 1)  THEN
3482          F1=P(2,4)*P(3,4)-P(2,1)*P(3,1)-P(2,2)*P(3,2)-P(2,3)*P(3,3)
3483          IF (MBST.EQ.1)  THEN
3484 C        WT = P0(5)*P(1,4)*F1
3485          WT = P0(5)*(P(1,4)+FLOAT(LA/L)*P(1,3))*F1
3486        ENDIF
3487          IF (MBST.EQ.0)  THEN
3488          WT=F1*(P(1,4)*P0(4)-P(1,1)*P0(1)-P(1,2)*P0(2)-P(1,3)*P0(3))
3489          WT= WT-FLOAT(LA/L)*(P0(4)*BETA*P(1,4)-P0(4)*P(1,3))*F1
3490        ENDIF
3491          WTMAX = P0(5)**4/8.
3492          IF(WT.LT.S_RNDM(0)*WTMAX)   GOTO 240
3493       ENDIF
3494 
3495 C...Boost back for rapidly moving particle
3496       IF (MBST .EQ. 1)   THEN
3497          DO 440 J=1,3
3498 440      BE(J)=P0(J)/P0(4)
3499          GA= P0(4)/P0(5)
3500          DO 460 I=1,ND
3501          BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
3502          DO 450 J=1,3
3503 450         P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
3504 460         P(I,4)=GA*(P(I,4)+BEP)
3505       ENDIF
3506 
3507 C...labels for antiparticle decay
3508       IF (LA .LT. 0 .AND. L .GT. 18)  THEN
3509            DO J=1,ND
3510             LL(J) = LBARP(LL(J))
3511          ENDDO
3512       ENDIF
3513 
3514       RETURN
3515       END
3516 
3517 
3518 
3519       BLOCK DATA DATDEC
3520 C-----------------------------------------------------------------------
3521 C...initialization of SIBYLL particle data
3522 C-----------------------------------------------------------------------
3523 
3524       COMMON /S_CSYDEC/ CBR(102), KDEC(612), LBARP(49), IDB(49)
3525       COMMON /S_MASS1/ AM(49), AM2(49)
3526       COMMON /S_CHP/ ICHP(49), ISTR(49), IBAR(49)
3527       COMMON /S_CNAM/ NAMP (0:49)
3528       CHARACTER NAMP*6
3529       SAVE
3530       DATA CBR /3*1.,0.,1.,1.,0.6351,0.8468,0.9027,0.9200,0.9518,1.,
3531      +   0.6351,0.8468,0.9027,0.9200,0.9518,1.,0.2160,0.3398,0.4748,
3532      +   0.6098,0.8049,1.,0.6861,1.,3*0.,0.5,1.,0.5,1.,
3533      +   0.3890,0.7080,0.9440,0.9930,1.,0.,0.4420,0.6470,0.9470,0.9770,
3534      +   0.9990,4*1.,0.6670,1.,9*0.,0.6670,1.,0.6670,1.,0.6670,1.,
3535      +   0.8880,0.9730,1.,0.4950,0.8390,0.9870,1.,0.5160,5*1.,0.6410,1.,
3536      +   1.,0.67,1.,0.33,1.,1.,0.88,0.94,1.,0.88,0.94,1.,0.88,0.94,1.,
3537      +   0.33,1.,0.67,1.,0.678,0.914,1./
3538       DATA AM / 0.,2*0.511E-3, 2*0.10566, 0.13497, 2*0.13957,
3539      +   2*0.49365, 2*0.49767, 0.93827, 0.93957, 4*0.,0.93827,
3540      +   0.93957, 2*0.49767, 0.54880,0.95750,2*0.76830,0.76860,
3541      +   2*0.89183,2*0.89610,0.78195,1.01941,1.18937,1.19255,
3542      +   1.19743,1.31490,1.32132,1.11563,1.23100,1.23500,
3543      +   1.23400,1.23300,1.38280,1.38370,1.38720,
3544      +   1.53180,1.53500,1.67243 /
3545       DATA AM2 /0.,2*2.61121E-07,2*0.011164,0.018217,0.019480,
3546      + 0.019480,0.243690,0.243690,0.247675,0.247675,0.880351,0.882792,
3547      + 0.000000,0.000000,0.000000,0.000000,0.880351,0.882792,0.247675,
3548      + 0.247675,0.301181,0.916806,0.590285,0.590285,0.590746,0.795361,
3549      + 0.795361,0.802995,0.802995,0.611446,1.039197,1.414601,1.422176,
3550      + 1.433839,1.728962,1.745887,1.244630,1.515361,1.525225,1.522765,
3551      + 1.520289,1.912136,1.914626,1.924324,2.346411,2.356225,2.797022/
3552       DATA IDB /
3553      +    0,0,0,1,2,3,5,6,7,13,19,25,8*0,30,32,34,40,46,47,48,49,60,62,
3554      +    64,66,69,73,75,76,77,78,79,81,82,84,86,87,90,93,96,98,100/
3555       DATA KDEC /
3556      + 3,1,15,2,18,0,3,1,16,3,17,0,2,0,1,1,8*0,2,0,4,17,0,0,2,0,5,18,0,
3557      + 0,2,0,4,17,0,0,2,0,7,6,0,0,3,0,7,7,8,0,3,0,7,6,6,0,3,1,17,4,6,0,
3558      + 3,1,15,2,6,0,2,0,5,18,0,0,2,0,8,6,0,0,3,0,8,8,7,0,3,0,8,6,6,0,3,
3559      + 1,18,5,6,0,3,1,16,3,6,0,3,0,6,6,6,0,3,0,7,8,6,0,3,1,18,5,7,0,3,
3560      + 1,17,4,8,0,3,1,16,3,7,0,3,1,15,2,8,0,2,0,7,8,0,0,2,0,6,6,20*0,1,
3561      + 0,11,3*0,1,0,12,0,0,0,1,0,11,0,0,0,1,0,12,0,0,0,2,0,1,1,0,0,3,0,
3562      + 6,6,6,0,3,0,7,8,6,0,3,0,1,7,8,0,3,0,1,3,2,7*0,3,0,7,8,23,0,3,0,6
3563      + ,6,23,0,2,0,1,27,0,0,2,0,1,32,0,0,2,0,1,1,0,0,3,0,6,6,6,0,2,0,7,
3564      + 6,0,0,2,0,8,6,0,0,2,0,7,8,0,0,2,0,21,7,0,0,2,0,9,6,0,0,54*0,2,0,
3565      + 22,8,0,0,2,0,10,6,0,0,2,0,9,8,0,0,2,0,21,6,0,0,2,0,10,7,0,0,
3566      + 2,0,22,6,0,0,3,0,7,8,6,0,2,0,1,6,0,0,2,0,7,8,0,0,2,0,9,10,0,
3567      + 0,2,0,11,12,0,0,3,0,7,
3568      + 8,6,0,2,0,1,23,0,0,2,0,13,6,0,0,2,0,14,7,0,0,2,0,39,1,0,0,2,
3569      + 0,14,8,0,0,2,0,39,6,0,0,2,0,39,8,0,0,2,0,13,8,0,0,2,0,
3570      + 14,6,0,0,2,0,13,7,0,0,2,0,13,6,
3571      + 0,0,2,0,14,7,0,0,2,0,13,8,0,0,2,0,14,6,0,0,2,0,14,8,0,0,2,0,
3572      + 39,7,0,0,2,0,34,6,0,0,2,0,35,7,0,0,2,0,39,6,0,0,2,0,34,8,0,0,
3573      + 2,0,36,7,0,0,2,0,39,8,0,0,2,
3574      + 0,35,8,0,0,2,0,36,6,0,0,2,0,37,6,0,0,2,0,38,7,0,0,2,0,
3575      + 37,8,0,0,2,0,38,6,0,0,2,0,39,10,0,0,2,0,37,8,0,0,2,0,38,6,0,0/
3576       DATA LBARP/1,3,2,5,4,6,8,7,10,9,11,12,-13,-14,16,15,18,17,13,14,
3577      +  22,21,23,24,26,25,27,29,28,31,30,32,33,-34,-35,-36,-37,-38,-39,
3578      +  -40,-41,-42,-43,-44,-45,-46,-47,-48,-49/
3579       DATA ICHP /0,1,-1,1,-1,0,1,-1,1,-1,0,0,1,0,4*0,-1,0,4*0,
3580      +    1,-1,0,1,-1,4*0,1,0,-1,0,-1,0,2,1,0,-1,1,0,-1,0,-1,-1/
3581       DATA ISTR /8*0,-1,+1,10,10,8*0,-1,+1,5*0,-1,+1,-1,+1,2*0,
3582      +           3*1,2*2,1,4*0,3*1,2*2,3 /
3583       DATA IBAR /12*0,2*1,4*0,2*-1,13*0,16*1/
3584       DATA NAMP /
3585      +     '     ','gam   ','e+','e-','mu+','mu-','pi0',
3586      +     'pi+','pi-','k+', 'k-', 'k0l','k0s',
3587      +     'p', 'n', 'nue', 'nueb', 'num', 'numb', 'pbar', 'nbar',
3588      +     'k0', 'k0b', 'eta', 'etap', 'rho+', 'rho-','rho0',
3589      +     'k*+','k*-','k*0','k*0b','omeg', 'phi', 'SIG+', 'SIG0',
3590      +     'SIG-','XI0','XI-','LAM','DELT++','DELT+','DELT0','DELT-',
3591      +     'SIG*+ ','SIG*0','SIG*-', 'XI*0', 'XI*-', 'OME*-'/
3592       END
3593 C->
3594       SUBROUTINE DECPR (LUN)
3595 C...Print on unit LUN the list of particles and decay channels
3596       COMMON /S_CSYDEC/ CBR(102), KDEC(612), LBARP(49), IDB(49)
3597       COMMON /S_MASS1/ AM(49), AM2(49)
3598       COMMON /S_CNAM/ NAMP (0:49)
3599       CHARACTER*6 NAMP
3600       DIMENSION LL(3)
3601       SAVE
3602 
3603       DO L=1,49
3604          IDC = IDB(L)-1
3605          NC = 0
3606          WRITE (LUN,10) L,NAMP(L), AM(L)
3607          IF(IDC+1 .GT. 0)  THEN
3608             CB = 0.
3609 110         IDC=IDC+1
3610             NC = NC+1
3611             CBOLD = CB
3612             CB = CBR(IDC)
3613             BR = CB-CBOLD
3614             KD = 6*(IDC-1)+1
3615             ND = KDEC(KD)
3616             MAT= KDEC(KD+1)
3617             DO J=1,ND
3618               LL(J) = KDEC(KD+1+J)
3619             ENDDO
3620             WRITE (LUN,15) NC,BR,ND,MAT, (NAMP(LL(J)),J=1,ND)
3621             IF (CB .LT. 1.)  GOTO 110
3622          ENDIF
3623       ENDDO
3624       RETURN
3625 10    FORMAT(1X,I3,2X,A6,3X,F10.4)
3626 15    FORMAT(5X,I2,2X,F9.4,I4,I4,2X,3(A6,2X))
3627       END
3628 
3629 
3630 
3631       SUBROUTINE DEC_DEBUG (L,P0, ND, LL, PD)
3632       COMMON /S_CNAM/ NAMP (0:49)
3633       CHARACTER*6 NAMP
3634       DIMENSION P0(5), LL(10), PD(10,5)
3635       SAVE
3636 
3637       ETOT = 0.
3638       DO J=1,ND
3639          ETOT = ETOT + PD(J,4)
3640       ENDDO
3641       WRITE(*,*)  NAMP(IABS(L)),' -> ', (NAMP(IABS(LL(J))),J=1,ND)
3642       WRITE(*,*)  ' Ei, Ef = ', P0(4), ETOT, ' L = ', L
3643       RETURN
3644       END
3645 
3646 
3647 
3648       SUBROUTINE SIB_SIG(Jint,SIB_SQS,SIB_PTmin,SIB_SIG_tot,
3649      &                   SIB_SIG_ine,SIB_diff,SIB_B_el,SIB_PJET)
3650 C-----------------------------------------------------------------------
3651 C
3652 C...SIBYLL 2.1 cross sections
3653 C
3654 C   input parameter: SIB_SQS   c.m.s. energy (GeV)
3655 C                    Jint      1 p-p cross sections
3656 C                              2 pi-p cross sections
3657 C
3658 C-----------------------------------------------------------------------
3659       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3660 
3661       PARAMETER (NS_max = 20, NH_max = 50)
3662       REAL SIB_PJET(0:NS_max,0:NH_max)
3663       REAL SIB_SQS,SIB_PTmin,
3664      &     SIB_SIG_ine,SIB_SIG_tot,SIB_diff(3),SIB_B_el
3665 
3666       COMMON /S_CFLAFR/ PAR(20), IPAR(10)
3667       REAL PAR
3668 
3669 
3670       COMMON /SIGMAS/SQS,SIGTOT,SIGEL,SIGINE,
3671      &               SIGSD1(2),SIGSD2(2),SIGDD(2),
3672      &               SLOPE,SLOPEc,RHO,PROB(0:NS_max,0:NH_max),SIGSUM
3673 C
3674       COMMON /PROFILE/XNUS2,XMUS2,XNUSPI2,
3675      &                XNUH2,XMUH2,XNUHPI2,
3676      &                ENHPP,ENHPIP,al1,be1,al2,be2
3677 C
3678       COMMON /S_CHDCNV/ABR(2,400),ABP(2,400),ABH(2,400),DB,NB
3679 
3680       DIMENSION XI(20)
3681 C
3682       DIMENSION SIG_BRN(3)
3683       DIMENSION SIG_dif_1(2),SIG_dif_2(2),SIG_dd(2)
3684 
3685       DIMENSION IHAR(2),SIGQCD(61,2)
3686       SAVE
3687 
3688       DATA (SIGQCD(K,1),K=    1,   61) /
3689      &8.4925E-02,1.8301E-01,3.4031E-01,5.7346E-01,9.0097E-01,1.3443E+00,
3690      &1.9284E+00,2.6830E+00,3.6426E+00,4.8474E+00,6.3440E+00,8.1867E+00,
3691      &1.0438E+01,1.3169E+01,1.6462E+01,2.0412E+01,2.5128E+01,3.0733E+01,
3692      &3.7368E+01,4.5196E+01,5.4401E+01,6.5191E+01,7.7807E+01,9.2520E+01,
3693      &1.0964E+02,1.2951E+02,1.5254E+02,1.7918E+02,2.0993E+02,2.4538E+02,
3694      &2.8618E+02,3.3307E+02,3.8689E+02,4.4859E+02,5.1923E+02,6.0003E+02,
3695      &6.9234E+02,7.9769E+02,9.1782E+02,1.0547E+03,1.2104E+03,1.3876E+03,
3696      &1.5888E+03,1.8174E+03,2.0767E+03,2.3707E+03,2.7038E+03,3.0810E+03,
3697      &3.5077E+03,3.9903E+03,4.5357E+03,5.1517E+03,5.8471E+03,6.6317E+03,
3698      &7.5163E+03,8.5134E+03,9.6365E+03,1.0901E+04,1.2324E+04,1.3925E+04,
3699      &1.5724E+04/
3700 
3701       DATA (SIGQCD(K,2),K=    1,   61) /
3702      &1.5891E-01,2.9085E-01,4.8190E-01,7.4586E-01,1.0985E+00,1.5582E+00,
3703      &2.1466E+00,2.8890E+00,3.8146E+00,4.9572E+00,6.3556E+00,8.0544E+00,
3704      &1.0104E+01,1.2564E+01,1.5500E+01,1.8987E+01,2.3112E+01,2.7972E+01,
3705      &3.3679E+01,4.0359E+01,4.8154E+01,5.7228E+01,6.7762E+01,7.9965E+01,
3706      &9.4071E+01,1.1035E+02,1.2909E+02,1.5063E+02,1.7536E+02,2.0370E+02,
3707      &2.3613E+02,2.7321E+02,3.1553E+02,3.6379E+02,4.1875E+02,4.8129E+02,
3708      &5.5238E+02,6.3311E+02,7.2470E+02,8.2854E+02,9.4614E+02,1.0793E+03,
3709      &1.2298E+03,1.3999E+03,1.5920E+03,1.8089E+03,2.0534E+03,2.3291E+03,
3710      &2.6396E+03,2.9892E+03,3.3825E+03,3.8248E+03,4.3219E+03,4.8803E+03,
3711      &5.5073E+03,6.2109E+03,7.0001E+03,7.8849E+03,8.8764E+03,9.9871E+03,
3712      &1.1231E+04/
3713 
3714 
3715       DATA CMBARN /0.389385/
3716       DATA PI /3.1415926/
3717       DATA INIT /0/
3718 
3719       IF(INIT.EQ.0) THEN
3720 *        CALL HAR_INI
3721         CALL FACT_INI
3722         IHAR(1) = 0
3723         IHAR(2) = 0
3724         INIT = 1
3725       ENDIF
3726 
3727       ECM = SIB_SQS
3728 
3729       IF(JINT.EQ.1) THEN
3730 
3731         XI( 1) =  4.989E+01
3732         XI( 2) =  8.203E-05
3733         XI( 3) =  2.449E-02
3734         XI( 4) = -4.000E-01
3735         XI( 5) =  2.000E-01
3736         XI( 6) =  5.000E-01
3737         XI( 7) =  0.000E+00
3738         XI( 8) =  6.000E-01
3739         XI( 9) =  9.000E-02
3740         XI(10) =  1.000E+00
3741         XI(11) =  2.000E+00
3742         XI(12) =  3.175E+00
3743         XI(13) =  2.500E-01
3744         XI(14) =  5.400E-01
3745         XI(15) =  7.700E-01
3746         XI(16) = -8.800E-01
3747         XI(17) =  5.400E-01
3748         XI(18) =  5.000E-01
3749         XI(19) =  9.000E-01
3750 
3751       ELSE IF(JINT.EQ.2) THEN
3752 
3753         XI( 1) = 2.653E+01
3754         XI( 2) = 2.707E+01
3755         XI( 3) = 2.449E-02
3756         XI( 4) =-4.000E-01
3757         XI( 5) = 2.000E-01
3758         XI( 6) = 5.000E-01
3759         XI( 7) = 0.000E+00
3760         XI( 8) = 6.000E-01
3761         XI( 9) = 9.000E-02
3762         XI(10) = 1.000E+00
3763         XI(11) = 2.000E+00
3764         XI(12) = 1.216E+00
3765         XI(13) = 2.500E-01
3766         XI(14) = 5.400E-01
3767         XI(15) = 7.700E-01
3768         XI(16) =-8.800E-01
3769         XI(17) = 5.400E-01
3770         XI(18) = 8.640E+00
3771         XI(19) = 9.000E-01
3772 
3773       ENDIF
3774 
3775       XNUS2   = XI(12)
3776       XMUS2   = XI(13)
3777       XNUSPI2 = XI(14)
3778 
3779       XNUH2   = XI(15)
3780       XMUH2   = XI(16)
3781       XNUHPI2 = XI(17)
3782 
3783       CALL HAD_CONV(IABS(JINT))
3784 
3785       PTCUT = XI(10)+0.065D0*EXP(0.9D0*SQRT(2.D0*LOG(ECM)))
3786       INDX = abs(JINT)
3787       IHAR(INDX) = IHAR(INDX)+1
3788       SIGHAR = SIGQCD(IHAR(INDX),INDX)
3789 
3790       S = ECM**2
3791 
3792       BREG =  ABS(XI(18)) + XI(19)*LOG(S)
3793       BPOM =  ABS(XI(12)) + XI(13)*LOG(S)
3794       IK = ABS(JINT)
3795       DO JB=1,NB
3796         B = DB*FLOAT(JB-1)
3797         ABR(IK,JB) = 2./(8.*PI*BREG)*EXP(-B**2/(4.*BREG))
3798         ABP(IK,JB) = 2./(8.*PI*BPOM)*EXP(-B**2/(4.*BPOM))
3799       ENDDO
3800 
3801 C  reggeon
3802       SIGSR = ABS(XI(2))*S**(-ABS(XI(4)))
3803       SIG_BRN(1) = SIGSR/CMBARN
3804 C  pomeron (soft part)
3805       SIGSP = ABS(XI(1))*S**ABS(XI(3))
3806       SIG_BRN(2) = SIGSP/CMBARN
3807 C  pomeron (hard part)
3808       SIG_BRN(3) = SIGHAR/CMBARN
3809 
3810 C  2x2 channel low-mass model and separate high-mass diffraction
3811 
3812       al1 = XI(5)
3813       be1 = XI(6)
3814       al2 = al1
3815       be2 = be1
3816       EnhPP  = XI(9)
3817       EnhPiP = EnhPP
3818 
3819       CALL SIG_JET_3 (SIG_brn,JINT,SIG_tot,SIG_ela,SIG_ine,SIG_sum,
3820      &                SIG_dif_1,SIG_dif_2,SIG_dd,B_el,PROB)
3821 
3822       SIGTOT = SIG_tot*CMBARN
3823       SIGINE = SIG_ine*CMBARN
3824       SIGSUM = SIG_sum*CMBARN
3825       SIGELc = SIGTOT-SIGINE
3826       SIGEL  = SIG_ela*CMBARN
3827       SIGSD1(1) = SIG_dif_1(1)*CMBARN
3828       SIGSD1(2) = SIG_dif_1(2)*CMBARN
3829       SIGSD2(1) = SIG_dif_2(1)*CMBARN
3830       SIGSD2(2) = SIG_dif_2(2)*CMBARN
3831       SIGDD(1)  = SIG_dd(1)*CMBARN
3832       SIGDD(2)  = SIG_dd(2)*CMBARN
3833       SLOPE  = B_EL
3834       SLOPEc = SIG_tot**2/(16.*Pi*SIG_ela)
3835 
3836       DE = ABS(SIGEL+SIGINE-SIGTOT)/SIGTOT
3837       IF(DE.GT.0.01) THEN
3838         print *,'SIBSIG:      Ecm: ',ECM
3839         print *,'          SIGTOT: ',SIGTOT
3840         print *,'        SIGEL1/2: ',SIGEL,SIGELc
3841         print *,'        SLOPE1/2: ',SLOPE,SLOPEc
3842         print *,'        SIGDIF 1: ',SIGSD1
3843         print *,'        SIGDIF 2: ',SIGSD2
3844         print *,'         SIGDDIF: ',SIGDD
3845         print *,'      SUM-SIGTOT: ',SIGEL+SIGINE-SIGTOT
3846       ENDIF
3847 
3848 C  SIBYLL interface to single precision
3849 
3850       SIB_PTmin   = PTCUT
3851       SIB_SIG_tot = SIGTOT
3852       SIB_SIG_ine = SIGINE
3853       SIB_diff(1) = SIGSD1(1)+SIGSD1(2)
3854       SIB_diff(2) = SIGSD2(1)+SIGSD2(2)
3855       SIB_diff(3) = SIGDD(1)+SIGDD(2)
3856       SIB_B_el    = SLOPE
3857       DO I=0,NS_max
3858         DO K=0,NH_max
3859           SIB_PJET(I,K) = PROB(I,K)
3860         ENDDO
3861       ENDDO
3862 
3863       RETURN
3864       END
3865 
3866 
3867       SUBROUTINE SIG_JET_3 (SIG_brn, JINT, SIG_TOT, SIG_ELA,
3868      &        SIG_INE, SIG_sum, SIG_DIF1, SIG_DIF2, SIG_DD, B_EL, P_int)
3869 C-----------------------------------------------------------------------
3870 C
3871 C...This subroutine  receives in INPUT:
3872 C.       SIG_brn (GeV-2)  Born graph cross sections
3873 C.       JINT (1 = pp interaction)    (2 pi-p interaction)
3874 C.       neg. value: without calculation of interaction probabilities
3875 C.
3876 C.  and returns as output:
3877 C.       SIG_???  , B_el
3878 C.       and P_int(0:NS_max,0:NH_max)   interaction probabilities
3879 C
3880 C   two x two -channel approximation for diffraction
3881 C
3882 C-----------------------------------------------------------------------
3883       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3884 
3885       DIMENSION SIG_brn(3)
3886       PARAMETER (NS_max = 20, NH_max = 50)
3887 
3888       COMMON /S_CFACT/ FACT(0:NH_max), CO_BIN(0:NH_max,0:NH_max)
3889       COMMON /S_CHDCNV/ABR(2,400),ABP(2,400),ABH(2,400),DB,NB
3890 
3891       COMMON /PROFILE/XNUS2,XMUS2,XNUSPI2,
3892      &                XNUH2,XMUH2,XNUHPI2,
3893      &                EnhPP,EnhPiP,al1,be1,al2,be2
3894 
3895       DIMENSION SIG_DIF1(2),SIG_DIF2(2),SIG_DD(2),
3896      &          P_int(0:NS_max,0:NH_max)
3897       SAVE
3898       DATA PI /3.1415926/
3899 
3900       DO J=0,NH_max
3901         DO I=0,NS_max
3902           P_int(I,J) = 0.
3903         ENDDO
3904       ENDDO
3905 
3906       ga1 = sqrt(al1*al1+be1*be1)
3907       ga2 = sqrt(al2*al2+be2*be2)
3908 
3909       fe_a_1  = (1.+al1/ga1)/2.
3910       fe_a_2  = (1.-al1/ga1)/2.
3911       fd_a_1  = sqrt(1.-(al1/ga1)**2)/2.
3912       fd_a_2  = -fd_a_1
3913 
3914       fe_b_1  = (1.+al2/ga2)/2.
3915       fe_b_2  = (1.-al2/ga2)/2.
3916       fd_b_1  = sqrt(1.-(al2/ga2)**2)/2.
3917       fd_b_2  = -fd_b_1
3918 
3919       fe_11 = fe_a_1*fe_b_1
3920       fe_22 = fe_a_2*fe_b_2
3921       fe_12 = fe_a_1*fe_b_2
3922       fe_21 = fe_a_2*fe_b_1
3923 
3924       fd_a_11 = fd_a_1*fe_b_1
3925       fd_a_22 = fd_a_2*fe_b_2
3926       fd_a_12 = fd_a_1*fe_b_2
3927       fd_a_21 = fd_a_2*fe_b_1
3928 
3929       fd_b_11 = fe_a_1*fd_b_1
3930       fd_b_22 = fe_a_2*fd_b_2
3931       fd_b_12 = fe_a_1*fd_b_2
3932       fd_b_21 = fe_a_2*fd_b_1
3933 
3934       fdd_11 = fd_a_1*fd_b_1
3935       fdd_22 = fd_a_2*fd_b_2
3936       fdd_12 = fd_a_1*fd_b_2
3937       fdd_21 = fd_a_2*fd_b_1
3938 
3939 
3940       sum_abs = 0.
3941       sum_tot = 0.
3942       sum_ela = 0.
3943       sum_sd_a = 0.
3944       sum_sd_b = 0.
3945       sum_dd  = 0.
3946       sum_B   = 0.
3947 
3948       IK = ABS(JINT)
3949       if(JINT.GT.0) then
3950         I0MAX = NS_max
3951         J0MAX = NH_max
3952       ELSE
3953         I0MAX = 1
3954         J0MAX = 1
3955       ENDIF
3956       SIG_REG = SIG_BRN(1)
3957       SIG_POM = SIG_BRN(2)
3958       SIG_HAR = SIG_BRN(3)
3959 
3960       DO JB=1,NB
3961 
3962          B = DB*FLOAT(JB-1)
3963 
3964          ABREG = ABR(IK,JB)
3965          ABPOM = ABP(IK,JB)
3966          ABHAR = ABH(IK,JB)
3967 
3968          chi2_soft = ABREG*SIG_REG+ABPOM*SIG_POM
3969          chi2_soft_11 = (1.-al1+ga1)*(1.-al2+ga2)*chi2_soft
3970          chi2_soft_22 = (1.-al1-ga1)*(1.-al2-ga2)*chi2_soft
3971          chi2_soft_12 = (1.-al1+ga1)*(1.-al2-ga2)*chi2_soft
3972          chi2_soft_21 = (1.-al1-ga1)*(1.-al2+ga2)*chi2_soft
3973 
3974          chi2_hard = ABHAR*SIG_HAR
3975          chi2_hard_11 = (1.-al1+ga1)*(1.-al2+ga2)*chi2_hard
3976          chi2_hard_22 = (1.-al1-ga1)*(1.-al2-ga2)*chi2_hard
3977          chi2_hard_12 = (1.-al1+ga1)*(1.-al2-ga2)*chi2_hard
3978          chi2_hard_21 = (1.-al1-ga1)*(1.-al2+ga2)*chi2_hard
3979 
3980 
3981          ef_11  = exp(-0.5*(chi2_soft_11+chi2_hard_11))
3982          ef_22  = exp(-0.5*(chi2_soft_22+chi2_hard_22))
3983          ef_12 = exp(-0.5*(chi2_soft_12+chi2_hard_12))
3984          ef_21 = exp(-0.5*(chi2_soft_21+chi2_hard_21))
3985 
3986          esf_11  = ef_11**2
3987          esf_22  = ef_22**2
3988          esf_12  = ef_12**2
3989          esf_21  = ef_21**2
3990 
3991          F_ine = B*(1. - fe_11*esf_11 - fe_12*esf_12
3992      &                 - fe_21*esf_21 - fe_22*esf_22)
3993          F_tot = 1. - fe_11*ef_11 - fe_12*ef_12
3994      &              - fe_21*ef_21 - fe_22*ef_22
3995          F_ela = B*F_tot**2
3996          F_tot = B*F_tot
3997 
3998          F_sd_a = B*(fd_a_11*ef_11 + fd_a_12*ef_12
3999      &              + fd_a_21*ef_21 + fd_a_22*ef_22)**2
4000          F_sd_b = B*(fd_b_11*ef_11 + fd_b_12*ef_12
4001      &              + fd_b_21*ef_21 + fd_b_22*ef_22)**2
4002          F_dd  = B*(fdd_11*ef_11 + fdd_12*ef_12
4003      &              + fdd_21*ef_21 + fdd_22*ef_22)**2
4004 
4005          sum_abs = sum_abs+F_ine
4006          sum_tot = sum_tot+F_tot
4007          sum_ela = sum_ela+F_ela
4008 
4009          sum_sd_a = sum_sd_a+F_sd_a
4010          sum_sd_b = sum_sd_b+F_sd_b
4011          sum_dd  = sum_dd +F_dd
4012 
4013          sum_B   = sum_b+B**2*F_tot
4014 
4015          fac_11 = B*esf_11
4016          fac_22 = B*esf_22
4017          fac_12 = B*esf_12
4018          fac_21 = B*esf_21
4019          soft_rec_11 = 1./chi2_soft_11
4020          soft_rec_22 = 1./chi2_soft_22
4021          soft_rec_12 = 1./chi2_soft_12
4022          soft_rec_21 = 1./chi2_soft_21
4023          chi2_hard_11 = max(chi2_hard_11,1.d-10)
4024          chi2_hard_22 = max(chi2_hard_22,1.d-10)
4025          chi2_hard_12 = max(chi2_hard_12,1.d-10)
4026          chi2_hard_21 = max(chi2_hard_21,1.d-10)
4027          DO I=0,I0MAX
4028            soft_rec_11 = soft_rec_11*chi2_soft_11
4029            soft_rec_22 = soft_rec_22*chi2_soft_22
4030            soft_rec_12 = soft_rec_12*chi2_soft_12
4031            soft_rec_21 = soft_rec_21*chi2_soft_21
4032            hard_rec_11 = 1./chi2_hard_11
4033            hard_rec_22 = 1./chi2_hard_22
4034            hard_rec_12 = 1./chi2_hard_12
4035            hard_rec_21 = 1./chi2_hard_21
4036            DO J=0,J0MAX
4037              hard_rec_11 = hard_rec_11*chi2_hard_11
4038              hard_rec_22 = hard_rec_22*chi2_hard_22
4039              hard_rec_12 = hard_rec_12*chi2_hard_12
4040              hard_rec_21 = hard_rec_21*chi2_hard_21
4041              P_int(I,J) = P_int(I,J)
4042      &                + fe_11*soft_rec_11*hard_rec_11*fac_11
4043      &                + fe_22*soft_rec_22*hard_rec_22*fac_22
4044      &                + fe_12*soft_rec_12*hard_rec_12*fac_12
4045      &                + fe_21*soft_rec_21*hard_rec_21*fac_21
4046            ENDDO
4047          ENDDO
4048 
4049       ENDDO
4050 
4051       SIG_abs  = SUM_abs*2.*PI*DB
4052       SIG_tot  = SUM_tot*4.*PI*DB
4053       SIG_ela  = SUM_ela*2.*PI*DB
4054       SIG_dif1(1) = SUM_sd_a*2.*PI*DB
4055       SIG_dif2(1) = SUM_sd_b*2.*PI*DB
4056       SIG_dd(1)   = SUM_dd*2.*PI*DB
4057       SIG_ine  = SIG_abs + SIG_dif1(1) + SIG_dif2(1) + SIG_dd(1)
4058       B_EL     = sum_B/SUM_tot/2.
4059 
4060       SA = 0.
4061       P_int(0,0) = 0.
4062       DO I=0,I0MAX
4063         DO J=0,J0MAX
4064           fac = FACT(I)*FACT(J)
4065           P_int(I,J) = P_int(I,J)/fac
4066           SA = SA + P_int(I,J)
4067         ENDDO
4068       ENDDO
4069 
4070       SIG_hmsd = EnhPP*(P_int(1,0)+P_int(0,1))*2.*PI*DB
4071       SIG_hmdd = be1**2*SIG_hmsd + be2**2*SIG_hmsd
4072      &          + EnhPP**2*P_int(1,1)*2.*PI*DB
4073 
4074       SIG_dif1(2) = SIG_hmsd
4075       SIG_dif2(2) = SIG_hmsd
4076       SIG_dd(2)   = SIG_hmdd
4077 
4078       SIG_sum = SA*2.*PI*DB
4079 
4080       DO I=0,I0MAX
4081         DO J=0,J0MAX
4082           P_int(I,J) = P_int(I,J)/SA
4083         ENDDO
4084       ENDDO
4085 
4086       RETURN
4087       END
4088 C
4089 C
4090       SUBROUTINE HAD_CONV(JINT)
4091 C-----------------------------------------------------------------------
4092 C
4093 C...Convolution of hadrons profile
4094 C.  [function A(b) of Durand and Pi]
4095 C.  precalculate and put  in COMMON block
4096 C
4097 C-----------------------------------------------------------------------
4098       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4099 C
4100       COMMON /S_CHDCNV/ABR(2,400),ABP(2,400),ABH(2,400),DB,NB
4101 
4102       DOUBLE PRECISION NU2, MU2, NUPI2, NU, MU, NUPI
4103       COMMON /S_CH0CNV/ NU2, MU2, NUPI2, NU, MU, NUPI
4104 
4105 C
4106       COMMON /PROFILE/XNUS2,XMUS2,XNUSPI2,
4107      &                XNUH2,XMUH2,XNUHPI2,
4108      &                ENHPP,ENHPIP,al1,be1,al2,be2
4109       SAVE
4110 
4111 C...integration constants
4112       BMAX = 50.
4113       NB  = 400
4114       DB = BMAX/FLOAT(NB)
4115 
4116 C  soft reggeon interactions
4117 
4118       NU2   = XNUS2
4119       MU2   = XMUS2
4120       NUPI2 = XNUSPI2
4121 
4122       NU = SQRT(NU2)
4123       MU = SQRT(ABS(MU2))
4124       NUPI = SQRT(NUPI2)
4125 
4126       DO JB=1,NB
4127          B = DB*FLOAT(JB-1)
4128          IF(JINT.EQ.1) THEN
4129            ABR(JINT,JB) = A_PP(B)
4130          ELSE
4131            ABR(JINT,JB) = A_PIP(B)
4132          ENDIF
4133       ENDDO
4134 
4135 C  soft pomeron interactions
4136 
4137       NU2   = XNUS2
4138       MU2   = XMUS2
4139       NUPI2 = XNUSPI2
4140 
4141       NU = SQRT(NU2)
4142       MU = SQRT(ABS(MU2))
4143       NUPI = SQRT(NUPI2)
4144 
4145       DO JB=1,NB
4146          B = DB*FLOAT(JB-1)
4147          IF(JINT.EQ.1) THEN
4148            ABP(JINT,JB) = A_PP(B)
4149          ELSE
4150            ABP(JINT,JB) = A_PIP(B)
4151          ENDIF
4152       ENDDO
4153 
4154 C  hard pomeron interactions
4155 
4156       NU2   = XNUH2
4157       MU2   = XMUH2
4158       NUPI2 = XNUHPI2
4159 
4160       NU = SQRT(NU2)
4161       MU = SQRT(ABS(MU2))
4162       NUPI = SQRT(NUPI2)
4163 
4164       DB = BMAX/FLOAT(NB)
4165       DO JB=1,NB
4166          B = DB*FLOAT(JB-1)
4167          IF(JINT.EQ.1) THEN
4168            ABH(JINT,JB) = A_PP(B)
4169          ELSE
4170            ABH(JINT,JB) = A_PIP(B)
4171          ENDIF
4172       ENDDO
4173 
4174       RETURN
4175       END
4176 C
4177 C
4178       DOUBLE PRECISION FUNCTION A_pp (b)
4179 C...Convolution of parton distribution for pp interaction
4180       IMPLICIT DOUBLE PRECISION (A-Z)
4181 C
4182       DOUBLE PRECISION NU2, MU2, NUPI2, NU, MU, NUPI
4183       COMMON /S_CH0CNV/ NU2, MU2, NUPI2, NU, MU, NUPI
4184       SAVE
4185       data pi / 3.1415926/
4186 
4187       ETA = NU2/MU2
4188 
4189       IF(ETA.LT.0.D0) THEN
4190 
4191         c = nu**5/(96.*pi)
4192         if (b .gt. 0.0001D0)  then
4193            A_pp = c*b**3 * bessk (3, b*nu)
4194         else
4195            A_pp = nu**2/(12.*pi)
4196         endif
4197 
4198       ELSE
4199 
4200         X = B*NU
4201         Y = B*MU
4202         C = NU2/(12.*PI)/(1.-ETA)**2
4203         IF(X.GT.0.0001D0) THEN
4204           A_PP = C*(1./8.*X**3*BESSK(3,X)
4205      &          -3./2.*ETA/(1.-ETA)*X**2*BESSK(2,X)
4206      &          +9*ETA**2/(1.-ETA)**2*X*BESSK1(X)
4207      &          -24*ETA**3/(1.-ETA)**3*(BESSK0(X)-BESSK0(Y))
4208      &          +3.*ETA**3/(1.-ETA)**2*Y*BESSK1(Y))
4209         ELSE
4210           A_PP = C*(1./8.*8.
4211      &          -3./2.*ETA/(1.-ETA)*2.
4212      &          +9*ETA**2/(1.-ETA)**2*1.
4213      &          -24*ETA**3/(1.-ETA)**3*LOG(MU/NU)
4214      &          +3.*ETA**3/(1.-ETA)**2*1.)
4215         ENDIF
4216 
4217       ENDIF
4218 
4219       RETURN
4220       END
4221 C
4222 C
4223       DOUBLE PRECISION FUNCTION A_pip (b)
4224 C...Convolution of parton distribution for pip interaction
4225       IMPLICIT DOUBLE PRECISION (A-Z)
4226 C
4227       DOUBLE PRECISION NU2, MU2, NUPI2, NU, MU, NUPI
4228       COMMON /S_CH0CNV/ NU2, MU2, NUPI2, NU, MU, NUPI
4229       SAVE
4230       data pi / 3.1415926/
4231 
4232       eta = nu2/nupi2
4233       c = nu2/(2.*pi) * 1./(1.-eta)
4234 
4235       if (b .gt. 0.0001D0)  then
4236          b1 = b*nu
4237          b2 = b*nupi
4238          f1 = 0.5*b1 * bessk1(b1)
4239          f2 = eta/(1.-eta)*(bessk0(b2)- bessk0(b1))
4240          A_pip = c*(f1+f2)
4241       else
4242          A_pip = c*(0.5 + eta/(1.-eta)*log(nu/nupi))
4243       endif
4244       return
4245       end
4246 C
4247 C
4248 C----------------------------------------------------------------------------
4249 C  Bessel functions
4250 C----------------------------------------------------------------------------
4251 C
4252       FUNCTION BESSK0(X)
4253       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4254 *     REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,
4255 *    *    Q1,Q2,Q3,Q4,Q5,Q6,Q7
4256       SAVE
4257 C
4258       DATA P1,P2,P3,P4,P5,P6,P7/-0.57721566D0,0.42278420D0,
4259      *    0.23069756D0,0.3488590D-1,0.262698D-2,0.10750D-3,0.74D-5/
4260       DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,-0.7832358D-1,
4261      * 0.2189568D-1,-0.1062446D-1,0.587872D-2,-0.251540D-2,0.53208D-3/
4262 
4263       IF (X.LE.2.0) THEN
4264         Y=X*X/4.0
4265         BESSK0=(-LOG(X/2.0)*BESSI0(X))+(P1+Y*(P2+Y*(P3+
4266      *        Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
4267       ELSE
4268         Y=(2.0/X)
4269         BESSK0=(EXP(-X)/SQRT(X))*(Q1+Y*(Q2+Y*(Q3+
4270      *        Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
4271       ENDIF
4272       RETURN
4273       END
4274 C
4275 C
4276       FUNCTION BESSK1(X)
4277       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4278 C
4279 *     REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,
4280 *    *    Q1,Q2,Q3,Q4,Q5,Q6,Q7
4281       SAVE
4282       DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,0.15443144D0,-0.67278579D0,
4283      *    -0.18156897D0,-0.1919402D-1,-0.110404D-2,-0.4686D-4/
4284       DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,0.23498619D0,
4285      *    -0.3655620D-1,0.1504268D-1,-0.780353D-2,0.325614D-2,
4286      *    -0.68245D-3/
4287 
4288       IF (X.LE.2.0) THEN
4289         Y=X*X/4.0
4290         BESSK1=(LOG(X/2.0)*BESSI1(X))+(1.0/X)*(P1+Y*(P2+
4291      *      Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
4292       ELSE
4293         Y=2.0/X
4294         BESSK1=(EXP(-X)/SQRT(X))*(Q1+Y*(Q2+Y*(Q3+
4295      *      Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
4296       ENDIF
4297       RETURN
4298       END
4299 C
4300 C
4301       FUNCTION BESSK(N,X)
4302       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4303       SAVE
4304 C
4305       IF (N.LT.2) STOP 'bad argument N in BESSK'
4306       TOX=2.0/X
4307       BKM=BESSK0(X)
4308       BK=BESSK1(X)
4309       DO 11 J=1,N-1
4310         BKP=BKM+J*TOX*BK
4311         BKM=BK
4312         BK=BKP
4313 11    CONTINUE
4314       BESSK=BK
4315       RETURN
4316       END
4317 C
4318 C
4319       FUNCTION BESSI0(X)
4320       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4321 C
4322 *     REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,
4323 *    *    Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9
4324       SAVE
4325       DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0,
4326      *    1.2067492D0,
4327      *    0.2659732D0,0.360768D-1,0.45813D-2/
4328       DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1,
4329      *    0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,
4330      *    0.2635537D-1,-0.1647633D-1,0.392377D-2/
4331 
4332       IF (ABS(X).LT.3.75) THEN
4333         Y=(X/3.75)**2
4334         BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
4335       ELSE
4336         AX=ABS(X)
4337         Y=3.75/AX
4338         BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
4339      *      +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
4340       ENDIF
4341       RETURN
4342       END
4343 C
4344 C
4345       FUNCTION BESSI1(X)
4346       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4347 C
4348 *     REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,
4349 *    *    Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9
4350       SAVE
4351       DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,
4352      *    0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
4353       DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1,
4354      *    -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,
4355      *    -0.2895312D-1,0.1787654D-1,-0.420059D-2/
4356 
4357       IF (ABS(X).LT.3.75) THEN
4358         Y=(X/3.75)**2
4359         BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
4360       ELSE
4361         AX=ABS(X)
4362         Y=3.75/AX
4363         BESSI1=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+
4364      *      Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
4365       ENDIF
4366       RETURN
4367       END
4368 
4369 
4370       SUBROUTINE FACT_INI
4371       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
4372 
4373       PARAMETER (NS_max = 20, NH_max = 50)
4374       COMMON /S_CFACT/ FACT(0:NH_max), CO_BIN(0:NH_max,0:NH_max)
4375       SAVE
4376 
4377       FACT(0) = 1.
4378       DO J=1,NH_max
4379          FACT(J) = FACT(J-1)*FLOAT(J)
4380       ENDDO
4381       DO J=0,NH_max
4382          DO K=0,J
4383             CO_BIN(J,K) = FACT(J)/(FACT(K)*FACT(J-K))
4384          ENDDO
4385       ENDDO
4386 
4387       RETURN
4388       END
4389 
4390 
4391       SUBROUTINE SIB_GAUSET(AX,BX,NX,Z,W)
4392 C-----------------------------------------------------------------------
4393 C
4394 C     N-point gauss zeros and weights for the interval (AX,BX) are
4395 C           stored in  arrays Z and W respectively.
4396 C
4397 C-----------------------------------------------------------------------
4398       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4399 C
4400       COMMON /GQCOM/A(273),X(273),KTAB(96)
4401       DIMENSION Z(NX),W(NX)
4402       SAVE
4403       DATA INIT/0/
4404 C
4405       ALPHA=0.5*(BX+AX)
4406       BETA=0.5*(BX-AX)
4407       N=NX
4408 *
4409 *  the N=1 case:
4410       IF(N.NE.1) GO TO 1
4411       Z(1)=ALPHA
4412       W(1)=BX-AX
4413       RETURN
4414 *
4415 *  the Gauss cases:
4416     1 IF((N.LE.16).AND.(N.GT.1)) GO TO 2
4417       IF(N.EQ.20) GO TO 2
4418       IF(N.EQ.24) GO TO 2
4419       IF(N.EQ.32) GO TO 2
4420       IF(N.EQ.40) GO TO 2
4421       IF(N.EQ.48) GO TO 2
4422       IF(N.EQ.64) GO TO 2
4423       IF(N.EQ.80) GO TO 2
4424       IF(N.EQ.96) GO TO 2
4425 *
4426 *  the extended Gauss cases:
4427       IF((N/96)*96.EQ.N) GO TO 3
4428 *
4429 C  jump to center of intervall intrgration:
4430       GO TO 100
4431 *
4432 C  get Gauss point array
4433 *
4434     2 CALL PO106BD
4435 C     -print out message
4436 *     IF(INIT.LE.20)THEN
4437 *       INIT=init+1
4438 *       WRITE (6,*) ' initialization of Gauss int. N=',N
4439 *     ENDIF
4440 C  extract real points
4441       K=KTAB(N)
4442       M=N/2
4443       DO 21 J=1,M
4444 C       extract values from big array
4445         JTAB=K-1+J
4446         WTEMP=BETA*A(JTAB)
4447         DELTA=BETA*X(JTAB)
4448 C       store them backward
4449         Z(J)=ALPHA-DELTA
4450         W(J)=WTEMP
4451 C       store them forward
4452         JP=N+1-J
4453         Z(JP)=ALPHA+DELTA
4454         W(JP)=WTEMP
4455    21 CONTINUE
4456 C     store central point (odd N)
4457       IF((N-M-M).EQ.0) RETURN
4458       Z(M+1)=ALPHA
4459       JMID=K+M
4460       W(M+1)=BETA*A(JMID)
4461       RETURN
4462 C
4463 C  get ND96 times chained 96 Gauss point array
4464 C
4465     3 CALL PO106BD
4466 C  print out message
4467       IF(INIT.LE.20)THEN
4468         INIT=init+1
4469         WRITE (6,*) ' initialization of extended Gauss int. N=',N
4470       ENDIF
4471 C     -extract real points
4472       K=KTAB(96)
4473       ND96=N/96
4474       DO 31 J=1,48
4475 C       extract values from big array
4476         JTAB=K-1+J
4477         WTEMP=BETA*A(JTAB)
4478         DELTA=BETA*X(JTAB)
4479         WTeMP=WTEMP/ND96
4480         DeLTA=DELTA/ND96
4481         DO 32 JD96=0,ND96-1
4482           ZCNTR= (ALPHA-BETA)+ BETA*FLOAT(2*JD96+1)/FLOAT(ND96)
4483 C         store them backward
4484           Z(J+JD96*96)=ZCNTR-DELTA
4485           W(J+JD96*96)=WTEMP
4486 C         store them forward
4487           JP=96+1-J
4488           Z(JP+JD96*96)=ZCNTR+DELTA
4489           W(JP+JD96*96)=WTEMP
4490    32   CONTINUE
4491    31 CONTINUE
4492       RETURN
4493 *
4494 C  the center of intervall cases:
4495   100 CONTINUE
4496 C  print out message
4497       IF(INIT.LE.20)THEN
4498         INIT=init+1
4499         WRITE (6,*) ' init. of center of intervall int. N=',N
4500       ENDIF
4501 C  put in constant weight and equally spaced central points
4502       N=IABS(N)
4503       DO 111 IN=1,N
4504         WIN=(BX-AX)/FLOAT(N)
4505         Z(IN)=AX  + (FLOAT(IN)-.5)*WIN
4506   111 W(IN)=WIN
4507       RETURN
4508       END
4509 C
4510 C
4511       SUBROUTINE PO106BD
4512 C-----------------------------------------------------------------------
4513 C
4514 C     store big arrays needed for Gauss integral, CERNLIB D106BD
4515 C     (arrays A,X,ITAB copied on B,Y,LTAB)
4516 C
4517 C-----------------------------------------------------------------------
4518       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4519 C
4520       COMMON /GQCOM/ B(273),Y(273),LTAB(96)
4521       DIMENSION      A(273),X(273),KTAB(96)
4522       SAVE
4523 C
4524 C-----TABLE OF INITIAL SUBSCRIPTS FOR N=2(1)16(4)96
4525       DATA KTAB(2)/1/
4526       DATA KTAB(3)/2/
4527       DATA KTAB(4)/4/
4528       DATA KTAB(5)/6/
4529       DATA KTAB(6)/9/
4530       DATA KTAB(7)/12/
4531       DATA KTAB(8)/16/
4532       DATA KTAB(9)/20/
4533       DATA KTAB(10)/25/
4534       DATA KTAB(11)/30/
4535       DATA KTAB(12)/36/
4536       DATA KTAB(13)/42/
4537       DATA KTAB(14)/49/
4538       DATA KTAB(15)/56/
4539       DATA KTAB(16)/64/
4540       DATA KTAB(20)/72/
4541       DATA KTAB(24)/82/
4542       DATA KTAB(28)/82/
4543       DATA KTAB(32)/94/
4544       DATA KTAB(36)/94/
4545       DATA KTAB(40)/110/
4546       DATA KTAB(44)/110/
4547       DATA KTAB(48)/130/
4548       DATA KTAB(52)/130/
4549       DATA KTAB(56)/130/
4550       DATA KTAB(60)/130/
4551       DATA KTAB(64)/154/
4552       DATA KTAB(68)/154/
4553       DATA KTAB(72)/154/
4554       DATA KTAB(76)/154/
4555       DATA KTAB(80)/186/
4556       DATA KTAB(84)/186/
4557       DATA KTAB(88)/186/
4558       DATA KTAB(92)/186/
4559       DATA KTAB(96)/226/
4560 C
4561 C-----TABLE OF ABSCISSAE (X) AND WEIGHTS (A) FOR INTERVAL (-1,+1).
4562 C
4563 C-----N=2
4564       DATA X(1)/0.577350269189626D0  /, A(1)/1.000000000000000D0  /
4565 C-----N=3
4566       DATA X(2)/0.774596669241483D0  /, A(2)/0.555555555555556D0  /
4567       DATA X(3)/0.000000000000000D0  /, A(3)/0.888888888888889D0  /
4568 C-----N=4
4569       DATA X(4)/0.861136311594053D0  /, A(4)/0.347854845137454D0  /
4570       DATA X(5)/0.339981043584856D0  /, A(5)/0.652145154862546D0  /
4571 C-----N=5
4572       DATA X(6)/0.906179845938664D0  /, A(6)/0.236926885056189D0  /
4573       DATA X(7)/0.538469310105683D0  /, A(7)/0.478628670499366D0  /
4574       DATA X(8)/0.000000000000000D0  /, A(8)/0.568888888888889D0  /
4575 C-----N=6
4576       DATA X(9)/0.932469514203152D0  /, A(9)/0.171324492379170D0  /
4577       DATA X(10)/0.661209386466265D0 /, A(10)/0.360761573048139D0 /
4578       DATA X(11)/0.238619186083197D0 /, A(11)/0.467913934572691D0 /
4579 C-----N=7
4580       DATA X(12)/0.949107912342759D0 /, A(12)/0.129484966168870D0 /
4581       DATA X(13)/0.741531185599394D0 /, A(13)/0.279705391489277D0 /
4582       DATA X(14)/0.405845151377397D0 /, A(14)/0.381830050505119D0 /
4583       DATA X(15)/0.000000000000000D0 /, A(15)/0.417959183673469D0 /
4584 C-----N=8
4585       DATA X(16)/0.960289856497536D0 /, A(16)/0.101228536290376D0 /
4586       DATA X(17)/0.796666477413627D0 /, A(17)/0.222381034453374D0 /
4587       DATA X(18)/0.525532409916329D0 /, A(18)/0.313706645877887D0 /
4588       DATA X(19)/0.183434642495650D0 /, A(19)/0.362683783378362D0 /
4589 C-----N=9
4590       DATA X(20)/0.968160239507626D0 /, A(20)/0.081274388361574D0 /
4591       DATA X(21)/0.836031107326636D0 /, A(21)/0.180648160694857D0 /
4592       DATA X(22)/0.613371432700590D0 /, A(22)/0.260610696402935D0 /
4593       DATA X(23)/0.324253423403809D0 /, A(23)/0.312347077040003D0 /
4594       DATA X(24)/0.000000000000000D0 /, A(24)/0.330239355001260D0 /
4595 C-----N=10
4596       DATA X(25)/0.973906528517172D0 /, A(25)/0.066671344308688D0 /
4597       DATA X(26)/0.865063366688985D0 /, A(26)/0.149451349150581D0 /
4598       DATA X(27)/0.679409568299024D0 /, A(27)/0.219086362515982D0 /
4599       DATA X(28)/0.433395394129247D0 /, A(28)/0.269266719309996D0 /
4600       DATA X(29)/0.148874338981631D0 /, A(29)/0.295524224714753D0 /
4601 C-----N=11
4602       DATA X(30)/0.978228658146057D0 /, A(30)/0.055668567116174D0 /
4603       DATA X(31)/0.887062599768095D0 /, A(31)/0.125580369464905D0 /
4604       DATA X(32)/0.730152005574049D0 /, A(32)/0.186290210927734D0 /
4605       DATA X(33)/0.519096129206812D0 /, A(33)/0.233193764591990D0 /
4606       DATA X(34)/0.269543155952345D0 /, A(34)/0.262804544510247D0 /
4607       DATA X(35)/0.000000000000000D0 /, A(35)/0.272925086777901D0 /
4608 C-----N=12
4609       DATA X(36)/0.981560634246719D0 /, A(36)/0.047175336386512D0 /
4610       DATA X(37)/0.904117256370475D0 /, A(37)/0.106939325995318D0 /
4611       DATA X(38)/0.769902674194305D0 /, A(38)/0.160078328543346D0 /
4612       DATA X(39)/0.587317954286617D0 /, A(39)/0.203167426723066D0 /
4613       DATA X(40)/0.367831498998180D0 /, A(40)/0.233492536538355D0 /
4614       DATA X(41)/0.125233408511469D0 /, A(41)/0.249147045813403D0 /
4615 C-----N=13
4616       DATA X(42)/0.984183054718588D0 /, A(42)/0.040484004765316D0 /
4617       DATA X(43)/0.917598399222978D0 /, A(43)/0.092121499837728D0 /
4618       DATA X(44)/0.801578090733310D0 /, A(44)/0.138873510219787D0 /
4619       DATA X(45)/0.642349339440340D0 /, A(45)/0.178145980761946D0 /
4620       DATA X(46)/0.448492751036447D0 /, A(46)/0.207816047536889D0 /
4621       DATA X(47)/0.230458315955135D0 /, A(47)/0.226283180262897D0 /
4622       DATA X(48)/0.000000000000000D0 /, A(48)/0.232551553230874D0 /
4623 C-----N=14
4624       DATA X(49)/0.986283808696812D0 /, A(49)/0.035119460331752D0 /
4625       DATA X(50)/0.928434883663574D0 /, A(50)/0.080158087159760D0 /
4626       DATA X(51)/0.827201315069765D0 /, A(51)/0.121518570687903D0 /
4627       DATA X(52)/0.687292904811685D0 /, A(52)/0.157203167158194D0 /
4628       DATA X(53)/0.515248636358154D0 /, A(53)/0.185538397477938D0 /
4629       DATA X(54)/0.319112368927890D0 /, A(54)/0.205198463721296D0 /
4630       DATA X(55)/0.108054948707344D0 /, A(55)/0.215263853463158D0 /
4631 C-----N=15
4632       DATA X(56)/0.987992518020485D0 /, A(56)/0.030753241996117D0 /
4633       DATA X(57)/0.937273392400706D0 /, A(57)/0.070366047488108D0 /
4634       DATA X(58)/0.848206583410427D0 /, A(58)/0.107159220467172D0 /
4635       DATA X(59)/0.724417731360170D0 /, A(59)/0.139570677926154D0 /
4636       DATA X(60)/0.570972172608539D0 /, A(60)/0.166269205816994D0 /
4637       DATA X(61)/0.394151347077563D0 /, A(61)/0.186161000015562D0 /
4638       DATA X(62)/0.201194093997435D0 /, A(62)/0.198431485327111D0 /
4639       DATA X(63)/0.000000000000000D0 /, A(63)/0.202578241925561D0 /
4640 C-----N=16
4641       DATA X(64)/0.989400934991650D0 /, A(64)/0.027152459411754D0 /
4642       DATA X(65)/0.944575023073233D0 /, A(65)/0.062253523938648D0 /
4643       DATA X(66)/0.865631202387832D0 /, A(66)/0.095158511682493D0 /
4644       DATA X(67)/0.755404408355003D0 /, A(67)/0.124628971255534D0 /
4645       DATA X(68)/0.617876244402644D0 /, A(68)/0.149595988816577D0 /
4646       DATA X(69)/0.458016777657227D0 /, A(69)/0.169156519395003D0 /
4647       DATA X(70)/0.281603550779259D0 /, A(70)/0.182603415044924D0 /
4648       DATA X(71)/0.095012509837637D0 /, A(71)/0.189450610455069D0 /
4649 C-----N=20
4650       DATA X(72)/0.993128599185094D0 /, A(72)/0.017614007139152D0 /
4651       DATA X(73)/0.963971927277913D0 /, A(73)/0.040601429800386D0 /
4652       DATA X(74)/0.912234428251325D0 /, A(74)/0.062672048334109D0 /
4653       DATA X(75)/0.839116971822218D0 /, A(75)/0.083276741576704D0 /
4654       DATA X(76)/0.746331906460150D0 /, A(76)/0.101930119817240D0 /
4655       DATA X(77)/0.636053680726515D0 /, A(77)/0.118194531961518D0 /
4656       DATA X(78)/0.510867001950827D0 /, A(78)/0.131688638449176D0 /
4657       DATA X(79)/0.373706088715419D0 /, A(79)/0.142096109318382D0 /
4658       DATA X(80)/0.227785851141645D0 /, A(80)/0.149172986472603D0 /
4659       DATA X(81)/0.076526521133497D0 /, A(81)/0.152753387130725D0 /
4660 C-----N=24
4661       DATA X(82)/0.995187219997021D0 /, A(82)/0.012341229799987D0 /
4662       DATA X(83)/0.974728555971309D0 /, A(83)/0.028531388628933D0 /
4663       DATA X(84)/0.938274552002732D0 /, A(84)/0.044277438817419D0 /
4664       DATA X(85)/0.886415527004401D0 /, A(85)/0.059298584915436D0 /
4665       DATA X(86)/0.820001985973902D0 /, A(86)/0.073346481411080D0 /
4666       DATA X(87)/0.740124191578554D0 /, A(87)/0.086190161531953D0 /
4667       DATA X(88)/0.648093651936975D0 /, A(88)/0.097618652104113D0 /
4668       DATA X(89)/0.545421471388839D0 /, A(89)/0.107444270115965D0 /
4669       DATA X(90)/0.433793507626045D0 /, A(90)/0.115505668053725D0 /
4670       DATA X(91)/0.315042679696163D0 /, A(91)/0.121670472927803D0 /
4671       DATA X(92)/0.191118867473616D0 /, A(92)/0.125837456346828D0 /
4672       DATA X(93)/0.064056892862605D0 /, A(93)/0.127938195346752D0 /
4673 C-----N=32
4674       DATA X(94)/0.997263861849481D0 /, A(94)/0.007018610009470D0 /
4675       DATA X(95)/0.985611511545268D0 /, A(95)/0.016274394730905D0 /
4676       DATA X(96)/0.964762255587506D0 /, A(96)/0.025392065309262D0 /
4677       DATA X(97)/0.934906075937739D0 /, A(97)/0.034273862913021D0 /
4678       DATA X(98)/0.896321155766052D0 /, A(98)/0.042835898022226D0 /
4679       DATA X(99)/0.849367613732569D0 /, A(99)/0.050998059262376D0 /
4680       DATA X(100)/0.794483795967942D0/, A(100)/0.058684093478535D0/
4681       DATA X(101)/0.732182118740289D0/, A(101)/0.065822222776361D0/
4682       DATA X(102)/0.663044266930215D0/, A(102)/0.072345794108848D0/
4683       DATA X(103)/0.587715757240762D0/, A(103)/0.078193895787070D0/
4684       DATA X(104)/0.506899908932229D0/, A(104)/0.083311924226946D0/
4685       DATA X(105)/0.421351276130635D0/, A(105)/0.087652093004403D0/
4686       DATA X(106)/0.331868602282127D0/, A(106)/0.091173878695763D0/
4687       DATA X(107)/0.239287362252137D0/, A(107)/0.093844399080804D0/
4688       DATA X(108)/0.144471961582796D0/, A(108)/0.095638720079274D0/
4689       DATA X(109)/0.048307665687738D0/, A(109)/0.096540088514727D0/
4690 C-----N=40
4691       DATA X(110)/0.998237709710559D0/, A(110)/0.004521277098533D0/
4692       DATA X(111)/0.990726238699457D0/, A(111)/0.010498284531152D0/
4693       DATA X(112)/0.977259949983774D0/, A(112)/0.016421058381907D0/
4694       DATA X(113)/0.957916819213791D0/, A(113)/0.022245849194166D0/
4695       DATA X(114)/0.932812808278676D0/, A(114)/0.027937006980023D0/
4696       DATA X(115)/0.902098806968874D0/, A(115)/0.033460195282547D0/
4697       DATA X(116)/0.865959503212259D0/, A(116)/0.038782167974472D0/
4698       DATA X(117)/0.824612230833311D0/, A(117)/0.043870908185673D0/
4699       DATA X(118)/0.778305651426519D0/, A(118)/0.048695807635072D0/
4700       DATA X(119)/0.727318255189927D0/, A(119)/0.053227846983936D0/
4701       DATA X(120)/0.671956684614179D0/, A(120)/0.057439769099391D0/
4702       DATA X(121)/0.612553889667980D0/, A(121)/0.061306242492928D0/
4703       DATA X(122)/0.549467125095128D0/, A(122)/0.064804013456601D0/
4704       DATA X(123)/0.483075801686178D0/, A(123)/0.067912045815233D0/
4705       DATA X(124)/0.413779204371605D0/, A(124)/0.070611647391286D0/
4706       DATA X(125)/0.341994090825758D0/, A(125)/0.072886582395804D0/
4707       DATA X(126)/0.268152185007253D0/, A(126)/0.074723169057968D0/
4708       DATA X(127)/0.192697580701371D0/, A(127)/0.076110361900626D0/
4709       DATA X(128)/0.116084070675255D0/, A(128)/0.077039818164247D0/
4710       DATA X(129)/0.038772417506050D0/, A(129)/0.077505947978424D0/
4711 C-----N=48
4712       DATA X(130)/0.998771007252426D0/, A(130)/0.003153346052305D0/
4713       DATA X(131)/0.993530172266350D0/, A(131)/0.007327553901276D0/
4714       DATA X(132)/0.984124583722826D0/, A(132)/0.011477234579234D0/
4715       DATA X(133)/0.970591592546247D0/, A(133)/0.015579315722943D0/
4716       DATA X(134)/0.952987703160430D0/, A(134)/0.019616160457355D0/
4717       DATA X(135)/0.931386690706554D0/, A(135)/0.023570760839324D0/
4718       DATA X(136)/0.905879136715569D0/, A(136)/0.027426509708356D0/
4719       DATA X(137)/0.876572020274247D0/, A(137)/0.031167227832798D0/
4720       DATA X(138)/0.843588261624393D0/, A(138)/0.034777222564770D0/
4721       DATA X(139)/0.807066204029442D0/, A(139)/0.038241351065830D0/
4722       DATA X(140)/0.767159032515740D0/, A(140)/0.041545082943464D0/
4723       DATA X(141)/0.724034130923814D0/, A(141)/0.044674560856694D0/
4724       DATA X(142)/0.677872379632663D0/, A(142)/0.047616658492490D0/
4725       DATA X(143)/0.628867396776513D0/, A(143)/0.050359035553854D0/
4726       DATA X(144)/0.577224726083972D0/, A(144)/0.052890189485193D0/
4727       DATA X(145)/0.523160974722233D0/, A(145)/0.055199503699984D0/
4728       DATA X(146)/0.466902904750958D0/, A(146)/0.057277292100403D0/
4729       DATA X(147)/0.408686481990716D0/, A(147)/0.059114839698395D0/
4730       DATA X(148)/0.348755886292160D0/, A(148)/0.060704439165893D0/
4731       DATA X(149)/0.287362487355455D0/, A(149)/0.062039423159892D0/
4732       DATA X(150)/0.224763790394689D0/, A(150)/0.063114192286254D0/
4733       DATA X(151)/0.161222356068891D0/, A(151)/0.063924238584648D0/
4734       DATA X(152)/0.097004699209462D0/, A(152)/0.064466164435950D0/
4735       DATA X(153)/0.032380170962869D0/, A(153)/0.064737696812683D0/
4736 C-----N=64
4737       DATA X(154)/0.999305041735772D0/, A(154)/0.001783280721696D0/
4738       DATA X(155)/0.996340116771955D0/, A(155)/0.004147033260562D0/
4739       DATA X(156)/0.991013371476744D0/, A(156)/0.006504457968978D0/
4740       DATA X(157)/0.983336253884625D0/, A(157)/0.008846759826363D0/
4741       DATA X(158)/0.973326827789910D0/, A(158)/0.011168139460131D0/
4742       DATA X(159)/0.961008799652053D0/, A(159)/0.013463047896718D0/
4743       DATA X(160)/0.946411374858402D0/, A(160)/0.015726030476024D0/
4744       DATA X(161)/0.929569172131939D0/, A(161)/0.017951715775697D0/
4745       DATA X(162)/0.910522137078502D0/, A(162)/0.020134823153530D0/
4746       DATA X(163)/0.889315445995114D0/, A(163)/0.022270173808383D0/
4747       DATA X(164)/0.865999398154092D0/, A(164)/0.024352702568710D0/
4748       DATA X(165)/0.840629296252580D0/, A(165)/0.026377469715054D0/
4749       DATA X(166)/0.813265315122797D0/, A(166)/0.028339672614259D0/
4750       DATA X(167)/0.783972358943341D0/, A(167)/0.030234657072402D0/
4751       DATA X(168)/0.752819907260531D0/, A(168)/0.032057928354851D0/
4752       DATA X(169)/0.719881850171610D0/, A(169)/0.033805161837141D0/
4753       DATA X(170)/0.685236313054233D0/, A(170)/0.035472213256882D0/
4754       DATA X(171)/0.648965471254657D0/, A(171)/0.037055128540240D0/
4755       DATA X(172)/0.611155355172393D0/, A(172)/0.038550153178615D0/
4756       DATA X(173)/0.571895646202634D0/, A(173)/0.039953741132720D0/
4757       DATA X(174)/0.531279464019894D0/, A(174)/0.041262563242623D0/
4758       DATA X(175)/0.489403145707052D0/, A(175)/0.042473515123653D0/
4759       DATA X(176)/0.446366017253464D0/, A(176)/0.043583724529323D0/
4760       DATA X(177)/0.402270157963991D0/, A(177)/0.044590558163756D0/
4761       DATA X(178)/0.357220158337668D0/, A(178)/0.045491627927418D0/
4762       DATA X(179)/0.311322871990210D0/, A(179)/0.046284796581314D0/
4763       DATA X(180)/0.264687162208767D0/, A(180)/0.046968182816210D0/
4764       DATA X(181)/0.217423643740007D0/, A(181)/0.047540165714830D0/
4765       DATA X(182)/0.169644420423992D0/, A(182)/0.047999388596458D0/
4766       DATA X(183)/0.121462819296120D0/, A(183)/0.048344762234802D0/
4767       DATA X(184)/0.072993121787799D0/, A(184)/0.048575467441503D0/
4768       DATA X(185)/0.024350292663424D0/, A(185)/0.048690957009139D0/
4769 C-----N=80
4770       DATA X(186)/0.999553822651630D0/, A(186)/0.001144950003186D0/
4771       DATA X(187)/0.997649864398237D0/, A(187)/0.002663533589512D0/
4772       DATA X(188)/0.994227540965688D0/, A(188)/0.004180313124694D0/
4773       DATA X(189)/0.989291302499755D0/, A(189)/0.005690922451403D0/
4774       DATA X(190)/0.982848572738629D0/, A(190)/0.007192904768117D0/
4775       DATA X(191)/0.974909140585727D0/, A(191)/0.008683945269260D0/
4776       DATA X(192)/0.965485089043799D0/, A(192)/0.010161766041103D0/
4777       DATA X(193)/0.954590766343634D0/, A(193)/0.011624114120797D0/
4778       DATA X(194)/0.942242761309872D0/, A(194)/0.013068761592401D0/
4779       DATA X(195)/0.928459877172445D0/, A(195)/0.014493508040509D0/
4780       DATA X(196)/0.913263102571757D0/, A(196)/0.015896183583725D0/
4781       DATA X(197)/0.896675579438770D0/, A(197)/0.017274652056269D0/
4782       DATA X(198)/0.878722567678213D0/, A(198)/0.018626814208299D0/
4783       DATA X(199)/0.859431406663111D0/, A(199)/0.019950610878141D0/
4784       DATA X(200)/0.838831473580255D0/, A(200)/0.021244026115782D0/
4785       DATA X(201)/0.816954138681463D0/, A(201)/0.022505090246332D0/
4786       DATA X(202)/0.793832717504605D0/, A(202)/0.023731882865930D0/
4787       DATA X(203)/0.769502420135041D0/, A(203)/0.024922535764115D0/
4788       DATA X(204)/0.744000297583597D0/, A(204)/0.026075235767565D0/
4789       DATA X(205)/0.717365185362099D0/, A(205)/0.027188227500486D0/
4790       DATA X(206)/0.689637644342027D0/, A(206)/0.028259816057276D0/
4791       DATA X(207)/0.660859898986119D0/, A(207)/0.029288369583267D0/
4792       DATA X(208)/0.631075773046871D0/, A(208)/0.030272321759557D0/
4793       DATA X(209)/0.600330622829751D0/, A(209)/0.031210174188114D0/
4794       DATA X(210)/0.568671268122709D0/, A(210)/0.032100498673487D0/
4795       DATA X(211)/0.536145920897131D0/, A(211)/0.032941939397645D0/
4796       DATA X(212)/0.502804111888784D0/, A(212)/0.033733214984611D0/
4797       DATA X(213)/0.468696615170544D0/, A(213)/0.034473120451753D0/
4798       DATA X(214)/0.433875370831756D0/, A(214)/0.035160529044747D0/
4799       DATA X(215)/0.398393405881969D0/, A(215)/0.035794393953416D0/
4800       DATA X(216)/0.362304753499487D0/, A(216)/0.036373749905835D0/
4801       DATA X(217)/0.325664370747701D0/, A(217)/0.036897714638276D0/
4802       DATA X(218)/0.288528054884511D0/, A(218)/0.037365490238730D0/
4803       DATA X(219)/0.250952358392272D0/, A(219)/0.037776364362001D0/
4804       DATA X(220)/0.212994502857666D0/, A(220)/0.038129711314477D0/
4805       DATA X(221)/0.174712291832646D0/, A(221)/0.038424993006959D0/
4806       DATA X(222)/0.136164022809143D0/, A(222)/0.038661759774076D0/
4807       DATA X(223)/0.097408398441584D0/, A(223)/0.038839651059051D0/
4808       DATA X(224)/0.058504437152420D0/, A(224)/0.038958395962769D0/
4809       DATA X(225)/0.019511383256793D0/, A(225)/0.039017813656306D0/
4810 C-----N=96
4811       DATA X(226)/0.999689503883230D0/, A(226)/0.000796792065552D0/
4812       DATA X(227)/0.998364375863181D0/, A(227)/0.001853960788946D0/
4813       DATA X(228)/0.995981842987209D0/, A(228)/0.002910731817934D0/
4814       DATA X(229)/0.992543900323762D0/, A(229)/0.003964554338444D0/
4815       DATA X(230)/0.988054126329623D0/, A(230)/0.005014202742927D0/
4816       DATA X(231)/0.982517263563014D0/, A(231)/0.006058545504235D0/
4817       DATA X(232)/0.975939174585136D0/, A(232)/0.007096470791153D0/
4818       DATA X(233)/0.968326828463264D0/, A(233)/0.008126876925698D0/
4819       DATA X(234)/0.959688291448742D0/, A(234)/0.009148671230783D0/
4820       DATA X(235)/0.950032717784437D0/, A(235)/0.010160770535008D0/
4821       DATA X(236)/0.939370339752755D0/, A(236)/0.011162102099838D0/
4822       DATA X(237)/0.927712456722308D0/, A(237)/0.012151604671088D0/
4823       DATA X(238)/0.915071423120898D0/, A(238)/0.013128229566961D0/
4824       DATA X(239)/0.901460635315852D0/, A(239)/0.014090941772314D0/
4825       DATA X(240)/0.886894517402420D0/, A(240)/0.015038721026994D0/
4826       DATA X(241)/0.871388505909296D0/, A(241)/0.015970562902562D0/
4827       DATA X(242)/0.854959033434601D0/, A(242)/0.016885479864245D0/
4828       DATA X(243)/0.837623511228187D0/, A(243)/0.017782502316045D0/
4829       DATA X(244)/0.819400310737931D0/, A(244)/0.018660679627411D0/
4830       DATA X(245)/0.800308744139140D0/, A(245)/0.019519081140145D0/
4831       DATA X(246)/0.780369043867433D0/, A(246)/0.020356797154333D0/
4832       DATA X(247)/0.759602341176647D0/, A(247)/0.021172939892191D0/
4833       DATA X(248)/0.738030643744400D0/, A(248)/0.021966644438744D0/
4834       DATA X(249)/0.715676812348967D0/, A(249)/0.022737069658329D0/
4835       DATA X(250)/0.692564536642171D0/, A(250)/0.023483399085926D0/
4836       DATA X(251)/0.668718310043916D0/, A(251)/0.024204841792364D0/
4837       DATA X(252)/0.644163403784967D0/, A(252)/0.024900633222483D0/
4838       DATA X(253)/0.618925840125468D0/, A(253)/0.025570036005349D0/
4839       DATA X(254)/0.593032364777572D0/, A(254)/0.026212340735672D0/
4840       DATA X(255)/0.566510418561397D0/, A(255)/0.026826866725591D0/
4841       DATA X(256)/0.539388108324357D0/, A(256)/0.027412962726029D0/
4842       DATA X(257)/0.511694177154667D0/, A(257)/0.027970007616848D0/
4843       DATA X(258)/0.483457973920596D0/, A(258)/0.028497411065085D0/
4844       DATA X(259)/0.454709422167743D0/, A(259)/0.028994614150555D0/
4845       DATA X(260)/0.425478988407300D0/, A(260)/0.029461089958167D0/
4846       DATA X(261)/0.395797649828908D0/, A(261)/0.029896344136328D0/
4847       DATA X(262)/0.365696861472313D0/, A(262)/0.030299915420827D0/
4848       DATA X(263)/0.335208522892625D0/, A(263)/0.030671376123669D0/
4849       DATA X(264)/0.304364944354496D0/, A(264)/0.031010332586313D0/
4850       DATA X(265)/0.273198812591049D0/, A(265)/0.031316425596861D0/
4851       DATA X(266)/0.241743156163840D0/, A(266)/0.031589330770727D0/
4852       DATA X(267)/0.210031310460567D0/, A(267)/0.031828758894411D0/
4853       DATA X(268)/0.178096882367618D0/, A(268)/0.032034456231992D0/
4854       DATA X(269)/0.145973714654896D0/, A(269)/0.032206204794030D0/
4855       DATA X(270)/0.113695850110665D0/, A(270)/0.032343822568575D0/
4856       DATA X(271)/0.081297495464425D0/, A(271)/0.032447163714064D0/
4857       DATA X(272)/0.048812985136049D0/, A(272)/0.032516118713868D0/
4858       DATA X(273)/0.016276744849602D0/, A(273)/0.032550614492363D0/
4859       DATA IBD/0/
4860       IF(IBD.NE.0) RETURN
4861       IBD=1
4862       DO 10 I=1,273
4863         B(I) = A(I)
4864 10      Y(I) = X(I)
4865       DO 20 I=1,96
4866 20      LTAB(I) = KTAB(I)
4867       RETURN
4868       END