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