PartonShowerVeto/src/ME2pythia.f

0001 C*******************************************************************
0002 C*          MadEvent - Pythia interface.                           *
0003 C*                                                                 *
0004 C*                                                                 *
0005 C*   Adapted version from ME2pythia 1.66 (J.Alwall)                *
0006 C*                                                                 *
0007 C*       Simon de Visscher August 2010   sdevissc@cern.ch          *
0008 C*        - Complete upgrade to 1.66                               *
0009 C*        - Addition of showerkt                                   *
0010 C*        - Addition of resonance exclusion for BSM matching       *
0011 C*                                                                 *
0012 C*       Christophe Saout 2009                                     *
0013 C*        - Improvement of JetMatching interface                   *
0014 C*                                                                 *
0015 C*       Dorian Kcira 2008.02.05                                         *
0016 C*        -First implementation of KtMLM in CMSSW                  *
0017 C*       - Improvement of matching routines                        *
0018 C*                                                                 *
0019 C*                                                                 *
0020 C*******************************************************************
0021
0022
0023 C*********************************************************************
0024 C...UPINIT
0025 C...Routine called by PYINIT to set up user-defined processes.
0026 C*********************************************************************
0027       SUBROUTINE MGINIT(npara,param,value)
0028
0029
0030       IMPLICIT NONE
0031
0032       integer npara
0033       character*20 param(*),value(*)
0034
0035
0036 c      CHARACTER*132 CHAR_READ
0037
0038 C...Pythia parameters.
0039       INTEGER MSTP,MSTI,MRPY
0040       DOUBLE PRECISION PARP,PARI,RRPY
0041       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0042       COMMON/PYDATR/MRPY(6),RRPY(100)
0043
0044 C...User process initialization commonblock.
0045       INTEGER MAXPUP
0046       PARAMETER (MAXPUP=100)
0047       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0048       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0049       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0050      &   IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0051      &   LPRUP(MAXPUP)
0052
0053 C...Extra commonblock to transfer run info.
0054       INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0055       COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0056       DATA LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE/77,6,1,0,0,1/
0057       SAVE /UPPRIV/
0058
0059 C...Inputs for the matching algorithm
0060       double precision etcjet,rclmax,etaclmax,qcut,clfact,showerkt
0061       integer maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres(30)
0062       integer nqmatch,nexcproc,iexcproc(MAXPUP),iexcval(MAXPUP)
0063       logical nosingrad,jetprocs
0064       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,showerkt,clfact,
0065      $   maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres,
0066      $   nqmatch,nexcproc,iexcproc,iexcval,nosingrad,jetprocs
0067
0068
0069 C...Parameter arrays (local)
0070 C      integer maxpara
0071 C      parameter (maxpara=1000)
0072 C      integer npara,iseed
0073 C      character*20 param(maxpara),value(maxpara)
0074
0075 C...Lines to read in assumed never longer than 200 characters.
0076 C      INTEGER MAXLEN,IBEG,IPR,I
0077 C      PARAMETER (MAXLEN=200)
0078 C      CHARACTER*(MAXLEN) STRING
0079
0080 C...Functions
0081 C      INTEGER iexclusive
0082 C      EXTERNAL iexclusive
0083
0084 C...Format for reading lines.
0085 C      CHARACTER*6 STRFMT
0086 C      STRFMT='(A000)'
0087 C      WRITE(STRFMT(3:5),'(I3)') MAXLEN
0088
0089 C...Extract the model parameter card and read it.
0090 C      CALL MODELPAR(LNHIN)
0091
0092 c...Read the <init> block information
0093
0094 C...Loop until finds line beginning with "<init>" or "<init ".
0095 C  100 READ(LNHIN,STRFMT,END=130,ERR=130) STRING
0096 C...Pick out random number seed and use for PYR initialization
0097 C      IF(INDEX(STRING,'iseed').NE.0)THEN
0098 C         READ(STRING,*) iseed
0099 C         IF(iseed.gt.0) THEN
0100 C            WRITE(LNHOUT,*) 'Initializing PYR with random seed ',iseed
0101 c            MRPY(1) = iseed
0102 C            MRPY(2) = 0
0103 C         ENDIF
0104 C      ENDIF
0105 C      IBEG=0
0106 C  110 IBEG=IBEG+1
0107 C...Allow indentation.
0108 C      IF(STRING(IBEG:IBEG).EQ.' '.AND.IBEG.LT.MAXLEN-5) GOTO 110
0109 C      IF(STRING(IBEG:IBEG+5).NE.'<init>'.AND.
0110 C     &STRING(IBEG:IBEG+5).NE.'<init ') GOTO 100
0111
0112 C...Read first line of initialization info.
0113 C      READ(LNHIN,*,END=130,ERR=130) IDBMUP(1),IDBMUP(2),EBMUP(1),
0114 C     &EBMUP(2),PDFGUP(1),PDFGUP(2),PDFSUP(1),PDFSUP(2),IDWTUP,NPRUP
0115
0116 C...Read NPRUP subsequent lines with information on each process.
0117 C      DO 120 IPR=1,NPRUP
0118 C        READ(LNHIN,*,END=130,ERR=130) XSECUP(IPR),XERRUP(IPR),
0119 C     &  XMAXUP(IPR),LPRUP(IPR)
0120 C  120 CONTINUE
0121
0122 C...Set PDFLIB or LHAPDF pdf number for Pythia
0123
0124 C      IF(PDFSUP(1).NE.19070.AND.(PDFSUP(1).NE.0.OR.PDFSUP(2).NE.0))THEN
0125 c     Not CTEQ5L, which is standard in Pythia
0126 C         CALL PYGIVE('MSTP(52)=2')
0127 c     The following works for both PDFLIB and LHAPDF (where PDFGUP(1)=0)
0128 c     But note that the MadEvent output uses the LHAPDF numbering scheme
0129 C        IF(PDFSUP(1).NE.0)THEN
0130 C           MSTP(51)=1000*PDFGUP(1)+PDFSUP(1)
0131 C        ELSE
0132 C           MSTP(51)=1000*PDFGUP(2)+PDFSUP(2)
0133 C        ENDIF
0134 C      ENDIF
0135
0136 C...Initialize widths and partial widths for resonances.
0137 C      CALL PYINRE
0138
0139 C...Calculate xsec reduction due to non-decayed resonances
0140 C...based on first event only!
0141
0142 C      CALL BRSUPP
0143
0144 C      REWIND(LNHIN)
0145
0146 C...Extract cuts and matching parameters
0147 C      CALL read_params(LNHIN,npara,param,value,maxpara)
0148
0149 C      call get_integer(npara,param,value," ickkw ",ickkw,0)
0150 C      if(ickkw.eq.1)then
0151 C         call get_integer(npara,param,value," ktscheme ",mektsc,1)
0152 C         write(*,*)'Running matching with ME ktscheme ',mektsc
0153 C      endif
0154 C
0155 C...Set kt clustering scheme (if not already set)
0156 C
0157       integer i
0158       call initpydata
0159       write(*,*)"MGINIT: ickkw is ",ickkw
0160       write(*,*)"MGINIT: ktscheme is ",mektsc
0161       write(*,*)"MGINIT: QCut is ",qcut
0162       write(*,*)"MGINIT: Showerkt is ",showerkt
0163           do 10 i = 1, nexcres
0164          write(*,*) 'EXCRES(', i,')=',EXCRES(i)
0165   10  continue
0166
0167       IF(ABS(IDBMUP(1)).EQ.11.AND.ABS(IDBMUP(2)).EQ.11.AND.
0168      $     IDBMUP(1).EQ.-IDBMUP(2).AND.ktsche.EQ.0)THEN
0169          ktsche=1
0170       ELSE IF(ktsche.EQ.0) THEN
0171          ktsche=4313
0172       ENDIF
0173
0174 C...Enhance primordial kt
0175 c      CALL PYGIVE('PARP(91)=2.5')
0176 c      CALL PYGIVE('PARP(93)=15')
0177
0178       IF(ickkw.gt.0) CALL set_matching(npara,param,value)
0179
0180 C...For photon initial states from protons: Set proton not to break up
0181       CALL PYGIVE('MSTP(98)=1')
0182
0183
0184 C      IF(ickkw.gt.0.and.(NPRUP.gt.1.or.iexclusive(LPRUP(1)).ne.-1))
0185 C     $     CALL set_matching(LNHIN,npara,param,value)
0186
0187 C...For photon initial states from protons: Set proton not to break up
0188 C      CALL PYGIVE('MSTP(98)=1')
0189
0190 C...Reset event numbering
0191 C      IEVNT=0
0192
0193       RETURN
0194
0195 C...Error exit: give up if initalization does not work.
0196 C  130 WRITE(*,*) ' Failed to read LHEF initialization information.'
0197 C      WRITE(*,*) ' Event generation will be stopped.'
0198 C      STOP
0199       END
0200
0201 C*********************************************************************
0202 C...UPEVNT
0203 C...Routine called by PYEVNT or PYEVNW to get user process event
0204 C*********************************************************************
0205       SUBROUTINE MGEVNT
0206
0207       IMPLICIT NONE
0208
0209 C...Pythia parameters.
0210       INTEGER MSTP,MSTI
0211       DOUBLE PRECISION PARP,PARI
0212       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0213
0214 C...User process initialization commonblock.
0215       INTEGER MAXPUP
0216       PARAMETER (MAXPUP=100)
0217       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0218       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0219       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0220      &   IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0221      &   LPRUP(MAXPUP)
0222 C...User process event common block.
0223       INTEGER MAXNUP
0224       PARAMETER (MAXNUP=500)
0225       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0226       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0227       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0228      &   ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0229      &   VTIMUP(MAXNUP),SPINUP(MAXNUP)
0230 C...Pythia common blocks
0231       INTEGER PYCOMP,KCHG,MINT,NPART,NPARTD,IPART,MAXNUR
0232       DOUBLE PRECISION PMAS,PARF,VCKM,VINT,PTPART
0233 C...Particle properties + some flavour parameters.
0234       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0235       COMMON/PYINT1/MINT(400),VINT(400)
0236       PARAMETER (MAXNUR=1000)
0237       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0238
0239 C...Extra commonblock to transfer run info.
0240       INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0241       COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0242
0243 C...Inputs for the matching algorithm
0244       double precision etcjet,rclmax,etaclmax,qcut,clfact,showerkt
0245       integer maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres(30)
0246       integer nqmatch,nexcproc,iexcproc(MAXPUP),iexcval(MAXPUP)
0247       logical nosingrad,jetprocs
0248       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,showerkt,clfact,
0249      $   maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres,
0250      $   nqmatch,nexcproc,iexcproc,iexcval,nosingrad,jetprocs
0251
0252 C...Commonblock to transfer event-by-event matching info
0253       INTEGER NLJETS,IEXC,Ifile
0254       DOUBLE PRECISION PTCLUS
0255       COMMON/MEMAEV/PTCLUS(20),NLJETS,IEXC,Ifile
0256
0257 C...Local variables
0258       INTEGER I,NEX,KP(MAXNUP),MOTH,NUPREAD,II,iexcl
0259       DOUBLE PRECISION PSUM,ESUM
0260 C...Lines to read in assumed never longer than 200 characters.
0261       INTEGER MAXLEN
0262       PARAMETER (MAXLEN=200)
0263
0264 C...Functions
0265       INTEGER iexclusive
0266       EXTERNAL iexclusive
0267 C...Format for reading lines.
0268 C      CHARACTER*6 STRFMT
0269 C      CHARACTER*1 CDUM
0270
0271 C      STRFMT='(A000)'
0272 C      WRITE(STRFMT(3:5),'(I3)') MAXLEN
0273
0274 C...Loop until finds line beginning with "<event>" or "<event ".
0275 C  100 READ(LNHIN,STRFMT,END=900,ERR=900) STRING
0276 C      IBEG=0
0277 C  110 IBEG=IBEG+1
0278 C...Allow indentation.
0279 C      IF(STRING(IBEG:IBEG).EQ.' '.AND.IBEG.LT.MAXLEN-6) GOTO 110
0280 C      IF(STRING(IBEG:IBEG+6).NE.'<event>'.AND.
0281 C     &STRING(IBEG:IBEG+6).NE.'<event ') GOTO 100
0282
0283 C...Read first line of event info.
0284 C      READ(LNHIN,*,END=900,ERR=900) NUPREAD,IDPRUP,XWGTUP,SCALUP,
0285 C     &AQEDUP,AQCDUP
0286
0287       NUPREAD=NUP
0288
0289 C...Read NUP subsequent lines with information on each particle.
0290       ESUM=0d0
0291       PSUM=0d0
0292       NEX=2
0293       NUP=1
0294 C      write(*,*)'SCALUP=',SCALUP
0295       DO 120 I=1,NUPREAD
0296 C      write(*,*)IDUP(NUP),' ',ISTUP(NUP),' ',MOTHUP(1,NUP)
0297 C     &  ,' ',MOTHUP(2,NUP),' ',ICOLUP(1,NUP),' ',ICOLUP(2,NUP)
0298 C     &  ,' ',(PUP(J,NUP),J=1,5),' ',VTIMUP(NUP),' ',SPINUP(NUP)
0299 c        READ(LNHIN,*,END=900,ERR=900) IDUP(NUP),ISTUP(NUP),
0300 c     &  MOTHUP(1,NUP),MOTHUP(2,NUP),ICOLUP(1,NUP),ICOLUP(2,NUP),
0301 c     &  (PUP(J,NUP),J=1,5),VTIMUP(NUP),SPINUP(NUP)
0302 C...Reset resonance momentum to prepare for mass shifts
0303         IF(ISTUP(NUP).EQ.2) PUP(3,NUP)=0
0304         IF(ISTUP(NUP).EQ.1)THEN
0305            NEX=NEX+1
0306 C...Mrenna:  only if 4 < pdgId < 21
0307            IF(PUP(5,NUP).EQ.0D0.AND.IABS(IDUP(NUP)).GT.3
0308      $         .AND.IDUP(NUP).LT.21) THEN
0309 C...Set massless particle masses to Pythia default. Adjust z-momentum.
0310               PUP(5,NUP)=PMAS(IABS(PYCOMP(IDUP(NUP))),1)
0311               PUP(3,NUP)=SIGN(SQRT(MAX(0d0,PUP(4,NUP)**2-PUP(5,NUP)**2-
0312      $           PUP(1,NUP)**2-PUP(2,NUP)**2)),PUP(3,NUP))
0313            ENDIF
0314            PSUM=PSUM+PUP(3,NUP)
0315 C...Set mother resonance momenta
0316            MOTH=MOTHUP(1,NUP)
0317            DO WHILE (MOTH.GT.2)
0318              PUP(3,MOTH)=PUP(3,MOTH)+PUP(3,NUP)
0319              MOTH=MOTHUP(1,MOTH)
0320            ENDDO
0321         ENDIF
0322         NUP=NUP+1
0323   120 CONTINUE
0324       NUP=NUP-1
0325
0326 C...Increment event number
0327 C      IEVNT=IEVNT+1
0328
0329 C..Adjust mass of resonances
0330       DO I=1,NUP
0331          IF(ISTUP(I).EQ.2)THEN
0332             PUP(5,I)=SQRT(PUP(4,I)**2-PUP(1,I)**2-PUP(2,I)**2-
0333      $             PUP(3,I)**2)
0334          ENDIF
0335       ENDDO
0336
0337 C...Adjust energy and momentum of incoming particles
0338 C...In massive case need to solve quadratic equation
0339 c      PM1=PUP(5,1)**2
0340 c      PM2=PUP(5,2)**2
0341 c      A1=4d0*(ESUM**2-PSUM**2)
0342 c      A2=ESUM**2-PSUM**2+PM2-PM1
0343 c      A3=2d0*PSUM*A2
0344 c      A4=A3/A1
0345 c      A5=(A2**2-4d0*ESUM**2*PM2)/A1
0346 c
0347 c      PUP(3,2)=A4+SIGN(SQRT(A4**2+A5),PUP(3,2))
0348 c      PUP(3,1)=PSUM-PUP(3,2)
0349 c      PUP(4,1)=SQRT(PUP(3,1)**2+PM1)
0350 c      PUP(4,2)=SQRT(PUP(3,2)**2+PM2)
0351
0352       ESUM=PUP(4,1)+PUP(4,2)
0353
0354 C...Assuming massless incoming particles - otherwise Pythia adjusts
0355 C...the momenta to make them massless
0356 C     IF(IDBMUP(1).GT.100.AND.IDBMUP(2).GT.100)THEN
0357 C       DO I=1,2
0358 C          PUP(3,I)=0.5d0*(PSUM+SIGN(ESUM,PUP(3,I)))
0359 C          PUP(5,I)=0d0
0360 C        ENDDO
0361 C        PUP(4,1)=ABS(PUP(3,1))
0362 C        PUP(4,2)=ESUM-PUP(4,1)
0363 C      ENDIF
0364
0365 C...If you want to use some other scale for parton showering then the
0366 C...factorisation scale given by MadEvent, please implement the function PYMASC
0367 C...(example function included below)
0368
0369       IF(ickkw.eq.0.AND.MSCAL.GT.0) CALL PYMASC(SCALUP)
0370 c      IF(MINT(35).eq.3.AND.ickkw.EQ.1) SCALUP=SQRT(PARP(67))*SCALUP
0371
0372 C...Read FSR scale for all FS particles (as comment in event file)
0373 C      IF(ickkw.eq.1)THEN
0374 C        READ(LNHIN,*,END=900,ERR=130) CDUM,(PTPART(I),I=1,NEX)
0375 C 130    CONTINUE
0376 C      ENDIF
0377
0378       IF(ickkw.gt.0) THEN
0379 c
0380 c   Set up number of jets
0381 c
0382 C         write(*,*)'Setting up the number of jets'
0383          NLJETS=0
0384          NPART=0
0385 C         write(*,*)'Cycling on 3->NUP'
0386          do i=3,NUP
0387 C            write(*,*)'Iteration: i=',i
0388             if(ISTUP(i).ne.1) cycle
0389 C            write(*,*)'Npart++'
0390             NPART=NPART+1
0391             IPART(NPART)=i
0392             if(iabs(IDUP(i)).gt.nqmatch.and.IDUP(i).ne.21) cycle
0393             if(MOTHUP(1,i).gt.2) cycle
0394 C     Remove final-state partons that combine to color singlets
0395             IF((ABS(IDBMUP(1)).NE.11.OR.IDBMUP(1).NE.-IDBMUP(2)).AND.
0396      $           nosingrad) THEN
0397                DO II=3,NUP
0398                   IF(II.NE.i.AND.ISTUP(II).EQ.1)THEN
0399                      IF((IDUP(II).EQ.-IDUP(i).OR.
0400      $                    IDUP(i).EQ.21.AND.IDUP(II).EQ.21).AND.
0401      $                    ICOLUP(1,II).EQ.ICOLUP(2,i).AND.
0402      $                    ICOLUP(2,II).EQ.ICOLUP(1,i))then
0403 c                        print *,'Found color singlet'
0404                         CALL PYLIST(7)
0405                         GOTO 140
0406                      endif
0407                   ENDIF
0408                ENDDO
0409             ENDIF
0410             NLJETS=NLJETS+1
0411 C            WRITE(*,*) ' NLJETS=',NLJETS
0412             PTCLUS(NLJETS)=PTPART(NPART)
0413 c            print *,'   Adding a jet and NLJETS=',NLJETS,' and
0414 c     $        PTCLUS(',NLJETS,')=',PTCLUS(NLJETS)
0415  140        continue
0416          enddo
0417          CALL ALPSOR(PTCLUS,nljets,KP,1)
0418
0419          if(jetprocs) IDPRUP=LPRUP(NLJETS-MINJETS+1)
0420
0421          IF(ickkw.eq.1) THEN
0422 c   ... and decide whether exclusive or inclusive
0423             iexcl=iexclusive(IDPRUP)
0424             if((IEXCFILE.EQ.0.and.NLJETS.eq.MAXJETS.or.
0425      $           iexcl.eq.0).and.
0426      $           iexcl.ne.1)then
0427                IEXC=0
0428             else if(iexcl.eq.-1)then
0429                IEXC=-1
0430             else
0431                IEXC=1
0432             endif
0433          ENDIF
0434       ENDIF
0435 c      write( *,*)'finishing MGEVNT'
0436       RETURN
0437
0438 C...Error exit, typically when no more events.
0439 C  900 WRITE(*,*) ' Failed to read LHEF event information,'
0440 C      WRITE(*,*) ' assume end of file has been reached.'
0441 C      NUP=0
0442 C      MINT(51)=2
0443 C      write( *,*)'finishing MGEVNT'
0444 C      RETURN
0445       END
0446
0447 C*********************************************************************
0448 C...UPVETO
0449 C...Subroutine to implement the MLM jet matching criterion
0450 C*********************************************************************
0451       SUBROUTINE MGVETO(IPVETO)
0452
0453       IMPLICIT NONE
0454
0455
0456
0457
0458
0459 C...Pythia common blocks
0460       INTEGER MINT
0461       DOUBLE PRECISION VINT
0462       COMMON/PYINT1/MINT(400),VINT(400)
0463       INTEGER MSTP,MSTI
0464       DOUBLE PRECISION PARP,PARI
0465       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0466
0467 C...GUP Event common block
0468       INTEGER MAXNUP
0469       PARAMETER (MAXNUP=500)
0470       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0471       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0472       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
0473      &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
0474      &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
0475      &              SPINUP(MAXNUP)
0476 C...User process initialization commonblock.
0477       INTEGER MAXPUP
0478       PARAMETER (MAXPUP=100)
0479       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0480       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0481       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0482      &   IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0483      &   LPRUP(MAXPUP)
0484 C...HEPEVT commonblock.
0485       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0486       PARAMETER (NMXHEP=4000)
0487       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0488      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0489       DOUBLE PRECISION PHEP,VHEP
0490       SAVE /HEPEVT/
0491       INTEGER IPVETO
0492 C...GETJET commonblocks
0493       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
0494       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
0495      &  SPHCAL,PCJET,ETJET
0496       PARAMETER (MNCY=200)
0497       PARAMETER (MNCPHI=200)
0498       COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
0499      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
0500      $YCMIN,YCMAX,NCY,NCPHI
0501       PARAMETER (NJMAX=500)
0502       COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
0503      $NCJET
0504       DOUBLE PRECISION PI
0505       PARAMETER (PI=3.141593D0)
0506 C
0507       DOUBLE PRECISION PSERAP
0508       INTEGER K(NJMAX),KP(NJMAX),kpj(njmax)
0509
0510 C...Variables for the kT-clustering
0511       INTEGER NMAX,NN,NSUB,JET,NJETM,IHARD,IP1,IP2
0512       DOUBLE PRECISION PP,PJET
0513       DOUBLE PRECISION ECUT,Y,YCUT
0514       PARAMETER (NMAX=512)
0515       DIMENSION JET(NMAX),Y(NMAX),PP(4,NMAX),PJET(4,NMAX)
0516       INTEGER NNM
0517       DOUBLE PRECISION YM(NMAX),PPM(4,NMAX)
0518
0519 C...kt clustering common block
0520       INTEGER NMAXKT,NUM,HIST
0521       PARAMETER (NMAXKT=512)
0522       DOUBLE PRECISION PPP,KT,ETOT,RSQ,KTP,KTS,KTLAST
0523       COMMON /KTCOMM/ETOT,RSQ,PPP(9,NMAXKT),KTP(NMAXKT,NMAXKT),
0524      $   KTS(NMAXKT),KT(NMAXKT),KTLAST(NMAXKT),HIST(NMAXKT),NUM
0525
0526 C...Extra commonblock to transfer run info.
0527       INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0528       COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
0529
0530 C...Inputs for the matching algorithm
0531 C   clfact determines how jet-to parton matching is done
0532 C   kt-jets: default=1
0533 C    clfact >= 0: Max mult. if within clfact*max(qcut,Q(partNmax)) from jet, others within clfact*qcut
0534 C    clfact < 0: Max mult. if within |clfact|*Q(jetNmax) from jet, other within |clfact|*qcut
0535 C   cone-jets: default=1.5
0536 C    Matching if within clfact*RCLMAX
0537
0538 C...Inputs for the matching algorithm
0539       double precision etcjet,rclmax,etaclmax,qcut,clfact,showerkt
0540       integer maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres(30)
0541       integer nqmatch,nexcproc,iexcproc(MAXPUP),iexcval(MAXPUP)
0542       logical nosingrad,jetprocs
0543       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,showerkt,clfact,
0544      $   maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres,
0545      $   nqmatch,nexcproc,iexcproc,iexcval,nosingrad,jetprocs
0546
0547 C...Commonblock to transfer event-by-event matching info
0548       INTEGER NLJETS,IEXC,Ifile
0549       DOUBLE PRECISION PTCLUS
0550       COMMON/MEMAEV/PTCLUS(20),NLJETS,IEXC,Ifile
0551
0552       INTEGER nvarev,nvar2
0553       PARAMETER (nvarev=57,nvar2=6)
0554
0555       REAL*4 varev(nvarev)
0556       COMMON/HISTDAT/varev
0557
0558 C...Pythia common blocks
0559       INTEGER NPART,NPARTD,IPART,MAXNUR
0560           DOUBLE PRECISION PTPART
0561       PARAMETER (MAXNUR=1000)
0562       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0563
0564           INTEGER flag
0565           COMMON/OUTTREE/flag
0566
0567           CHARACTER*8 htit(nvarev),htit2(nvar2)
0568       DATA htit/'Npart','Qjet1','Qjet2','Qjet3','Qjet4',
0569      $   'Ptcjet1','Ptcjet2','Ptcjet3','Ptcjet4',
0570      $   'Etacjet1','Etacjet2','Etacjet3','Etacjet4',
0571      $   'Phicjet1','Phicjet2','Phicjet3','Phicjet4',
0572      $   'Ptjet1','Ptjet2','Ptjet3','Ptjet4',
0573      $   'Etajet1','Etajet2','Etajet3','Etajet4',
0574      $   'Phijet1','Phijet2','Phijet3','Phijet4',
0575      $   'Idres1','Ptres1','Etares1','Phires1',
0576      $   'Idres2','Ptres2','Etares2','Phires2',
0577      $   'Ptlep1','Etmiss','Htjets',
0578      $   'Ptb','Etab','Ptbbar','Etabbar','Ptbj','Etabj',
0579      $   'Qpar1','Qpar2','Qpar3','Qpar4',
0580      $   'Ptpar1','Ptpar2','Ptpar3','Ptpar4',
0581      $   'Ncjets','Njets','Nfile'/
0582       DATA htit2/'Npart','Qjet1','Qjet2','Qjet3','Qjet4','Nfile'/
0583
0584
0585
0586
0587
0588 C   local variables
0589       double precision tiny
0590       parameter (tiny=1d-3)
0591       integer idbg
0592       data idbg/0/
0593
0594       integer i,j,ihep,nmatch,jrmin,KPT(MAXNUP),nres,ii
0595       double precision etajet,phijet,delr,dphi,delrmin,ptjet
0596       double precision p(4,10),pt(10),eta(10),phi(10)
0597       INTEGER IMO
0598       logical norad(20)
0599       REAL*4 var2(nvar2)
0600
0601 c      if(NLJETS.GT.0)then
0602 c        idbg=1
0603 c      else
0604 c        idbg=0
0605 c      endif
0606
0607
0608 c      write(*,*)'Entering MGVETO'
0609 c      write(*,*)'qcut is ',qcut,' and showerkt is ',showerkt
0610       IPVETO=0
0611       YCUT=-1.0
0612 c     Return if not MLM matching (or non-matched subprocess)
0613
0614       IF(ICKKW.LE.0.OR.IEXC.eq.-1) RETURN
0615
0616       IF(NLJETS.LT.MINJETS.OR.NLJETS.GT.MAXJETS)THEN
0617         if(idbg.eq.1)
0618      $     WRITE(LNHOUT,*) 'Failed due to NLJETS ',NLJETS,' < ',MINJETS,
0619      $        ' or > ',MAXJETS
0620          GOTO 999
0621       ENDIF
0622
0623 C      write(*,*)'Throw event if it contains an excluded resonance'
0624       NRES=0
0625       DO I=1,NUP
0626 c            write(*,*)'cycling on particles'
0627         IF(ISTUP(I).EQ.2)THEN
0628 c                       write(*,*)'found S2, now comparin with the ',nexcres,' excres'
0629            DO J=1,nexcres
0630 c                       write(*,*)'comparing ',IDUP(I),' and ',EXCRES(J)
0631               IF(IDUP(I).EQ.EXCRES(J)) NRES=NRES+1
0632            ENDDO
0633         ENDIF
0634       ENDDO
0635       IF(NRES.GT.0)THEN
0636 c               write(*,*)'Event',IEVNT,
0637 c     &  ' thrown because of ',NRES,'e r'
0638 c     CALL PYLIST(7)
0639          GOTO 999
0640       ENDIF
0641
0642 c init uninit variables
0643       jrmin = 0
0644
0645 C   Set up vetoed mothers
0646 c      DO I=1,MAXNUP
0647 c        INORAD(I)=0
0648 c      ENDDO
0649 c      DO IHEP=1,NUP-2
0650 c        if(ISTHEP(ihep).gt.1.and.iabs(IDHEP(ihep)).gt.8) then
0651 c        if(iabs(IDHEP(ihep)).gt.5.and.IDHEP(ihep).ne.21) then
0652 c          INORAD(ihep)=1
0653 c        endif
0654 c      ENDDO
0655
0656 c
0657 c     reconstruct parton-level event
0658 c     Set norad for daughters of decayed particles, to not include
0659 c     radiation from these in matched jets
0660 c
0661       if(idbg.eq.1) then
0662         write(LNHOUT,*) ' '
0663         write(LNHOUT,*) 'new event '
0664 c        CALL PYLIST(1)
0665         CALL PYLIST(7)
0666         CALL PYLIST(5)
0667         write(LNHOUT,*) 'PARTONS'
0668       endif
0669       i=0
0670       do ihep=3,nup
0671          NORAD(ihep)=.false.
0672         if((ABS(IDBMUP(1)).NE.11.OR.IDBMUP(1).NE.-IDBMUP(2)).AND.
0673      $        MOTHUP(1,ihep).gt.2) goto 100
0674         if(ISTUP(ihep).ne.1.or.
0675      $     (iabs(IDUP(ihep)).gt.nqmatch.and.IDUP(ihep).ne.21)) cycle
0676 c     If quark or gluon making singlet system with other final-state parton
0677 c     remove (since unseen singlet resonance) unless e+e- collision
0678         IF((ABS(IDBMUP(1)).NE.11.OR.IDBMUP(1).NE.-IDBMUP(2)).AND.
0679      $       nosingrad)THEN
0680            DO II=3,NUP
0681               IF(II.NE.ihep.AND.ISTUP(II).EQ.1)THEN
0682                  IF((IDUP(II).EQ.-IDUP(ihep).OR.
0683      $                IDUP(ihep).EQ.21.AND.IDUP(II).EQ.21).AND.
0684      $                ICOLUP(1,II).EQ.ICOLUP(2,ihep).AND.
0685      $                ICOLUP(2,II).EQ.ICOLUP(1,ihep))
0686      $                GOTO 100
0687               ENDIF
0688            ENDDO
0689         ENDIF
0690         i=i+1
0691         do j=1,4
0692           p(j,i)=pup(j,ihep)
0693         enddo
0694         pt(i)=sqrt(p(1,i)**2+p(2,i)**2)
0695         if(i.LE.4) varev(50+i)=pt(i)
0696         eta(i)=-log(tan(0.5d0*atan2(pt(i)+tiny,p(3,i))))
0697         phi(i)=atan2(p(2,i),p(1,i))
0698         if(idbg.eq.1) then
0699           write(LNHOUT,*) pt(i),eta(i),phi(i)
0700         endif
0701         cycle
0702  100    norad(ihep)=.true.
0703       enddo
0704       if(i.ne.NLJETS)then
0705         print *,'Error in UPVETO: Wrong number of jets found ',i,NLJETS
0706         CALL PYLIST(7)
0707         CALL PYLIST(2)
0708         stop
0709       endif
0710 C Bubble-sort PTs in descending order
0711       DO I=1,3
0712          DO J=4,I+1,-1
0713             IF(varev(50+J).GT.varev(50+I))THEN
0714                PTJET=varev(50+J)
0715                varev(50+J)=varev(50+I)
0716                varev(50+I)=PTJET
0717             ENDIF
0718          ENDDO
0719       ENDDO
0720 C     Set status for non-clustering partons to 2
0721       DO ihep=1,NHEP
0722 c         ISTORG(ihep)=ISTHEP(ihep)
0723          IF(ISTHEP(ihep).EQ.1.AND.iabs(IDHEP(ihep)).GT.5.AND.
0724      $        IDHEP(ihep).NE.21) THEN
0725             ISTHEP(ihep)=2
0726          ELSEIF(ISTHEP(ihep).EQ.1.AND.JMOHEP(1,ihep).GT.0) then
0727             IMO=JMOHEP(1,ihep)
0728             DO WHILE(IMO.GT.0)
0729 c           Trace mothers, if non-radiating => daughter is decay - remove
0730               IF(IMO.le.NUP-2.and.norad(IMO+2)) GOTO 105
0731               IMO=JMOHEP(1,IMO)
0732             ENDDO
0733             cycle
0734  105        ISTHEP(ihep)=2
0735          ENDIF
0736       ENDDO
0737
0738 c      DO ihep=1,NHEP
0739 c            print *,'Part ',ihep,' status=',ISTHEP(ihep),'
0740 c     $   PID=',iabs(IDHEP(ihep)),' mother number=',
0741 c     $  JMOHEP(1,ihep),' status=',ISTHEP(JMOHEP(1,
0742 c     $     ihep)),' PID=',IDHEP(JMOHEP(1,ihep))
0743 c      ENDDO
0744
0745       DO ihep=1,NHEP
0746
0747         if ( jmohep(1,ihep) .gt. 0 ) then
0748 c
0749 c         If valid mother and status is 2 and a mother of 6<PID>nqmatch =>reject from particle list
0750 c
0751          IF(ISTHEP(JMOHEP(1,ihep)).EQ.2
0752      $    .AND.iabs(IDHEP(JMOHEP(1,ihep))).GT.nqmatch.AND.
0753      $    iabs(IDHEP(JMOHEP(1,ihep))).LT.6) THEN
0754 c            print *,'Have found: part ',ihep,' status=',ISTHEP(ihep),
0755 c     $      'PID=',iabs(IDHEP(ihep)),' mother number=',
0756 c     $      JMOHEP(1,ihep),' status=',ISTHEP(JMOHEP(1,
0757 c     $      ihep)),' PID=',IDHEP(JMOHEP(1,ihep))
0758           ISTHEP(ihep)=2
0759          ENDIF
0760          IF(ISTHEP(ihep).eq.1.AND.iabs(IDHEP(ihep)).GT.
0761      $     nqmatch.AND.iabs(IDHEP(ihep)).LT.6.AND.
0762      $     ISTHEP(JMOHEP(1,ihep)).EQ.2.AND.iabs(IDHEP(JMOHEP(1,ihep)))
0763      $     .EQ.21) goto 999
0764 c
0765         endif
0766
0767       ENDDO
0768 c
0769 c      DO ihep=1,NHEP
0770 c          IF(ISTHEP(ihep).EQ.1)print *,'After selection:  Part ',
0771 c     $ ihep,' status=',ISTHEP(ihep),'PID=',iabs(IDHEP(ihep))
0772 c     $   ,' mother number=',JMOHEP(1,ihep),' status=',
0773 c     $   ISTHEP(JMOHEP(1,ihep)),' PID=',IDHEP(JMOHEP(1,ihep))
0774 c       ENDDO
0775
0776
0777
0778 C     Prepare histogram filling
0779         DO I=1,4
0780           var2(1+I)=-1
0781           varev(46+I)=-1
0782           varev(50+I)=-1
0783         ENDDO
0784
0785       I=0
0786       if(idbg.eq.1) then
0787         do i=1,nhep
0788           write(LNHOUT,1000)i,isthep(i),idhep(i),jmohep(1,i),jmohep(2,i)
0789      $         ,phep(1,i),phep(2,i),phep(3,i)
0790         enddo
0791  1000   format(5(i4,1x),3(f12.5,1x))
0792       endif
0793
0794       IF(ICKKW.EQ.2) GOTO 150
0795       IF(MSTP(61).eq.0..and.MSTP(71).eq.0)then
0796 c      write(*,*)'No showering - just print out event'
0797       ELSE IF(qcut.le.0d0)then
0798 c      write(*,*)'qcut<0'
0799
0800       IF(clfact.EQ.0d0) clfact=1.5d0
0801
0802 c      CALL PYLIST(7)
0803 c      CALL PYLIST(2)
0804 c      CALL PYLIST(5)
0805 c     Start from the partonic system
0806       IF(NLJETS.GT.0) CALL ALPSOR(pt,nljets,KP,2)
0807 c     reconstruct showered jets
0808 c
0809       YCMAX=ETACLMAX+RCLMAX
0810       YCMIN=-YCMAX
0811       CALL CALINIM
0812       CALL CALDELM(1,1)
0813       CALL GETJETM(RCLMAX,ETCJET,ETACLMAX)
0814 c     analyse only events with at least nljets-reconstructed jets
0815       IF(NCJET.GT.0) CALL ALPSOR(ETJET,NCJET,K,2)
0816       if(idbg.eq.1) then
0817         write(LNHOUT,*) 'JETS'
0818         do i=1,ncjet
0819           j=k(ncjet+1-i)
0820           ETAJET=PSERAP(PCJET(1,j))
0821           PHIJET=ATAN2(PCJET(2,j),PCJET(1,j))
0822           write(LNHOUT,*) etjet(j),etajet,phijet
0823         enddo
0824       endif
0825       IF(NCJET.LT.NLJETS) THEN
0826         if(idbg.eq.1)
0827      $     WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' < ',NLJETS
0828         GOTO 999
0829       endif
0830 c     associate partons and jets, using min(delr) as criterion
0831       NMATCH=0
0832       DO I=1,NCJET
0833         KPJ(I)=0
0834       ENDDO
0835       DO I=1,NLJETS
0836         DELRMIN=1D5
0837         DO 110 J=1,NCJET
0838           IF(KPJ(J).NE.0) GO TO 110
0839           ETAJET=PSERAP(PCJET(1,J))
0840           PHIJET=ATAN2(PCJET(2,J),PCJET(1,J))
0841           DPHI=ABS(PHI(KP(NLJETS-I+1))-PHIJET)
0842           IF(DPHI.GT.PI) DPHI=2.*PI-DPHI
0843           DELR=SQRT((ETA(KP(NLJETS-I+1))-ETAJET)**2+(DPHI)**2)
0844           IF(DELR.LT.DELRMIN) THEN
0845             DELRMIN=DELR
0846             JRMIN=J
0847           ENDIF
0848  110    CONTINUE
0849         IF(DELRMIN.LT.clfact*RCLMAX) THEN
0850           NMATCH=NMATCH+1
0851           KPJ(JRMIN)=I
0852         ENDIF
0853 C     WRITE(*,*) 'PARTON-JET',I,' best match:',k(ncjet+1-jrmin)
0854 c     $           ,delrmin
0855       ENDDO
0856       IF(NMATCH.LT.NLJETS)  THEN
0857         if(idbg.eq.1)
0858      $     WRITE(LNHOUT,*) 'Failed due to NMATCH ',NMATCH,' < ',NLJETS
0859         GOTO 999
0860       endif
0861 C REJECT EVENTS WITH LARGER JET MULTIPLICITY FROM EXCLUSIVE SAMPLE
0862       IF(NCJET.GT.NLJETS.AND.IEXC.EQ.1)  THEN
0863         if(idbg.eq.1)
0864      $     WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' > ',NLJETS
0865         GOTO 999
0866       endif
0867 C     VETO EVENTS WHERE MATCHED JETS ARE SOFTER THAN NON-MATCHED ONES
0868       IF(IEXC.NE.1) THEN
0869         J=NCJET
0870         DO I=1,NLJETS
0871           IF(KPJ(K(J)).EQ.0) GOTO 999
0872           J=J-1
0873         ENDDO
0874       ENDIF
0875
0876       else                      ! qcut.gt.0
0877       if(showerkt.eq.1.0) then
0878 c      write(*,*)"qcut>=0 and showerkt=1 ==> Veto events where
0879 c     & first shower emission has kt > YCUT"
0880
0881
0882         IF(NLJETS.EQ.0)THEN
0883            VINT(358)=0
0884         ENDIF
0885
0886         IF(idbg.eq.1) THEN
0887 C           PRINT *,'Using shower emission pt method'
0888 C           write(*,*)'Using shower emission pt method'
0889 C           write(*,*)'qcut, ptclus(1), vint(357),vint(358),vint(360): ',
0890 C     $          qcut,ptclus(1),vint(357),vint(358),vint(360)
0891
0892 C      PRINT *,'qcut, ptclus(1), vint(357),vint(358),vint(360): ',
0893 C     $          qcut,ptclus(1),vint(357),vint(358),vint(360)
0894         ENDIF
0895         YCUT=qcut**2
0896
0897         IF(NLJETS.GT.0.AND.PTCLUS(1)**2.LT.YCUT) THEN
0898           if(idbg.eq.1)
0899      $       WRITE(LNHOUT,*) 'Failed due to KT ',
0900      $       PTCLUS(1),' < ',SQRT(YCUT)
0901           GOTO 999
0902         ENDIF
0903 c        PRINT *,'Y,VINT:',SQRT(Y(NLJETS+1)),SQRT(VINT(390))
0904 C        write(*,*)'Y,VINT:',SQRT(Y(NLJETS+1)),SQRT(VINT(390))
0905 C        write(*,*)'mektsc 357, 358: ',mektsc,' ',VINT(357),' ',VINT(358)
0906         IF(IEXC.EQ.1.AND.
0907      $       ((mektsc.eq.1.and.MAX(VINT(357),VINT(358)).GT.SQRT(YCUT))
0908      $       .OR.
0909      $       (mektsc.eq.2.and.MAX(VINT(360),VINT(358)).GT.SQRT(YCUT))))
0910      $       THEN
0911 C            write(*,*)'rejection'
0912           if(idbg.eq.1)
0913      $       WRITE(LNHOUT,*)
0914      $       'Failed due to ',max(VINT(357),VINT(358)),' > ',SQRT(YCUT)
0915           GOTO 999
0916         ENDIF
0917 C        PRINT *,NLJETS,IEXC,SQRT(VINT(390)),PTCLUS(1),SQRT(YCUT)
0918 C        write(*,*)'NLJets, iexc, VINT, ptclus(1), sqrt(ycut)',NLJETS
0919 c     &,IEXC,SQRT(VINT(390)),PTCLUS(1),SQRT(YCUT)
0920 c     Highest multiplicity case
0921         IF(IEXC.EQ.0.AND.NLJETS.GT.0.AND.
0922      $       ((mektsc.eq.1.and.MAX(VINT(357),VINT(358)).GT.PTCLUS(1))
0923      $       .OR.
0924      $       (mektsc.eq.2.and.MAX(VINT(360),VINT(358)).GT.PTCLUS(1))))
0925      $       THEN
0926 c     $     VINT(390).GT.PTCLUS(1)**2)THEN
0927           if(idbg.eq.1)
0928      $       WRITE(LNHOUT,*)
0929      $       'Failed due to ',max(VINT(357),VINT(358)),' > ',PTCLUS(1)
0930           GOTO 999
0931         ENDIF
0932 c
0933       else                      ! not shower kt method
0934
0935         IF(clfact.EQ.0d0) clfact=1d0
0936
0937 C---FIND FINAL STATE COLOURED PARTICLES
0938         NN=0
0939         DO IHEP=1,NHEP
0940           IF (ISTHEP(IHEP).EQ.1
0941      $       .AND.(ABS(IDHEP(IHEP)).LE.5.OR.IDHEP(IHEP).EQ.21)) THEN
0942             PTJET=sqrt(PHEP(1,IHEP)**2+PHEP(2,IHEP)**2)
0943             ETAJET=ABS(LOG(MIN((SQRT(PTJET**2+PHEP(3,IHEP)**2)+
0944      $       ABS(PHEP(3,IHEP)))/PTJET,1d5)))
0945             IF(ETAJET.GT.etaclmax) cycle
0946             NN=NN+1
0947             IF (NN.GT.NMAX) then
0948               CALL PYLIST(2)
0949               PRINT *, 'Too many particles: ', NN
0950               NN=NN-1
0951               GOTO 120
0952             endif
0953             DO I=1,4
0954               PP(I,NN)=PHEP(I,IHEP)
0955             ENDDO
0956           ELSE if(idbg.eq.1)THEN
0957             PRINT *,'Skipping particle ',IHEP,ISTHEP(IHEP),IDHEP(IHEP)
0958           ENDIF
0959         ENDDO
0960
0961
0962 C...Cluster event to find values of Y including jet matching but not veto of too many jets
0963 C...Only used to fill the beforeveto Root tree
0964  120    ECUT=1
0965         IF (NN.GT.1) then
0966           CALL KTCLUS(KTSCHE,PP,NN,ECUT,Y,*999)
0967           if(idbg.eq.1)
0968      $       WRITE(LNHOUT,*) 'Clustering values:',
0969      $       (SQRT(Y(i)),i=1,MIN(NN,3))
0970
0971 C       Print out values in the case where all jets are matched at the
0972 C       value of the NLJETS:th clustering
0973         var2(1)=NLJETS
0974         var2(6)= Ifile
0975
0976         if(NLJETS.GT.MINJETS)then
0977           YCUT=Y(NLJETS)
0978           CALL KTRECO(MOD(KTSCHE,10),PP,NN,ECUT,YCUT,YCUT,PJET,JET,
0979      $       NCJET,NSUB,*999)
0980
0981 C     Cluster jets with first hard parton
0982           DO I=1,NLJETS
0983             DO J=1,4
0984               PPM(J,I)=PJET(J,I)
0985             ENDDO
0986           ENDDO
0987
0988           NJETM=NLJETS
0989           DO IHARD=1,NLJETS
0990             NNM=NJETM+1
0991             DO J=1,4
0992               PPM(J,NNM)=p(J,IHARD)
0993             ENDDO
0994             CALL KTCLUS(KTSCHE,PPM,NNM,ECUT,YM,*999)
0995             IF(YM(NNM).GT.YCUT) THEN
0996 C       Parton not clustered
0997               GOTO 130
0998             ENDIF
0999
1000 C       Find jet clustered with parton
1001
1002             IP1=HIST(NNM)/NMAXKT
1003             IP2=MOD(HIST(NNM),NMAXKT)
1004             IF(IP2.NE.NNM.OR.IP1.LE.0)THEN
1005               GOTO 130
1006             ENDIF
1007             DO I=IP1,NJETM-1
1008               DO J=1,4
1009                 PPM(J,I)=PPM(J,I+1)
1010               ENDDO
1011             ENDDO
1012             NJETM=NJETM-1
1013           ENDDO                 ! IHARD=1,NLJETS
1014         endif                   ! NLJETS.GT.MINJETS
1015
1016         DO I=1,MIN(NN,4)
1017           var2(1+I)=SQRT(Y(I))
1018         ENDDO
1019         WRITE(15,4001) (var2(I),I=1,nvar2)
1020
1021  130    CONTINUE
1022 C   Now perform jet clustering at the value chosen in qcut
1023
1024         CALL KTCLUS(KTSCHE,PP,NN,ECUT,Y,*999)
1025
1026         YCUT=qcut**2
1027         NCJET=0
1028
1029 C     Reconstruct jet momenta
1030           CALL KTRECO(MOD(KTSCHE,10),PP,NN,ECUT,YCUT,YCUT,PJET,JET,
1031      $       NCJET,NSUB,*999)
1032
1033         ELSE IF (NN.EQ.1) THEN
1034
1035           Y(1)=PP(1,1)**2+PP(2,1)**2
1036           IF(Y(1).GT.YCUT)THEN
1037             NCJET=1
1038             DO I=1,4
1039               PJET(I,1)=PP(I,1)
1040             ENDDO
1041           ENDIF
1042         endif
1043
1044         if(idbg.eq.1) then
1045           write(LNHOUT,*) 'JETS'
1046           do i=1,ncjet
1047             PTJET =SQRT(PJET(1,i)**2+PJET(2,i)**2)
1048             ETAJET=PSERAP(PJET(1,i))
1049             PHIJET=ATAN2(PJET(2,i),PJET(1,i))
1050             write(LNHOUT,*) ptjet,etajet,phijet
1051           enddo
1052         endif
1053
1054         IF(NCJET.LT.NLJETS) THEN
1055           if(idbg.eq.1)
1056      $       WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' < ',NLJETS
1057           GOTO 999
1058         endif
1059
1060 C...Right number of jets - but the right jets?
1061 C     For max. multiplicity case, count jets only to the NHARD:th jet
1062         IF(IEXC.EQ.0)THEN
1063            IF(NLJETS.GT.0)THEN
1064               YCUT=Y(NLJETS)
1065               CALL KTRECO(MOD(KTSCHE,10),PP,NN,ECUT,YCUT,YCUT,PJET,JET,
1066      $             NCJET,NSUB,*999)
1067               IF(clfact.GE.0d0) THEN
1068                  CALL ALPSOR(PTCLUS,nljets,KPT,2)
1069                  YCUT=MAX(qcut,PTCLUS(KPT(1)))**2
1070               ENDIF
1071            ENDIF
1072         ELSE IF(NCJET.GT.NLJETS) THEN
1073            if(idbg.eq.1)
1074      $       WRITE(LNHOUT,*) 'Failed due to NCJET ',NCJET,' > ',NLJETS
1075            GOTO 999
1076         ENDIF
1077 C     Cluster jets with hard partons, one at a time
1078         DO I=1,NLJETS
1079           DO J=1,4
1080             PPM(J,I)=PJET(J,I)
1081           ENDDO
1082         ENDDO
1083
1084         NJETM=NLJETS
1085         IF(clfact.NE.0) YCUT=clfact**2*YCUT
1086 c        YCUT=qcut**2
1087 c        YCUT=(1.5*qcut)**2
1088
1089         DO 140 IHARD=1,NLJETS
1090           NN=NJETM+1
1091           DO J=1,4
1092             PPM(J,NN)=p(J,IHARD)
1093           ENDDO
1094           CALL KTCLUS(KTSCHE,PPM,NN,ECUT,Y,*999)
1095
1096           IF(Y(NN).GT.YCUT) THEN
1097 C       Parton not clustered
1098           if(idbg.eq.1)
1099      $       WRITE(LNHOUT,*) 'Failed due to parton ',IHARD,
1100      $         ' not clustered: ',Y(NN)
1101             GOTO 999
1102           ENDIF
1103
1104 C       Find jet clustered with parton
1105
1106           IP1=HIST(NN)/NMAXKT
1107           IP2=MOD(HIST(NN),NMAXKT)
1108           IF(IP2.NE.NN.OR.IP1.LE.0)THEN
1109           if(idbg.eq.1)
1110      $       WRITE(LNHOUT,*) 'Failed due to parton ',IHARD,
1111      $         ' not clustered: ',IP1,IP2,NN,HIST(NN)
1112             GOTO 999
1113           ENDIF
1114 C     Remove jet clustered with parton
1115           DO I=IP1,NJETM-1
1116             DO J=1,4
1117               PPM(J,I)=PPM(J,I+1)
1118             ENDDO
1119           ENDDO
1120           NJETM=NJETM-1
1121  140   CONTINUE
1122
1123       endif                     ! pt-ordered showers
1124       endif                     ! qcut.gt.0
1125 C...Cluster particles with |eta| < etaclmax for histograms
1126  150  NN=0
1127       DO IHEP=1,NHEP
1128          IF (ISTHEP(IHEP).EQ.1
1129      $        .AND.(ABS(IDHEP(IHEP)).LE.5.OR.IDHEP(IHEP).EQ.21)) THEN
1130             PTJET=sqrt(PHEP(1,IHEP)**2+PHEP(2,IHEP)**2)
1131             ETAJET=ABS(LOG(MIN((SQRT(PTJET**2+PHEP(3,IHEP)**2)+
1132      $           ABS(PHEP(3,IHEP)))/PTJET,1d5)))
1133             IF(ETAJET.GT.etaclmax) cycle
1134             NN=NN+1
1135             IF (NN.GT.NMAX) then
1136                CALL PYLIST(2)
1137                PRINT *, 'Too many particles: ', NN
1138                NN=NN-1
1139                GOTO 160
1140             ENDIF
1141             DO I=1,4
1142                PP(I,NN)=PHEP(I,IHEP)
1143             ENDDO
1144          ELSE if(idbg.eq.1)THEN
1145             PRINT *,'Skipping particle ',IHEP,ISTHEP(IHEP),IDHEP(IHEP)
1146          ENDIF
1147       ENDDO
1148
1149  160  ECUT=1
1150       IF (NN.GT.1) THEN
1151          CALL KTCLUS(KTSCHE,PP,NN,ECUT,Y,*999)
1152       ELSE IF(NN.EQ.1) THEN
1153          Y(1)=PP(1,NN)**2+PP(2,NN)**2
1154       ENDIF
1155
1156       DO I=1,MIN(NN,4)
1157          varev(46+I)=SQRT(Y(I))
1158       ENDDO
1159
1160 c      write(*,*)' finishing up mgveto, with ipveto= ', ipveto
1161         OPEN (10, FILE='events.tree')
1162 c       WRITE(10,'(a)') '# File with ntuple events with the variables:'
1163 c      WRITE(10,CGIVE0) (htit(I)(1:len_trim(htit(I))),I=1,nvarev)
1164 c      write (*,'(a)'),(htit(I)(1:len_trim(htit(I))),I=47,50)
1165 c      write (*,4001),(varev(I),I=47,50)
1166        if (flag.eq.1) then
1167           varev(1)=NLJETS
1168           WRITE(10,4001) varev(1),(varev(I),I=47,50)
1169 c          WRITE(*,4001) varev(1),(varev(I),I=47,50)
1170       endif
1171
1172       RETURN
1173  4001 FORMAT(50E15.6)
1174 c HERWIG/PYTHIA TERMINATION:
1175
1176
1177
1178  999  IPVETO=1
1179 c      write(*,*)' finishing up mgveto, with ipveto= ', ipveto
1180       END
1181
1182 C*********************************************************************
1183 C   PYMASC
1184 C   Implementation of scale used in Pythia parton showers
1185 C*********************************************************************
1186       SUBROUTINE PYMASC(scale)
1187       IMPLICIT NONE
1188
1189 C...Arguments
1190       REAL*8 scale
1191
1192 C...Functions
1193       REAL*8 SMDOT5
1194
1195 C...User process initialization commonblock.
1196       INTEGER MAXPUP
1197       PARAMETER (MAXPUP=100)
1198       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
1199       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
1200       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
1201      &   IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
1202      &   LPRUP(MAXPUP)
1203 C...User process event common block.
1204       INTEGER MAXNUP
1205       PARAMETER (MAXNUP=500)
1206       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
1207       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
1208       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
1209      &   ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
1210      &   VTIMUP(MAXNUP),SPINUP(MAXNUP)
1211
1212 C...Extra commonblock to transfer run info.
1213       INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
1214       COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
1215
1216 C...Local variables
1217       INTEGER ICC1,ICC2,IJ,IDC1,IDC2,IC,IC1,IC2
1218       REAL*8 QMIN,QTMP
1219
1220 C   Just use the scale read off the event record
1221       scale=SCALUP
1222
1223 C   Alternatively:
1224
1225 C...  Guesses for the correct scale
1226 C     Assumptions:
1227 C     (1) if the initial state is a color singlet, then
1228 C     use s-hat for the scale
1229 C
1230 C     (2) if color flow to the final state, use the minimum
1231 C     of the dot products of color connected pairs
1232 C     (times two for consistency with above)
1233
1234         QMIN=SMDOT5(PUP(1,1),PUP(1,2))
1235         ICC1=1
1236         ICC2=2
1237 C
1238 C     For now, there is no generic way to guarantee the "right"
1239 C     scale choice.  Here, we take the HERWIG pt. of view and
1240 C     choose the dot product of the colored connected "primary"
1241 C     pairs.
1242 C
1243
1244         DO 101 IJ=1,NUP
1245           IF(MOTHUP(2,IJ).GT.2) GOTO 101
1246           IDC1=ICOLUP(1,IJ)
1247           IDC2=ICOLUP(2,IJ)
1248           IF(IDC1.EQ.0) IDC1=-1
1249           IF(IDC2.EQ.0) IDC2=-2
1250
1251           DO 201 IC=IJ+1,NUP
1252             IF(MOTHUP(2,IC).GT.2) GOTO 201
1253             IC1=ICOLUP(1,IC)
1254             IC2=ICOLUP(2,IC)
1255             IF(ISTUP(IC)*ISTUP(IJ).GE.1) THEN
1256               IF(IDC1.EQ.IC2.OR.IDC2.EQ.IC1) THEN
1257                 QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC))
1258                 IF(QTMP.LT.QMIN) THEN
1259                   QMIN=QTMP
1260                   ICC1=IJ
1261                   ICC2=IC
1262                 ENDIF
1263               ENDIF
1264             ELSEIF(ISTUP(IC)*ISTUP(IJ).LE.-1) THEN
1265               IF(IDC1.EQ.IC1.OR.IDC2.EQ.IC2) THEN
1266                 QTMP=SMDOT5(PUP(1,IJ),PUP(1,IC))
1267                 IF(QTMP.LT.QMIN) THEN
1268                   QMIN=QTMP
1269                   ICC1=IJ
1270                   ICC2=IC
1271                 ENDIF
1272               ENDIF
1273             ENDIF
1274  201      CONTINUE
1275  101    CONTINUE
1276
1277         scale=QMIN
1278
1279       RETURN
1280       END
1281
1282 C...SMDOT5
1283 C   Helper function
1284
1285       FUNCTION SMDOT5(V1,V2)
1286       IMPLICIT NONE
1287       REAL*8 SMDOT5,TEMP
1288       REAL*8 V1(5),V2(5)
1289       INTEGER I
1290
1291       SMDOT5=0D0
1292       TEMP=V1(4)*V2(4)
1293       DO I=1,3
1294         TEMP=TEMP-V1(I)*V2(I)
1295       ENDDO
1296
1297       SMDOT5=SQRT(ABS(TEMP))
1298
1299       RETURN
1300       END
1301
1302 C*********************************************************************
1303 C...set_matching
1304 C...Sets parameters for the matching, i.e. cuts and jet multiplicities
1305 C*********************************************************************
1306
1307       SUBROUTINE set_matching(npara,param,value)
1308       implicit none
1309 c
1310 c   arguments
1311 c
1312       integer npara
1313       character*20 param(*),value(*)
1314
1315 C...Pythia parameters.
1316       INTEGER MSTP,MSTI
1317       DOUBLE PRECISION PARP,PARI
1318       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
1319
1320 C...User process initialization commonblock.
1321       INTEGER MAXPUP
1322       PARAMETER (MAXPUP=100)
1323       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
1324       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
1325       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
1326      &   IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
1327      &   LPRUP(MAXPUP)
1328
1329 C...User process event common block.
1330       INTEGER MAXNUP
1331       PARAMETER (MAXNUP=500)
1332       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
1333       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
1334       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
1335      &   ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
1336      &   VTIMUP(MAXNUP),SPINUP(MAXNUP)
1337
1338 C...Extra commonblock to transfer run info.
1339       INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
1340       COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
1341
1342 C...Inputs for the matching algorithm
1343       double precision etcjet,rclmax,etaclmax,qcut,clfact,showerkt
1344       integer maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres(30)
1345       integer nqmatch,nexcproc,iexcproc(MAXPUP),iexcval(MAXPUP)
1346       logical nosingrad,jetprocs
1347       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,showerkt,clfact,
1348      $   maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres,
1349      $   nqmatch,nexcproc,iexcproc,iexcval,nosingrad,jetprocs
1350
1351 c      DATA ktsche,maxjets,minjets,nexcres/0,-1,-1,0/
1352 c      DATA ktsche,nexcres/0,0/
1353 c      DATA qcut,clfact,showerkt/0d0,0d0,0d0/
1354
1355 C...Commonblock to transfer event-by-event matching info
1356       INTEGER NLJETS,IEXC,Ifile
1357       DOUBLE PRECISION PTCLUS
1358       COMMON/MEMAEV/PTCLUS(20),NLJETS,IEXC,Ifile
1359
1360 C...Local variables
1361       INTEGER I,MAXNJ,NREAD,MINJ,MAXJ
1362       parameter(MAXNJ=6)
1363       DOUBLE PRECISION XSTOT(MAXNJ),XSECTOT
1364       DOUBLE PRECISION ptjmin,etajmax,drjmin,ptbmin,etabmax,xqcut
1365
1366       integer icount
1367
1368 C...Functions
1369       INTEGER iexclusive
1370       EXTERNAL iexclusive
1371
1372 C...Initialize the icount counter to detect infinite loops
1373       icount=0
1374
1375 C...Need lower scale for final state radiation in e+e-
1376       IF(IABS(IDBMUP(1)).EQ.11.AND.IABS(IDBMUP(2)).EQ.11) then
1377         CALL PYGIVE('PARP(71)=1')
1378       ENDIF
1379
1380 C...CRUCIAL FOR JET-PARTON MATCHING: CALL UPVETO, ALLOW JET-PARTON MATCHING
1381 C      call pygive('MSTP(143)=1')
1382
1383 C
1384 C...Check jet multiplicities and set processes
1385 C
1386       DO I=1,MAXNJ
1387         XSTOT(I)=0D0
1388       ENDDO
1389       MINJ=MAXNJ
1390       MAXJ=0
1391       NREAD=0
1392       NUP=0
1393       DO WHILE(.true.)
1394 C         write(LNHOUT,*)'Launching MGEVNT'
1395         CALL MGEVNT()
1396                 write(LNHOUT,*)'NLJETS=',NLJETS
1397
1398         icount = icount+1
1399         if (icount.gt.10) then
1400           write (LNHOUT,*)
1401      &      'GeneratorInterface/PartonShowerVeto ME2phythia:'
1402      &      //' Aborting, loop in set_matching above ',icount,' cycles'
1403           write (LNHOUT,*) 'NUP = ',NUP,' IEXC = ',IEXC
1404           stop
1405         endif
1406
1407         IF(NUP.eq.0) goto 20
1408         IF(IEXC.EQ.-1) cycle
1409
1410         if(NLJETS.GT.MAXJ) MAXJ=NLJETS
1411         if(NLJETS.LT.MINJ) MINJ=NLJETS
1412 c        XSTOT(NLJETS+1)=XSTOT(NLJETS+1)+XWGTUP
1413         XSTOT(NLJETS+1)=XSTOT(NLJETS+1)+1
1414         NREAD=NREAD+1
1415       ENDDO
1416
1417  20   continue
1418
1419 C      REWIND(iunit)
1420
1421       write(LNHOUT,*) 'Minimum number of jets in file: ',MINJ
1422       write(LNHOUT,*) 'Maximum number of jets in file: ',MAXJ
1423
1424       XSECTOT=0d0
1425       DO I=1,NPRUP
1426          XSECTOT=XSECTOT+XSECUP(I)
1427       ENDDO
1428                 write(LNHOUT,*)'NPRUP=',NPRUP
1429       IF(NPRUP.eq.1.AND.MINJ.lt.MAXJ)THEN
1430 C...If different process ids not set by user, set by jet number
1431
1432          jetprocs=.true.
1433          IF(IEXCFILE.eq.0.AND.iexclusive(LPRUP(1)).ne.1) THEN
1434             nexcproc=1
1435             IEXCPROC(1)=MAXJ-MINJ
1436             IEXCVAL(1)=0
1437          ENDIF
1438          NPRUP=1+MAXJ-MINJ
1439          DO I=MINJ,MAXJ
1440             XSECUP(1+I-MINJ) = XSECTOT*XSTOT(I+1)/NREAD
1441             XMAXUP(1+I-MINJ) = XMAXUP(1)
1442             LPRUP(1+I-MINJ)  = I-MINJ
1443          ENDDO
1444       ELSE IF(IEXCFILE.EQ.0) THEN
1445 C...Check if any IEXCPROC set, then set IEXCFILE=1
1446          DO I=1,NPRUP
1447             IF(iexclusive(LPRUP(I)).EQ.0) IEXCFILE=1
1448          ENDDO
1449       ENDIF
1450
1451       WRITE(LNHOUT,*) ' Number of Events Read:: ',NREAD
1452       WRITE(LNHOUT,*) ' Total cross section (pb):: ',XSECTOT
1453       WRITE(LNHOUT,*) ' Process   Cross Section (pb):: '
1454       DO I=1,NPRUP
1455         WRITE(LNHOUT,'(I5,E23.5)') I,XSECUP(I)
1456       ENDDO
1457
1458       IF(MINJETS.EQ.-1) MINJETS=MINJ
1459       IF(MAXJETS.EQ.-1) MAXJETS=MAXJ
1460       write(LNHOUT,*) 'Minimum number of jets allowed: ',MINJETS
1461       write(LNHOUT,*) 'Maximum number of jets allowed: ',MAXJETS
1462       write(LNHOUT,*) 'IEXCFILE = ',IEXCFILE
1463       write(LNHOUT,*) 'jetprocs = ',jetprocs
1464       DO I=1,NPRUP
1465          write(LNHOUT,*) 'IEXCPROC(',LPRUP(I),') = ',
1466      $        iexclusive(LPRUP(I))
1467       ENDDO
1468
1469 C      CALL FLUSH()
1470
1471 C...Run PYPTFS instead of PYSHOW
1472 c        CALL PYGIVE("MSTJ(41)=12")
1473
1474 c***********************************************************************
1475 c   Read in jet cuts
1476 c***********************************************************************
1477
1478         call get_real   (npara,param,value," ptj " ,ptjmin,7d3)
1479         call get_real   (npara,param,value," etaj " ,etajmax,7d3)
1480         call get_real   (npara,param,value," ptb " ,ptbmin,7d3)
1481         call get_real   (npara,param,value," etab " ,etabmax,7d3)
1482         call get_real   (npara,param,value," drjj " ,drjmin,7d3)
1483         call get_real   (npara,param,value," xqcut " ,xqcut,0d0)
1484
1485         if(qcut.lt.xqcut) then
1486            if(showerkt.eq.1) then
1487               qcut=xqcut
1488            else
1489               qcut=max(xqcut*1.2,xqcut+5)
1490            endif
1491         endif
1492         if(xqcut.le.0)then
1493            write(*,*) 'Warning! ME generation QCUT = 0. QCUT set to 0!'
1494            qcut=0
1495         endif
1496
1497 c        etajmax=min(etajmax,etabmax)
1498 c        ptjmin=max(ptjmin,ptbmin)
1499
1500 c      IF(ICKKW.EQ.1) THEN
1501 c        WRITE(*,*) ' '
1502 c        WRITE(*,*) 'INPUT 0 FOR INCLUSIVE JET SAMPLE, 1 FOR EXCLUSIVE'
1503 c        WRITE(*,*) '(SELECT 0 FOR HIGHEST PARTON MULTIPLICITY SAMPLE)'
1504 c        WRITE(*,*) '(SELECT 1 OTHERWISE)'
1505 c        READ(*,*) IEXCFILE
1506 c      ENDIF
1507
1508 C     INPUT PARAMETERS FOR CONE ALGORITHM
1509
1510         IF(ETCJET.LE.PTJMIN)THEN
1511            ETCJET=MAX(PTJMIN+5,1.2*PTJMIN)
1512         ENDIF
1513
1514         RCLMAX=DRJMIN
1515         ETACLMAX=ETAJMAX
1516         IF(qcut.le.0)THEN
1517           WRITE(*,*) 'JET CONE PARAMETERS FOR MATCHING:'
1518           WRITE(*,*) 'ET>',ETCJET,' R=',RCLMAX
1519           WRITE(*,*) 'DR(PARTON-JET)<',1.5*RCLMAX
1520           WRITE(*,*) 'ETA(JET)<',ETACLMAX
1521       ELSE IF(ickkw.eq.1) THEN
1522         WRITE(*,*) 'KT JET PARAMETERS FOR MATCHING:'
1523         WRITE(*,*) 'QCUT=',qcut
1524         WRITE(*,*) 'ETA(JET)<',ETACLMAX
1525         WRITE(*,*) 'Note that in ME generation, qcut = ',xqcut
1526         write(*,*)'the showerkt param is ',showerkt
1527         if(showerkt.eq.1.0)THEN
1528 C             WRITE(*,*) 'shower kt is activated'
1529         endif
1530         if(showerkt.eq.1.0.and.MSTP(81).LT.20)THEN
1531           WRITE(*,*)'WARNING: "shower kt" needs pT-ordered showers'
1532           WRITE(*,*)'         Setting MSTP(81)=',20+MOD(MSTP(81),10)
1533           MSTP(81)=20+MOD(MSTP(81),10)
1534        endif
1535       else if(ickkw.eq.2)then
1536 c     Turn off color coherence suppressions (leave this to ME)
1537         CALL PYGIVE('MSTP(62)=2')
1538         CALL PYGIVE('MSTP(67)=0')
1539         if(MSTP(81).LT.20)THEN
1540           WRITE(*,*)'WARNING: Must run CKKW with pt-ordered showers'
1541           WRITE(*,*)'         Setting MSTP(81)=',20+MOD(MSTP(81),10)
1542           MSTP(81)=20+MOD(MSTP(81),10)
1543         endif
1544       endif
1545       return
1546       end
1547
1548       subroutine get_real(npara,param,value,name,var,def_value)
1549 c----------------------------------------------------------------------------------
1550 c   finds the parameter named "name" in param and associate to "value" in value
1551 c----------------------------------------------------------------------------------
1552       implicit none
1553
1554 c
1555 c   arguments
1556 c
1557       integer npara
1558       character*20 param(*),value(*)
1559       character*(*)  name
1560       real*8 var,def_value
1561 c
1562 c   local
1563 c
1564       logical found
1565       integer i
1566 c
1567 c   start
1568 c
1569 c      write(*,*)'entered get_real subroutine, looking for ',name,
1570 c     &' there are ',npara,' parameters'
1571       i=1
1572       found=.false.
1573       do while(.not.found.and.i.le.npara)
1574 c         write(*,*)'trying ',param(i)
1575         found = (index(param(i),name).ne.0)
1576         if (found) read(value(i),*) var
1577 c     if (found) write (*,*) name,var
1578         i=i+1
1579       enddo
1580       if (.not.found) then
1581         write (*,*) "Warning: parameter ",name," not found"
1582         write (*,*) "         setting it to default value ",def_value
1583         var=def_value
1584       else
1585         write(*,*) 'Found parameter ',name,var
1586       endif
1587       return
1588
1589       end
1590 c
1591
1592       subroutine get_integer(npara,param,value,name,var,def_value)
1593 c----------------------------------------------------------------------------------
1594 c   finds the parameter named "name" in param and associate to "value" in value
1595 c----------------------------------------------------------------------------------
1596       implicit none
1597 c
1598 c   arguments
1599 c
1600       integer npara
1601       character*20 param(*),value(*)
1602       character*(*)  name
1603       integer var,def_value
1604 c
1605 c   local
1606 c
1607       logical found
1608       integer i
1609 c
1610 c   start
1611 c
1612       i=1
1613       found=.false.
1614       do while(.not.found.and.i.le.npara)
1615         found = (index(param(i),name).ne.0)
1616         if (found) read(value(i),*) var
1617 c     if (found) write (*,*) name,var
1618         i=i+1
1619       enddo
1620       if (.not.found) then
1621         write (*,*) "Warning: parameter ",name," not found"
1622         write (*,*) "         setting it to default value ",def_value
1623         var=def_value
1624       else
1625         write(*,*)'Found parameter ',name,var
1626       endif
1627       return
1628
1629       end
1630
1631 C-----------------------------------------------------------------------
1632       SUBROUTINE ALPSOR(A,N,K,IOPT)
1633 C-----------------------------------------------------------------------
1634 C     Sort A(N) into ascending order
1635 C     IOPT = 1 : return sorted A and index array K
1636 C     IOPT = 2 : return index array K only
1637 C-----------------------------------------------------------------------
1638       DOUBLE PRECISION A(N),B(5000)
1639       INTEGER N,I,J,IOPT,K(N),IL(5000),IR(5000)
1640       IF (N.GT.5000) then
1641         write(*,*) 'Too many entries to sort in alpsrt, stop'
1642         stop
1643       endif
1644       if(n.le.0) return
1645       IL(1)=0
1646       IR(1)=0
1647       DO 10 I=2,N
1648       IL(I)=0
1649       IR(I)=0
1650       J=1
1651    2  IF(A(I).GT.A(J)) GOTO 5
1652       IF(IL(J).EQ.0) GOTO 4
1653       J=IL(J)
1654       GOTO 2
1655    4  IR(I)=-J
1656       IL(J)=I
1657       GOTO 10
1658    5  IF(IR(J).LE.0) GOTO 6
1659       J=IR(J)
1660       GOTO 2
1661    6  IR(I)=IR(J)
1662       IR(J)=I
1663   10  CONTINUE
1664       I=1
1665       J=1
1666       GOTO 8
1667   20  J=IL(J)
1668    8  IF(IL(J).GT.0) GOTO 20
1669    9  K(I)=J
1670       B(I)=A(J)
1671       I=I+1
1672       IF(IR(J)) 12,30,13
1673   13  J=IR(J)
1674       GOTO 8
1675   12  J=-IR(J)
1676       GOTO 9
1677   30  IF(IOPT.EQ.2) RETURN
1678       DO 31 I=1,N
1679   31  A(I)=B(I)
1680       END
1681
1682 C-----------------------------------------------------------------------
1683 C----Calorimeter simulation obtained from Frank Paige 23 March 1988-----
1684 C
1685 C          USE
1686 C
1687 C     CALL CALINIM
1688 C     CALL CALSIMM
1689 C
1690 C          THEN TO FIND JETS WITH A SIMPLIFIED VERSION OF THE UA1 JET
1691 C          ALGORITHM WITH JET RADIUS RJET AND MINIMUM SCALAR TRANSVERSE
1692 C          ENERGY EJCUT
1693 C            (RJET=1., EJCUT=5. FOR UA1)
1694 C          USE
1695 C
1696 C     CALL GETJETM(RJET,EJCUT)
1697 C
1698 C
1699 C-----------------------------------------------------------------------
1700 C
1701 C          ADDED BY MIKE SEYMOUR: PARTON-LEVEL CALORIMETER. ALL PARTONS
1702 C          ARE CONSIDERED TO BE HADRONS, SO IN FACT RESEM IS IGNORED
1703 C
1704 C     CALL CALPARM
1705 C
1706 C          HARD PARTICLE CALORIMETER. ONLY USES THOSE PARTICLES WHICH
1707 C          CAME FROM THE HARD PROCESS, AND NOT THE UNDERLYING EVENT
1708 C
1709 C     CALL CALHARM
1710 C
1711 C-----------------------------------------------------------------------
1712       SUBROUTINE CALINIM
1713 C
1714 C          INITIALIZE CALORIMETER FOR CALSIMM AND GETJETM.  NOTE THAT
1715 C          BECAUSE THE INITIALIZATION IS SEPARATE, CALSIMM CAN BE
1716 C          CALLED MORE THAN ONCE TO SIMULATE PILEUP OF SEVERAL EVENTS.
1717 C
1718       IMPLICIT NONE
1719 C...GETJET commonblocks
1720       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
1721       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
1722      &  SPHCAL,PCJET,ETJET
1723       PARAMETER (MNCY=200)
1724       PARAMETER (MNCPHI=200)
1725       COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
1726      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
1727      $YCMIN,YCMAX,NCY,NCPHI
1728       PARAMETER (NJMAX=500)
1729       COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
1730      $     NCJET
1731
1732       INTEGER IPHI,IY
1733       DOUBLE PRECISION PI,PHIX,YX,THX
1734       PARAMETER (PI=3.141593D0)
1735       LOGICAL FSTCAL
1736       DATA FSTCAL/.TRUE./
1737 C
1738 C          INITIALIZE ET ARRAY.
1739       DO 100 IPHI=1,NCPHI
1740       DO 100 IY=1,NCY
1741 100   ET(IY,IPHI)=0.
1742 C
1743       IF (FSTCAL) THEN
1744 C          CALCULATE TRIG. FUNCTIONS.
1745         DELPHI=2.*PI/FLOAT(NCPHI)
1746         DO 200 IPHI=1,NCPHI
1747         PHIX=DELPHI*(IPHI-.5)
1748         CPHCAL(IPHI)=COS(PHIX)
1749         SPHCAL(IPHI)=SIN(PHIX)
1750 200     CONTINUE
1751         DELY=(YCMAX-YCMIN)/FLOAT(NCY)
1752         DO 300 IY=1,NCY
1753         YX=DELY*(IY-.5)+YCMIN
1754         THX=2.*ATAN(EXP(-YX))
1755         CTHCAL(IY)=COS(THX)
1756         STHCAL(IY)=SIN(THX)
1757 300     CONTINUE
1758         FSTCAL=.FALSE.
1759       ENDIF
1760       END
1761 C
1762       SUBROUTINE CALSIMM
1763 C
1764 C          SIMPLE CALORIMETER SIMULATION.  ASSUME UNIFORM Y AND PHI
1765 C          BINS
1766 C...HEPEVT commonblock.
1767       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
1768       PARAMETER (NMXHEP=4000)
1769       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1770      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1771       DOUBLE PRECISION PHEP,VHEP
1772       SAVE /HEPEVT/
1773
1774 C...GETJET commonblocks
1775       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
1776       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
1777      &  SPHCAL,PCJET,ETJET
1778       PARAMETER (MNCY=200)
1779       PARAMETER (MNCPHI=200)
1780       COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
1781      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
1782      $YCMIN,YCMAX,NCY,NCPHI
1783       PARAMETER (NJMAX=500)
1784       COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
1785      $     NCJET
1786
1787       INTEGER IHEP,ID,IY,IPHI
1788       DOUBLE PRECISION PI,YIP,PSERAP,PHIIP,EIP
1789       PARAMETER (PI=3.141593D0)
1790 C
1791 C          FILL CALORIMETER
1792 C
1793       DO 200 IHEP=1,NHEP
1794       IF (ISTHEP(IHEP).EQ.1) THEN
1795         YIP=PSERAP(PHEP(1,IHEP))
1796         IF(YIP.LT.YCMIN.OR.YIP.GT.YCMAX) GOTO 200
1797         ID=ABS(IDHEP(IHEP))
1798 C---EXCLUDE TOP QUARK, LEPTONS, PROMPT PHOTONS
1799         IF ((ID.GE.11.AND.ID.LE.16).OR.ID.EQ.6.OR.ID.EQ.22) GOTO 200
1800 C
1801         PHIIP=ATAN2(PHEP(2,IHEP),PHEP(1,IHEP))
1802         IF(PHIIP.LT.0.) PHIIP=PHIIP+2.*PI
1803         IY=INT((YIP-YCMIN)/DELY)+1
1804         IPHI=INT(PHIIP/DELPHI)+1
1805         EIP=PHEP(4,IHEP)
1806 C            WEIGHT BY SIN(THETA)
1807         ET(IY,IPHI)=ET(IY,IPHI)+EIP*STHCAL(IY)
1808       ENDIF
1809   200 CONTINUE
1810       END
1811       SUBROUTINE GETJETM(RJET,EJCUT,ETAJCUT)
1812 C
1813 C          SIMPLE JET-FINDING ALGORITHM (SIMILAR TO UA1).
1814 C
1815 C     FIND HIGHEST REMAINING CELL > ETSTOP AND SUM SURROUNDING
1816 C          CELLS WITH--
1817 C            DELTA(Y)**2+DELTA(PHI)**2<RJET**2
1818 C            ET>ECCUT.
1819 C          KEEP JETS WITH ET>EJCUT AND ABS(ETA)<ETAJCUT
1820 C          THE UA1 PARAMETERS ARE RJET=1.0 AND EJCUT=5.0
1821 C
1822       IMPLICIT NONE
1823 C...GETJET commonblocks
1824       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
1825       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
1826      &  SPHCAL,PCJET,ETJET
1827       PARAMETER (MNCY=200)
1828       PARAMETER (MNCPHI=200)
1829       COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
1830      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
1831      $YCMIN,YCMAX,NCY,NCPHI
1832       PARAMETER (NJMAX=500)
1833       COMMON/GETCOMM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),
1834      $     NCJET
1835
1836       INTEGER IPHI,IY,J,K,NPHI1,NPHI2,NY1,
1837      &  NY2,IPASS,IYMX,IPHIMX,ITLIS,IPHI1,IPHIX,IY1,IYX
1838       DOUBLE PRECISION PI,RJET,
1839      &  ETMAX,ETSTOP,RR,ECCUT,PX,EJCUT
1840       PARAMETER (PI=3.141593D0)
1841       DOUBLE PRECISION ETAJCUT,PSERAP
1842 C
1843 C          PARAMETERS
1844       DATA ECCUT/0.1D0/
1845       DATA ETSTOP/1.5D0/
1846       DATA ITLIS/6/
1847 C
1848 C          INITIALIZE
1849 C
1850       DO 100 IPHI=1,NCPHI
1851       DO 100 IY=1,NCY
1852 100   JETNO(IY,IPHI)=0
1853       DO 110 J=1,NJMAX
1854       ETJET(J)=0.
1855       DO 110 K=1,4
1856 110   PCJET(K,J)=0.
1857       NCJET=0
1858       NPHI1=RJET/DELPHI
1859       NPHI2=2*NPHI1+1
1860       NY1=RJET/DELY
1861       NY2=2*NY1+1
1862       IPASS=0
1863 C-ap  initialize these two too to avoid compiler warnings
1864       iymx = 0
1865       iphimx = 0
1866 C-ap  end
1867 C
1868 C          FIND HIGHEST CELL REMAINING
1869 C
1870 1     ETMAX=0.
1871       DO 200 IPHI=1,NCPHI
1872       DO 210 IY=1,NCY
1873       IF(ET(IY,IPHI).LT.ETMAX) GOTO 210
1874       IF(JETNO(IY,IPHI).NE.0) GOTO 210
1875       ETMAX=ET(IY,IPHI)
1876       IYMX=IY
1877       IPHIMX=IPHI
1878 210   CONTINUE
1879 200   CONTINUE
1880       IF(ETMAX.LT.ETSTOP) RETURN
1881 C
1882 C          SUM CELLS
1883 C
1884       IPASS=IPASS+1
1885       IF(IPASS.GT.NCY*NCPHI) THEN
1886         WRITE(ITLIS,8888) IPASS
1887 8888    FORMAT(//' ERROR IN GETJETM...IPASS > ',I6)
1888         RETURN
1889       ENDIF
1890       NCJET=NCJET+1
1891       IF(NCJET.GT.NJMAX) THEN
1892         WRITE(ITLIS,9999) NCJET
1893 9999    FORMAT(//' ERROR IN GETJETM...NCJET > ',I5)
1894         RETURN
1895       ENDIF
1896       DO 300 IPHI1=1,NPHI2
1897       IPHIX=IPHIMX-NPHI1-1+IPHI1
1898       IF(IPHIX.LE.0) IPHIX=IPHIX+NCPHI
1899       IF(IPHIX.GT.NCPHI) IPHIX=IPHIX-NCPHI
1900       DO 310 IY1=1,NY2
1901       IYX=IYMX-NY1-1+IY1
1902       IF(IYX.LE.0) GOTO 310
1903       IF(IYX.GT.NCY) GOTO 310
1904       IF(JETNO(IYX,IPHIX).NE.0) GOTO 310
1905       RR=(DELY*(IY1-NY1-1))**2+(DELPHI*(IPHI1-NPHI1-1))**2
1906       IF(RR.GT.RJET**2) GOTO 310
1907       IF(ET(IYX,IPHIX).LT.ECCUT) GOTO 310
1908       PX=ET(IYX,IPHIX)/STHCAL(IYX)
1909 C          ADD CELL TO JET
1910       PCJET(1,NCJET)=PCJET(1,NCJET)+PX*STHCAL(IYX)*CPHCAL(IPHIX)
1911       PCJET(2,NCJET)=PCJET(2,NCJET)+PX*STHCAL(IYX)*SPHCAL(IPHIX)
1912       PCJET(3,NCJET)=PCJET(3,NCJET)+PX*CTHCAL(IYX)
1913       PCJET(4,NCJET)=PCJET(4,NCJET)+PX
1914       ETJET(NCJET)=ETJET(NCJET)+ET(IYX,IPHIX)
1915       JETNO(IYX,IPHIX)=NCJET
1916 310   CONTINUE
1917 300   CONTINUE
1918 C
1919 C          DISCARD JET IF ET < EJCUT.
1920 C
1921       IF(ETJET(NCJET).GT.EJCUT.AND.ABS(PSERAP(PCJET(1,NCJET))).LT
1922      $     .ETAJCUT) GOTO 1
1923       ETJET(NCJET)=0.
1924       DO 400 K=1,4
1925 400   PCJET(K,NCJET)=0.
1926       NCJET=NCJET-1
1927       GOTO 1
1928       END
1929 C-----------------------------------------------------------------------
1930       SUBROUTINE CALDELM(ISTLO,ISTHI)
1931 C     LABEL ALL PARTICLES WITH STATUS BETWEEN ISTLO AND ISTHI (UNTIL A
1932 C     PARTICLE WITH STATUS ISTOP IS FOUND) AS FINAL-STATE, CALL CALSIMM
1933 C     AND THEN PUT LABELS BACK TO NORMAL
1934 C-----------------------------------------------------------------------
1935       IMPLICIT NONE
1936       INTEGER MAXNUP
1937       PARAMETER(MAXNUP=500)
1938 C...HEPEVT commonblock.
1939       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
1940       PARAMETER (NMXHEP=4000)
1941       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1942      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1943       DOUBLE PRECISION PHEP,VHEP
1944       SAVE /HEPEVT/
1945       INTEGER ISTLO,ISTHI
1946
1947 C...Avoid gcc 4.3.4 warnings.
1948       ISTLO=ISTLO
1949       ISTHI=ISTHI
1950       CALL CALSIMM
1951       END
1952
1953 C****************************************************
1954 C iexclusive returns whether exclusive process or not
1955 C****************************************************
1956
1957       integer function iexclusive(iproc)
1958       implicit none
1959
1960       integer iproc, i
1961       INTEGER MAXPUP
1962       PARAMETER (MAXPUP=100)
1963
1964 C...Inputs for the matching algorithm
1965       double precision etcjet,rclmax,etaclmax,qcut,clfact,showerkt
1966       integer maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres(30)
1967       integer nqmatch,nexcproc,iexcproc(MAXPUP),iexcval(MAXPUP)
1968       logical nosingrad,jetprocs
1969       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,showerkt,clfact,
1970      $   maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres,
1971      $   nqmatch,nexcproc,iexcproc,iexcval,nosingrad,jetprocs
1972
1973
1974       iexclusive=-2
1975       do i=1,nexcproc
1976          if(iproc.eq.iexcproc(i)) then
1977             iexclusive=iexcval(i)
1978             return
1979          endif
1980       enddo
1981
1982       return
1983       end
1984
1985 C*********************************************************************
1986
1987       subroutine initpydata
1988
1989       INTEGER KCHG
1990       DOUBLE PRECISION PMAS,PARF,VCKM
1991       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
1992       INTEGER MSTP,MSTI
1993       DOUBLE PRECISION PARP,PARI
1994       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
1995
1996       if(PMAS(1,1).lt.0.1D0.or.PMAS(1,1).gt.1.) THEN
1997
1998         PMAS(1,1) = 0.33D0
1999         PMAS(2,1) = 0.33D0
2000         PMAS(3,1) = 0.5D0
2001         PMAS(4,1) = 1.5D0
2002         PMAS(5,1) = 4.8D0
2003         PMAS(6,1) = 175D0
2004         PMas(7,1) = 400D0
2005         PMas(8,1) = 400D0
2006         PMas(9,1) = 0D0
2007         PMas(10,1) = 0D0
2008         PMas(11,1) = 0.0005D0
2009         PMas(12,1) = 0D0
2010         PMas(13,1) = 0.10566D0
2011         PMas(14,1) = 0D0
2012         PMas(15,1) = 1.777D0
2013         PMas(16,1) = 0D0
2014         PMas(17,1) = 400D0
2015         PMas(18,1) = 0D0
2016         PMas(19,1) = 0D0
2017         PMas(20,1) = 0D0
2018         PMAS(21,1) = 0D0
2019
2020         MSTP(61) = 2
2021         MSTP(71) = 1
2022         MSTP(183)= 2013
2023
2024       endif
2025
2026       return
2027       end
2028
2029 C***********************************
2030 C Common block initialization block
2031 C***********************************
2032
2033       BLOCK DATA MEPYDAT
2034
2035       INTEGER MAXPUP
2036       PARAMETER (MAXPUP=100)
2037 C...Inputs for the matching algorithm
2038       double precision etcjet,rclmax,etaclmax,qcut,clfact,showerkt
2039       integer maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres(30)
2040       integer nqmatch,nexcproc,iexcproc(MAXPUP),iexcval(MAXPUP)
2041       logical nosingrad,jetprocs
2042       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,showerkt,clfact,
2043      $   maxjets,minjets,iexcfile,ktsche,mektsc,nexcres,excres,
2044      $   nqmatch,nexcproc,iexcproc,iexcval,nosingrad,jetprocs
2045
2046 C...GETJET commonblocks
2047       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
2048       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
2049      &  SPHCAL,PCJET,ETJET
2050       PARAMETER (MNCY=200)
2051       PARAMETER (MNCPHI=200)
2052       COMMON/CALORM/DELY,DELPHI,ET(MNCY,MNCPHI),
2053      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
2054      $YCMIN,YCMAX,NCY,NCPHI
2055
2056 C...Extra commonblock to transfer run info.
2057       INTEGER LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
2058       COMMON/UPPRIV/LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE
2059
2060 C...Initialization statements
2061       DATA showerkt/0.0/
2062       DATA qcut,clfact,etcjet/0d0,0d0,0d0/
2063       DATA ktsche,mektsc,maxjets,minjets,nexcres/0,1,-1,-1,0/
2064       DATA nqmatch/5/
2065       DATA nexcproc/0/
2066       DATA iexcproc/MAXPUP*-1/
2067       DATA iexcval/MAXPUP*-2/
2068 C      DATA nosingrad,showerkt,jetprocs/.false.,.false.,.false./
2069
2070       DATA NCY,NCPHI/50,60/
2071
2072       DATA LNHIN,LNHOUT,MSCAL,IEVNT,ICKKW,ISCALE/77,6,0,0,0,1/
2073
2074 C      nosingrad = .false.
2075 C      showerkt = .false.
2076 C      jetprocs = .false.
2077
2078       END
2079