PartonShowerVeto/src/getjet.f

0001 C-----------------------------------------------------------------------
0002 C----Calorimeter simulation obtained from Frank Paige 23 March 1988-----
0003 C
0004 C          USE
0005 C
0006 C     CALL CALINI
0007 C     CALL CALSIM
0008 C
0009 C          THEN TO FIND JETS WITH A SIMPLIFIED VERSION OF THE UA1 JET
0010 C          ALGORITHM WITH JET RADIUS RJET AND MINIMUM SCALAR TRANSVERSE
0011 C          ENERGY EJCUT
0012 C            (RJET=1., EJCUT=5. FOR UA1)
0013 C          USE
0014 C
0015 C     CALL GETJET(RJET,EJCUT)
0016 C
0017 C
0018 C-----------------------------------------------------------------------
0019 C
0020 C          ADDED BY MIKE SEYMOUR: PARTON-LEVEL CALORIMETER. ALL PARTONS
0021 C          ARE CONSIDERED TO BE HADRONS, SO IN FACT RESEM IS IGNORED
0022 C
0023 C     CALL CALPAR
0024 C
0025 C          HARD PARTICLE CALORIMETER. ONLY USES THOSE PARTICLES WHICH
0026 C          CAME FROM THE HARD PROCESS, AND NOT THE UNDERLYING EVENT
0027 C
0028 C     CALL CALHAR
0029 C
0030 C-----------------------------------------------------------------------
0031       SUBROUTINE CALINI
0032 C
0033 C          INITIALIZE CALORIMETER FOR CALSIM AND GETJET.  NOTE THAT
0034 C          BECAUSE THE INITIALIZATION IS SEPARATE, CALSIM CAN BE
0035 C          CALLED MORE THAN ONCE TO SIMULATE PILEUP OF SEVERAL EVENTS.
0036 C
0037       IMPLICIT NONE
0038 C...GETJET commonblocks
0039       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
0040       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
0041      &  SPHCAL,PCJET,ETJET
0042       PARAMETER (MNCY=200)
0043       PARAMETER (MNCPHI=200)
0044       COMMON/CALOR/DELY,DELPHI,ET(MNCY,MNCPHI),
0045      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
0046      $YCMIN,YCMAX,NCY,NCPHI
0047       PARAMETER (NJMAX=500)
0048       COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),NCJET
0049
0050       INTEGER IPHI,IY
0051       DOUBLE PRECISION PI,PHIX,YX,THX
0052       PARAMETER (PI=3.141593D0)
0053       LOGICAL FSTCAL
0054       DATA FSTCAL/.TRUE./
0055 C
0056 C          INITIALIZE ET ARRAY.
0057       DO 100 IPHI=1,NCPHI
0058       DO 100 IY=1,NCY
0059 100   ET(IY,IPHI)=0.
0060 C
0061       IF (FSTCAL) THEN
0062 C          CALCULATE TRIG. FUNCTIONS.
0063         DELPHI=2.*PI/FLOAT(NCPHI)
0064         DO 200 IPHI=1,NCPHI
0065         PHIX=DELPHI*(IPHI-.5)
0066         CPHCAL(IPHI)=COS(PHIX)
0067         SPHCAL(IPHI)=SIN(PHIX)
0068 200     CONTINUE
0069         DELY=(YCMAX-YCMIN)/FLOAT(NCY)
0070         DO 300 IY=1,NCY
0071         YX=DELY*(IY-.5)+YCMIN
0072         THX=2.*ATAN(EXP(-YX))
0073         CTHCAL(IY)=COS(THX)
0074         STHCAL(IY)=SIN(THX)
0075 300     CONTINUE
0076         FSTCAL=.FALSE.
0077       ENDIF
0078       END
0079 C
0080       SUBROUTINE CALSIM
0081 C
0082 C          SIMPLE CALORIMETER SIMULATION.  ASSUME UNIFORM Y AND PHI
0083 C          BINS
0084 C...HEPEVT commonblock.
0085       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0086       PARAMETER (NMXHEP=4000)
0087       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0088      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0089       DOUBLE PRECISION PHEP,VHEP
0090       SAVE /HEPEVT/
0091
0092 C...GETJET commonblocks
0093       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
0094       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
0095      &  SPHCAL,PCJET,ETJET
0096       PARAMETER (MNCY=200)
0097       PARAMETER (MNCPHI=200)
0098       COMMON/CALOR/DELY,DELPHI,ET(MNCY,MNCPHI),
0099      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
0100      $YCMIN,YCMAX,NCY,NCPHI
0101       PARAMETER (NJMAX=500)
0102       COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),NCJET
0103
0104       INTEGER IHEP,ID,IY,IPHI
0105       DOUBLE PRECISION PI,YIP,PSERAP,PHIIP,EIP
0106       PARAMETER (PI=3.141593D0)
0107 C
0108 C          FILL CALORIMETER
0109 C
0110       DO 200 IHEP=1,NHEP
0111       IF (ISTHEP(IHEP).EQ.1) THEN
0112         YIP=PSERAP(PHEP(1,IHEP))
0113         IF(YIP.LT.YCMIN.OR.YIP.GT.YCMAX) GOTO 200
0114         ID=ABS(IDHEP(IHEP))
0115 C---EXCLUDE TOP QUARK, LEPTONS, PROMPT PHOTONS
0116         IF ((ID.GE.11.AND.ID.LE.16).OR.ID.EQ.6.OR.ID.EQ.22) GOTO 200
0117 C
0118         PHIIP=ATAN2(PHEP(2,IHEP),PHEP(1,IHEP))
0119         IF(PHIIP.LT.0.) PHIIP=PHIIP+2.*PI
0120         IY=INT((YIP-YCMIN)/DELY)+1
0121         IPHI=INT(PHIIP/DELPHI)+1
0122         EIP=PHEP(4,IHEP)
0123 C            WEIGHT BY SIN(THETA)
0124         ET(IY,IPHI)=ET(IY,IPHI)+EIP*STHCAL(IY)
0125       ENDIF
0126   200 CONTINUE
0127       END
0128       SUBROUTINE GETJET(RJET,EJCUT,ETAJCUT)
0129 C
0130 C          SIMPLE JET-FINDING ALGORITHM (SIMILAR TO UA1).
0131 C
0132 C     FIND HIGHEST REMAINING CELL > ETSTOP AND SUM SURROUNDING
0133 C          CELLS WITH--
0134 C            DELTA(Y)**2+DELTA(PHI)**2<RJET**2
0135 C            ET>ECCUT.
0136 C          KEEP JETS WITH ET>EJCUT AND ABS(ETA)<ETAJCUT
0137 C          THE UA1 PARAMETERS ARE RJET=1.0 AND EJCUT=5.0
0138 C
0139       IMPLICIT NONE
0140 C...GETJET commonblocks
0141       INTEGER MNCY,MNCPHI,NCY,NCPHI,NJMAX,JETNO,NCJET
0142       DOUBLE PRECISION YCMIN,YCMAX,DELY,DELPHI,ET,STHCAL,CTHCAL,CPHCAL,
0143      &  SPHCAL,PCJET,ETJET
0144       PARAMETER (MNCY=200)
0145       PARAMETER (MNCPHI=200)
0146       COMMON/CALOR/DELY,DELPHI,ET(MNCY,MNCPHI),
0147      $CTHCAL(MNCY),STHCAL(MNCY),CPHCAL(MNCPHI),SPHCAL(MNCPHI),
0148      $YCMIN,YCMAX,NCY,NCPHI
0149       PARAMETER (NJMAX=500)
0150       COMMON/GETCOM/PCJET(4,NJMAX),ETJET(NJMAX),JETNO(MNCY,MNCPHI),NCJET
0151
0152       INTEGER IPHI,IY,J,K,NPHI1,NPHI2,NY1,
0153      &  NY2,IPASS,IYMX,IPHIMX,ITLIS,IPHI1,IPHIX,IY1,IYX
0154       DOUBLE PRECISION PI,RJET,
0155      &  ETMAX,ETSTOP,RR,ECCUT,PX,EJCUT
0156       PARAMETER (PI=3.141593D0)
0157       DOUBLE PRECISION ETAJCUT,PSERAP
0158 C
0159 C          PARAMETERS
0160       DATA ECCUT/0.1D0/
0161       DATA ETSTOP/1.5D0/
0162       DATA ITLIS/6/
0163 C
0164 C          INITIALIZE
0165 C
0166
0167 c  init uninit variables
0168       iymx = 0
0169       iphimx = 0
0170 c
0171       DO 100 IPHI=1,NCPHI
0172       DO 100 IY=1,NCY
0173 100   JETNO(IY,IPHI)=0
0174       DO 110 J=1,NJMAX
0175       ETJET(J)=0.
0176       DO 110 K=1,4
0177 110   PCJET(K,J)=0.
0178       NCJET=0
0179       NPHI1=RJET/DELPHI
0180       NPHI2=2*NPHI1+1
0181       NY1=RJET/DELY
0182       NY2=2*NY1+1
0183       IPASS=0
0184 C
0185 C          FIND HIGHEST CELL REMAINING
0186 C
0187 1     ETMAX=0.
0188       DO 200 IPHI=1,NCPHI
0189       DO 210 IY=1,NCY
0190       IF(ET(IY,IPHI).LT.ETMAX) GOTO 210
0191       IF(JETNO(IY,IPHI).NE.0) GOTO 210
0192       ETMAX=ET(IY,IPHI)
0193       IYMX=IY
0194       IPHIMX=IPHI
0195 210   CONTINUE
0196 200   CONTINUE
0197       IF(ETMAX.LT.ETSTOP) RETURN
0198 C
0199 C          SUM CELLS
0200 C
0201       IPASS=IPASS+1
0202       IF(IPASS.GT.NCY*NCPHI) THEN
0203         WRITE(ITLIS,8888) IPASS
0204 8888    FORMAT(//' ERROR IN GETJET...IPASS > ',I6)
0205         RETURN
0206       ENDIF
0207       NCJET=NCJET+1
0208       IF(NCJET.GT.NJMAX) THEN
0209         WRITE(ITLIS,9999) NCJET
0210 9999    FORMAT(//' ERROR IN GETJET...NCJET > ',I5)
0211         RETURN
0212       ENDIF
0213       DO 300 IPHI1=1,NPHI2
0214       IPHIX=IPHIMX-NPHI1-1+IPHI1
0215       IF(IPHIX.LE.0) IPHIX=IPHIX+NCPHI
0216       IF(IPHIX.GT.NCPHI) IPHIX=IPHIX-NCPHI
0217       DO 310 IY1=1,NY2
0218       IYX=IYMX-NY1-1+IY1
0219       IF(IYX.LE.0) GOTO 310
0220       IF(IYX.GT.NCY) GOTO 310
0221       IF(JETNO(IYX,IPHIX).NE.0) GOTO 310
0222       RR=(DELY*(IY1-NY1-1))**2+(DELPHI*(IPHI1-NPHI1-1))**2
0223       IF(RR.GT.RJET**2) GOTO 310
0224       IF(ET(IYX,IPHIX).LT.ECCUT) GOTO 310
0225       PX=ET(IYX,IPHIX)/STHCAL(IYX)
0226 C          ADD CELL TO JET
0227       PCJET(1,NCJET)=PCJET(1,NCJET)+PX*STHCAL(IYX)*CPHCAL(IPHIX)
0228       PCJET(2,NCJET)=PCJET(2,NCJET)+PX*STHCAL(IYX)*SPHCAL(IPHIX)
0229       PCJET(3,NCJET)=PCJET(3,NCJET)+PX*CTHCAL(IYX)
0230       PCJET(4,NCJET)=PCJET(4,NCJET)+PX
0231       ETJET(NCJET)=ETJET(NCJET)+ET(IYX,IPHIX)
0232       JETNO(IYX,IPHIX)=NCJET
0233 310   CONTINUE
0234 300   CONTINUE
0235 C
0236 C          DISCARD JET IF ET < EJCUT.
0237 C
0238       IF(ETJET(NCJET).GT.EJCUT.AND.ABS(PSERAP(PCJET(1,NCJET))).LT
0239      $     .ETAJCUT) GOTO 1
0240       ETJET(NCJET)=0.
0241       DO 400 K=1,4
0242 400   PCJET(K,NCJET)=0.
0243       NCJET=NCJET-1
0244       GOTO 1
0245       END
0246 C-----------------------------------------------------------------------
0247       SUBROUTINE CALDEL(ISTLO,ISTHI)
0248 C     LABEL ALL PARTICLES WITH STATUS BETWEEN ISTLO AND ISTHI (UNTIL A
0249 C     PARTICLE WITH STATUS ISTOP IS FOUND) AS FINAL-STATE, CALL CALSIM
0250 C     AND THEN PUT LABELS BACK TO NORMAL
0251 C-----------------------------------------------------------------------
0252       IMPLICIT NONE
0253       INTEGER MAXNUP
0254       PARAMETER(MAXNUP=500)
0255 C...HEPEVT commonblock.
0256       INTEGER NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0257       PARAMETER (NMXHEP=4000)
0258       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0259      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0260       DOUBLE PRECISION PHEP,VHEP
0261       SAVE /HEPEVT/
0262       INTEGER ISTLO,ISTHI
0263 C...Inputs for the matching algorithm
0264       double precision etcjet,rclmax,etaclmax,qcut,qfact
0265       integer maxjets,minjets,iexcfile,ktsche
0266       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,qfact,
0267      $   maxjets,minjets,iexcfile,ktsche
0268
0269 C...Silence gcc 4.3.4 warnings.
0270       ISTLO=ISTLO
0271       ISTHI=ISTHI
0272       CALL CALSIM
0273       END
0274 C-----------------------------------------------------------------------
0275       FUNCTION PSERAP(P)
0276 C     PSEUDO-RAPIDITY (-LOG TAN THETA/2)
0277 C-----------------------------------------------------------------------
0278       DOUBLE PRECISION PSERAP,P(3),PT,PL,TINY,THETA
0279       PARAMETER (TINY=1D-3)
0280       PT=SQRT(P(1)**2+P(2)**2)+TINY
0281       PL=P(3)
0282       THETA=ATAN2(PT,PL)
0283       PSERAP=-LOG(TAN(0.5*THETA))
0284       END
0285 C-----------------------------------------------------------------------