AlpgenInterface/plugins/alpg_match.F

0001 *-- AUTHOR :    MICHELANGELO MANGANO
0002 C----------------------------------------------------------------------
0003       SUBROUTINE ALVETO(IPVETO)
0004 C----------------------------------------------------------------------
0005 C     SUBROUTINE TO IMPLEMENT THE MLM JET MATCHING CRITERION
0006 C----------------------------------------------------------------------
0007       IMPLICIT NONE
0008 C...GUP EVENT COMMON BLOCK
0009       INTEGER MAXNUP
0010       PARAMETER (MAXNUP=500)
0011       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0012       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0013       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
0014      &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
0015      &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
0016      &              SPINUP(MAXNUP)
0017 C...HEPEVT COMMONBLOCK.
0018       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0019       PARAMETER (NMXHEP=4000)
0020       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0021      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0022       DOUBLE PRECISION PHEP,VHEP
0023       SAVE /HEPEVT/
0024 c process code
0025       integer ihrd
0026       integer itopprc
0027 c--event file data
0028       integer Nunit,NunitOut,NunitIni
0029       character*512 filename,tmpstr
0030 c total number of input events
0031       common/AHio/Nunit,NunitOut,NunitIni,filename
0032 c inputs for the matching algorithm
0033       integer iexc,npfst,nplst,nljets,njstart,njlast
0034      $     ,ickkw
0035       double precision etclus,rclus,etaclmax
0036       common/AHopts/etclus,rclus,etaclmax,iexc,npfst
0037      $     ,nplst,nljets,njstart,njlast,ickkw
0038 c process and particles parameters
0039       double precision mc,mb,mt,mw,mz,mh
0040       double precision ebeam
0041       integer ndns,ih1,ih2
0042       integer nw,nz,nh,nph
0043       integer ihvy,ihvy2
0044 c pdf set type
0045       character pdftyp*25
0046 c total number of partons
0047       integer npart
0048       common/AHppara/mc,mb,mt,mw,mz,mh,
0049      & ebeam,ih1,ih2,
0050      & ihrd,itopprc,
0051      & nw,nz,nh,nph,
0052      & ihvy,ihvy2,
0053      & npart,ndns,pdftyp
0054 c weight information
0055       real *8 maxwgt,avgwgt,errwgt,totlum
0056       integer unwev
0057       common/AHwgts/maxwgt,avgwgt,errwgt,totlum,unwev
0058 c general parameters
0059       integer nparam
0060       parameter (nparam=200)
0061       integer parlen,partyp
0062       character chpar*8,chpdes*70
0063       double precision parval
0064       common/AHpars/parval(nparam),chpar(nparam),chpdes(nparam)
0065      $     ,parlen(nparam),partyp(nparam)
0066 c global event cuts
0067       double precision ptjmin,ptjmax,etajmax,drjmin,
0068      +     ptbmin,ptbmax,etabmax,drbmin,
0069      +     ptcmin,ptcmax,etacmax,drcmin,
0070      +     ptphmin,etaphmax,drphjmin,drphmin,drphlmin,
0071      +     ptlmin,etalmax,drlmin,metmin,mllmin,mllmax
0072       common/AHcuts/ptjmin,ptjmax,etajmax,drjmin,
0073      +     ptbmin,ptbmax,etabmax,drbmin,
0074      +        ptcmin,ptcmax,etacmax,drcmin,
0075      +        ptphmin,etaphmax,drphjmin,drphmin,drphlmin,
0076      +        ptlmin,etalmax,drlmin,metmin,mllmin,mllmax
0077       INTEGER IPVETO
0078 C     CALSIM AND JET VARIABLES
0079       INTEGER NCY,NCPHI,NJMAX,JETNO,NCJET
0080       DOUBLE PRECISION YCMIN,YCMAX,PI,ET,DELPHI,CPHCAL,SPHCAL,DELY,
0081      &     CTHCAL,STHCAL,PCJET,ETJET
0082       PARAMETER (NCY=100)
0083       PARAMETER (NCPHI=60,PI=3.141593D0)
0084       COMMON/CALOR/DELY,DELPHI,ET(NCY,NCPHI),
0085      $     CTHCAL(NCY),STHCAL(NCY),CPHCAL(NCPHI),SPHCAL(NCPHI),YCMIN
0086      $     ,YCMAX
0087       PARAMETER (NJMAX=500)
0088       COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(NCY,NCPHI),NCJET
0089 C
0090       DOUBLE PRECISION PSERAP
0091       INTEGER K(NJMAX),KP(NJMAX),KPJ(NJMAX)
0092 C LOCAL VARIABLES
0093       INTEGER I,J,IHEP,NMATCH,JRMIN
0094       DOUBLE PRECISION ETAJET,PHIJET,DELR,DPHI,DELRMIN
0095       DOUBLE PRECISION P(4,10),PT(10),ETA(10),PHI(10)
0096 C HEAVY QUARK MATCHING
0097       INTEGER IHVQ(10),NHVQ,NMJET,ID
0098       DOUBLE PRECISION ETAHVQ(10),PHIHVQ(10)
0099       INTEGER IEND,INORAD
0100       COMMON/SHVETO/IEND,INORAD(MAXNUP)
0101
0102 C DEBUGGING OPTIONS
0103       INTEGER IDBG
0104       PARAMETER (IDBG=0)
0105       DOUBLE PRECISION PTPART,PTJETS,ETAPART,ETAJETS
0106       INTEGER NMAX
0107       COMMON/MTCHDBG/PTPART(10),PTJETS(10),ETAPART(10),ETAJETS(10),NMAX
0108 C
0109       DOUBLE PRECISION ETMIN,ETMAX
0110       DOUBLE PRECISION TINY
0111       PARAMETER (TINY=1D-3)
0112       INTEGER ICOUNT
0113       DATA ICOUNT/0/
0114 C
0115       IPVETO=0
0116       IF(ICKKW.EQ.0) RETURN
0117       IF(IHRD.EQ.7.OR.IHRD.EQ.8.OR.IHRD.EQ.13) THEN
0118         WRITE(*,*) 'JET MATCHING FOR HARD PROCESS ',IHRD
0119      $       ,' NOT IMPLEMENTED, STOP'
0120         STOP
0121       ENDIF
0122       IF(NLJETS.EQ.0.AND.IEXC.EQ.0) RETURN
0123 C CHECK FOR EVENT ERROR OR ZERO WGT
0124       I=0
0125 C     HERWIG/PYTHIA SPECIFIC
0126 C     CALL ALSHER(I)
0127       IF(I.EQ.1) RETURN
0128 C
0129 C     INITIALIZE DEBUG IF NEEDED
0130       IF(IDBG.EQ.1) THEN
0131         WRITE(1,*) ' '
0132         WRITE(1,*) 'NEW EVENT '
0133         WRITE(1,*) 'PARTONS'
0134       ENDIF
0135       IF(IDBG.EQ.2) THEN
0136         DO I=1,10
0137           PTPART(I)=0D0
0138           ETAPART(I)=0D0
0139           PTJETS(I)=0D0
0140           ETAJETS(I)=0D0
0141         ENDDO
0142         NMAX=0
0143       ENDIF
0144 C
0145 C     RECONSTRUCT PARTON-LEVEL EVENT
0146 C     START FROM THE PARTONIC SYSTEM
0147 C     NOTICE THAT THIS INFO COMES FROM THE HEPEUP
0148       DO I=1,NLJETS
0149         IHEP=I+NJSTART
0150         DO J=1,4
0151           P(J,I)=PUP(J,IHEP)
0152         ENDDO
0153         PT(I)=SQRT(P(1,I)**2+P(2,I)**2)
0154         ETA(I)=-LOG(TAN(0.5D0*ATAN2(PT(I)+TINY,P(3,I))))
0155         PHI(I)=ATAN2(P(2,I),P(1,I))
0156         IF(IDBG.EQ.1) THEN
0157           WRITE(1,*) PT(I),ETA(I),PHI(I)
0158         ENDIF
0159       ENDDO
0160       IF(NLJETS.GT.0) CALL ALPSOR(PT,NLJETS,KP,2)
0161 C
0162 C     DISPLAY EVENT SEEN BY UPVETO:
0163 C     NOW THIS INFO COMES FROM THE HEPEVT
0164       IF(IDBG.EQ.1) THEN
0165         DO I=1,NHEP
0166           WRITE(1,111) I,ISTHEP(I),IDHEP(I),JMOHEP(1,I),JMOHEP(2,I)
0167      $         ,PHEP(1,I),PHEP(2,I),PHEP(3,I)
0168         ENDDO
0169  111    FORMAT(5(I4,1X),3(F12.5,1X))
0170       ENDIF
0171
0172 C     PRINT THE INORAD STATUS OF THE PARTICLES.
0173 C     (UNCOMMENT TO MAKE IT WORK)
0174 C      CALL DBINRD
0175
0176 C     DISPLAY PYTHIA EVENT:
0177 C      CALL PYLIST(7)    ! PYTHIA USER PROCESS EVENT DUMP
0178 C      CALL PYLIST(2)    ! PYTHIA FULL EVENT DUMP
0179 C      CALL PYLIST(5)    ! PYTHIA HEPEVT DUMP
0180 C
0181 C     RECONSTRUCT SHOWERED JETS:
0182       CALL CALINI
0183       CALL CALDEL(NPFST,NPLST,NJLAST)
0184       CALL GETJET(RCLUS,ETCLUS,ETACLMAX)
0185       IF(NCJET.GT.0) CALL ALPSOR(ETJET,NCJET,K,2)
0186       IF(IDBG.EQ.1) THEN
0187         WRITE(1,*) 'JETS'
0188         DO I=1,NCJET
0189           J=K(NCJET+1-I)
0190           ETAJET=PSERAP(PCJET(1,J))
0191           PHIJET=ATAN2(PCJET(2,J),PCJET(1,J))
0192           WRITE(1,*) ETJET(J),ETAJET,PHIJET
0193         ENDDO
0194       ENDIF
0195 C     ANALYSE ONLY EVENTS WITH AT LEAST NLJETS-RECONSTRUCTED JETS
0196       IF(NCJET.LT.NLJETS) GOTO 999
0197 C     ASSOCIATE PARTONS AND JETS, USING MIN(DELR) AS CRITERION
0198       NMATCH=0
0199       DO I=1,NCJET
0200         KPJ(I)=0
0201       ENDDO
0202       DO I=1,NLJETS
0203         DELRMIN=1D5
0204         DO 110 J=1,NCJET
0205           IF(KPJ(J).NE.0) GO TO 110
0206           ETAJET=PSERAP(PCJET(1,J))
0207           PHIJET=ATAN2(PCJET(2,J),PCJET(1,J))
0208           DPHI=ABS(PHI(KP(NLJETS-I+1))-PHIJET)
0209           IF(DPHI.GT.PI) DPHI=2.*PI-DPHI
0210           DELR=SQRT((ETA(KP(NLJETS-I+1))-ETAJET)**2+(DPHI)**2)
0211           IF(DELR.LT.DELRMIN) THEN
0212             DELRMIN=DELR
0213             JRMIN=J
0214           ENDIF
0215  110    CONTINUE
0216         ETMIN=1D10
0217         IF(DELRMIN.LT.1.5*RCLUS) THEN
0218           NMATCH=NMATCH+1
0219           KPJ(JRMIN)=I
0220           ETMIN=MIN(ETMIN,ETJET(JRMIN))
0221 C     ASSOCIATE PARTONS AND MATCHED JETS:
0222           IF(IDBG.EQ.2) THEN
0223             PTPART(I)=PT(KP(NLJETS-I+1))
0224             ETAPART(I)=ETA(KP(NLJETS-I+1))
0225             PTJETS(I)=ETJET(JRMIN)
0226             ETAJETS(I)=PSERAP(PCJET(1,JRMIN))
0227             NMAX=NCJET
0228           ENDIF
0229 C          WRITE(*,*) 'PARTON-JET',I,' BEST MATCH:',K(NCJET+1-JRMIN)
0230 C     $           ,DELRMIN
0231         ENDIF
0232       ENDDO
0233       IF(NMATCH.LT.NLJETS) GOTO 999
0234 C     REJECT EVENTS WITH LARGER JET MULTIPLICITY FROM EXCLUSIVE SAMPLE
0235       IF(NCJET.GT.NLJETS.AND.IEXC.EQ.1) GOTO 999
0236 C     VETO EVENTS WHERE MATCHED JETS ARE SOFTER THAN NON-MATCHED ONES
0237       IF(IEXC.NE.1) THEN
0238         J=NCJET
0239         DO I=1,NLJETS
0240           IF(KPJ(K(J)).EQ.0) GOTO 999
0241           J=J-1
0242         ENDDO
0243       ENDIF
0244 C
0245 C     ADDITIONAL TREATMENT FOR HVQ EVENTS: VETO/ACCEPT JETS EMITTED BY HVQ'S
0246       IF(IHRD.LE.2.OR.IHRD.EQ.6.OR.IHRD.EQ.10.OR.IHRD.EQ.15) THEN
0247         CALL CALINI
0248 C     RECOSTRUCT POSSIBLE JETS FROM RADIATION OFF THE TOP QUARKS
0249         CALL CALDEL_HVQ(NPFST,NPLST,NJLAST)
0250         CALL GETJET(RCLUS,ETCLUS,ETACLMAX)
0251 C     IF NO EXTRA JET: ACCEPT EVENT
0252         IF(NCJET.EQ.0) RETURN
0253 C     IF EXTRA JETS, REMOVE THOSE LYING WITHIN DRJMIN OF A B/C QUARK, TO
0254 C     ALLOW THE SHOWER TO GOVERN THE DEVELOPMENT OF A B/C JET.
0255 C     START BY FLAGGING VETOED Q AND QBAR OBJECTS:
0256         NHVQ=0
0257         DO I=1,NHEP
0258           ID=IDHEP(I)
0259           IF(INORAD(I).EQ.1.AND.ABS(ID).LE.5.AND.ABS(ID)
0260      $         .GE.4) THEN
0261             NHVQ=NHVQ+1
0262             IHVQ(NHVQ)=I
0263             ETAHVQ(NHVQ)=PSERAP(PHEP(1,I))
0264             PHIHVQ(NHVQ)=ATAN2(PHEP(2,I),PHEP(1,I))
0265           ENDIF
0266         ENDDO
0267         NMJET=NCJET
0268         DO I=1,NCJET
0269           ETAJET=PSERAP(PCJET(1,I))
0270           PHIJET=ATAN2(PCJET(2,I),PCJET(1,I))
0271           DO J=1,NHVQ
0272             DPHI=ABS(PHIHVQ(J)-PHIJET)
0273             IF(DPHI.GT.PI) DPHI=ABS(DPHI-2*PI)
0274             DELR=SQRT(DPHI**2+(ETAJET-ETAHVQ(J))**2)
0275             IF(DELR.LT.DRJMIN) THEN
0276               NMJET=NMJET-1
0277               ETJET(I)=0D0
0278             ENDIF
0279           ENDDO
0280         ENDDO
0281 C     IF NO UNMATCHED JET: ACCEPT EVENT
0282         IF(NMJET.EQ.0) RETURN
0283 C     IF JETS AND IEXC=1: REJECT EVENT
0284         IF(IEXC.EQ.1) GOTO 999
0285 C     IF JETS AND IEXC=0: CHECK THAT JETS ARE SOFTER THAN MATCHED ONES
0286         ETMAX=0D0
0287         DO I=1,NCJET
0288           ETMAX=MAX(ETMAX,ETJET(I))
0289         ENDDO
0290         IF(ETMAX.GT.ETMIN) GOTO 999
0291       ENDIF
0292       RETURN
0293 C     HERWIG/PYTHIA TERMINATION:
0294  999  CALL ALSHEN
0295       IPVETO=1
0296       END
0297
0298
0299 C----------------------------------------------------------------------
0300       SUBROUTINE ALSETP
0301 C----------------------------------------------------------------------
0302 C     EVENT INITIALITION ROUTINE THAT STARTS OFF FILLED COMMON BLOCKS
0303 C----------------------------------------------------------------------
0304       IMPLICIT NONE
0305 c process code
0306       integer ihrd
0307       integer itopprc
0308 c inputs for the matching algorithm
0309       integer iexc,npfst,nplst,nljets,njstart,njlast
0310      $     ,ickkw
0311       double precision etclus,rclus,etaclmax
0312       common/AHopts/etclus,rclus,etaclmax,iexc,npfst
0313      $     ,nplst,nljets,njstart,njlast,ickkw
0314 c process and particles parameters
0315       double precision mc,mb,mt,mw,mz,mh
0316       double precision ebeam
0317       integer ndns,ih1,ih2
0318       integer nw,nz,nh,nph
0319       integer ihvy,ihvy2
0320 c pdf set type
0321       character pdftyp*25
0322 c total number of partons
0323       integer npart
0324       common/AHppara/mc,mb,mt,mw,mz,mh,
0325      & ebeam,ih1,ih2,
0326      & ihrd,itopprc,
0327      & nw,nz,nh,nph,
0328      & ihvy,ihvy2,
0329      & npart,ndns,pdftyp
0330 c general parameters
0331       integer nparam
0332       parameter (nparam=200)
0333       integer parlen,partyp
0334       character chpar*8,chpdes*70
0335       double precision parval
0336       common/AHpars/parval(nparam),chpar(nparam),chpdes(nparam)
0337      $     ,parlen(nparam),partyp(nparam)
0338 c global event cuts
0339       double precision ptjmin,ptjmax,etajmax,drjmin,
0340      +     ptbmin,ptbmax,etabmax,drbmin,
0341      +     ptcmin,ptcmax,etacmax,drcmin,
0342      +     ptphmin,etaphmax,drphjmin,drphmin,drphlmin,
0343      +     ptlmin,etalmax,drlmin,metmin,mllmin,mllmax
0344       common/AHcuts/ptjmin,ptjmax,etajmax,drjmin,
0345      +     ptbmin,ptbmax,etabmax,drbmin,
0346      +        ptcmin,ptcmax,etacmax,drcmin,
0347      +        ptphmin,etaphmax,drphjmin,drphmin,drphlmin,
0348      +        ptlmin,etalmax,drlmin,metmin,mllmin,mllmax
0349       CHARACTER *3 CSHO
0350 C     CALSIM AND JET VARIABLES
0351       INTEGER NCY,NCPHI,NJMAX,JETNO,NCJET
0352       DOUBLE PRECISION YCMIN,YCMAX,PI,ET,DELPHI,CPHCAL,SPHCAL,DELY,
0353      &     CTHCAL,STHCAL,PCJET,ETJET
0354       PARAMETER (NCY=100)
0355       PARAMETER (NCPHI=60,PI=3.141593D0)
0356       COMMON/CALOR/DELY,DELPHI,ET(NCY,NCPHI),
0357      $     CTHCAL(NCY),STHCAL(NCY),CPHCAL(NCPHI),SPHCAL(NCPHI),YCMIN
0358      $     ,YCMAX
0359 C
0360       INTEGER MAXNUP
0361       PARAMETER (MAXNUP=500)
0362       INTEGER IEND,INORAD
0363       COMMON/SHVETO/IEND,INORAD(MAXNUP)
0364 c
0365 C     LOCAL VARIABLES
0366       CHARACTER*70 STDUMMY
0367       INTEGER I,NTMP,IHEPMIN
0368       DOUBLE PRECISION TMP
0369 C     USER ACCESS TO MATCHING PARAMETERS
0370       INTEGER IUSRMAT
0371       PARAMETER (IUSRMAT=1)
0372 C      character*512                           CXpar
0373 C      integer         Ixpar
0374 C      double precision            RXpar
0375 C      common /EXGPAR/ IXpar(100), RXpar(100), CXpar(100)
0376
0377 C     WRITE PARAMETER VALUES
0378       CALL AHSPAR
0379       IH1=1
0380 C     CONVERT PDF TYPES
0381       CALL PDFCONVH(NDNS,NTMP,PDFTYP)
0382 C     DEFINE RANGE FOR PARTONS TO BE USED IN MATCHING
0383       NLJETS=PARVAL(10)     !  NJETS
0384       DO I=1,MAXNUP
0385         INORAD(I)=0
0386       ENDDO
0387       CALL ALSHCD(CSHO)
0388       IF(CSHO.EQ.'HER') THEN
0389         NPFST=149
0390         NPLST=149
0391 C     HERWIG: ALL SHOWERS ORIGINATE FROM IHEP=6
0392         IEND=6
0393 C     HERWIG: HEPEVT EVENT RECORD FOR FINAL STATE STARTS AT 7=6+1
0394         IHEPMIN=6
0395       ELSE
0396         NPFST=1
0397         NPLST=1
0398 C     PYTHIA: ALL SHOWERS ORIGINATE FROM IHEP=O
0399         IEND=0
0400 C     PYTHIA: HEPEVT EVENT RECORD FOR FINAL STATE STARTS AT 7=6+1
0401         IHEPMIN=6
0402       ENDIF
0403       IF(IHRD.LE.2) THEN
0404 C        NLJETS=NPART-6
0405         NJSTART=4
0406         NJLAST=155
0407 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM THE W
0408         INORAD(IHEPMIN+2+NLJETS+1)=1
0409 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM HEAVY QUARK
0410 C     PAIR
0411         INORAD(IHEPMIN+1)=1
0412         INORAD(IHEPMIN+2)=1
0413       ELSEIF(IHRD.LE.4) THEN
0414 C        NLJETS=NPART-4
0415         NJSTART=2
0416         NJLAST=155
0417 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM THE W
0418         IF(CSHO.EQ.'HER'.AND.NLJETS.EQ.0) THEN
0419 C     ALLOW RADIATION FROM IHEPMIN+1 IN HERWIG WHEN NJET=0
0420           INORAD(IHEPMIN+NLJETS+1)=0
0421         ELSE
0422           INORAD(IHEPMIN+NLJETS+1)=1
0423         ENDIF
0424       ELSEIF(IHRD.EQ.5) THEN
0425 C        NLJETS=NPART-3*(NW+NZ)-NH-2
0426         NJSTART=2
0427         NJLAST=155
0428 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM THE W
0429         IF(CSHO.EQ.'HER'.AND.NLJETS.EQ.0.AND.(NW+NZ+NH+NPH).EQ.1) THEN
0430 C     ALLOW RADIATION FROM IHEPMIN+1 IN HERWIG WHEN NJET=0 AND NW+NZ+NH
0431 C     +NPH=1
0432           INORAD(IHEPMIN+NLJETS+1)=0
0433         ELSE
0434           DO I=1,NW+NZ+NH+NPH
0435             INORAD(IHEPMIN+NLJETS+I)=1
0436           ENDDO
0437         ENDIF
0438       ELSEIF(IHRD.EQ.6) THEN
0439 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM HEAVY QUARK
0440 C     PAIR
0441 C     NLJETS=NPART-8  (IHVY.EQ.6)     NLJETS=NPART-4  (IHVY.LT.6)
0442         INORAD(IHEPMIN+1)=1
0443         INORAD(IHEPMIN+2)=1
0444         NJSTART=4
0445         NJLAST=155
0446       ELSEIF(IHRD.EQ.9) THEN
0447 C     NLJETS=NPART-2
0448         NJSTART=2
0449         NJLAST=155
0450       ELSEIF(IHRD.EQ.10) THEN
0451 C        NLJETS=NPART-4
0452         NJSTART=3
0453         NJLAST=155
0454 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM THE W
0455         INORAD(IHEPMIN+NLJETS+1+1)=1
0456 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM THE CHARM
0457         INORAD(IHEPMIN+1)=1
0458       ELSEIF(IHRD.EQ.11) THEN
0459 C        NLJETS=NPART-2-NPH
0460         NJSTART=2
0461         NJLAST=155
0462 C     DO NOT INCLUDE IN MATCHING THE HARD PHOTONS
0463         DO I=1,NPH
0464           INORAD(IHEPMIN+NLJETS+I)=1
0465         ENDDO
0466       ELSEIF(IHRD.EQ.12) THEN
0467 C        NLJETS=NPART-2-NH
0468         NJSTART=2
0469         NJLAST=155
0470 C     DO NOT INCLUDE IN MATCHING THE HIGGS DECAY PRODUCTS
0471         IF(NLJETS+NH.GT.1) THEN
0472           DO I=1,NH
0473             INORAD(IHEPMIN+NLJETS+I)=1
0474           ENDDO
0475         ELSE
0476           IF(CSHO.EQ.'HER'.AND.NLJETS.EQ.0) THEN
0477 C     ALLOW RADIATION FROM IHEPMIN+1 IN HERWIG WHEN NJET=0
0478             INORAD(IHEPMIN+1)=0
0479           ELSE
0480             INORAD(IHEPMIN+NLJETS+1)=1
0481           ENDIF
0482         ENDIF
0483       ELSEIF(IHRD.EQ.14) THEN
0484 C        NLJETS=NPART-4-NPH
0485         NJSTART=2
0486         NJLAST=155
0487 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM THE W
0488         IF(CSHO.EQ.'HER'.AND.NLJETS.EQ.0.AND.NPH.EQ.0) THEN
0489 C     ALLOW RADIATION FROM IHEPMIN+1 IN HERWIG WHEN NJET=0
0490           INORAD(IHEPMIN+NLJETS+1)=0
0491         ELSE
0492           INORAD(IHEPMIN+NLJETS+1)=1
0493         ENDIF
0494       ELSEIF(IHRD.EQ.15) THEN
0495 C        NLJETS=NPART-6-NPH
0496         NJSTART=2
0497         NJLAST=155
0498 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM THE W
0499         INORAD(IHEPMIN+2+NLJETS+NPH+1)=1
0500 C     DO NOT INCLUDE IN MATCHING RADIATION ORIGINATING FROM HEAVY QUARK
0501 C     PAIR
0502 C     NOTICE THAT HERE THE LIGHT JETS PRECEDE THE QQ PAIR
0503         INORAD(IHEPMIN+NLJETS+1)=1
0504         INORAD(IHEPMIN+NLJETS+2)=1
0505       ENDIF
0506 C
0507 C     INPUT JET MATCHING CRITERIA
0508 C
0509       IF(ICKKW.EQ.1) THEN
0510 C       WRITE(*,*) ' '
0511 C       WRITE(*,*) 'INPUT 0 FOR INCLUSIVE JET SAMPLE, 1 FOR EXCLUSIVE'
0512 C       WRITE(*,*) '(SELECT 0 FOR HIGHEST PARTON MULTIPLICITY SAMPLE)'
0513 C       WRITE(*,*) '(SELECT 1 OTHERWISE)'
0514 C       READ(*,*) IEXC
0515 C       IEXC = IXpar(2)
0516         WRITE(*,*) ' IEXC = ',IEXC
0517         IF(NLJETS.GT.0) THEN
0518           WRITE(*,*) 'INPUT ET(CLUS), R(CLUS)'
0519 C         ETCLUS=MAX(PTJMIN+5,1.2*PTJMIN)
0520 C         RCLUS=DRJMIN
0521 C         WRITE(*,*) '(SUGGESTED VALUES:',ETCLUS,RCLUS,')'
0522 C         READ(*,*) ETCLUS,RCLUS,ETACLMAX
0523 C         ETCLUS = RXpar(1)
0524 C         RCLUS  = RXpar(2)
0525           ETACLMAX = ETAJMAX
0526         ELSEIF(NLJETS.EQ.0) THEN
0527           WRITE(*,*) 'INPUT ET(CLUS), R(CLUS)'
0528           WRITE(*,*) '(MUST MATCH VALUES USED IN PROCESSING',
0529      +         ' OF NJET>0 EVENTS; THESE DEFAULT TO:'
0530           WRITE(*,*)
0531      +         'ETCLUS=MAX(PTJMIN+5,1.2*PTJMIN) RCLUS=DRJMIN, ',
0532      +         ' ETACLMAX=ETAJMAX)'
0533 C         READ(*,*) ETCLUS,RCLUS,ETACLMAX
0534 C         ETCLUS = RXpar(1)
0535 C         RCLUS  = RXpar(2)
0536           ETACLMAX = ETAJMAX
0537         ENDIF
0538         WRITE(*,*) ' '
0539         WRITE(*,*) 'JET PARAMETERS FOR MATCHING:'
0540         WRITE(*,*) 'ET >',ETCLUS,' ETACLUS <',ETACLMAX,' R =',RCLUS
0541         WRITE(*,*) 'DR(PARTON-JET) <',1.5*RCLUS
0542       ENDIF
0543 C
0544 C     CALORIMETER ETA RANGE
0545 C     YCMAX=ETAJMAX+DRJMIN
0546       YCMAX=ETACLMAX+RCLUS      ! MLM August 24
0547       YCMIN=-YCMAX
0548       END
0549
0550 c-------------------------------------------------------------------
0551       subroutine AHspar
0552 c     set list of parameters types and assign default values
0553 c-------------------------------------------------------------------
0554       implicit none
0555 c process code
0556       integer ihrd
0557       integer itopprc
0558 c--event file data
0559       integer Nunit,NunitOut,NunitIni
0560       character*512 filename,tmpstr
0561 c total number of input events
0562       common/AHio/Nunit,NunitOut,NunitIni,filename
0563 c inputs for the matching algorithm
0564       integer iexc,npfst,nplst,nljets,njstart,njlast
0565      $     ,ickkw
0566       double precision etclus,rclus,etaclmax
0567       common/AHopts/etclus,rclus,etaclmax,iexc,npfst
0568      $     ,nplst,nljets,njstart,njlast,ickkw
0569 c process and particles parameters
0570       double precision mc,mb,mt,mw,mz,mh
0571       double precision ebeam
0572       integer ndns,ih1,ih2
0573       integer nw,nz,nh,nph
0574       integer ihvy,ihvy2
0575 c pdf set type
0576       character pdftyp*25
0577 c total number of partons
0578       integer npart
0579       common/AHppara/mc,mb,mt,mw,mz,mh,
0580      & ebeam,ih1,ih2,
0581      & ihrd,itopprc,
0582      & nw,nz,nh,nph,
0583      & ihvy,ihvy2,
0584      & npart,ndns,pdftyp
0585 c weight information
0586       real *8 maxwgt,avgwgt,errwgt,totlum
0587       integer unwev
0588       common/AHwgts/maxwgt,avgwgt,errwgt,totlum,unwev
0589 c general parameters
0590       integer nparam
0591       parameter (nparam=200)
0592       integer parlen,partyp
0593       character chpar*8,chpdes*70
0594       double precision parval
0595       common/AHpars/parval(nparam),chpar(nparam),chpdes(nparam)
0596      $     ,parlen(nparam),partyp(nparam)
0597 c global event cuts
0598       double precision ptjmin,ptjmax,etajmax,drjmin,
0599      +     ptbmin,ptbmax,etabmax,drbmin,
0600      +     ptcmin,ptcmax,etacmax,drcmin,
0601      +     ptphmin,etaphmax,drphjmin,drphmin,drphlmin,
0602      +     ptlmin,etalmax,drlmin,metmin,mllmin,mllmax
0603       common/AHcuts/ptjmin,ptjmax,etajmax,drjmin,
0604      +     ptbmin,ptbmax,etabmax,drbmin,
0605      +        ptcmin,ptcmax,etacmax,drcmin,
0606      +        ptphmin,etaphmax,drphjmin,drphmin,drphlmin,
0607      +        ptlmin,etalmax,drlmin,metmin,mllmin,mllmax
0608 c
0609       ih2=parval(2)
0610       ebeam=parval(3)
0611       ndns=parval(4)
0612       ickkw=parval(7)
0613       ihvy=parval(11)
0614       ihvy2=parval(12)
0615       nw=parval(13)
0616       nz=parval(14)
0617       nh=parval(15)
0618       nph=parval(16)
0619       ptjmin=parval(30)
0620       ptbmin=parval(31)
0621       ptcmin=parval(32)
0622       ptlmin=parval(33)
0623       metmin=parval(34)
0624       ptphmin=parval(35)
0625       etajmax=parval(40)
0626       etabmax=parval(41)
0627       etacmax=parval(42)
0628       etalmax=parval(43)
0629       etaphmax=parval(44)
0630       drjmin=parval(50)
0631       drbmin=parval(51)
0632       drcmin=parval(52)
0633       drlmin=parval(55)
0634       drphjmin=parval(56)
0635       drphlmin=parval(57)
0636       drphmin=parval(58)
0637       mllmin=parval(61)
0638       mllmax=parval(62)
0639       itopprc=parval(102)
0640 c
0641       end
0642
0643 c-------------------------------------------------------------------
0644       subroutine pdfconvH(nin,nout,type)
0645 c-------------------------------------------------------------------
0646 c converts ALPHA convention for PDF namings to hvqpdf conventions
0647       implicit none
0648       integer nin,nout
0649       character*25 type
0650       character*25 pdftyp(20,2)
0651       data pdftyp/
0652 c cteq sets
0653      $     'CTEQ4M ','CTEQ4L ','CTEQ4HJ',
0654      $     'CTEQ5M ','CTEQ5L ','CTEQ5HJ',
0655      $     'CTEQ6M ','CTEQ6L ',12*' ',
0656 C MRST SETS
0657      $     'MRST99 ',
0658      $     'MRST01; as=0.119','MRST01; as=0.117','MRST01; as=0.121'
0659      $     ,'MRST01J; as=0.121','MRST02LO',14*' '/
0660       integer pdfmap(20,2)
0661       data pdfmap/
0662      $   81,83,88,   101,103, 104,   131,133, 12*0,
0663      $  111,  185,186,187,188,189,   14*0/
0664 c
0665       nout=pdfmap(mod(nin ,100),1+nin /100)
0666       type=pdftyp(mod(nin ,100),1+nin /100)
0667
0668       end
0669
0670 C-----------------------------------------------------------------------
0671 C----Calorimeter simulation obtained from Frank Paige 23 March 1988-----
0672 C
0673 C          USE
0674 C
0675 C     CALL CALINI
0676 C     CALL CALSIM
0677 C
0678 C          THEN TO FIND JETS WITH A SIMPLIFIED VERSION OF THE UA1 JET
0679 C          ALGORITHM WITH JET RADIUS RJET AND MINIMUM SCALAR TRANSVERSE
0680 C          ENERGY EJCUT
0681 C            (RJET=1., EJCUT=5. FOR UA1)
0682 C          USE
0683 C
0684 C     CALL GETJET(RJET,EJCUT)
0685 C
0686 C
0687 C-----------------------------------------------------------------------
0688 C
0689 C          ADDED BY MIKE SEYMOUR: PARTON-LEVEL CALORIMETER. ALL PARTONS
0690 C          ARE CONSIDERED TO BE HADRONS, SO IN FACT RESEM IS IGNORED
0691 C
0692 C     CALL CALPAR
0693 C
0694 C          HARD PARTICLE CALORIMETER. ONLY USES THOSE PARTICLES WHICH
0695 C          CAME FROM THE HARD PROCESS, AND NOT THE UNDERLYING EVENT
0696 C
0697 C     CALL CALHAR
0698 C
0699 C-----------------------------------------------------------------------
0700       SUBROUTINE CALINI
0701 C
0702 C          INITIALIZE CALORIMETER FOR CALSIM AND GETJET.  NOTE THAT
0703 C          BECAUSE THE INITIALIZATION IS SEPARATE, CALSIM CAN BE
0704 C          CALLED MORE THAN ONCE TO SIMULATE PILEUP OF SEVERAL EVENTS.
0705 C
0706       IMPLICIT NONE
0707       INTEGER NCY,NCPHI,NJMAX,IPHI,IY,JETNO,NCJET
0708       DOUBLE PRECISION YCMIN,YCMAX,PI,ET,DELPHI,PHIX,CPHCAL,SPHCAL,DELY,
0709      &  YX,THX,CTHCAL,STHCAL,PCJET,ETJET
0710       PARAMETER (NCY=100)
0711       PARAMETER (NCPHI=60,PI=3.141593D0)
0712       COMMON/CALOR/DELY,DELPHI,ET(NCY,NCPHI),
0713      $CTHCAL(NCY),STHCAL(NCY),CPHCAL(NCPHI),SPHCAL(NCPHI),YCMIN,YCMAX
0714       PARAMETER (NJMAX=500)
0715       COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(NCY,NCPHI),NCJET
0716       LOGICAL FSTCAL
0717       DATA FSTCAL/.TRUE./
0718 C
0719 C          INITIALIZE ET ARRAY.
0720       DO 100 IPHI=1,NCPHI
0721       DO 100 IY=1,NCY
0722 100   ET(IY,IPHI)=0.
0723 C
0724       IF (FSTCAL) THEN
0725 C          CALCULATE TRIG. FUNCTIONS.
0726         DELPHI=2.*PI/FLOAT(NCPHI)
0727         DO 200 IPHI=1,NCPHI
0728         PHIX=DELPHI*(IPHI-.5)
0729         CPHCAL(IPHI)=COS(PHIX)
0730         SPHCAL(IPHI)=SIN(PHIX)
0731 200     CONTINUE
0732         DELY=(YCMAX-YCMIN)/FLOAT(NCY)
0733         DO 300 IY=1,NCY
0734         YX=DELY*(IY-.5)+YCMIN
0735         THX=2.*ATAN(EXP(-YX))
0736         CTHCAL(IY)=COS(THX)
0737         STHCAL(IY)=SIN(THX)
0738 300     CONTINUE
0739         FSTCAL=.FALSE.
0740       ENDIF
0741       END
0742 C
0743       SUBROUTINE CALSIM
0744 C
0745 C          SIMPLE CALORIMETER SIMULATION.  ASSUME UNIFORM Y AND PHI
0746 C          BINS
0747 C...HEPEVT commonblock.
0748       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0749       PARAMETER (NMXHEP=4000)
0750       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0751      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0752       DOUBLE PRECISION PHEP,VHEP
0753       SAVE /HEPEVT/
0754       INTEGER NCY,NCPHI,NJMAX,IHEP,ID,IY,IPHI,JETNO,NCJET
0755       DOUBLE PRECISION YCMIN,YCMAX,PI,YIP,PSERAP,
0756      &  PHIIP,DELY,DELPHI,EIP,ET,STHCAL,CTHCAL,CPHCAL,SPHCAL,
0757      &  PCJET,ETJET
0758       PARAMETER (NCY=100)
0759       PARAMETER (NCPHI=60,PI=3.141593D0)
0760       COMMON/CALOR/DELY,DELPHI,ET(NCY,NCPHI),
0761      $CTHCAL(NCY),STHCAL(NCY),CPHCAL(NCPHI),SPHCAL(NCPHI),YCMIN,YCMAX
0762       PARAMETER (NJMAX=500)
0763       COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(NCY,NCPHI),NCJET
0764 C
0765 C          FILL CALORIMETER
0766 C
0767       DO 200 IHEP=1,NHEP
0768       IF (ISTHEP(IHEP).EQ.1) THEN
0769         YIP=PSERAP(PHEP(1,IHEP))
0770         IF(YIP.LT.YCMIN.OR.YIP.GT.YCMAX) GOTO 200
0771         ID=ABS(IDHEP(IHEP))
0772 C---EXCLUDE TOP QUARK, LEPTONS, PROMPT PHOTONS
0773         IF ((ID.GE.11.AND.ID.LE.16).OR.ID.EQ.6.OR.ID.EQ.22) GOTO 200
0774 C
0775         PHIIP=ATAN2(PHEP(2,IHEP),PHEP(1,IHEP))
0776         IF(PHIIP.LT.0.) PHIIP=PHIIP+2.*PI
0777         IY=INT((YIP-YCMIN)/DELY)+1
0778         IPHI=INT(PHIIP/DELPHI)+1
0779         EIP=PHEP(4,IHEP)
0780 C            WEIGHT BY SIN(THETA)
0781         ET(IY,IPHI)=ET(IY,IPHI)+EIP*STHCAL(IY)
0782       ENDIF
0783   200 CONTINUE
0784   999 END
0785       SUBROUTINE GETJET(RJET,EJCUT,ETAJCUT)
0786 C
0787 C          SIMPLE JET-FINDING ALGORITHM (SIMILAR TO UA1).
0788 C
0789 C     FIND HIGHEST REMAINING CELL > ETSTOP AND SUM SURROUNDING
0790 C          CELLS WITH--
0791 C            DELTA(Y)**2+DELTA(PHI)**2<RJET**2
0792 C            ET>ECCUT.
0793 C          KEEP JETS WITH ET>EJCUT AND ABS(ETA)<ETAJCUT
0794 C          THE UA1 PARAMETERS ARE RJET=1.0 AND EJCUT=5.0
0795 C
0796       IMPLICIT NONE
0797       INTEGER NCY,NCPHI,NJMAX,IPHI,IY,JETNO,J,K,NCJET,NPHI1,NPHI2,NY1,
0798      &  NY2,IPASS,IYMX,IPHIMX,ITLIS,IPHI1,IPHIX,IY1,IYX
0799       DOUBLE PRECISION YCMIN,YCMAX,PI,ETJET,PCJET,RJET,DELPHI,DELY,
0800      &  ETMAX,ET,ETSTOP,RR,ECCUT,PX,STHCAL,CPHCAL,SPHCAL,CTHCAL,EJCUT
0801       PARAMETER (NCY=100)
0802       PARAMETER (NCPHI=60,PI=3.141593D0)
0803       COMMON/CALOR/DELY,DELPHI,ET(NCY,NCPHI),
0804      &CTHCAL(NCY),STHCAL(NCY),CPHCAL(NCPHI),SPHCAL(NCPHI),YCMIN,YCMAX
0805       PARAMETER (NJMAX=500)
0806       COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(NCY,NCPHI),NCJET
0807       DOUBLE PRECISION ETAJCUT,PSERAP
0808 C
0809 C          PARAMETERS
0810       DATA ECCUT/0.1D0/
0811       DATA ETSTOP/1.5D0/
0812       DATA ITLIS/6/
0813 C
0814 C          INITIALIZE
0815 C
0816       DO 100 IPHI=1,NCPHI
0817       DO 100 IY=1,NCY
0818 100   JETNO(IY,IPHI)=0
0819       DO 110 J=1,NJMAX
0820       ETJET(J)=0.
0821       DO 110 K=1,4
0822 110   PCJET(K,J)=0.
0823       NCJET=0
0824       NPHI1=RJET/DELPHI
0825       NPHI2=2*NPHI1+1
0826       NY1=RJET/DELY
0827       NY2=2*NY1+1
0828       IPASS=0
0829 C
0830 C          FIND HIGHEST CELL REMAINING
0831 C
0832 1     ETMAX=0.
0833       DO 200 IPHI=1,NCPHI
0834       DO 210 IY=1,NCY
0835       IF(ET(IY,IPHI).LT.ETMAX) GOTO 210
0836       IF(JETNO(IY,IPHI).NE.0) GOTO 210
0837       ETMAX=ET(IY,IPHI)
0838       IYMX=IY
0839       IPHIMX=IPHI
0840 210   CONTINUE
0841 200   CONTINUE
0842       IF(ETMAX.LT.ETSTOP) RETURN
0843 C
0844 C          SUM CELLS
0845 C
0846       IPASS=IPASS+1
0847       IF(IPASS.GT.NCY*NCPHI) THEN
0848         WRITE(ITLIS,8888) IPASS
0849 8888    FORMAT(//' ERROR IN GETJET...IPASS > ',I6)
0850         RETURN
0851       ENDIF
0852       NCJET=NCJET+1
0853       IF(NCJET.GT.NJMAX) THEN
0854         WRITE(ITLIS,9999) NCJET
0855 9999    FORMAT(//' ERROR IN GETJET...NCJET > ',I5)
0856         RETURN
0857       ENDIF
0858       DO 300 IPHI1=1,NPHI2
0859       IPHIX=IPHIMX-NPHI1-1+IPHI1
0860       IF(IPHIX.LE.0) IPHIX=IPHIX+NCPHI
0861       IF(IPHIX.GT.NCPHI) IPHIX=IPHIX-NCPHI
0862       DO 310 IY1=1,NY2
0863       IYX=IYMX-NY1-1+IY1
0864       IF(IYX.LE.0) GOTO 310
0865       IF(IYX.GT.NCY) GOTO 310
0866       IF(JETNO(IYX,IPHIX).NE.0) GOTO 310
0867       RR=(DELY*(IY1-NY1-1))**2+(DELPHI*(IPHI1-NPHI1-1))**2
0868       IF(RR.GT.RJET**2) GOTO 310
0869       IF(ET(IYX,IPHIX).LT.ECCUT) GOTO 310
0870       PX=ET(IYX,IPHIX)/STHCAL(IYX)
0871 C          ADD CELL TO JET
0872       PCJET(1,NCJET)=PCJET(1,NCJET)+PX*STHCAL(IYX)*CPHCAL(IPHIX)
0873       PCJET(2,NCJET)=PCJET(2,NCJET)+PX*STHCAL(IYX)*SPHCAL(IPHIX)
0874       PCJET(3,NCJET)=PCJET(3,NCJET)+PX*CTHCAL(IYX)
0875       PCJET(4,NCJET)=PCJET(4,NCJET)+PX
0876       ETJET(NCJET)=ETJET(NCJET)+ET(IYX,IPHIX)
0877       JETNO(IYX,IPHIX)=NCJET
0878 310   CONTINUE
0879 300   CONTINUE
0880 C
0881 C          DISCARD JET IF ET < EJCUT.
0882 C
0883       IF(ETJET(NCJET).GT.EJCUT.AND.ABS(PSERAP(PCJET(1,NCJET))).LT
0884      $     .ETAJCUT) GOTO 1
0885       ETJET(NCJET)=0.
0886       DO 400 K=1,4
0887 400   PCJET(K,NCJET)=0.
0888       NCJET=NCJET-1
0889       GOTO 1
0890       END
0891 C-----------------------------------------------------------------------
0892       SUBROUTINE CALDEL(ISTLO,ISTHI,ISTOP)
0893 C     LABEL ALL PARTICLES WITH STATUS BETWEEN ISTLO AND ISTHI (UNTIL A
0894 C     PARTICLE WITH STATUS ISTOP IS FOUND) AS FINAL-STATE, CALL CALSIM
0895 C     AND THEN PUT LABELS BACK TO NORMAL
0896 C
0897 C     THIS VERSION LEAVES OUT PARTICLES THAT POINT BACK TO MOTHERS TO BE
0898 C     LEFT OUT OF MATCHING
0899 C-----------------------------------------------------------------------
0900       IMPLICIT NONE
0901       INTEGER MAXNUP
0902       PARAMETER(MAXNUP=500)
0903       INTEGER IEND,INORAD
0904       COMMON/SHVETO/IEND,INORAD(MAXNUP)
0905 C...HEPEVT commonblock.
0906       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0907       PARAMETER (NMXHEP=4000)
0908       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0909      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0910       DOUBLE PRECISION PHEP,VHEP
0911       SAVE /HEPEVT/
0912       INTEGER ISTOLD(NMXHEP),IHEP,IST,ISTLO,ISTHI,ISTOP,IMO
0913       LOGICAL FOUND
0914 c      write(3,*) 'new event',nevhep
0915       FOUND=.FALSE.
0916       DO 10 IHEP=1,NHEP
0917         IST=ISTHEP(IHEP)
0918         ISTOLD(IHEP)=IST
0919         IF (IST.EQ.ISTOP) FOUND=.TRUE.
0920         IF (IST.GE.ISTLO.AND.IST.LE.ISTHI.AND..NOT.FOUND) THEN
0921 C     FOUND A RADIATED PARTON, CHECK MOTHER
0922           IMO=IHEP
0923  1        IMO=JMOHEP(1,IMO)
0924           IF(IMO.EQ.IEND) THEN
0925 C     PARENTHOOD OK
0926             IST=1
0927 c            write(3,*) ihep,ist
0928             GOTO 9
0929           ENDIF
0930           IF(INORAD(IMO).EQ.1) THEN
0931 C     PARTON COMES FROM A VETOED MOTHER
0932             IST=0
0933             GOTO 9
0934           ELSE
0935 C     CHECK GRANDMOTHER
0936             GOTO 1
0937           ENDIF
0938         ELSE
0939           IST=0
0940         ENDIF
0941  9      ISTHEP(IHEP)=IST
0942  10   CONTINUE
0943
0944 C     PRINT PARTICLE STATUS FOR THE CALORIMETER SIMULATION.
0945 C     (UNCOMMENT TO MAKE IT WORK)
0946 C      CALL DBCAL
0947
0948       CALL CALSIM
0949       DO 20 IHEP=1,NHEP
0950         ISTHEP(IHEP)=ISTOLD(IHEP)
0951  20   CONTINUE
0952       END
0953 C-----------------------------------------------------------------------
0954       SUBROUTINE CALDEL_HVQ(ISTLO,ISTHI,ISTOP)
0955 C     LABEL ALL PARTICLES WITH STATUS BETWEEN ISTLO AND ISTHI (UNTIL A
0956 C     PARTICLE WITH STATUS ISTOP IS FOUND) AS FINAL-STATE, CALL CALSIM
0957 C     AND THEN PUT LABELS BACK TO NORMAL
0958 C
0959 C     THIS VERSION KEEPS ONLY ALL IST=1 PARTICLES REJECTED BY CALDEL AS
0960 C     DAUGHTERS OF VETOED HEAVY-QUARK MOTHERS: JETS COMPLEMENTARY TO
0961 C     THOSE RECONSTRUCTED BY CALDEL
0962 C-----------------------------------------------------------------------
0963       IMPLICIT NONE
0964       INTEGER MAXNUP
0965       PARAMETER(MAXNUP=500)
0966       INTEGER IEND,INORAD
0967       COMMON/SHVETO/IEND,INORAD(MAXNUP)
0968 C...HEPEVT commonblock.
0969       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0970       PARAMETER (NMXHEP=4000)
0971       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0972      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0973       DOUBLE PRECISION PHEP,VHEP
0974       SAVE /HEPEVT/
0975       INTEGER ISTOLD(NMXHEP),IHEP,IST,ISTLO,ISTHI,ISTOP,IMO
0976       INTEGER IDMOTH,IDDAUG
0977       LOGICAL FOUND
0978       FOUND=.FALSE.
0979 C TAG W-B FROM TOP DECAYS
0980       DO IHEP=1,NHEP
0981         IMO=JMOHEP(1,IHEP)
0982         IDMOTH=ABS(IDHEP(IMO))
0983         IDDAUG=ABS(IDHEP(IHEP))
0984         IF(IDMOTH.EQ.6.AND.(IDDAUG.EQ.5.OR.IDDAUG.EQ.24)) INORAD(IHEP)=2
0985       ENDDO
0986 C
0987       DO 10 IHEP=1,NHEP
0988         IST=ISTHEP(IHEP)
0989         ISTOLD(IHEP)=IST
0990         IF (IST.EQ.ISTOP) FOUND=.TRUE.
0991         IF (IST.GE.ISTLO.AND.IST.LE.ISTHI.AND..NOT.FOUND) THEN
0992 C     FOUND A RADIATED PARTON, CHECK MOTHER
0993           IMO=IHEP
0994  1        IMO=JMOHEP(1,IMO)
0995           IF(IMO.EQ.IEND) THEN
0996 C     PARENTHOOD OK, VETO
0997             IST=0
0998             GOTO 9
0999           ENDIF
1000           IF(INORAD(IMO).EQ.1) THEN
1001             IDMOTH=ABS(IDHEP(IMO))
1002             IDDAUG=ABS(IDHEP(IHEP))
1003 C     VERIFY IT'S A HEAVY QUARK -- LEAVE OUT GAUGE BOSON DECAYS
1004             IF(IDMOTH.GE.4.AND.IDMOTH.LE.6) THEN
1005 C     PARTON COMES FROM A VETOED MOTHER, KEEP UNLESS IT IS THE EVOLVED
1006 C     MOTHER
1007               IF(IDMOTH.NE.IDDAUG) THEN
1008                 IST=1
1009                 GOTO 9
1010               ELSE
1011                 IST=0
1012               ENDIF
1013             ELSE
1014 C     NOT A HEAVY QUARK MOTHER, LEAVE OUT OF JET RECONSUTRCTION
1015               IST=0
1016             ENDIF
1017           ELSEIF(INORAD(IMO).EQ.2) THEN
1018 C     IT'S A TOP DECAY PRODUCT
1019             IST=0
1020           ELSE
1021 C     GO CHECK GRANDMOTHER
1022             GOTO 1
1023           ENDIF
1024         ELSE
1025           IST=0
1026         ENDIF
1027  9      ISTHEP(IHEP)=IST
1028  10   CONTINUE
1029       CALL CALSIM
1030       DO 20 IHEP=1,NHEP
1031         ISTHEP(IHEP)=ISTOLD(IHEP)
1032  20   CONTINUE
1033       END
1034 C-----------------------------------------------------------------------
1035       FUNCTION PSERAP(P)
1036 C     PSEUDO-RAPIDITY (-LOG TAN THETA/2)
1037 C-----------------------------------------------------------------------
1038       DOUBLE PRECISION PSERAP,P(3),PT,PL,TINY,THETA
1039       PARAMETER (TINY=1D-3)
1040       PT=SQRT(P(1)**2+P(2)**2)+TINY
1041       PL=P(3)
1042       THETA=ATAN2(PT,PL)
1043       PSERAP=-LOG(TAN(0.5*THETA))
1044       END
1045 C-----------------------------------------------------------------------
1046 C-----------------------------------------------------------------------
1047       SUBROUTINE ALPSOR(A,N,K,IOPT)
1048 C-----------------------------------------------------------------------
1049 C     Sort A(N) into ascending order
1050 C     IOPT = 1 : return sorted A and index array K
1051 C     IOPT = 2 : return index array K only
1052 C-----------------------------------------------------------------------
1053       DOUBLE PRECISION A(N),B(5000)
1054       INTEGER N,I,J,IOPT,K(N),IL(5000),IR(5000)
1055       IF (N.GT.5000) then
1056         write(*,*) 'Too many entries to sort in alpsrt, stop'
1057         stop
1058       endif
1059       if(n.le.0) return
1060       IL(1)=0
1061       IR(1)=0
1062       DO 10 I=2,N
1063       IL(I)=0
1064       IR(I)=0
1065       J=1
1066    2  IF(A(I).GT.A(J)) GOTO 5
1067    3  IF(IL(J).EQ.0) GOTO 4
1068       J=IL(J)
1069       GOTO 2
1070    4  IR(I)=-J
1071       IL(J)=I
1072       GOTO 10
1073    5  IF(IR(J).LE.0) GOTO 6
1074       J=IR(J)
1075       GOTO 2
1076    6  IR(I)=IR(J)
1077       IR(J)=I
1078   10  CONTINUE
1079       I=1
1080       J=1
1081       GOTO 8
1082   20  J=IL(J)
1083    8  IF(IL(J).GT.0) GOTO 20
1084    9  K(I)=J
1085       B(I)=A(J)
1086       I=I+1
1087
1088 C Compatibility with gfortran (CMSSW 3_X)
1089 C      IF(IR(J)) 12,30,13
1090 C  13  J=IR(J)
1091 C      GOTO 8
1092 C  12  J=-IR(J)
1093 C      GOTO 9
1094 C  30  IF(IOPT.EQ.2) RETURN
1095
1096       IF (IR(J).GT.0) THEN
1097          J=IR(J)
1098          GOTO 8
1099       ELSEIF (IR(J).LT.0) then
1100          J=-IR(J)
1101          GOTO 9
1102       ENDIF
1103       IF (IOPT.EQ.2) RETURN
1104
1105       DO 31 I=1,N
1106   31  A(I)=B(I)
1107  999  END
1108
1109 C-----------------------------------------------------------------------
1110       subroutine getunit(n)
1111       implicit none
1112       integer n,i
1113       logical yes
1114       WRITE(*,*)'CALLED GETUNIT'
1115       do i=10,100
1116         inquire(unit=i,opened=yes)
1117         if(.not.yes) goto 10
1118       enddo
1119       write(*,*) 'no free units to write to available, stop'
1120       stop
1121  10   n=i
1122       end